|
| 1 | +from numpy import * |
| 2 | +from numpy.linalg import * |
| 3 | +from scipy.integrate import solve_ivp |
| 4 | +from matplotlib.pyplot import * |
| 5 | +# Python 3.x Standard Library |
| 6 | +import gc |
| 7 | +import os |
| 8 | + |
| 9 | +# Third-Party Packages |
| 10 | +import numpy as np; np.seterr(all="ignore") |
| 11 | +import numpy.linalg as la |
| 12 | +import scipy.misc |
| 13 | +import matplotlib as mpl; mpl.use("Agg") |
| 14 | +import matplotlib.pyplot as pp |
| 15 | +import matplotlib.axes as ax |
| 16 | +import matplotlib.patches as pa |
| 17 | + |
| 18 | + |
| 19 | +# |
| 20 | +# Matplotlib Configuration & Helper Functions |
| 21 | +# -------------------------------------------------------------------------- |
| 22 | + |
| 23 | +# TODO: also reconsider line width and markersize stuff "for the web |
| 24 | +# settings". |
| 25 | +fontsize = 10 |
| 26 | + |
| 27 | +width = 345 / 72.27 |
| 28 | +height = width / (16/9) |
| 29 | + |
| 30 | +rc = { |
| 31 | + "text.usetex": True, |
| 32 | + "pgf.preamble": r"\usepackage{amsmath,amsfonts,amssymb}", |
| 33 | + #"font.family": "serif", |
| 34 | + "font.serif": [], |
| 35 | + #"font.sans-serif": [], |
| 36 | + "legend.fontsize": fontsize, |
| 37 | + "axes.titlesize": fontsize, |
| 38 | + "axes.labelsize": fontsize, |
| 39 | + "xtick.labelsize": fontsize, |
| 40 | + "ytick.labelsize": fontsize, |
| 41 | + "figure.max_open_warning": 100, |
| 42 | + #"savefig.dpi": 300, |
| 43 | + #"figure.dpi": 300, |
| 44 | + "figure.figsize": [width, height], |
| 45 | + "lines.linewidth": 1.0, |
| 46 | +} |
| 47 | +mpl.rcParams.update(rc) |
| 48 | + |
| 49 | +# Web target: 160 / 9 inches (that's ~45 cm, this is huge) at 90 dpi |
| 50 | +# (the "standard" dpi for Web computations) gives 1600 px. |
| 51 | +width_in = 160 / 9 |
| 52 | + |
| 53 | +def save(name, **options): |
| 54 | + cwd = os.getcwd() |
| 55 | + root = os.path.dirname(os.path.realpath(__file__)) |
| 56 | + os.chdir(root) |
| 57 | + pp.savefig(name + ".svg", **options) |
| 58 | + os.chdir(cwd) |
| 59 | + |
| 60 | +def set_ratio(ratio=1.0, bottom=0.1, top=0.1, left=0.1, right=0.1): |
| 61 | + height_in = (1.0 - left - right)/(1.0 - bottom - top) * width_in / ratio |
| 62 | + pp.gcf().set_size_inches((width_in, height_in)) |
| 63 | + pp.gcf().subplots_adjust(bottom=bottom, top=1.0-top, left=left, right=1.0-right) |
| 64 | + |
| 65 | +width |
| 66 | +def Q(f, xs, ys): |
| 67 | + X, Y = meshgrid(xs, ys) |
| 68 | + fx = vectorize(lambda x, y: f([x, y])[0]) |
| 69 | + fy = vectorize(lambda x, y: f([x, y])[1]) |
| 70 | + return X, Y, fx(X, Y), fy(X, Y) |
| 71 | +def fun(t, y): |
| 72 | + return y * y |
| 73 | +t0, tf, y0 = 0.0, 3.0, array([1.0]) |
| 74 | +result = solve_ivp(fun, t_span=[t0, tf], y0=y0) |
| 75 | +figure() |
| 76 | +plot(result["t"], result["y"][0], "k") |
| 77 | +xlim(t0, tf); xlabel("$t$"); ylabel("$x(t)$") |
| 78 | +tight_layout() |
| 79 | +save("images/finite-time-blowup") |
| 80 | +tf = 1.0 |
| 81 | +r = solve_ivp(fun, [t0, tf], y0, |
| 82 | + dense_output=True) |
| 83 | +figure() |
| 84 | +t = linspace(t0, tf, 1000) |
| 85 | +plot(t, r["sol"](t)[0], "k") |
| 86 | +ylim(0.0, 10.0); grid(); |
| 87 | +xlabel("$t$"); ylabel("$x(t)$") |
| 88 | +tight_layout() |
| 89 | +save("images/finite-time-blowup-2") |
| 90 | +def f(x1x2): |
| 91 | + x1, x2 = x1x2 |
| 92 | + dx1 = 1.0 if x1 < 0.0 else -1.0 |
| 93 | + return array([dx1, 0.0]) |
| 94 | +figure() |
| 95 | +x1 = x2 = linspace(-1.0, 1.0, 20) |
| 96 | +gca().set_aspect(1.0) |
| 97 | +quiver(*Q(f, x1, x2), color="k") |
| 98 | +tight_layout() |
| 99 | +save("images/discont") |
| 100 | +def sigma(x): |
| 101 | + return 1 / (1 + exp(-x)) |
| 102 | +figure() |
| 103 | +x = linspace(-7.0, 7.0, 1000) |
| 104 | +plot(x, sigma(x), label="$y=\sigma(x)$") |
| 105 | +grid(True) |
| 106 | +xlim(-5, 5) |
| 107 | +xticks([-5.0, 0.0, 5.0]) |
| 108 | +yticks([0.0, 0.5, 1.0]) |
| 109 | +xlabel("$x$") |
| 110 | +ylabel("$y$") |
| 111 | +legend() |
| 112 | +pp.gcf().subplots_adjust(bottom=0.2) |
| 113 | +save("images/sigmoid") |
| 114 | +alpha = 2 / 3; beta = 4 / 3; delta = gamma = 1.0 |
| 115 | + |
| 116 | +def fun(t, y): |
| 117 | + x, y = y |
| 118 | + u = alpha * x - beta * x * y |
| 119 | + v = delta * x * y - gamma * y |
| 120 | + return array([u, v]) |
| 121 | +tf = 3.0 |
| 122 | +result = solve_ivp( |
| 123 | + fun, |
| 124 | + t_span=(0.0, tf), |
| 125 | + y0=[1.5, 1.5], |
| 126 | + max_step=0.01) |
| 127 | +x, y = result["y"][0], result["y"][1] |
| 128 | +def display_streamplot(): |
| 129 | + ax = gca() |
| 130 | + xr = yr = linspace(0.0, 2.0, 1000) |
| 131 | + def f(y): |
| 132 | + return fun(0, y) |
| 133 | + streamplot(*Q(f, xr, yr), color="grey") |
| 134 | +def display_reference_solution(): |
| 135 | + for xy in zip(x, y): |
| 136 | + x_, y_ = xy |
| 137 | + gca().add_artist(Circle((x_, y_), |
| 138 | + 0.2, color="#d3d3d3")) |
| 139 | + gca().add_artist(Circle((x[0], y[0]), 0.1, |
| 140 | + color="#808080")) |
| 141 | + plot(x, y, "k") |
| 142 | +def display_alternate_solution(): |
| 143 | + result = solve_ivp(fun, |
| 144 | + t_span=[0.0, tf], |
| 145 | + y0=[1.5, 1.575], |
| 146 | + max_step=0.01) |
| 147 | + x, y = result["y"][0], result["y"][1] |
| 148 | + plot(x, y, "k--") |
| 149 | +figure() |
| 150 | +display_streamplot() |
| 151 | +display_reference_solution() |
| 152 | +display_alternate_solution() |
| 153 | +axis([0,2,0,2]); axis("square") |
| 154 | +save("images/continuity") |
| 155 | +def fun(t, y): |
| 156 | + x = y[0] |
| 157 | + dx = sqrt(abs(y)) |
| 158 | + return [dx] |
| 159 | +tspan = [0.0, 3.0] |
| 160 | +t = linspace(tspan[0], tspan[1], 1000) |
| 161 | +figure() |
| 162 | +for x0 in [0.1, 0.01, 0.001, 0.0001, 0.0]: |
| 163 | + r = solve_ivp(fun, tspan, [x0], |
| 164 | + dense_output=True) |
| 165 | + plot(t, r["sol"](t)[0], |
| 166 | + label=f"$x_0 = {x0}$") |
| 167 | +xlabel("$t$"); ylabel("$x(t)$") |
| 168 | +legend() |
| 169 | +pp.gcf().subplots_adjust(bottom=0.2) |
| 170 | +save("images/eps") |
0 commit comments