% Asymptotic Behavior % 👤 Sébastien Boisgérault

Control Engineering with Python

📖 Documents (GitHub)
©️ License CC BY 4.0
🏦 Mines ParisTech, PSL University

Symbols


🐍	Code	🔍	Worked Example
📈	Graph	🧩	Exercise
🏷️	Definition	💻	Numerical Method
💎	Theorem	🧮	Analytical Method
📝	Remark	🧠	Theory
ℹ️	Information	🗝️	Hint
⚠️	Warning	🔓	Solution

🐍 Imports

from numpy import *
from numpy.linalg import *
from scipy.linalg import *
from matplotlib.pyplot import *
from mpl_toolkits.mplot3d import *
from scipy.integrate import solve_ivp

::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

# Python 3.x Standard Library
import gc
import os

# Third-Party Packages
import numpy as np; np.seterr(all="ignore")
import numpy.linalg as la
import scipy.misc
import matplotlib as mpl; mpl.use("Agg")
import matplotlib.pyplot as pp
import matplotlib.axes as ax
import matplotlib.patches as pa


#
# Matplotlib Configuration & Helper Functions
# --------------------------------------------------------------------------

# TODO: also reconsider line width and markersize stuff "for the web
#       settings".
fontsize = 10

width = 345 / 72.27
height = width / (16/9)

rc = {
    "text.usetex": True,
    "pgf.preamble": r"\usepackage{amsmath,amsfonts,amssymb}",
    #"font.family": "serif",
    "font.serif": [],
    #"font.sans-serif": [],
    "legend.fontsize": fontsize,
    "axes.titlesize":  fontsize,
    "axes.labelsize":  fontsize,
    "xtick.labelsize": fontsize,
    "ytick.labelsize": fontsize,
    "figure.max_open_warning": 100,
    #"savefig.dpi": 300,
    #"figure.dpi": 300,
    "figure.figsize": [width, height],
    "lines.linewidth": 1.0,
}
mpl.rcParams.update(rc)

# Web target: 160 / 9 inches (that's ~45 cm, this is huge) at 90 dpi
# (the "standard" dpi for Web computations) gives 1600 px.
width_in = 160 / 9

def save(name, **options):
    cwd = os.getcwd()
    root = os.path.dirname(os.path.realpath(__file__))
    os.chdir(root)
    pp.savefig(name + ".svg", **options)
    os.chdir(cwd)

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)

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

🐍 Streamplot Helper

def Q(f, xs, ys):
    X, Y = meshgrid(xs, ys)
    v = vectorize
    fx = v(lambda x, y: f([x, y])[0])
    fy = v(lambda x, y: f([x, y])[1])
    return X, Y, fx(X, Y), fy(X, Y)

ℹ️ Assumption

From now on, we only deal with well-posed systems.

🏷️ Asymptotic

Asymptotic = Long-Term: when $t \to + \infty$

::: notes ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

So far, we have given little attention to the properties that characterize the long-term behavior of dynamical systems. Well-posedness ensures merely the local (in time) existence of uniqueness of solutions to an IVP and continuity with respect to the initial condition is only applicable to solutions defined on a finite (compact) time intervals.

We will now turn to the study of such properties specifically.

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

⚠️

Even simple dynamical systems may exhibit

complex asymptotic behaviors.

Lorenz System

$$ \begin{array}{lll} \dot{x} &=& \sigma (y - x) \\ \dot{y} &=& x (\rho - z) \\ \dot{z} &=& xy - \beta z \end{array} $$

::: notes ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

But first, we need to realize that dynamical systems, even when they are governed by a simple, low-dimensional system of equations, may exhibit very complex asymptotic patterns.

The Lorenz system is a classical example of such system.

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

Visualized with Fibre

::: notes ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

The solutions $x(t)$ of this system, which are global, have no limit when $t \to + \infty$ and don't blow up either but instead oscillates forever between two regions of the state space. This long-term behavior is quantitatively very sensitive to the choice of the initial condition, making long-term predictions about the system practically impossible: the system is chaotic.

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

Hadley System

$$ \begin{array}{lll} \dot{x} &=& -y^2 - z^2 - ax + af\\ \dot{y} &=& xy - b xz - y + g \\ \dot{z} &=& bxy + xz - z \end{array} $$

Visualized with Fibre

::: notes ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

Asymptotic patterns of chaotic system are by no means restricted to this "switching" behavior. The Hadley system is another -- even "messier" -- example of chaotic behavior.

Fortunately, many systems behaved more predictabily, and/or can be controlled so that their asymptotic behavior is more acceptable.

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

🏷️ Equilibrium

An equilibrium of system $\dot{x} = f(x)$ is a state $x_e$ such that the maximal solution $x(t)$ such that $x(0) = x_e$

is global and,
is $x(t) = x_e$ for any $t > 0$.

💎 Equilibrium

The state $x_e$ is an equilibrium of $\dot{x} = f(x)$

$$\Longleftrightarrow$$

$$f(x_e) = 0.$$

::: notes ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

If $f(x_e) = 0$, then $d x_e / dt = f(x_e) = 0$ and thus the (constant) function $t \in \left[0, +\infty \right) \mapsto x_e$ is a global solution of the IVP $\dot{x} =f(x)$ and $x_e = 0$.

Conversely, if $x(t):= x_e$, $t\geq 0$ is a global solution of $\dot{x} = f(x)$, then $0 = \dot{x}(0) = f(x(0)) = f(x_e)$.

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

Stability

About the long-term behavior of solutions.

"Stability" subtle concept,
"Asymptotic Stability" simpler (and stronger),
"Attractivity" simpler yet, (but often too weak).

Attractivity

Context: system $\dot{x} = f(x)$ with equilibrium $x_e$.

🏷️ Global Attractivity

The equilibrium $x_e$ is globally attractive if for every $x_0,$ the maximal solution $x(t)$ such that $x(0)=x_0$

is global and,
$x(t) \to x_e$ when $t \to +\infty$.

🏷️ Local Attractivity

The equilibrium $x_e$ is locally attractive if for every $x_0$ close enough to $x_e$, the maximal solution $x(t)$ such that $x(0)=x_0$

is global and,
$x(t) \to x_e$ when $t \to +\infty$.

::: notes ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

Technically, an equilibrium $x_e$ is locally attractive if there is a $r>0$ such that the maximal solution to the IVP $\dot{x} = f(x)$ and $x(0) = x_0$ is global and satisfies $x(t) \to +\infty$ when $t \to +\infty$ whenever $|x_0 - x_e| \leq r$.

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

🔍 Global Attractivity

The system

$$ \begin{array}{cc} \dot{x} &=& -2x + y \\ \dot{y} &=& -2y + x \end{array} $$

is well-posed,
has an equilibrium at $(0, 0)$.

🐍 Vector field

def f(xy):
    x, y = xy
    dx = -2*x + y
    dy = -2*y + x
    return array([dx, dy])

📈 Stream plot

figure()
x = y = linspace(-5.0, 5.0, 1000)
streamplot(*Q(f, x, y), color="k")
plot([0], [0], "k.", ms=20.0)
axis("square")
axis("off")

::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

tight_layout()
save("images/globally-attractive")

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

{.section data-background="images/globally-attractive.svg" data-background-size="contain"}

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

import matplotlib.animation as ani
from matplotlib.colors import to_rgb
from tqdm import tqdm


neutral = grey_4 = to_rgb("#ced4da")
#grey_5 = to_rgb("#adb5bd")
#grey_8 = to_rgb("#343a40")
good = to_rgb("#51cf66")
bad = to_rgb("#ff6b6b")

ft = lambda t, y: f(y)
fps = df = 60.0
dt = 1.0 / df
t_span = t_i, t_f = (0.0, 3.0)
t = np.arange(t_i, t_f + dt, dt)

y0s = [[-4.0, 1.0], [4.0, -1.0]]
colors = [good, good]
xys = []
for y0 in tqdm(y0s):
    r = solve_ivp(fun=ft, y0=y0, t_span=t_span, t_eval=t)
    xys.append(r.y)


fig = figure()
x = y = linspace(-5.0, 5.0, 1000)
streamplot(*Q(f, x, y), color=grey_4)
plot([0], [0], lw=3.0, marker="o", ms=10.0, markevery=[-1],
        markeredgecolor="white", color=neutral)
axis("square")
axis("off")

lines = []
for x, y in xys:
    line = plot(
        [x[0]], [y[0]],
        lw=3.0,
        ms=10.0,
        color=neutral,
        marker="o", markevery=[-1],
        markeredgecolor="white")[0]
    lines.append(line)
tight_layout()

num_frames = len(t) * len(lines)

def gamma(x):
    return pow(x, 0.5)

def update(i):
    j, k = divmod(i, len(t))
    x, y = xys[j]
    line = lines[j]
    line.set_data(x[:k+1], y[:k+1])
    alpha = gamma(k / (len(t)-1))
    final_color = colors[j]
    line.set_color(
      tuple((1-alpha)*array(neutral) + alpha*array(final_color))
    )

animation = ani.FuncAnimation(fig, func=update, frames=num_frames)
writer = ani.FFMpegWriter(fps=fps)
bar = tqdm(total=num_frames)
animation.save("videos/globally-attractive.mp4", writer=writer, dpi=300,
progress_callback = lambda i, n: bar.update(1))
bar.close()

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

<video controls style="width:100vw;">
  <source src="videos/globally-attractive.mp4" type="video/mp4">
  Your browser does not support the video tag.
</video>

🔍 Local Attractivity

The system

$$ \begin{array}{cc} \dot{x} &=& -2x + y^3 \\ \dot{y} &=& -2y + x^3 \end{array} $$

is well-posed,
has an equilibrium at $(0, 0)$.

🐍 Vector field

def f(xy):
    x, y = xy
    dx = -2*x + y**3
    dy = -2*y + x**3
    return array([dx, dy])

📈 Stream plot

figure()
x = y = linspace(-5.0, 5.0, 1000)
streamplot(*Q(f, x, y), color="k")
plot([0], [0], "k.", ms=10.0)
axis("square")
axis("off")

::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

tight_layout()
save("images/locally-attractive")

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

{.section data-background="images/locally-attractive.svg" data-background-size="contain"}

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

import matplotlib.animation as ani
from matplotlib.colors import to_rgb
from tqdm import tqdm


neutral = grey_4 = to_rgb("#ced4da")
#grey_5 = to_rgb("#adb5bd")
#grey_8 = to_rgb("#343a40")
good = to_rgb("#51cf66")
bad = to_rgb("#ff6b6b")

ft = lambda t, y: f(y)
fps = df = 60.0
dt = 1.0 / df
t_span = t_i, t_f = (0.0, 3.0)
t = np.arange(t_i, t_f + dt, dt)

y0s = [[-2.4, 0.0], [-1.8, 0.0], [-1.2, 0.0]]
colors = [bad, good, good]
xys = []
for y0 in tqdm(y0s):
    r = solve_ivp(fun=ft, y0=y0, t_span=t_span, t_eval=t)
    xys.append(r.y)


fig = figure()
x = y = linspace(-5.0, 5.0, 1000)
streamplot(*Q(f, x, y), color=grey_4)
plot([0], [0], lw=3.0, marker="o", ms=10.0, markevery=[-1],
        markeredgecolor="white", color=neutral)
axis("square")
axis("off")

lines = []
for x, y in xys:
    line = plot(
        [x[0]], [y[0]],
        lw=3.0,
        ms=10.0,
        color=neutral,
        marker="o", markevery=[-1],
        markeredgecolor="white")[0]
    lines.append(line)
tight_layout()

num_frames = len(t) * len(lines)

def gamma(x):
    return pow(x, 0.5)

def update(i):
    j, k = divmod(i, len(t))
    x, y = xys[j]
    line = lines[j]
    line.set_data(x[:k+1], y[:k+1])
    alpha = gamma(k / (len(t)-1))
    final_color = colors[j]
    line.set_color(
      tuple((1-alpha)*array(neutral) + alpha*array(final_color))
    )

animation = ani.FuncAnimation(fig, func=update, frames=num_frames)
writer = ani.FFMpegWriter(fps=fps)
bar = tqdm(total=num_frames)
animation.save("videos/locally-attractive.mp4", writer=writer, dpi=300,
progress_callback = lambda i, n: bar.update(1))
bar.close()

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

<video controls style="width:100vw;">
  <source src="videos/locally-attractive.mp4" type="video/mp4">
  Your browser does not support the video tag.
</video>

🔍 No Attractivity

The system

$$ \begin{array}{cr} \dot{x} &=& -2x + y \\ \dot{y} &=& 2y - x \end{array} $$

is well-posed,
has a (unique) equilibrium at $(0, 0)$.

🐍 Vector field

def f(xy):
    x, y = xy
    dx = -2*x + y
    dy =  2*y - x
    return array([dx, dy])

📈 Stream plot

figure()
x = y = linspace(-5.0, 5.0, 1000)
streamplot(*Q(f, x, y), color="k")
plot([0], [0], "k.", ms=10.0)
axis("square")
axis("off")

::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

tight_layout()
save("images/not-attractive")

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

{.section data-background="images/not-attractive.svg" data-background-size="contain"}

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

import matplotlib.animation as ani
from matplotlib.colors import to_rgb
from tqdm import tqdm


neutral = grey_4 = to_rgb("#ced4da")
#grey_5 = to_rgb("#adb5bd")
#grey_8 = to_rgb("#343a40")
good = to_rgb("#51cf66")
bad = to_rgb("#ff6b6b")

ft = lambda t, y: f(y)
fps = df = 60.0
dt = 1.0 / df
t_span = t_i, t_f = (0.0, 3.0)
t = np.arange(t_i, t_f + dt, dt)

y0s = [[-4*0.9659258262890683, -4*0.2588190451025208],
       [-4*0.9659258262890683, -4*0.2588190451025208 + 0.2],
       [0.2*0.9659258262890683, 0.2*0.2588190451025208 - 0.2]]
colors = [good, bad, bad]
xys = []
for y0 in tqdm(y0s):
    r = solve_ivp(fun=ft, y0=y0, t_span=t_span, t_eval=t)
    xys.append(r.y)


fig = figure()
x = y = linspace(-5.0, 5.0, 1000)
streamplot(*Q(f, x, y), color=grey_4)
plot([0], [0], lw=3.0, marker="o", ms=10.0, markevery=[-1],
        markeredgecolor="white", color=neutral)
axis("square")
axis("off")

lines = []
for x, y in xys:
    line = plot(
        [x[0]], [y[0]],
        lw=3.0,
        ms=10.0,
        color=neutral,
        marker="o", markevery=[-1],
        markeredgecolor="white")[0]
    lines.append(line)
tight_layout()

num_frames = len(t) * len(lines)

def gamma(x):
    return pow(x, 0.5)

def update(i):
    j, k = divmod(i, len(t))
    x, y = xys[j]
    line = lines[j]
    line.set_data(x[:k+1], y[:k+1])
    alpha = gamma(k / (len(t)-1))
    final_color = colors[j]
    line.set_color(
      tuple((1-alpha)*array(neutral) + alpha*array(final_color))
    )

animation = ani.FuncAnimation(fig, func=update, frames=num_frames)
writer = ani.FFMpegWriter(fps=fps)
bar = tqdm(total=num_frames)
animation.save("videos/not-attractive.mp4", writer=writer, dpi=300,
progress_callback = lambda i, n: bar.update(1))
bar.close()

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

<video controls style="width:100vw;">
  <source src="videos/not-attractive.mp4" type="video/mp4">
  Your browser does not support the video tag.
</video>

🧩 Pendulum

The pendulum is governed by the equation

$$ m \ell^2 \ddot{\theta} + b \dot{\theta} + m g \ell \sin \theta = 0 $$

where $m>0$, $\ell>0$, $g>0$ and $b\geq0$.

1. 🧮

Compute the equilibria of this system.

2. 🧠

Can any of these equilibria be globally attractive?

3. 📈

Assume that $m=1$, $\ell=1$, $g=9.81$ and $b=1$.

Make a stream plot of the system.

4. 🔬

Determine which equilibria are locally attractive.

5. 📈 Frictionless Pendulum

Assume now that $b=0$.

Make a stream plot of the system.

6. 🧮 🧠

Prove that the equilibrium at $(0,0)$ is not locally attractive.

🗝️ Hint. Study how the total mechanical energy $E$

$$ E(\theta,\dot{\theta}) := m\ell^2 \dot{\theta}^2 / 2 - m g\ell \cos \theta $$

evolves in time.

🔓 Pendulum

1. 🔓

The 2nd-order differential equations of the pendulum are equivalent to the first order system

$$ \left| \begin{array}{rcl} \dot{\theta} &=& \omega \\ \dot{\omega} &=& (- b / m \ell^2) \omega - (g /\ell) \sin \theta \\ \end{array} \right. $$

Thus, the system state is $x =(\theta, \omega)$ and is governed by $\dot{x} = f(x)$ with

$$ f(\theta, \omega) = (\omega, (- b / m \ell^2) \omega - (g /\ell) \sin \theta). $$

Hence, the state $(\theta, \omega)$ is a solution to $f(\theta, \omega) = 0$ if and only if $\omega = 0$ and $\sin \theta = 0$. In other words, the equilibria of the system are characterized by $\theta = k \pi$ for some $k \in \mathbb{Z}$ and $\omega (= \dot{\theta}) = 0$.

2. 🔓

Since there are several equilibria, none of them can be globally attractive.

Indeed let $x_1$ be a globally attractive equilibrium and assume that $x_2$ is any other equilibrium. By definition, the maximal solution $x(t)$ such that $x(0) = x_2$ is $x(t) = x_2$ for every $t\geq0$. On the other hand, since $x_1$ is globally attractive, it also satisfies $x(t) \to x_1$ when $t\to +\infty$, hence there is a contradiction.

Thus, $x_1$ is the only possible equilibrium.

3. 🔓

m = l = b = 1; g=9.81

def f(theta_omega):
    theta, omega = theta_omega
    d_theta = omega
    d_omega = - b / (m * l * l) * omega
    d_omega -=  (g / l) * sin(theta)
    return (d_theta, d_omega)

figure()
theta = linspace(-2*pi*(1.2), 2*pi*(1.2), 1000)
d_theta = linspace(-5.0, 5.0, 1000)
streamplot(*Q(f, theta, d_theta), color="k")
plot([-2*pi, -pi, 0 ,pi, 2*pi], 5*[0.0], "k.")
xticks([-2*pi, -pi, 0 ,pi, 2*pi],
[r"$-2\pi$", r"$\pi$", r"$0$", r"$\pi$", r"$2\pi$"])
grid(True)

::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

tight_layout()
save("images/pendulum-friction")

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

{.section data-background="images/pendulum-friction.svg" data-background-size="contain"}

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

4. 🔓

From the streamplot, we see that the equilibria

$$(\theta, \dot{\theta}) = (2k\pi, 0), ; k\in \mathbb{Z}$$

are asymptotically stable, but that the equilibria

$$(\theta, \dot{\theta}) = (2(k+1) \pi, 0),; k\in \mathbb{Z}$$

are not (they are not locally attractive).

5. 🔓

b = 0
figure()
streamplot(*Q(f, theta, d_theta), color="k")
plot([-2*pi, -pi, 0 ,pi, 2*pi], 5*[0.0], "k.")
xticks([-2*pi, -pi, 0 ,pi, 2*pi],
[r"$-2\pi$", r"$\pi$", r"$0$", r"$\pi$", r"$2\pi$"])
grid(True)

::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

tight_layout()
save("images/pendulum-no-friction")

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

{.section data-background="images/pendulum-no-friction.svg" data-background-size="contain"}

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

6. 🔓

$$ \begin{split} \dot{E} &= \frac{d}{dt} \left(m\ell^2 \dot{\theta}^2 / 2 - m g\ell \cos \theta\right) \\ &= m\ell^2 \ddot{\theta}\dot{\theta} + m g \ell (\sin \theta) \dot{\theta} \\ &= \left(m\ell^2 \ddot{\theta} + m g \ell \sin \theta \right)\dot{\theta} \\ &= \left(- b \dot{\theta}\right)\dot{\theta} \\ &= 0 \end{split} $$

Therefore, $E(t)$ is constant.

On the other hand,

$$ \min , {E(\theta, \dot{\theta}) ; | ; (\theta,\dot{\theta}) \in \mathbb{R}^2} = E(0, 0) = -mgl. $$

Moreover, this minimum is locally strict. Precisely, for any $0< |\theta| < \pi$, $$ E(0, 0) < E(\theta,\dot{\theta}). $$

If the origin was locally attractive, for any $\theta(0)$ and $\dot{\theta}(0)$ small enough, we would have $$ E(\theta(t), \dot{\theta}(t)) \to E(0, 0) ; \mbox{ when } ; t\to +\infty $$ (by continuity). But if $0 < |\theta(0)| < \pi$, we have

$$E(\theta(0), \dot{\theta}(0)) > E(0, 0)$$

and that would contradict that $E(t)$ is constant.

Hence the origin is not locally attractive.

💎 Attractivity (Low-level)

The equilibrium $x_e$ is globally attractive iff:

for any state $x_0$ and for any $\epsilon > 0$ there is a $\tau \geq 0$,
such that the maximal solution $x(t)$ such that $x(0) = x_0$ is global and,
satisfies:

$$ |x(t) - x_e| \leq \epsilon ; \mbox{ when } ; t \geq \tau. $$

⚠️ Warning

Very close values of $x_0$ could theoretically lead to very different "speed of convergence" of $x(t)$ towards the equilibrium.
This is not contradictory with the well-posedness assumption: continuity w.r.t. the initial condition only works with finite time spans.

🔍 (Pathological) Example

$$ \begin{array}{lll} \dot{x} &=& x + xy - (x + y)\sqrt{x^2 + y^2} \\ \dot{y} &=& y - x^2 + (x - y) \sqrt{x^2 + y^2} \end{array} $$

Equivalently, in polar coordinates:

$$ \begin{array}{lll} \dot{r} &=& r (1 - r) \\ \dot{\theta} &=& r (1 - \cos \theta) \end{array} $$

🐍 Vector Field

def f(xy):
    x, y = xy
    r = sqrt(x*x + y*y)
    dx = x + x * y - (x + y) * r
    dy = y - x * x + (x - y) * r
    return array([dx, dy])

📈 Stream Plot

figure()
x = y = linspace(-2.0, 2.0, 1000)
streamplot(*Q(f, x, y), color="k")
plot([1], [0], "k.", ms=20.0)
axis("square")
axis("off")

::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

tight_layout()
save("images/attractive2")

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

{.section data-background="images/attractive2.svg" data-background-size="contain"}

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

import matplotlib.animation as ani
from matplotlib.colors import to_rgb
from tqdm import tqdm

neutral = grey_4 = to_rgb("#ced4da")
#grey_5 = to_rgb("#adb5bd")
#grey_8 = to_rgb("#343a40")
good = to_rgb("#51cf66")
bad = to_rgb("#ff6b6b")

ft = lambda t, y: f(y)
fps = df = 60.0
dt = 1.0 / df
t_span = t_i, t_f = (0.0, 10.0)
t = np.arange(t_i, t_f + 0.1*dt, dt)

y0s = [[1.5*cos(theta), 1.5*sin(theta)] for theta in linspace(0, (11/12)*2*pi, 12)]
# colors = [good] * len(y0s)
xys = []
for y0 in tqdm(y0s):
    r = solve_ivp(fun=ft, y0=y0, t_span=t_span, t_eval=t)
    xys.append(r.y)


fig = figure()
x = y = linspace(-2.0, 2.0, 1000)
streamplot(*Q(f, x, y), color=grey_4)
plot([1], [0], lw=3.0, marker="o", ms=10.0, markevery=[-1],
        markeredgecolor="white", color=neutral)
axis("square")
axis("off")

lines = []
for x, y in xys:
    line = plot(
        [x[0]], [y[0]],
        lw=3.0,
        ms=10.0,
        #color=neutral,
        marker="o", markevery=[-1],
        markeredgecolor="white")[0]
    lines.append(line)
tight_layout()

num_frames = len(t)

def gamma(x):
    return pow(x, 0.5)

def update(i):
  for line, (x, y) in zip(lines, xys):
      line.set_data(x[:i+1], y[:i+1])
    # j, k = divmod(i, len(t))
    # x, y = xys[j]
    # line = lines[j]
    # line.set_data(x[:k+1], y[:k+1])
    # alpha = gamma(k / (len(t)-1))
    # final_color = colors[j]
    # line.set_color(
    #   tuple((1-alpha)*array(neutral) + alpha*array(final_color))
    # )

animation = ani.FuncAnimation(fig, func=update, frames=num_frames)
writer = ani.FFMpegWriter(fps=fps)
bar = tqdm(total=num_frames)
animation.save("videos/pathological.mp4", writer=writer, dpi=300,
progress_callback = lambda i, n: bar.update(1))
bar.close()

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

<video controls style="width:100vw;">
  <source src="videos/pathological.mp4" type="video/mp4">
  Your browser does not support the video tag.
</video>

Asymptotic Stability

Asymptotic stability is a stronger version of attractivity which is by definition robust with respect to the choice of the initial state.

🏷️ Global Asympt. Stability

The equilibrium $x_e$ is globally asympt. stable iff:

for any state $x_0$ and for any $\epsilon > 0$ there is a $\tau \geq 0$,
and there is a $r > 0$ such that if $|x_0' - x_0| \leq r$,
such that the maximal solution $x(t)$ such that $x(0) = x_0'$ is global and,
satisfies:

$$ |x(t) - x_e| \leq \epsilon ; \mbox{ when } ; t \geq \tau. $$

::: notes ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

Global asymptotic stability is requirement which is very similar to attractivity. The only difference is that it requires the time $\tau$ that we should wait to be sure that the distance between $x(t)$ and the equilibrium is less that $\varepsilon$ should be valid not only for the initial condition $x(0) = x_0$ but also for any other initial condition in an (arbitrary small) neighbourhood of $x_0$. There is a common profile for the "speed of convergence" towards the equilibrium between an initial condition and its neighbors ; this condition is not always met if we merely have an attractive equilibrium.

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

Set of Initial Conditions

Let $f: \mathbb{R}^n \to \mathbb{R}^n$ and $X_0 \subset \mathbb{R}^n$.

Let $X(t)$ be the image of $X_0$ by the flow at time $t$:

$$ X(t) := {x(t) ; | ; \dot{x} = f(x), ; x(0)= x_0, ; x_0 \in X_0 \ }. $$

💎 Global Asympt. Stability

An equilibrium $x_e$ is globally asympt. stable iff

for every bounded set $X_0$ and any $x_0 \in X_0$ the associated maximal solution $x(t)$ is global and,
$X(t) \to {x_e}$ when $t\to +\infty$.

🏷️ Limits of Sets

$$ X(t) \to {x_e} $$

to be interpreted as

$$ \sup_{x(t) \in X(t)} |x(t) - x_e| \to 0. $$

🏷️ Hausdorff Distance

$$ \sup_{x(t) \in X(t)} |x(t) - x_e| = d_H(X(t), {x_e}) $$

where $d_H$ is the Hausdorff distance between sets:

$$ d_H(A, B) := \max \left{ \sup_{a \in A} d(a, B), \sup_{b\in B} d(A, b) \right}. $$

::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

# Third-Party Libraries
import numpy as np
import matplotlib.pyplot as plt

# Local Library
import mivp

# ------------------------------------------------------------------------------

# Vector field
def fun(t, xy):
    x, y = xy
    dx = -2*x + y
    dy = -2*y + x
    return array([dx, dy])

# Streamplot
fig = figure()
x = y = linspace(-5.0, 5.0, 1000)
streamplot(*Q(lambda xy: fun(0, xy), x, y), color=grey_4)
plot([0], [0], lw=3.0, marker="o", ms=10.0, markevery=[-1],
        markeredgecolor="white", color=neutral)
axis("square")
axis("off")
tight_layout()

# Time span & frame rate
t_span = (0.0, 10.0)

df = 60.0
dt = 1.0 / df
t = np.arange(t_span[0], t_span[1], dt)
t = np.r_[t, t_span[1]]

# Initial set boundary
y0 = [2.5, 0.0]
radius = 2.0
xc, yc = y0

def boundary(t):  # we assume that t is a 1-dim array
    return np.array(
        [
            [xc + radius * np.cos(theta), yc + radius * np.sin(theta)]
            for theta in 2 * np.pi * t
        ]
    )

# Precision
rtol = 1e-6  # default: 1e-3
atol = 1e-12  # default: 1e-6

# ------------------------------------------------------------------------------

data = mivp.solve_alt(
    fun=fun,
    t_eval=t,
    boundary=boundary,
    boundary_rtol=0.0,
    boundary_atol=0.1,
    rtol=rtol,
    atol=atol,
    method="LSODA",
)

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,
)

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

<video controls style="width:100vw;">
  <source src="videos/gas.mp4" type="video/mp4">
  Your browser does not support the video tag.
</video>

::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

# Third-Party Libraries
import numpy as np
import matplotlib.pyplot as plt

# Local Library
import mivp

# ------------------------------------------------------------------------------

# Vector field
def fun(t, xy):
    x, y = xy
    r = np.sqrt(x * x + y * y)
    dx = x + x * y - (x + y) * r
    dy = y - x * x + (x - y) * r
    return [dx, dy]

# Streamplot
fig = figure()
x = y = linspace(-2.0, 2.0, 1000)
streamplot(*Q(f, x, y), color=grey_4)
plot([1], [0], lw=3.0, marker="o", ms=10.0, markevery=[-1],
        markeredgecolor="white", color=neutral)
axis("square")
axis("off")
tight_layout()

# Time span & frame rate
t_span = (0.0, 10.0)

df = 60.0
dt = 1.0 / df
t = np.arange(t_span[0], t_span[1], dt)
t = np.r_[t, t_span[1]]

# Initial set boundary
y0 = [1.0, 0.0]
radius = 0.5
n = 10000
xc, yc = y0

def boundary(t):  # we assume that t is a 1-dim array
    return np.array(
        [
            [xc + radius * np.cos(theta), yc + radius * np.sin(theta)]
            for theta in 2 * np.pi * t
        ]
    )

# Precision
rtol = 1e-6  # default: 1e-3
atol = 1e-12  # default: 1e-6

# ------------------------------------------------------------------------------

data = mivp.solve_alt(
    fun=fun,
    t_eval=t,
    boundary=boundary,
    boundary_rtol=0.0,
    boundary_atol=0.1,
    rtol=rtol,
    atol=atol,
    method="LSODA",
)

bad = to_rgb("#ff6b6b")
mivp.generate_movie(data, filename="videos/movie.mp4", fps=df,
    axes=gca(), zorder=1000, color=bad, linewidth=3.0,
)

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

<video controls style="width:100vw;">
  <source src="videos/movie.mp4" type="video/mp4">
  Your browser does not support the video tag.
</video>

🏷️ Local Asymptotic Stability

The equilibrium $x_e$ is locally asympt. stable iff:

there is a $r>0$ such that for any $\epsilon > 0$,
there is a $\tau \geq 0$ such that,
if $|x_0 - x_e| \leq r$, the maximal solution $x(t)$ such that $x(0) = x_0$ is global and satisfies:

$$ |x(t) - x_e| \leq \epsilon ; \mbox{ when } ; t \geq \tau. $$

💎 Local Asympt. Stability {#LAS-ppty}

An equilibrium $x_e$ is locally asympt. stable iff:

There is a $r>0$ such that for every set $X_0$ such that

$$ X_0 \subset {x ; | ; |x - x_e| \leq r }, $$

and for any $x_0 \in X_0$, the associated maximal solution $x(t)$ is global and

$$ X(t) \to {x_e} \mbox{ when } t\to +\infty. $$

🏷️ Stability

An equilibrium $x_e$ is stable iff:

for any $r>0$,
there is a $\rho \leq r$ such that if $|x(0) - x_e| \leq \rho$, then
the solution $x(t)$ is global,
for any $t\geq 0$, $|x(t) - x_e| \leq r$.

🧩 Vinograd System

Consider the system:

$$ \begin{array}{rcl} \dot{x} &=& (x^2 (y-x) +y^5) / (x^2 + y^2 (1 + (x^2 + y^2)^2 )) \\ \dot{y} &=& y^2 (y - 2x) / (x^2 + y^2 (1 + (x^2 + y^2)^2 )) \end{array} $$

🐍 Vector field

def f(xy):
    x, y = xy
    q = x**2 + y**2 * (1 + (x**2 + y**2)**2)
    dx = (x**2 * (y - x) + y**5) / q
    dy = y**2 * (y - 2*x) / q
    return array([dx, dy])

📈 Stream plot

figure()
x = y = linspace(-1.0, 1.0, 1000)
streamplot(*Q(f, x, y), color="k")
xticks([-1, 0, 1])
plot([0], [0], "k.", ms=10.0)
axis("square")
axis("off")

::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

tight_layout()
save("images/vinograd")

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

{.section data-background="images/vinograd.svg" data-background-size="contain"}

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

1. 🧮

Show that the origin $(0, 0)$ is the unique equilibrium.

2. 📈🔬

Does this equilibrium seem to be attractive (graphically) ?

3. 🧠

Show that for any equilibrium of a well-posed system:

💎 (locally) asymptotically stable $\Rightarrow$ stable

4. 🧪📈

Does the origin seem to be stable (experimentally ?)

Conclude accordingly.

🔓 Vinograd System

1. 🔓

$(x,y)$ is an equilibrium of the Vinograd system iff

$$ \begin{array}{rcl} (x^2 (y-x) +y^5) / (x^2 + y^2 (1 + (x^2 + y^2)^2 )) &=& 0 \\ y^2 (y - 2x) / (x^2 + y^2 (1 + (x^2 + y^2)^2 )) &=& 0 \end{array} $$

or equivalently

$$ y^2 (y - 2x) = 0 ; \mbox{ and } ; x^2 (y-x) +y^5 = 0. $$

If we assume that $y=0$, then:

$y^2 (y - 2x) = 0$ is satisfied and
$x^2 (y-x) +y^5 = 0 , \Leftrightarrow , -x^3 = 0 , \Leftrightarrow , x=0.$

Under this assumption, $(0,0)$ is the only equilibrium.

Otherwise (if $y\neq 0$),

$y^2 (y - 2x) = 0$ yields $y = 2x$,
The substitution of $y$ by $2x$ in $x^2 (y-x) +y^5 = 0$ yields $$x^3(1+32x^2)=0$$ and therefore $x=0$.
$y^2 (y - 2x) = 0$ becomes $y^3=0$ and thus $y=0$.

The initial assumption cannot hold.

Conclusion:

The Vinograd system has a single equilibrium: $(0, 0)$.

2. 🔓

Yes, the origin seems to be (globally) attractive.

As far as we can tell, the streamplot displays trajectories that ultimately all converge towards the origin.

3. 🔓

Let's assume that $x_e$ is a (locally) asymptotically stable of a well-posed system.

Let $r>0$ such that this property is satisfied and let

$$ B := {x \in \mathbb{R}^n ; | ; |x - x_e| \leq r } \subset \mathrm{dom} , f. $$

The set $x(t, B)$ is defined for any $t\geq 0$ and since $B$ is a neighbourhood of $x_e$, there is $\tau \geq 0$ such that for any $t \geq \tau$, the image of $B$ by $x(t, \cdot)$ is included in $B$.

$$ t\geq \tau , \Rightarrow , x(t, B) \subset B. $$

Additionally, the system is well-posed.

Hence there is a $r' > 0$ such that for any $x_0$ in the closed ball $B'$ of radius $r'$ centered at $x_e$ and any $t \in [0, \tau]$, we have

$$ |x(x_0, t) - x(x_e, t) | \leq r. $$

Since $x_e$ is an equilibrium, $x(t, x_e) = x_e$, thus $|x(x_0, t) - x_e | \leq r.$ Equivalently,

$$ 0\leq t \leq \tau , \Rightarrow , x(t, B') \subset B. $$

Note that since $x(0, B') = B'$, this inclusion yields $B' \subset B$. Thus, for any $t \geq 0$, either $t\in [0, \tau]$ and $x(t, B') \subset B$ or $t\geq \tau$ and since $B' \subset B$,

$$x(t, B') \subset x(t, B) \subset B.$$

Conclusion: we have established that there is a $r > 0$ such that $B \subset \mathrm{dom} , f$ and a $r'>0$ such that $r' \leq r$ and

$$ t\geq 0 , \Rightarrow x(t, B') \subset B. $$

In other words, the system is stable! 🎉

4. 🔓

No! We can pick initial states $(0, \varepsilon)$, with $\varepsilon >0$ which are just above the origin and still the distance of their trajectory to the origin will exceed $1.0$ at some point:

def fun(t, xy):
    return f(xy)
eps = 1e-10; xy0 = (0, eps)
sol = solve_ivp(
  fun=fun,
  y0=xy0,
  t_span=(0.0, 100.0),
  dense_output=True)["sol"]

t = linspace(0.0, 100.0, 10000)
xt, yt = sol(t)
figure()
x = y = linspace(-1.0, 1.0, 1000)
streamplot(*Q(f, x, y), color="#ced4da")
xticks([-1, 0, 1])
plot([0], [0], "k.", ms=10.0)
plot(xt, yt, color="C0")

::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

axis("square")
axis("off")
save("images/unstable")

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

{.section data-background="images/unstable.svg" data-background-size="contain"}

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

Files

asymptotic.md

Latest commit

History

asymptotic.md

File metadata and controls

Control Engineering with Python

Symbols

🐍 Imports

🐍 Streamplot Helper

ℹ️ Assumption

🏷️ Asymptotic

⚠️

Lorenz System

Hadley System

🏷️ Equilibrium

💎 Equilibrium

Stability

Attractivity

🏷️ Global Attractivity

🏷️ Local Attractivity

🔍 Global Attractivity

🐍 Vector field

📈 Stream plot

{.section data-background="images/globally-attractive.svg" data-background-size="contain"}

🔍 Local Attractivity

🐍 Vector field

📈 Stream plot

{.section data-background="images/locally-attractive.svg" data-background-size="contain"}

🔍 No Attractivity

🐍 Vector field

📈 Stream plot

{.section data-background="images/not-attractive.svg" data-background-size="contain"}

🧩 Pendulum

1. 🧮

2. 🧠

3. 📈

4. 🔬

5. 📈 Frictionless Pendulum

6. 🧮 🧠

🔓 Pendulum

1. 🔓

2. 🔓

3. 🔓

{.section data-background="images/pendulum-friction.svg" data-background-size="contain"}

4. 🔓

5. 🔓

{.section data-background="images/pendulum-no-friction.svg" data-background-size="contain"}

6. 🔓

💎 Attractivity (Low-level)

⚠️ Warning

🔍 (Pathological) Example

🐍 Vector Field

📈 Stream Plot

{.section data-background="images/attractive2.svg" data-background-size="contain"}

Asymptotic Stability

🏷️ Global Asympt. Stability

Set of Initial Conditions

💎 Global Asympt. Stability

🏷️ Limits of Sets

🏷️ Hausdorff Distance

🏷️ Local Asymptotic Stability

💎 Local Asympt. Stability {#LAS-ppty}

🏷️ Stability

🧩 Vinograd System

🐍 Vector field

📈 Stream plot

{.section data-background="images/vinograd.svg" data-background-size="contain"}

1. 🧮

2. 📈🔬

3. 🧠

4. 🧪📈

🔓 Vinograd System

1. 🔓

2. 🔓

3. 🔓

4. 🔓

{.section data-background="images/unstable.svg" data-background-size="contain"}