from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())
On comprend bien que, en générale, il ne suffit pas qu'un schéma numérique soit convergent pour qu'il donne des bons résultats sur n'importe quelle équation différentielle. Encore faut-il que le problème soit
On dit qu'un problème de Cauchy est mathématiquement mal posé s'il existe pluesieurs solutions.
Par exemple, on cherche une fonction $y\colon t\in\R^+\mapsto y(t)\in\R$ qui satisfait
$$
\begin{cases}
y'(t) = \sqrt[3]{y(t)}, &\forall t>0,\\
y(0) = 0.
\end{cases}
$$
Il est simple de vérifier que les fonctions $y_1(t)=0$ et $y_{2,3}(t)=\pm\sqrt{8t^3/27}$, pour tout $t\ge0$, sont toutes les trois solution du problème de Cauchy donné.
Dans ce cas, lorsqu'on utilise une méthode numérique, on ne sait pas quelle solution ce schéma approche et différents schémas peuvent approcher différentes solutions.
Correction
La fonction $\varphi(t,y)=\sqrt{y}$ n'est pas lipschitzienne par rapport à $y$, donc le théorème d'existence et unicité locale n'est pas valable au voisinage de $(0,0)$.
L'EDO est à variables séparables, on peut donc expliciter toutes les solutions du problème de Cauchy.
En imposant la CI on trouve que, pour tout $b\in\R^+$, les fonctions $$ y_b(t)= \begin{cases} 0 , &\text{si }0\le t\le b, \\ \frac{1}{4}(t-b)^{2} , &\text{si }t\ge b, \end{cases} $$ sont de classe $\mathcal{C}^1(\R^+)$ et sont solution du problème de Cauchy donné.
Traçons quelques curbes:
%matplotlib inline
from matplotlib.pylab import *
b=1
xx=linspace(-1,2,101)
yyr=[0 for x in xx]
def f(x,b):
if x<=b:
return 0
else:
return (x-b)**2/4
yyb=[f(x,b) for x in xx]
plot(xx,yyr,'r-',xx,yyb,'b:');
legend([r'$y(t)=0$ $\forall$ $t$', '$y_{b=1}(t)$']);
La méthode d'Euler explicite construit la suite \begin{cases} u_0=y_0=0,\\ u_{n+1}=u_n+h \varphi(t_n,u_n)=u_n+hu_n^{1/2},& n=0,1,2,\dots N_h-1 \end{cases} par conséquent $u_n=0$ pour tout $n$. La méthode d'Euler explicite approche la solution constante $y(t)=$ pour tout $t\in\R^+$.
La méthode d'Euler implicite construit la suite \begin{cases} u_0=y_0=0,\\ u_{n+1}=u_n+h \varphi(t_{n+1},u_{n+1})=u_n+hu_{n+1}^{1/2},& n=0,1,2,\dots N_h-1. \end{cases} par conséquent $u_0=0$ mais $u_1$ dépend de la méthode de résolution de l'équation implicite $x=0+h\sqrt{x}$. Bien sur $x=0$ est une solution mais $x=h^{2}$ est aussi solution. Si le schéma choisit $u_1=h^{2}$, alors $u_n>0$ pour tout $n\in\N^*$.
Notons que le problème de Cauchy avec une CI $y(0)=y_0>0$ admet une et une seule solution, la fonction $y(t)=\frac{1}{4}(t+2\sqrt{y_0})^{2}$. Dans ce cas, les deux schémas approchent forcement la même solution.
Une fois calculée la solution numérique $\{u_n\}_{n=1}^{N_h}$ d'un problème de Cauchy, il est légitime de chercher à savoir dans quelle mesure l'erreur $|y(t_n)-u_n|$ est petite pour $n=1,2,\dots$. Cela dépend bien sur du schéma choisi mais aussi du problème en question. Essayons de voir cela sur un exemple.
On se donne $\varphi(t,y(t))=3t-3y(t)$ et $y_0=\alpha$ (un nombre quelconque). On cherche une fonction $y\colon t\in\R\mapsto y(t)\in\R$ qui satisfait $$ \begin{cases} y'(t) = 3t-3y(t), &\forall t\in\R,\\ y(0) = \alpha. \end{cases} $$ Sa solution, définie sur $\R$, est donnée par $y(t)=(\alpha+1/3)e^{-3t}+t-1/3$. En effet on a bien \begin{align*} y(0)&=(\alpha+1/3)e^0+0-1/3=\alpha, \\ y'(t)&=-3(\alpha+1/3)e^{-3t}+1=-3(\alpha+1/3)e^{-3t}+1-3t+3t=-3y(t)+3t. \end{align*} Si nous cherchons à résoudre le problème de Cauchy jusqu'à $t=10$ avec $\alpha=1/3$, nous obtenons $y(10)=10+1/3=31/3$. Par contre, si nous faisons le calcul avec l'approximation $\alpha=0.333333$ au lieu de $1/3$, nous avons $y(10)=(0.333333-1/3)e^{30}+10+1/3=-e^{30}/3000000+31/3$ ce qui représente une différence avec la précédente valeur de $e^{30}/3000000\approx10^7/3$. Cet exemple nous apprend qu'une petite erreur sur la condition initiale (erreur relative d'ordre $10^{-6}$) peut provoquer une très grande erreur sur $y(10)$ (erreur relative d'ordre $10^{6}$). Ainsi, si le calculateur mis à notre disposition ne calcul qu'avec $6$ chiffres significatifs (en virgule flottante), alors $\alpha=1/3$ devient $\alpha=0.333333$ et il est inutile d'essayer d'inventer une méthode numérique pour calculer $y(10)$. En effet, la seule erreur sur la condition initiale provoque déjà une erreur inadmissible sur la solution. Nous sommes en présence ici d'un problème numériquement mal posé
Voici d'autres exemples :
Considérons le problème de Cauchy \begin{cases} y'(t) = -y(t), &\forall t>0,\\ y(0) = y_0+\varepsilon. \end{cases} La solution est $y(t)=y_0e^{-t}+\varepsilon e^{-t}$: l'effet de la perturbation $\varepsilon$ s’atténue lorsque $t\to+\infty$ puisque $\varepsilon e^{-t}\xrightarrow[t\to+\infty]{}0$. Cela suggère que si une erreur est faite dans une étape d'une méthode d’itération, l'effet de cette erreur s’atténue au cours du temps: le problème est numériquement bien posé.
%matplotlib inline
from matplotlib.pylab import *
tt=linspace(0,5,101)
yye=[(1+0.1)*exp(-t) for t in tt]
yy0=[(1)*exp(-t) for t in tt]
axis([0,5,0,1.1])
plot(tt,yye,tt,yy0);
legend(["$y_0=1$","$y_0=1.01$"]);
Considérons maintenant le problème de Cauchy \begin{cases} y'(t) = y(t), &\forall t>0,\\ y(0) = y_0+\varepsilon. \end{cases} La solution est $y(t)=y_0e^{t}+\varepsilon e^{t}$: l'effet de la perturbation $\varepsilon$ s’amplifie lorsque $t\to+\infty$ puisque $\varepsilon e^{t}\xrightarrow[t\to+\infty]{}+\infty$. Cela suggère que si une erreur est faite dans une étape d'une méthode d’itération, l'effet de cette erreur s’amplifie au cours du temps: le problème est numériquement mal posé.
%matplotlib inline
from matplotlib.pylab import *
tt=linspace(0,5,101)
yye=[(0+0.1)*exp(t) for t in tt]
yy0=[(0)*exp(t) for t in tt]
axis([0,5,0,1.1])
plot(tt,yye,lw='2')
plot(tt,yy0,lw='5');
legend(["$y_0=0$","$y_0=0.01$"]);
Considérons le problème de Cauchy \begin{cases} y'(t) = y(t)\sin(t), &\forall t>0,\\ y(0) = -1+\varepsilon. \end{cases} La solution est $y(t)=-e^{1-\cos(t)}+\varepsilon e^{1-\cos(t)}$: l'effet de la perturbation $\varepsilon$ reste borné lorsque $t\to+\infty$ puisque $|\varepsilon e^{1-\cos(t)}|\le \varepsilon e^2$. Cela suggère que si une erreur est faite dans une étape d'une méthode d’itération, l'effet de cette erreur reste borné au cours du temps: le problème est numériquement bien posé.
%matplotlib inline
from matplotlib.pylab import *
tt=linspace(0,10,101)
yye=[(-1+0.1)*exp(1-cos(t)) for t in tt]
yy0=[(-1)*exp(1-cos(t)) for t in tt]
# axis([0,5,-5,0])
plot(tt,yye,tt,yy0);
legend(["$y_0=-0.9$","$y_0=-1$"]);
L'objectif de cet exercice est de bien mettre en évidence l'influence d'une petite perturbation de la condition initiale sur la solution d'un problème numériquement mal posé et d'appliquer la méthode d'Euler explicite pour en approcher la solution exacte.
Considérons le problème de Cauchy \begin{cases} y'(t)=3\dfrac{y(t)}{t}-\dfrac{5}{t^3},\\ y(1)=a \end{cases}
Correction
Calculons les solutions de l'edo pour $t>0$. On pose $u(t)=\frac{y(t)}{t}$ ainsi $y(t)=tu(t)$ et $y'(t)=u(t)+tu'(t)$. On obtient l'edo linéaire d'ordre 1 suivante: $$ tu'(t)-2u(t)=-5t^{-3}. $$ On a $a(t)=t$, $b(t)=-2$ et $g(t)=-5t^{-3}$ donc
En imposant la condition $a=y(1)$ on trouve l'unique solution du problème de Cauchy donné: $$ y(t)=(a-1)t^3+\dfrac{1}{t^2}. $$
%matplotlib inline
from math import *
from matplotlib.pylab import *
sol_exacte = lambda t,y0 : (y0-1.)*t**3+1./t**2
# INITIALISATION
t0 = 1.
tfinal = 5.
c=1.
tt=linspace(t0,tfinal,101)
Y0=[c+1.e-2,c+1.e-3,c,c-1.e-3,c-1.e-2]
for y0 in Y0:
yy=[sol_exacte(t,y0) for t in tt]
plot(tt,yy,label=('$y_0=$'+str(y0)))
legend()
title('Solution exacte - differents valeurs de $y_0$')
axis([t0,tfinal,-c,c]);
On constate qu'une petite perturbation de la condition initiale entraîne une grosse perturbation de la solution, en effet
$$
\lim_{t\to+\infty}y(t)=
\begin{cases}
+\infty &\text{si } a>1,\\
0^+ &\text{si } a=1,\\
-\infty &\text{si } a<1.
\end{cases}
$$
Si $a=1$, la solution exacte est $y(t)=t^{-2}$ mais le problème est numériquement mal posé car toute (petite) erreur de calcul a le même effet qu'une perturbation de la condition initiale: on "réveil" le terme dormant $t^{3}$.
from scipy.optimize import fsolve
def euler_progressif(phi,tt):
h=tt[1]-tt[0]
uu = [y0]
for i in range(len(tt)-1):
uu.append(uu[i]+h*phi(tt[i],uu[i]))
return uu
def euler_regressif(phi,tt):
h=tt[1]-tt[0]
uu = [y0]
for i in range(len(tt)-1):
uu.append(fsolve(lambda x : -x +uu[i]+h*phi(tt[i+1],x),uu[i]))
return uu
phi = lambda t,y : 3*y/t-5./t**3
y0=Y0[2]
NN=[20,100,300,5000]
figure(1, figsize=(15, 5))
subplot(1,2,1)
for N in NN:
tt=linspace(t0,tfinal,N)
uu=euler_progressif(phi,tt)
plot(tt,uu,'.-')
plot(tt,[sol_exacte(t,y0) for t in tt],'r-',lw=2)
legend(['N=20','N=100','N=300','N=5000','Exacte'])
title('$y_0=1$ - Exacte vs EE')
axis([t0,tfinal,-c,c]);
subplot(1,2,2)
for N in NN:
tt=linspace(t0,tfinal,N)
uu=euler_regressif(phi,tt)
plot(tt,uu,'.-')
plot(tt,[sol_exacte(t,y0) for t in tt],'r-',lw=2)
legend(['N=20','N=100','N=300','N=5000','Exacte'])
title('$y_0=1$ - Exacte vs EI')
axis([t0,tfinal,-c,c]);
Considérons un problème de Cauchy et supposons que l'on ait montré l'existence d'une solution $y$ (donc il est mathématiquement bien posé) et qu'il est numériquement bien posé. Deux questions naturelles se posent:
Dans les deux cas le nombre de nœuds tend vers l'infini mais dans le premier cas on s'intéresse à l'erreur en chaque point, dans le deuxième cas il s'agit du comportement asymptotique de la solution et de son approximation.
Ces deux questions font intervenir deux notions de stabilité:
Soit $\beta>0$ un nombre réel positif et considérons le problème de Cauchy \begin{cases} y'(t)=-\beta y(t), &\text{pour }t>0,\\ y(0)=y_0 \end{cases} où $y_0\neq0$ est une valeur donnée. Sa solution est $y(t)=y_0e^{-\beta t}$ donc $$\lim_{t\to+\infty}y(t)=0.$$
Soit $h>0$ un pas de temps donné, $t_n=nh$ pour $n\in\N$ et notons $u_n\approx y(t_n)$ une approximation de la solution $y$ au temps~$t_n$.
Si, sous d'éventuelles conditions sur $h$, on a $$ \lim_{n\to+\infty} u_n =0, $$ alors on dit que le schéma est A-stable.
On peut tirer des conclusions analogues quand $\beta$ est un complexe ou une fonction positive de $t$.
D'autre part, en général il n'y a aucune raison d'exiger qu'une méthode numérique soit absolument stable quand on l'applique à un autre problème.
Cependant, on peut montrer que quand une méthode absolument stable sur le problème modèle est utilisée pour un problème modèle généralisé, l'erreur de perturbation (qui est la valeur absolue de la différence entre la solution perturbée et la solution non perturbée) est bornée uniformément (par rapport à $h$).
En bref, on peut dire que les méthodes absolument stables permettent de contrôler les perturbations.
Remarque: ce résultat peut être généralisé à un système de $n$ EDO de la forme
$$
\mathbf{y}'(t)=-\mathbb{A}\mathbf{y}(t)
$$
où $\mathbb{A}$ est une matrice constante ayant $n$ valeurs propres positives $\lambda_i$, $i=1,2,...,n$.
La solution s'écrit
$$
\mathbf{y}(t)=\sum_i c_i \mathbf{v}_i e^{-\lambda_i t}
$$
où $\mathbf{v}_i$ est le $i$-ème vecteur propre associé à la valeure propre $\lambda_i$ donc
$$
\lim_{t\to+\infty}y_i(t)=0 \quad \forall i=1,2,...,n.
$$
Le schéma d'Euler progressif devient $$ u_{n+1}=(1-\beta h)u_n, \qquad n=0,1,2,\dots $$ et par suite $$ u_{n}=(1-\beta h)^nu_0, \qquad n=0,1,2,\dots $$ Par conséquente, $\lim\limits_{n\to+\infty}u_n=0$ si et seulement si $$ \left|1-\beta h\right|<1, $$ ce qui a pour effet de limiter $h$ à $$ h<\frac{2}{\beta}. $$ Cette condition de A-stabilité limite le pas $h$ d'avance en $t$ lorsqu'on utilise le schéma d'Euler progressif. Notons que si $1-\beta h>1$ alors $u_n$ tend vers $+\infty$ lorsque $t$ tend vers l'infini et si $1-\beta h<-1$ alors $u_n$ tend vers l'infini en alternant de signe lorsque $t$ tend vers l'infini. Nous dirons dans ces cas que le schéma d'Euler progressif est instable.
Remarque: dans le cas du système, le schéma d'Euler progressif est A-stable ssi $$ h<\frac{2}{\lambda_\text{max}}. $$
Le schéma d'Euler rétrograde devient dans le cadre de notre exemple $$ (1+\beta h)u_{n+1}=u_n, \qquad n=0,1,2,\dots $$ et par suite $$ u_{n}=\frac{1}{(1+\beta h)^{n}}u_0, \qquad n=0,1,2,\dots $$ Dans ce cas nous voyons que pour tout $h>0$ nous avons $\lim_{n\to\infty}u_n=0$, le schéma d'Euler rétrograde est donc toujours stable, sans limitations sur $h$.
Le schéma d'Euler modifié pour notre exemple devient $$ u_{n+1}=u_n+h\left( -\beta\left( u_n+\frac{h}{2}(-\beta u_n) \right) \right) =\left(1-\beta h +\frac{(\beta h )^2}{2} \right) u_{n} $$ Par induction on obtient $$ u_{n}=\left(1-\beta h +\frac{(\beta h )^2}{2} \right)^nu_0. $$ Par conséquent, $\lim\limits_{n\to+\infty}u_n=0$ si et seulement si $$ \left|1-\beta h +\beta^2\frac{( h )^2}{2}\right|<1. $$ Notons $x$ le produit $\beta h $ et $q$ le polynôme $q(x)=\frac{1}{2}x^2-x+1$ dont le graphe est représenté en figure:
%matplotlib inline
from matplotlib.pylab import *
xx=linspace(0,3,101)
yy=[1-x+0.5*x**2 for x in xx]
plot(xx,yy,[0,2,2],[1,1,0],'m:')
axis([0,3,0,3])
xlabel(r'$h\beta$')
ylabel(r'$q(h\beta)$');
Nous avons $|q(x)|<1$ si et seulement si $0<x<2$. La relation $\lim\limits_{n\to+\infty}u_n=0$ est donc satisfaite si et seulement si $$ h <\frac{2}{\beta}. $$ Cette condition de stabilité limite le pas $h$ d'avance en $t$ lorsqu'on utilise le schéma d'Euler modifié.
Le schéma de Crank-Nicolson appliqué à notre exemple s'écrit
$$
\left(1+\beta\frac{h}{2}\right)u_{n+1} = \left(1-\beta\frac{h}{2}\right) u_{n}
$$
et par suite
$$
u_{n}=\left( \frac{2-\beta h}{2+\beta h} \right)^{n}u_0, \qquad n=0,1,2,\dots
$$
Par conséquent, $\lim\limits_{n\to+\infty}u_n=0$ si et seulement si
$$
\left|\frac{2-\beta h}{2+\beta h}\right|<1.
$$
Notons $x$ le produit $\beta h>0$ et $q$ la fonction $q(x)=\frac{2-x}{2+x}=1-2\frac{x}{2+x}$. Nous avons $0<\frac{x}{2+x}<1$ pour tout $x\in\R_+^*$, donc $|q(x)|<1$ pour tout $x\in\R_+^*$.
La relation $\lim\limits_{n\to+\infty}u_n=0$ est donc satisfaite pour tout $h>0$: le schéma de Crank-Nicolson est donc toujours stable, sans limitations sur $h$.
Cependant, si on cherche $0<q<1$ on a (pour $x>0$)
$$
1-2\frac{x}{2+x}>0
\iff
x<2
$$
%matplotlib inline
from matplotlib.pylab import *
xx=linspace(0,10,101)
yy=[1-2*x/(2+x) for x in xx]
plot(xx,yy)
plot([0,10],[1,1],'r--',[0,10],[-1,-1],'r--')
plot([0,10],[0,0],'m:',[2,2],[-1.5,1.5],'m:')
axis([0,10,-1.5,1.5])
xlabel(r'$h\beta$')
ylabel(r'$q(h\beta)$');
Le schéma de Heun pour notre exemple devient $$ u_{n+1} =u_n+\frac{h}{2}\left( -\beta u_n-\beta\left( u_n+h(-\beta u_n) \right) \right) =\left(1-\beta h +\frac{(\beta h )^2}{2} \right) u_{n} $$ ce qui coïncide avec la méthode d'Euler modifiée. Par conséquent, $\lim\limits_{n\to+\infty}u_n=0$ si et seulement si $$ h <\frac{2}{\beta}. $$ Cette condition de stabilité limite le pas $h$ d'avance en $t$ lorsqu'on utilise le schéma de Heun.
Le schéma de Simpson implicite appliqué à notre exemple s'écrit $$ \left(1+\beta\frac{h}{6}\right)u_{n+1} = \left(1-\frac{5}{6}\beta h+\frac{1}{3}(\beta h)^2\right) u_{n} $$ et par suite $$ u_{n}=\left( \frac{ 1-\frac{5}{6}\beta h+\frac{1}{3}(\beta h)^2 }{ 1+\beta\frac{h}{6} } \right)^{n}u_0, \qquad n=0,1,2,\dots $$ Par conséquent, $\lim\limits_{n\to+\infty}u_n=0$ si et seulement si $$ \left|\frac{ 6-5\beta h+2(\beta h)^2 }{ 6+\beta h }\right|<1. $$ Notons $x$ le produit $\beta h>0$ et $q$ la fonction $q(x)=\frac{6-5x+2x^2}{6+x}$ dont le graphe est représenté en figure.
%matplotlib inline
from matplotlib.pylab import *
xx=linspace(0,4,101)
yy=[(6-5*x+2*x**2)/(6+x) for x in xx]
plot(xx,yy,[0,3,3],[1,1,0],'m:')
axis([0,4,0,3])
xlabel(r'$h\beta$')
ylabel(r'$q(h\beta)$');
Nous avons $q(x)>0$ pour tout $x\in\R_+^*$ et $q(x)<1$ si et seulement si $x<3$, donc $|q(x)|<1$ pour tout $x\in]0;3[$. Par conséquent, $\lim\limits_{n\to+\infty}u_n=0$ si et seulement si $$ h <\frac{3}{\beta}. $$ Cette condition de stabilité limite le pas $h$ d'avance en $t$.
Le schéma de Simpson explicite pour notre exemple devient $$ u_{n+1} = u_n+\frac{h}{6}\left( -\beta u_n -4 \beta\left( u_n+\frac{h}{2}(-\beta u_n) \right) -\beta\left( u_n+h(-\beta u_n) \right) \right) =\left(1-\beta h +\frac{3}{2}(\beta h )^2 \right) u_{n} $$ Par induction on obtient $$ u_{n}=\left(1-\beta h +\frac{3}{2}(\beta h )^2 \right)^nu_0. $$ Par conséquent, $\lim\limits_{n\to+\infty}u_n=0$ si et seulement si $$ \left|1-\beta h +\frac{3}{2}(\beta h )^2\right|<1. $$ Notons $x$ le produit $\beta h $ et $q$ le polynôme $q(x)=\frac{3}{2}x^2-x+1$ dont le graphe est représenté en figure.
%matplotlib inline
from matplotlib.pylab import *
xx=linspace(0,3,101)
yy=[1-x+1.5*x**2 for x in xx]
plot(xx,yy,[0,2/3,2/3],[1,1,0],'m:')
axis([0,3,0,3])
xlabel(r'$h\beta$')
ylabel(r'$q(h\beta)$');
Nous avons $|q(x)|<1$ si et seulement si $0<x<\frac{2}{3}$. La relation $\lim\limits_{n\to+\infty}u_n=0$ est donc satisfaite si et seulement si $$ h <\frac{2}{3\beta}. $$ Cette condition de stabilité limite le pas $h$ d'avance en $t$.
On considère le problème de Cauchy $$ \begin{cases} y'(t)=-y(t),\\ y(0)=1, \end{cases} $$ sur l'intervalle $[0;12]$.
Il s'agit d'une EDO à variables séparables. L'unique solution constante de l'EDO est la fonction $y(t)\equiv0$, toutes les autres solutions sont du type $y(t)=C e^{-t}$. Donc l'unique solution du problème de Cauchy est la fonction $y(t)=e^{-t}$ définie pour tout $t$ $\in$ $\R$.
La méthode d'Euler explicite pour cette EDO s'écrit $$ u_{n+1}=(1-h)u_n. $$ En procédant par récurrence sur $n$, on obtient $$ u_{n+1}=(1-h)^{n+1}. $$ De la formule $u_{n+1}=(1-h)^{n+1}$ on déduit que
Remarque: la suite obtenue est une suite géométrique de raison $q=1-h$. On sait qu'une telle suite
Voyons graphiquement ce que cela donne avec différentes valeurs de $h>0$:
%matplotlib inline
from matplotlib.pylab import *
# SCHEMA
def euler_progressif(phi,tt):
h=tt[1]-tt[0]
uu = [y0]
for i in range(len(tt)-1):
uu.append(uu[i]+h*phi(tt[i],uu[i]))
return uu
phi = lambda t,y : -y
sol_exacte = lambda t : exp(-t)
# INITIALISATION
t0, y0, tfinal = 0 , 1 , 12
# CALCUL
H = [ 2**(k-2) for k in range(5) ]
tt = [] # liste de liste
uu = [] # liste de liste
for h in H:
Nh = int((tfinal-t0)/h)
tt.append( [ t0+i*h for i in range(Nh+1) ] )
uu.append(euler_progressif(phi,tt[-1]))
# AFFICHAGE
figure(1, figsize=(15, 10))
yy = [sol_exacte(t) for t in tt[0]] # affichage de la sol exacte sur la grille la plus fine
axis([t0, tfinal, -8, 8])
plot(tt[0],yy,'-',label=('$y(t)=e^{-t}$'))
for k in range(5):
plot(tt[k],uu[k],'-o',label=('h='+str(H[k])))
legend()
grid();
La méthode d'Euler implicite pour cette EDO s'écrit $$ u_{n+1}=\frac{1}{1+h}u_n. $$ En procédant par récurrence sur $n$, on obtient $$ u_{n+1}=\frac{1}{(1+h)^{n+1}}. $$ De la formule $u_{n+1}=(1+h)^{-(n+1)}$ on déduit que la solution numérique est stable et convergente pour tout $h$. En effet, la méthode est inconditionnellement A-stable.
Remarque: la suite obtenue est une suite géométrique de raison $q=1/(1+h)\in]0;1[$.
%matplotlib inline
from matplotlib.pylab import *
from scipy.optimize import fsolve
# SCHEMA
def euler_regressif(phi,tt):
h=tt[1]-tt[0]
uu = [y0]
for i in range(len(tt)-1):
uu.append(fsolve(lambda x:-x+uu[i]+h*phi(tt[i+1],x),uu[i]))
return uu
phi = lambda t,y : -y
sol_exacte = lambda t : exp(-t)
# INITIALISATION
t0, y0, tfinal = 0 , 1 , 12
# CALCUL
H = [ 2**(k-2) for k in range(5) ]
tt = [] # liste de liste
uu = [] # liste de liste
for h in H:
Nh = int((tfinal-t0)/h)
tt.append( [ t0+i*h for i in range(Nh+1) ] )
uu.append(euler_regressif(phi,tt[-1]))
# AFFICHAGE
figure(1, figsize=(15, 10))
yy = [sol_exacte(t) for t in tt[0]] # affichage de la sol exacte sur la grille la plus fine
axis([t0, tfinal, -1, 1])
plot(tt[0],yy,'-',label=('$y(t)=e^{-t}$'))
for k in range(5):
plot(tt[k],uu[k],'-o',label=('h='+str(H[k])))
legend()
grid();
Soient $b,g>0$ et considérons le problème de Cauchy \begin{cases} y'(t)=-by(t)+g,& t\in[0;1]\\ y(0)=y_0. \end{cases}
Correction
Solution exacte
Il s'agit d'une edo linéaire d'ordre 1 dont la solution est
$$
y(t)=\left( y_0-\frac{g}{b} \right)e^{-bt}+\frac{g}{b}
$$
et
$$
\lim_{t\to+\infty}y(t)=\frac{g}{b} \text{ pour tout } y_0\in\R
$$
le problème est numériquement bien posé, c'est-à-dire que l'erreur sur la solution sera au pire de l'ordre de l'erreur sur la condition initiale.
Lorsque $b$ devient de plus en plus grand, la solution devient de plus en plus raide:
%reset -f
%matplotlib inline
from matplotlib.pylab import *
exacte = lambda t,b,g,y0 : (y0-g/b)*exp(-b*t)+g/b
t0 = 0.
tfinal = 1.0
tt=linspace(t0,tfinal,10001)
figure(1, figsize=(15, 5))
subplot(1,2,1)
g=10.
b=10.
Y0=[g/b,g/b+1.e-8]
for y0 in Y0:
yy=[exacte(t,b,g,y0) for t in tt]
plot(tt,yy,label=('$y_0=$'+str(y0)));
title("$g=b=$"+str(g))
axis([0,1,g/b-1.e-8,g/b+1.e-8])
legend()
grid()
subplot(1,2,2)
g=40.
b=40.
Y0=[g/b,g/b+1.e-8]
for y0 in Y0:
yy=[exacte(t,b,g,y0) for t in tt]
plot(tt,yy,label=(r'$y_0=$'+str(y0)));
title("$g=b=$"+str(g))
axis([0,1,g/b-1.e-8,g/b+1.e-8])
legend()
grid()
Méthode d'Euler explicite
u_n-\frac{g}{b}=(1-hb)^n\left( y_0 -\frac{g}{b}\right)
$$%reset -f
%matplotlib inline
from matplotlib.pylab import *
phi = lambda t,y,b,g : -b*y+g
def euler_progressif(phi,tt,y0,b,g):
h=tt[1]-tt[0]
uu = [y0]
for i in range(len(tt)-1):
uu.append(uu[i]+h*phi(tt[i],uu[i],b,g))
return uu
t0 = 0.
tfinal = 1.0
g=40.
b=40.
Y0=[g/b,g/b+1.e-8]
NN=[19,23,45,100]
figure(1, figsize=(15, 5))
subplot(1,2,1)
y0=Y0[0]
for N in NN:
tt=linspace(t0,tfinal,N+1)
yy=euler_progressif(phi,tt,y0,b,g)
plot(tt,yy,label=('$N=$'+str(N)+' noeuds'));
title("Euler explicite, $g=b=$"+str(g)+ ", $y_0$="+str(y0))
axis([0,1,g/b-1.e-8,g/b+1.e-8])
legend()
grid()
subplot(1,2,2)
y0=Y0[1]
for N in NN:
tt=linspace(t0,tfinal,N+1)
yy=euler_progressif(phi,tt,y0,b,g)
plot(tt,yy,label=('$N=$'+str(N)+' noeuds'));
title("Euler explicite, $g=b=$"+str(g)+ ", $y_0$="+str(y0))
axis([0,1,g/b-1.e-8,g/b+1.e-8])
legend()
grid()
Méthode d'Euler implicite
%reset -f
%matplotlib inline
from matplotlib.pylab import *
from scipy.optimize import fsolve
phi = lambda t,y,b,g : -b*y+g
def euler_regressif(phi,tt,y0,b,g):
h=tt[1]-tt[0]
uu = [y0]
for i in range(len(tt)-1):
uu.append(g/b + (y0-g/b)/(1.+h*b)**i)
#uu.append(fsolve(lambda x: x-uu[i]-h*phi(tt[i+1],x,b,g), uu[i]))
return uu
t0 = 0.
tfinal = 1.0
g=40.
b=40.
Y0=[g/b,g/b+1.e-8]
NN=[19,23,45,100]
figure(1, figsize=(15, 5))
subplot(1,2,1)
y0=Y0[0]
for N in NN:
tt=linspace(t0,tfinal,N+1)
yy=euler_regressif(phi,tt,y0,b,g)
plot(tt,yy,label=('$N=$'+str(N)+' noeuds'));
title("Euler implicite, $g=b=$"+str(g)+ ", $y_0$="+str(y0))
axis([0,1,g/b-1.e-8,g/b+1.e-8])
legend()
grid()
subplot(1,2,2)
y0=Y0[1]
for N in NN:
tt=linspace(t0,tfinal,N+1)
yy=euler_regressif(phi,tt,y0,b,g)
plot(tt,yy,'-*',label=('$N=$'+str(N)+' noeuds'));
title("Euler implicite, $g=b=$"+str(g)+ ", $y_0$="+str(y0))
axis([0,1,g/b-1.e-8,g/b+1.e-8])
legend()
grid()
Soit le problème de Cauchy: \begin{cases} y'(t)+10y(t)=0, & \forall t \in \R,\\ y(0)=y_0>0. \end{cases}
C'est un problème deCauchy du type \begin{cases} y'(t)=\varphi(t,y(t)), & \forall t \in \R,\\ y(0)=y_0>0, \end{cases} avec $\varphi(t,y)=g(y)=-10y$.
On montre d'abord qu'il existe une et une seule solution locale (ie sur $[-T;T]$) de classe $\mathscr{C}^1([-T,T],\R)$.
On montre ensuite que cette solution est de classe $\mathscr{C}^\infty([-T,T],\R)$.
On montre enfin que la solution admet un prolongement sur $\R$.
Pour en calculer la solution, on remarque qu'il s'agit d'une EDO à variables séparables. L'unique solution constante est $y(t)\equiv0$, toutes les autres solutions sont du type $y(t)=C e^{-10t}$. En prenant en compte la condition initiale on conclut que l'unique solution du problème de Cauchy est $$ y(t)=y_0e^{-10t} $$ définie pour tout $t$ $\in$ $\R$.
y0=1
exacte = lambda t : y0*exp(-10*t)
tt=linspace(0,3,101)
yy=[exacte(t) for t in tt]
plot(tt,yy);
Soit le schéma numérique de Crank-Nicolson défini par la suite $(u_n)_{n\in \N}$ vérifiant $$ \frac{u_{n+1}-u_n}{h}+5(u_{n+1}+u_n)=0, \qquad \forall n \in \N, $$ pour $h>0$ fixé. On obtient une formule de récurrence rendue explicite par un calcul élémentaire: $$ u_{n+1}=-5hu_{n+1}-5hu_n+u_n $$ d'où $$ u_{n+1}=\frac{1-5h}{1+5h}u_n. $$ Par récurrence on obtient $$ u_{n}=\left(\frac{1-5h}{1+5h}\right)^nu_0. $$ Il s'agit d'une suite géométrique de raison $$ r=\frac{1-5h}{1+5h}. $$
Le schéma génère une suite positive ssi $$ 1-\frac{10h}{1+5h}>0 $$ c'est-à-dire ssi $$ h<\frac{1}{5}. $$
%matplotlib inline
from matplotlib.pylab import *
N=6
tt=linspace(0,1,N+1)
y0=1
exacte = lambda t : y0*exp(-10*t)
yy=[exacte(t) for t in tt]
phi = lambda t,y : -y**5
def schema(phi,tt):
h=tt[1]-tt[0]
uu=[y0]
for i in range(len(tt)-1):
uu.append((1-5*h)/(1+5*h)*uu[i])
return uu
plot(tt,yy,label=("Exacte"))
plot(tt,schema(phi,tt),"-.",label=("Schéma"))
legend();
Soit le problème de Cauchy: \begin{cases} y'(t)+\dfrac{\sqrt{y(t)}}{2}=0, & \forall t \in \R^+,\\ y(0)=y_0>0. \end{cases}
C'est un problème de Cauchy du type \begin{cases} y'(t)=\varphi(t,y(t)), & \forall t \in \R,\\ y(0)=y_0>0, \end{cases} avec $\varphi(t,y)=g(y)=-\frac{\sqrt{y}}{2}$.
Pour $h>0$ fixé on obtient une formule de récurrence rendue explicite par un calcul élémentaire: $$ u_{n+1}=\frac{u_n}{1 + \frac{h}{2\sqrt{u_n}}}=\frac{2(u_n)^{3/2}}{2\sqrt{u_n}+h}, \qquad \forall n \in \N. $$
On étudie la suite \begin{cases} u_0>0,\\ u_{n+1}=\frac{2(u_n)^{3/2}}{2\sqrt{u_n}+h}, \qquad \forall n \in \N, \end{cases} pour $h>0$ fixé.
Par récurrence on montre que si $u_0>0$ alors $u_n>0$ pour tout $n$ $\in$ $\N$. De plus, $\frac{u_{n+1}}{u_{n}}=\frac{1}{1 + \frac{h}{2\sqrt{u_n}}}<1$ pour tout $n$ $\in$ $\N$: la suite est monotone décroissante. On a alors une suite décroissante et bornée par zéro, donc elle converge. Soit $\ell$ la limite de cette suite, alors $\ell\ge0$ et $\ell= \frac{2\ell^{3/2}}{2\sqrt{\ell}+h}$ donc $\ell=0$.
%matplotlib inline
from matplotlib.pylab import *
y0=1
phi = lambda t,y : -0.5*sqrt(y)
exacte = lambda t: 0 if t>=4*sqrt(y0) else (sqrt(y0)-t/4)**2
N=10
tt=linspace(0,10,N+1)
def schema(phi,tt):
h=tt[1]-tt[0]
uu=[y0]
for i in range(len(tt)-1):
uu.append(2*(uu[i])**(3/2)/(h+2*sqrt(uu[i])))
return uu
yy=[exacte(t) for t in tt]
plot(tt,yy,'r-',label=('Exacte'))
plot(tt,schema(phi,tt),'b*',label=("Schéma"))
legend();
Soit le problème de Cauchy: \begin{cases} y'(t)+y^5(t)=0, & \forall t \in \R^+,\\ y(0)=y_0>0. \end{cases}
C'est un problème deCauchy du type \begin{cases} y'(t)=\varphi(t,y(t)), & \forall t \in \R,\\ y(0)=y_0>0, \end{cases} avec $\varphi(t,y)=g(y)=-y^5$.
On montre d'abord qu'il existe une et une seule solution locale (ie sur $[0;T]$) de classe $\mathscr{C}^1([0,T],\R)$.
On montre ensuite que cette solution est de classe $\mathscr{C}^\infty([0,T],\R)$.
On montre enfin que la solution admet un prolongement sur $\R$.
Pour en calculer la solution, on remarque qu'il s'agit d'une EDO à variables séparables. L'unique solution du problème de Cauchy est $$ y(t)=y_0(4ty_0^4+1)^{-1/4} $$
y0=1
exacte = lambda t : y0*(4*t*y0**4+1)**(-0.25)
tt=linspace(0,10,101)
yy=[exacte(t) for t in tt]
plot(tt,yy);
Pour $h>0$ fixé on obtient une formule de récurrence rendue explicite par un calcul élémentaire: $$ u_{n+1}=\frac{u_n}{1 + u_n^4 h}, \qquad \forall n \in \N. $$
On étudie la suite \begin{cases} u_0=y_0>0,\\ u_{n+1}=\frac{u_n}{1 + u_n^4 h}, \qquad \forall n \in \N, \end{cases} pour $h>0$ fixé.
Par récurrence on montre que si $u_0>0$ alors $u_n>0$ pour tout $n$ $\in$ $\N$. De plus, $\frac{u_{n+1}}{u_{n}}=\frac{1}{1 + u_n^4 h}<1$ pour tout $n$ $\in$ $\N$: la suite est monotone décroissante. On a alors une suite décroissante et bornée par zéro, donc elle converge. Soit $\ell$ la limite de cette suite, alors $\ell\ge0$ et $\ell= \frac{\ell}{1 + \ell^4 h}$ donc $\ell=0$. Ce schéma est donc inconditionnellement A-stable.
N=50
tt=linspace(0,10,N+1)
y0=1
exacte = lambda t : y0*(4*t*y0**4+1)**(-0.25)
yy=[exacte(t) for t in tt]
phi = lambda t,y : -y**5
def schema(phi,tt):
h=tt[1]-tt[0]
uu=[y0]
for i in range(len(tt)-1):
uu.append(uu[i]/(1+h*uu[i]**4))
return uu
plot(tt,yy,label=("Exacte"),lw=2)
plot(tt,schema(phi,tt),'o',label=("Schéma"))
legend();
Soit le problème de Cauchy: \begin{cases} y'(t)+\sin(y(t))=0, & \forall t \in \R,\\ y(0)=y_0>0. \end{cases}
C'est un problème deCauchy du type \begin{cases} y'(t)=\varphi(t,y(t)), & \forall t \in \R,\\ y(0)=y_0>0, \end{cases} avec $\varphi(t,y)=g(y)=-\sin(y)$.
On montre d'abord qu'il existe une et une seule solution locale (ie sur $[0;T]$) de classe $\mathscr{C}^1([0,T],\R)$.
On montre ensuite que cette solution est de classe $\mathscr{C}^\infty([0,T],\R)$.
On montre enfin que la solution admet un prolongement sur $\R$.
Soit $h>0$ fixé et $t_n=nh$ pour tout $n\in\Z$. Le schéma d'Euler explicite pour l'EDO donnée construit la suite $$ u_{n+1}=u_n - h \sin(u_n), \qquad \forall n \in \N. $$
Comme $|a+b|\le |a|+|b|$ et comme $|-\sin(x)|\le1$ pour tout $x\in\R$, on conclut que $$ |u_{n+1}|=|u_n - h \sin(u_n)| \le |u_n| + |h \sin(u_n)| \le |u_n| + h $$ pour $h>0$ fixé.
Par récurrence: $|u_{n+1}|\le |u_n| + h \le |u_{n-1}| + 2h \le\dots\le |u_0| + (n+1)h$.
N=100
tt=linspace(0,10,N+1)
phi = lambda t,y : -sin(y)
def schema(phi,tt):
h=tt[1]-tt[0]
uu=[y0]
for i in range(len(tt)-1):
uu.append(uu[i]-h*sin(uu[i]))
return uu
plot(tt,schema(phi,tt),'o',label=("Schéma"))
legend();
Un modèle pour la diffusion d'une épidémie se base sur l'hypothèse que sa vitesse de propagation est proportionnelle au nombre d'individus infectés et au nombre d'individus sains.
Si on note $I(t)\ge0$ le nombre d'individus infectés à l'instant $t\ge0$ et $A>0$ le nombre d'individus total, il existe une constante $k$ $\in$ $\R^+$ telle que $I'(t)=k I(t)(A-I(t))$.
C'est un problème deCauchy du type \begin{cases} I'(t)=\varphi(t,I(t)), & \forall t \in \R^+,\\ I(0)=I_0>0, \end{cases} avec $\varphi(t,I(t))=g(I(t))=kI(t)(A-I(t))$.
Soit $0<I_0<A$. On considère le schéma semi-implicite $$ \frac{I_{n+1}-I_n}{h}=kI_n(A-I_{n+1}) $$ pour $h>0$ fixé. On obtient une formule de récurrence rendue explicite par un calcul élémentaire: $$ I_{n+1}=\frac{1+kAh}{1+kI_nh}I_n. $$ Si $0<I_0<A$ alors
donc $I_n\to A$ lorsque $n\to+\infty$ indépendamment du pas $h>0$ choisi.
On a déjà observé qu'il y a deux solutions constantes de l'EDO: la fonction $I(t)\equiv0$ et la fonction $I(t)\equiv A$.
Pour chercher toutes les solutions non constantes on remarque qu'il s'agit d'une EDO à variables séparables donc on a \begin{align*} % I'(t)=k I(t)(5000-I(t))\\ % \frac{I'(t)}{I(t)(5000-I(t))}&=k\\ % \frac{\mathrm{d}\,I}{I(5000-I)}&=k\mathrm{d}\,t\\ % \int\frac{1}{I(5000-I)}\mathrm{d}\,I&=k\int\mathrm{d}\,t\\ % \int\frac{1}{I}\mathrm{d}\,I-\int\frac{1}{5000-I}\mathrm{d}\,I&=5000k\int\mathrm{d}\,t\\ % \ln(I)+\ln(5000-I)&=5000kt+c\\ % \ln\frac{I}{5000-I}&=5000kt+c\\ % \frac{I}{5000-I}&=De^{5000kt}\\ % I(t)&=\frac{5000De^{5000 k t}}{1+De^{5000 k t}}\\ I(t)&=\frac{A}{De^{-A k t}+1} \end{align*}
La valeur numérique de la constante d'intégration $D$ est obtenue grâce à la CI: \begin{align*} % I(0)&=\frac{A}{De^{0}+1}\\ D&=\frac{A-I_0}{I_0} \end{align*}
Exemple avec $A=5000$, $I_0=160$, $k=\frac{\ln(363/38)}{35000}$ et $h=1$:
A=5000
I0=160
k=log(363/38)/35000
D=(A-I0)/I0
exacte = lambda t : A/(D*exp(-A*k*t)+1)
h=1
tt=arange(0,30,h)
yy=[exacte(t) for t in tt]
phi = lambda t,I : k*I*(A-I)
def schema(phi,tt):
uu=[I0]
for i in range(len(tt)-1):
uu.append((1+k*A*h)/(1+k*uu[-1]*h)*uu[-1])
return uu
plot(tt,yy,label=("Exacte"))
plot(tt,schema(phi,tt),label=("Schéma"))
legend();
Soit le problème de Cauchy: \begin{cases} y''(t)+1001y'(t)+1000y(t)=0, & \forall t \in \R,\\ y(0)=1,\\ y'(0)=0. \end{cases}
Soit $\mathbf{z}(t)=(y(t),y'(t))^T$ alors $y''(t)+1001y'(t)+1000y(t)=0$ se réécrit $$ \mathbf{z}'(t)= -\begin{pmatrix} 0&-1\\ 1000&1001 \end{pmatrix} \mathbf{z}(t) $$ Cette matrice a pour valeurs propres $\lambda_1=1$ et $\lambda_2=1000$ et pour vecteurs propres associés $\mathbf{v}_1=(-1,1)^T$ et $\mathbf{v}_2=(0,1000)^T$ (ci dessous le calcul formel des valeurs et vecteurs propes). Ainsi $$ \mathbf{z}(t)=c_1\mathbf{v}_1e^{-\lambda_1t}+c_2\mathbf{v}_2e^{-\lambda_2t} $$ et $$ y(t)=-c_1e^{-t}-c_2e^{-1001t}. $$ Il s'agit d'un problème stiff!
%reset -f
import sympy as symb
symb.init_printing()
symb.var('t')
M = symb.Matrix(((0,-1), (1000,1001)))
M
p=M.charpoly(t)
symb.factor(p)
P, D = M.diagonalize()
P,D
symb.simplify(P*D*P**-1)
Soit $h>0$ fixé et $t_n=nh$ pour tout $n\in\Z$. Le schéma d'Euler explicite pour le système donné construit la suite $$ \begin{align} \mathbf{u}_{n+1} &= \mathbf{u}_n -h \begin{pmatrix} 0&-1\\ 1000&1001 \end{pmatrix} \mathbf{u}_n \\ &= \begin{pmatrix} 1 & h\\ -1000h & 1-1001h \end{pmatrix} \mathbf{u}_n \\ &= \begin{pmatrix} 1 & h\\ -1000h & 1-1001h \end{pmatrix}^{n+1} \mathbf{u}_0,\\ &= \begin{pmatrix} -1 & -1\\ 1000 & 1 \end{pmatrix} \begin{pmatrix} 1-1000h & 0\\ 0 & 1-h \end{pmatrix}^{n+1} \begin{pmatrix} -1 & -1\\ 1000 & 1 \end{pmatrix}^{-1} \mathbf{u}_0\\ &= \begin{pmatrix} -1 & -1\\ 1000 & 1 \end{pmatrix} \begin{pmatrix} (1-1000h)^{n+1} & 0\\ 0 & (1-h)^{n+1} \end{pmatrix} \begin{pmatrix} -1 & -1\\ 1000 & 1 \end{pmatrix}^{-1} \mathbf{u}_0, \qquad \forall n \in \N. \end{align} $$
symb.var('h')
M = symb.Matrix(((1,h), (-1000*h,1-1001*h)))
M
P, D = M.diagonalize()
P , D
Le schéma est A-stable si $\lim_{n\to+\infty}\mathbf{u}_{n+1}=\mathbf{0}$. Il faut donc que \begin{cases} |1-1000h|<1\\ |1-h|<1 \end{cases} d'où la condition $h<\frac{2}{1000}$.
Pour calculer la solution approchée considérons la $\vartheta$-méthode suivante: \begin{cases} u_0=y_0,\\ u_{n+1}=u_{n}+h\left(\vartheta\varphi(t_{n+1},u_{n+1})+(1-\vartheta)\varphi(t_{n},u_{n})\right). \end{cases} avec $\vartheta \in [0;1]$.
Soit $\varphi(t,y(t))=-\beta y(t)$ avec $\beta>0$. Donner une condition sur $h$ pour que la $\vartheta$-méthode soit A-stable. Pour quelles valeurs de $\vartheta$ la méthode est inconditionnellement A-stable?
Remarque: pour $\vartheta=0$ on retrouve la méthode d’Euler explicite, pour $\vartheta=1$ la méthode d’Euler implicite, pour $\vartheta=\frac{1}{2}$ la méthode de Crank-Nicholson.
Par induction $$ u_n=\left(1-\frac{\beta h }{1+\beta h \vartheta}\right)^n y_0 $$ Il s'agit d'une suite géométrique de raison $q=1-\frac{\beta h }{1+\beta h \vartheta}$.
La suite tends à zéro ssi $|q|<1$ (la convergence est monotone ssi $0\le q<1$):
$$
-1<1-\frac{\beta h }{1+\beta h \vartheta}<1
\iff
0<\frac{\beta h }{1+\beta h \vartheta}<2
\iff
(1-2\vartheta)\beta h < 2
$$
Pour $\vartheta\ge \frac{1}{2}$ la méthode est inconditionnellement A-stable,
pour $\vartheta< \frac{1}{2}$ la méthode est A-stable ssi $h<\frac{2}{(1-2\vartheta)\beta}$.