% Well-Posedness % 👤 Sébastien Boisgérault
| 🐍 | Code | 🔍 | Worked Example |
| 📈 | Graph | 🧩 | Exercise |
| 🏷️ | Definition | 💻 | Numerical Method |
| 💎 | Theorem | 🧮 | Analytical Method |
| 📝 | Remark | 🧠 | Theory |
| ℹ️ | Information | 🗝️ | Hint |
| Warning | 🔓 | Solution |
from numpy import *
from numpy.linalg import *
from scipy.integrate import solve_ivp
from matplotlib.pyplot import *::: notebook :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
rcParams['figure.dpi'] = 200::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::: 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)
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)Make sure that a system is "sane" (not "pathological"):
Well-Posedness:
- Existence +
- Uniqueness +
- Continuity.
We will define and study each one in the sequel.
So far, we have mostly dealt with global solutions
This concept is sometimes too stringent.
Consider the IVP
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)$")::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
tight_layout()
save("images/finite-time-blowup")
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
{.section data-background-color="#f3f0ff" data-background="images/finite-time-blowup.svg" data-background-size="contain"}
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
🤕 Ouch.
There is actually no global solution.
However there is a local solution
-
defined for
$t \in \left[t_0, \tau\right[$ -
for some
$\tau > t_0$ .
Indeed, the function
$$
\dot{x}(t) = \frac{d}{dt} x(t) = -\frac{-1}{(1 - t)^2}
= (x(t))^2
$$
and
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)$")::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
tight_layout()
save("images/finite-time-blowup-2")
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
{.section data-background="images/finite-time-blowup-2.svg" data-background-size="contain" data-background-color="#f3f0ff"}
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
This local solution is also maximal:
You cannot extend this solution beyond
A solution
A solution
A (local) solution
-
defined on
$[0, \tau'[$ with$\tau' > \tau$ , -
whose restriction to
$[0, \tau[$ is$x$ .
Consider the IVP
Find a closed-formed local solution
🗝️ Hint: assume that
Make sure that your solutions are maximal.
As long as
$$ \frac{d}{dt} \frac{1}{x(t)} =
- \frac{\dot{x}(t)}{x(t)^2} = 1. $$
By integration, this leads to
and thus provides
which is indeed a solution as long as the denominator is not zero.
-
If
$x_0 < 0$ , this solution is valid for all$t\geq 0$ and thus maximal. -
If
$x_0 > 0$ , the solution is defined until$t=1/x(0)$ where it blows up. Thus, this solution is also maximal.
Sometimes things get worse than simply having no global solution.
Consider the scalar IVP with initial value
def f(x1x2):
x1, x2 = x1x2
dx1 = 1.0 if x1 < 0.0 else -1.0
return array([dx1, 0.0])
figure()
x1 = x2 = linspace(-1.0, 1.0, 20)
gca().set_aspect(1.0)
quiver(*Q(f, x1, x2), color="k")::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
tight_layout()
save("images/discont")
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
This system has no solution, not even a local one, when
-
Assume that
$x: [0, \tau[ \to \mathbb{R}$ is a local solution. -
Since
$\dot{x}(0) = -1 < 0$ , for some small enough$0 < \epsilon < \tau$ and any$t \in \left]0, \epsilon\right]$ , we have$x(t) < 0$ . -
Consequently,
$\dot{x}(t) = +1$ and thus by integration$$ x(\epsilon) = x(0) + \int_0^{\epsilon} \dot{x}(t) , dt = \epsilon > 0, $$
which is a contradiction.
However, a local solution exists under very mild assumptions.
If
-
There is a (at least one) local solution to the IVP
$\dot{x} = f(x)$ and$x(t_0) = x_0$ . -
Any local solution on some
$\left[t_0, \tau \right[$ can be extended to a (at least one) maximal one on some$\left[t_0, t_{\infty}\right[$ .
📝 Note: a maximal solution is global iff
A solution on
-
$\tau = +\infty$ : the solution is global, or -
$\tau < +\infty$ and$\displaystyle \lim_{t \to \tau} |x(t)| = +\infty.$
In plain words : a non-global solution cannot be extended further in time if and only if it "blows up".
Let's assume that a local maximal solution exists.
You wonder if this solution is defined in
For example, you wonder if a solution is global
(if
Task. Show that any solution which defined on some
sub-interval
Then, no solution can be maximal on any
such
Consider the dynamical system
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)::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
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/sigmoid")::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
data-background-color="#f3f0ff"}
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
Show that there is a (at least one) maximal solution to each initial condition.
Show that any such solution is global.
The sigmoid function
Consequently, [💎 Existence] proves the existence of a (at least one) maximal solution.
Let
and by integration,
Thus, it cannot blow-up in finite time; by [💎 Maximal Solutions], it is global.
Consider the pendulum, subject to a torque
We assume that the torque provides a bounded power:
Show that for any initial state,
there is a global solution
🗝️ Hint. Compute the derivative with respect to
Since the system vector field
is continuous, [💎 Existence] yields the existence of a (at least one) maximal solution.
Additionally,
By integration
Hence, since
Thus,
By [💎 Maximal Solutions], any maximal solution is global.
Let
Consider the dynamical system
Show that
is differentiable and satisfies
for some
Let
Compute
Prove that for any initial state
By definition of
Let
For any vector
By the triangle inequality and the Cauchy-Schwarz inequality, we obtain
and thus
Since
Since
By integration
hence
The vector field
$$
x \in \mathbb{R}^n \to A x
$$
is continuous, thus by [💎 Existence] there is a maximal solution
Moreover,
Hence there is no finite-time blow-up and the maximal solution is global.
In the current context, uniqueness means uniqueness of the maximal solution to an IVP.
Uniqueness of solutions, even the maximal ones, is not granted either.
The IVP
has several maximal (global) solutions.
For any
However, uniqueness of maximal solution holds under mild assumptions.
Jacobian matrix of
If
An infinitely small error in the initial value could result in a finite error in the solution, even in finite time.
That would severely undermine the utility of any approximation method.
Instead of denoting
Continuity w.r.t. the initial state means that
if
and that this convergence is uniform w.r.t.
However, continuity w.r.t. the initial value holds under mild assumptions.
Assume that
Then the dynamical system is continous w.r.t. the initial state.
Let
with
alpha = 2 / 3; beta = 4 / 3; delta = gamma = 1.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,
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")::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
save("images/continuity")
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
Let
Let
Find the largest
What is the behavior of
The solution
Thus, the smallest
For any
Consider the IVP
Solve numerically this IVP for
Then, solve it again for
Does the solution seem to be continuous with respect to the initial value?
Explain this experimental result.
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()
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()::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
pp.gcf().subplots_adjust(bottom=0.2)
save("images/eps")::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
The solution does not seem to be continuous with respect to the initial value
since the graph of the solution seems to have a limit when
The jacobian matrix of the vector field is not defined when
Consider the system
where
Prove that the system is well-posed.
Prove that all maximal solutions such that
Hint 🗝️. Compute the ODE satisfied by
The jacobian matrix of the system vector field $$ f(x, y)= (\alpha x - \beta xy, \delta x y - \gamma y) $$ is defined and continuous: $$ \frac{\partial f}{\partial (x, y)}
\left[ \begin{array}{rr} \alpha -\beta y & - \beta x \ \delta y & \delta x - \gamma \ \end{array} \right] $$ thus the sytem is well-posed.
The (continuously differentiable) change of variable
$$
F: (x, y) \mapsto (u, v) := (\ln x, \ln y)
$$
is a bijection between
Since $$ \frac{d}{dt} \ln x = \frac{\dot{x}}{x}, ; \frac{d}{dt} \ln y = \frac{\dot{y}}{y} $$ the prey-predator ODE is equivalent to $$ \begin{array}{rcl} \dot{u} &=& \alpha - \beta e^v \ \dot{v} &=& \delta e^u - \gamma \ \end{array} $$
Accordingly, $$ \begin{split} \frac{d}{dt}{V} &= \delta e^u \dot{u} - \gamma \dot{u} +\beta e^v \dot{v} - \alpha \dot{v} \ &= (\delta e^u - \gamma) \dot{u} + (\beta e^v - \alpha \dot{v}) \ &= (\delta e^u - \gamma) (\alpha - \beta e^v) + (\beta e^v - \alpha) (\delta e^u - \gamma) \ & =0 \end{split} $$
Therefore
Now, the function
$$
\phi(u) := \delta e^u - \gamma u, ;
\psi(v) := \beta e^v - \alpha v
$$
are continuous and
$$
\lim_{|u| \to +\infty} \phi(u) = +\infty, ;
\lim_{|v| \to +\infty} \phi(v) = +\infty.
$$
As
Consequently, since
Therefore the solution
Since