Deploying to gh-pages from @ c6c2f35 🚀

Author: boisgera

.jupyterlite.doit.db

@@ -1,9 +1,7 @@
 from numpy import *
 from numpy.linalg import *
 from scipy.integrate import solve_ivp
-from scipy.linalg import solve_continuous_are
 from matplotlib.pyplot import *
-from numpy.testing import *
 # Python 3.x Standard Library
 import gc
 import os
@@ -63,82 +61,110 @@ def set_ratio(ratio=1.0, bottom=0.1, top=0.1, left=0.1, right=0.1):
     height_in = (1.0 - left - right)/(1.0 - bottom - top) * width_in / ratio
     pp.gcf().set_size_inches((width_in, height_in))
     pp.gcf().subplots_adjust(bottom=bottom, top=1.0-top, left=left, right=1.0-right)
+
+width
 def Q(f, xs, ys):
     X, Y = meshgrid(xs, ys)
     fx = vectorize(lambda x, y: f([x, y])[0])
     fy = vectorize(lambda x, y: f([x, y])[1])
     return X, Y, fx(X, Y), fy(X, Y)
-from scipy.signal import place_poles
-A = array([[0, 1], [0, 0]])
-C = array([[1, 0]])
-poles = [-1, -2]
-K = place_poles(A.T, C.T, poles).gain_matrix
-L = K.T
-assert_almost_equal(K, [[3.0, 2.0]])
-def fun(t, X_Xhat):
-    x, x_hat = X_Xhat[0:2], X_Xhat[2:4]
-    y, y_hat = C.dot(x), C.dot(x_hat)
-    dx = A.dot(x)
-    dx_hat = A.dot(x_hat) - L.dot(y_hat - y)
-    return r_[dx, dx_hat]
-y0 = [-2.0, 1.0, 0.0, 0.0]
-result = solve_ivp(
-    fun=fun, 
-    t_span=[0.0, 5.0], 
-    y0=y0, 
-    max_step=0.1
-)
+def fun(t, y):
+    return y * y
+t0, tf, y0 = 0.0, 3.0, array([1.0])
+result = solve_ivp(fun, t_span=[t0, tf], y0=y0)
+figure()
+plot(result["t"], result["y"][0], "k")
+xlim(t0, tf); xlabel("$t$"); ylabel("$x(t)$")
+tight_layout()
+save("images/finite-time-blowup")
+tf = 1.0
+r = solve_ivp(fun, [t0, tf], y0,
+              dense_output=True)
+figure()
+t = linspace(t0, tf, 1000)
+plot(t, r["sol"](t)[0], "k")
+ylim(0.0, 10.0); grid();
+xlabel("$t$"); ylabel("$x(t)$")
+tight_layout()
+save("images/finite-time-blowup-2")
+def f(x1x2):
+    x1, x2 = x1x2
+    dx1 = 1.0 if x1 < 0.0 else -1.0
+    return array([dx1, 0.0])
 figure()
-t = result["t"]
-y = result["y"]
-plot(t, y[0], "C0", label="$x_1$")
-plot(t, y[2], "C0--", label=r"$\hat{x}_1$")
-plot(t, y[1], "C1", label="$x_2$")
-plot(t, y[3], "C1--", label=r"$\hat{x}_2$")
-xlabel("$t$"); grid(); legend()
+x1 = x2 = linspace(-1.0, 1.0, 20)
+gca().set_aspect(1.0)
+quiver(*Q(f, x1, x2), color="k")
+tight_layout()
+save("images/discont")
+def sigma(x):
+  return 1 / (1 + exp(-x))
+figure()
+x = linspace(-7.0, 7.0, 1000)
+plot(x, sigma(x), label="$y=\sigma(x)$")
+grid(True)
+xlim(-5, 5)
+xticks([-5.0, 0.0, 5.0])
+yticks([0.0, 0.5, 1.0])
+xlabel("$x$")
+ylabel("$y$")
+legend()
 pp.gcf().subplots_adjust(bottom=0.2)
-save("images/observer-trajectories")
-A = array([[0, 1], [0, 0]])
-B = array([[0], [1]])
-Q = array([[1, 0], [0, 1]])
-R = array([[1]])
-sca = solve_continuous_are
-Sigma = sca(A.T, C.T, inv(Q), inv(R))
-L = Sigma @ C.T @ R
+save("images/sigmoid")
+alpha = 2 / 3; beta = 4 / 3; delta = gamma = 1.0
 
-eigenvalues, _ = eig(A - L @ C)
-assert all([real(s) < 0 for s in eigenvalues])
-figure()
-x = [real(s) for s in eigenvalues]
-y = [imag(s) for s in eigenvalues]
-plot(x, y, "kx"); xlim(-2, 2); ylim(-2, 2)
-plot([0, 0], [-2, 2], "k"); 
-plot([-2, 2], [0, 0], "k")   
-grid(True); axis("square")
-title("Eigenvalues")
-axis("square")
-axis([-2, 2, -2, 2])
-save("images/poles-Kalman")
-def fun(t, X_Xhat):
-    x, x_hat = X_Xhat[0:2], X_Xhat[2:4]
-    y, y_hat = C.dot(x), C.dot(x_hat)
-    dx = A.dot(x)
-    dx_hat = A.dot(x_hat) - L.dot(y_hat - y)
-    return r_[dx, dx_hat]
-y0 = [-2.0, 1.0, 0.0, 0.0]
+def fun(t, y):
+    x, y = y
+    u = alpha * x - beta * x * y
+    v = delta * x * y - gamma * y
+    return array([u, v])
+tf = 3.0
 result = solve_ivp(
-    fun=fun, 
-    t_span=[0.0, 5.0], 
-    y0=y0, 
-    max_step=0.1
-)
+  fun, 
+  t_span=(0.0, tf), 
+  y0=[1.5, 1.5], 
+  max_step=0.01)
+x, y = result["y"][0], result["y"][1]
+def display_streamplot():
+    ax = gca()
+    xr = yr = linspace(0.0, 2.0, 1000)
+    def f(y):
+        return fun(0, y)
+    streamplot(*Q(f, xr, yr), color="grey")
+def display_reference_solution():
+    for xy in zip(x, y):
+        x_, y_ = xy
+        gca().add_artist(Circle((x_, y_), 
+                         0.2, color="#d3d3d3"))
+    gca().add_artist(Circle((x[0], y[0]), 0.1, 
+                     color="#808080"))
+    plot(x, y, "k")
+def display_alternate_solution():
+    result = solve_ivp(fun, 
+                       t_span=[0.0, tf],
+                       y0=[1.5, 1.575], 
+                       max_step=0.01)
+    x, y = result["y"][0], result["y"][1]
+    plot(x, y, "k--")
+figure()
+display_streamplot()
+display_reference_solution()
+display_alternate_solution()
+axis([0,2,0,2]); axis("square")
+save("images/continuity")
+def fun(t, y):
+  x = y[0]
+  dx = sqrt(abs(y))
+  return [dx]
+tspan = [0.0, 3.0]
+t = linspace(tspan[0], tspan[1], 1000)
 figure()
-t = result["t"]; y = result["y"]
-plot(t, y[0], "C0", label="$x_1$")
-plot(t, y[2], "C0--", label=r"$\hat{x}_1$")
-plot(t, y[1], "C1", label="$x_2$")
-plot(t, y[3], "C1--", label=r"$\hat{x}_2$")
-xlabel("$t$")
-grid(); legend()
+for x0 in [0.1, 0.01, 0.001, 0.0001, 0.0]:
+    r = solve_ivp(fun, tspan, [x0], 
+        dense_output=True)
+    plot(t, r["sol"](t)[0], 
+         label=f"$x_0 = {x0}$")
+xlabel("$t$"); ylabel("$x(t)$")
+legend()
 pp.gcf().subplots_adjust(bottom=0.2)
-save("images/observer-Kalman-trajectories")
+save("images/eps")

.tmp.py

@@ -1260,7 +1260,7 @@ data = mivp.solve_alt(
 good = to_rgb("#51cf66")
 bad = to_rgb("#ff6b6b")
 mivp.generate_movie(data, filename="videos/gas.mp4", fps=df,
-    axes=gca(), zorder=1000, color=good, linewidth=3.0,
+    axes=gca(), zorder=1000, color=good, linewidth=1.0,
 )
 ```
 
@@ -1286,8 +1286,6 @@ import matplotlib.pyplot as plt
 # Local Library
 import mivp
 
-# ------------------------------------------------------------------------------
-
 # Vector field
 def fun(t, xy):
     x, y = xy
@@ -1347,7 +1345,7 @@ data = mivp.solve_alt(
 
 bad = to_rgb("#ff6b6b")
 mivp.generate_movie(data, filename="videos/movie.mp4", fps=df,
-    axes=gca(), zorder=1000, color=bad, linewidth=3.0,
+    axes=gca(), zorder=1000, color=bad, linewidth=1.0,
 )
 ```
 

__pycache__/mivp.cpython-39.pyc

@@ -2,11 +2,11 @@
   "content": [
     {
       "content": null,
-      "created": "2025-02-07T18:42:32.403648Z",
+      "created": "2025-02-08T16:51:28.644619Z",
       "format": null,
       "hash": null,
       "hash_algorithm": null,
-      "last_modified": "2025-02-07T18:26:32.576940Z",
+      "last_modified": "2025-02-08T16:33:58.925278Z",
       "mimetype": null,
       "name": "build",
       "path": "JOSE/build",
@@ -16,11 +16,11 @@
     },
     {
       "content": null,
-      "created": "2025-02-07T18:42:32.403648Z",
+      "created": "2025-02-08T16:51:28.645618Z",
       "format": null,
       "hash": null,
       "hash_algorithm": null,
-      "last_modified": "2025-02-07T18:26:32.576940Z",
+      "last_modified": "2025-02-08T16:33:58.925278Z",
       "mimetype": "text/x-bibtex",
       "name": "paper.bib",
       "path": "JOSE/paper.bib",
@@ -30,11 +30,11 @@
     },
     {
       "content": null,
-      "created": "2025-02-07T18:42:32.404648Z",
+      "created": "2025-02-08T16:51:28.645618Z",
       "format": null,
       "hash": null,
       "hash_algorithm": null,
-      "last_modified": "2025-02-07T18:26:32.576940Z",
+      "last_modified": "2025-02-08T16:33:58.925278Z",
       "mimetype": "text/markdown",
       "name": "paper.md",
       "path": "JOSE/paper.md",
@@ -43,11 +43,11 @@
       "writable": true
     }
   ],
-  "created": "2025-02-07T18:42:32.403648Z",
+  "created": "2025-02-08T16:51:28.645618Z",
   "format": "json",
   "hash": null,
   "hash_algorithm": null,
-  "last_modified": "2025-02-07T18:42:32.403648Z",
+  "last_modified": "2025-02-08T16:51:28.645618Z",
   "mimetype": null,
   "name": "JOSE",
   "path": "JOSE",