|
1 | 1 | from numpy import * |
2 | 2 | from numpy.linalg import * |
3 | | -from scipy.linalg import * |
4 | | -from matplotlib.pyplot import * |
5 | | -from mpl_toolkits.mplot3d import * |
6 | 3 | from scipy.integrate import solve_ivp |
| 4 | +from matplotlib.pyplot import * |
7 | 5 | # Python 3.x Standard Library |
8 | 6 | import gc |
9 | 7 | import os |
@@ -63,206 +61,110 @@ def set_ratio(ratio=1.0, bottom=0.1, top=0.1, left=0.1, right=0.1): |
63 | 61 | height_in = (1.0 - left - right)/(1.0 - bottom - top) * width_in / ratio |
64 | 62 | pp.gcf().set_size_inches((width_in, height_in)) |
65 | 63 | pp.gcf().subplots_adjust(bottom=bottom, top=1.0-top, left=left, right=1.0-right) |
| 64 | + |
| 65 | +width |
66 | 66 | def Q(f, xs, ys): |
67 | 67 | X, Y = meshgrid(xs, ys) |
68 | | - v = vectorize |
69 | | - fx = v(lambda x, y: f([x, y])[0]) |
70 | | - fy = v(lambda x, y: f([x, y])[1]) |
| 68 | + fx = vectorize(lambda x, y: f([x, y])[0]) |
| 69 | + fy = vectorize(lambda x, y: f([x, y])[1]) |
71 | 70 | return X, Y, fx(X, Y), fy(X, Y) |
72 | | -a = 2.0; x0 = 1.0 |
73 | | -figure() |
74 | | -t = linspace(0.0, 3.0, 1000) |
75 | | -plot(t, exp(a*t)*x0, "k") |
76 | | -xlabel("$t$"); ylabel("$x(t)$"); title(f"$a={a}$") |
77 | | -grid(); axis([0.0, 2.0, 0.0, 10.0]) |
78 | | -tight_layout() |
79 | | -save("images/scalar-LTI-2") |
80 | | -figure() |
81 | | -plot(real(a), imag(a), "x", color="k") |
82 | | -gca().set_aspect(1.0) |
83 | | -xlim(-3,3); ylim(-3,3); |
84 | | -plot([-3,3], [0,0], "k") |
85 | | -plot([0, 0], [-3, 3], "k") |
86 | | -xticks([-2,-1,0,1,2]); yticks([-2,-1,0,1,2]) |
87 | | -title(f"$a={a}$"); grid(True) |
88 | | -tight_layout() |
89 | | -save("images/scalar-LTI-2-poles") |
90 | | -a = 1.0; x0 = 1.0 |
91 | | -figure() |
92 | | -t = linspace(0.0, 3.0, 1000) |
93 | | -plot(t, exp(a*t)*x0, "k") |
94 | | -xlabel("$t$"); ylabel("$x(t)$"); title(f"$a={a}$") |
95 | | -grid(); axis([0.0, 2.0, 0.0, 10.0]) |
96 | | -tight_layout() |
97 | | -save("images/scalar-LTI-1") |
98 | | -figure() |
99 | | -plot(real(a), imag(a), "x", color="k") |
100 | | -gca().set_aspect(1.0) |
101 | | -xlim(-3,3); ylim(-3,3); |
102 | | -plot([-3,3], [0,0], "k") |
103 | | -plot([0, 0], [-3, 3], "k") |
104 | | -xticks([-2,-1,0,1,2]); yticks([-2,-1,0,1,2]) |
105 | | -title(f"$a={a}$"); grid(True) |
106 | | -tight_layout() |
107 | | -save("images/scalar-LTI-1-poles") |
108 | | -a = 0.0; x0 = 1.0 |
109 | | -figure() |
110 | | -t = linspace(0.0, 3.0, 1000) |
111 | | -plot(t, exp(a*t)*x0, "k") |
112 | | -xlabel("$t$"); ylabel("$x(t)$"); title(f"$a={a}$") |
113 | | -grid(); axis([0.0, 2.0, 0.0, 10.0]) |
114 | | -tight_layout() |
115 | | -save("images/scalar-LTI-0") |
116 | | -figure() |
117 | | -plot(real(a), imag(a), "x", color="k") |
118 | | -gca().set_aspect(1.0) |
119 | | -xlim(-3,3); ylim(-3,3); |
120 | | -plot([-3,3], [0,0], "k") |
121 | | -plot([0, 0], [-3, 3], "k") |
122 | | -xticks([-2,-1,0,1,2]); yticks([-2,-1,0,1,2]) |
123 | | -title(f"$a={a}$"); grid(True) |
124 | | -tight_layout() |
125 | | -save("images/scalar-LTI-0-poles") |
126 | | -a = -1.0; x0 = 1.0 |
127 | | -figure() |
128 | | -t = linspace(0.0, 3.0, 1000) |
129 | | -plot(t, exp(a*t)*x0, "k") |
130 | | -xlabel("$t$"); ylabel("$x(t)$"); title(f"$a={a}$") |
131 | | -grid(); axis([0.0, 2.0, 0.0, 10.0]) |
132 | | -tight_layout() |
133 | | -save("images/scalar-LTI-m1") |
134 | | -figure() |
135 | | -plot(real(a), imag(a), "x", color="k") |
136 | | -gca().set_aspect(1.0) |
137 | | -xlim(-3,3); ylim(-3,3); |
138 | | -plot([-3,3], [0,0], "k") |
139 | | -plot([0, 0], [-3, 3], "k") |
140 | | -xticks([-2,-1,0,1,2]); yticks([-2,-1,0,1,2]) |
141 | | -title(f"$a={a}$"); grid(True) |
142 | | -tight_layout() |
143 | | -save("images/scalar-LTI-m1-poles") |
144 | | -a = -2.0; x0 = 1.0 |
145 | | -figure() |
146 | | -t = linspace(0.0, 3.0, 1000) |
147 | | -plot(t, exp(a*t)*x0, "k") |
148 | | -xlabel("$t$"); ylabel("$x(t)$"); title(f"$a={a}$") |
149 | | -grid(); axis([0.0, 2.0, 0.0, 10.0]) |
150 | | -tight_layout() |
151 | | -save("images/scalar-LTI-m2") |
152 | | -figure() |
153 | | -plot(real(a), imag(a), "x", color="k") |
| 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) |
154 | 96 | gca().set_aspect(1.0) |
155 | | -xlim(-3,3); ylim(-3,3); |
156 | | -plot([-3,3], [0,0], "k") |
157 | | -plot([0, 0], [-3, 3], "k") |
158 | | -xticks([-2,-1,0,1,2]); yticks([-2,-1,0,1,2]) |
159 | | -title(f"$a={a}$"); grid(True) |
| 97 | +quiver(*Q(f, x1, x2), color="k") |
160 | 98 | tight_layout() |
161 | | -save("images/scalar-LTI-m2-poles") |
162 | | -a1 = -1.0; a2 = 2.0; x10 = x20 = 1.0 |
| 99 | +save("images/discont") |
| 100 | +def sigma(x): |
| 101 | + return 1 / (1 + exp(-x)) |
163 | 102 | figure() |
164 | | -t = linspace(0.0, 3.0, 1000) |
165 | | -x1 = exp(a1*t)*x10; x2 = exp(a2*t)*x20 |
166 | | -xn = sqrt(x1**2 + x2**2) |
167 | | -plot(t, xn , "k") |
168 | | -plot(t, x1, "k--") |
169 | | -plot(t, x2 , "k--") |
170 | | -xlabel("$t$"); ylabel("$\|x(t)\|$"); title(f"$a_1={a1}, \; a_2={a2}$") |
171 | | -grid(); axis([0.0, 2.0, 0.0, 10.0]) |
172 | | -tight_layout() |
173 | | -save("images/scalar-LTI-m1p2") |
174 | | -figure() |
175 | | -plot(real(a1), imag(a1), "x", color="k") |
176 | | -plot(real(a2), imag(a2), "x", color="k") |
177 | | -gca().set_aspect(1.0) |
178 | | -xlim(-3,3); ylim(-3,3); |
179 | | -plot([-3,3], [0,0], "k") |
180 | | -plot([0, 0], [-3, 3], "k") |
181 | | -xticks([-2,-1,0,1,2]); yticks([-2,-1,0,1,2]) |
182 | | -title(f"$a_1={a1}, \; a_2={a2}$") |
183 | | -grid(True) |
184 | | -tight_layout() |
185 | | -save("images/scalar-LTI-m1p2-poles") |
186 | | -a1 = -1.0; a2 = -2.0; x10 = x20 = 1.0 |
187 | | -figure() |
188 | | -t = linspace(0.0, 3.0, 1000) |
189 | | -x1 = exp(a1*t)*x10; x2 = exp(a2*t)*x20 |
190 | | -xn = sqrt(x1**2 + x2**2) |
191 | | -plot(t, xn , "k") |
192 | | -plot(t, x1, "k--") |
193 | | -plot(t, x2 , "k--") |
194 | | -xlabel("$t$"); ylabel("$\|x(t)\|$"); title(f"$a_1={a1}, \; a_2={a2}$") |
195 | | -grid(); axis([0.0, 2.0, 0.0, 10.0]) |
196 | | -tight_layout() |
197 | | -save("images/scalar-LTI-m1m2") |
198 | | -figure() |
199 | | -plot(real(a1), imag(a1), "x", color="k") |
200 | | -plot(real(a2), imag(a2), "x", color="k") |
201 | | -gca().set_aspect(1.0) |
202 | | -xlim(-3,3); ylim(-3,3); |
203 | | -plot([-3,3], [0,0], "k") |
204 | | -plot([0, 0], [-3, 3], "k") |
205 | | -xticks([-2,-1,0,1,2]); yticks([-2,-1,0,1,2]) |
206 | | -title(f"$a_1={a1}, \; a_2={a2}$") |
207 | | -grid(True) |
208 | | -tight_layout() |
209 | | -save("images/scalar-LTI-m1m2-poles") |
210 | | -a = 1.0j; x0=1.0 |
211 | | -figure() |
212 | | -t = linspace(0.0, 20.0, 1000) |
213 | | -plot(t, real(exp(a*t)*x0), label="$\Re(x(t))$") |
214 | | -plot(t, imag(exp(a*t)*x0), label="$\mathrm{Im}(x(t))$") |
215 | | -xlabel("$t$") |
216 | | -legend(); grid() |
217 | | -tight_layout() |
218 | | -save("images/scalar-LTI-alt-1") |
219 | | -fig = figure() |
220 | | -ax = fig.add_subplot(111, projection="3d") |
221 | | -zticks = ax.set_zticks |
222 | | -ax.plot(t, real(exp(a*t)*x0), imag(exp(a*t)*x0)) |
223 | | -xticks([0.0, 20.0]); yticks([]); zticks([]) |
224 | | -ax.set_xlabel("$t$") |
225 | | -ax.set_ylabel("$\Re(x(t))$") |
226 | | -ax.set_zlabel("$\mathrm{Im}(x(t))$") |
227 | | -tight_layout() |
228 | | -save("images/scalar-LTI-3d") |
229 | | -figure() |
230 | | -plot(real(a), imag(a), "x", color="k") |
231 | | -gca().set_aspect(1.0) |
232 | | -xlim(-3,3); ylim(-3,3); |
233 | | -plot([-3,3], [0,0], "k") |
234 | | -plot([0, 0], [-3, 3], "k") |
235 | | -xticks([-2,-1,0,1,2]); yticks([-2,-1,0,1,2]) |
236 | | -title(f"$a={a}$"); grid(True) |
237 | | -tight_layout() |
238 | | -save("images/scalar-LTI-1j-poles") |
239 | | -a = -0.5 + 1.0j; x0=1.0 |
240 | | -figure() |
241 | | -t = linspace(0.0, 20.0, 1000) |
242 | | -plot(t, real(exp(a*t)*x0), label="$\Re(x(t))$") |
243 | | -plot(t, imag(exp(a*t)*x0), label="$\mathrm{Im}(x(t))$") |
244 | | -xlabel("$t$") |
245 | | -legend(); grid() |
246 | | -tight_layout() |
247 | | -save("images/scalar-LTI-alt-2") |
248 | | -fig = figure() |
249 | | -ax = fig.add_subplot(111, projection="3d") |
250 | | -zticks = ax.set_zticks |
251 | | -ax.plot(t, real(exp(a*t)*x0), imag(exp(a*t)*x0)) |
252 | | -xticks([0.0, 20.0]); yticks([]); zticks([]) |
253 | | -ax.set_xlabel("$t$") |
254 | | -ax.set_ylabel("$\Re(x(t))$") |
255 | | -ax.set_zlabel("$\mathrm{Im}(x(t))$") |
256 | | -tight_layout() |
257 | | -save("images/scalar-LTI-3d-2") |
258 | | -figure() |
259 | | -plot(real(a), imag(a), "x", color="k") |
260 | | -gca().set_aspect(1.0) |
261 | | -xlim(-3,3); ylim(-3,3); |
262 | | -plot([-3,3], [0,0], "k") |
263 | | -plot([0, 0], [-3, 3], "k") |
264 | | -xticks([-2,-1,0,1,2]); yticks([-2,-1,0,1,2]) |
265 | | -title(f"$a={a}$") |
| 103 | +x = linspace(-7.0, 7.0, 1000) |
| 104 | +plot(x, sigma(x), label="$y=\sigma(x)$") |
266 | 105 | grid(True) |
267 | | -tight_layout() |
268 | | -save("images/scalar-LTI-m11j-poles") |
| 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