from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())

import sys #only needed to determine Python version number
print('Python version ' + sys.version)

Python version 3.8.10 (default, Nov 26 2021, 20:14:08) 
[GCC 9.3.0]

M62_CM3 : La boîte de Pandore des approximations : problèmes bien posés (mathématiquement et numériquement), schémas A-stables¶

En général, il ne suffit pas qu'un schéma numérique soit convergent pour qu'il donne de bons résultats sur n'importe quelle équation différentielle. Dès qu'on tente d’approcher une solution, de nombreux facteurs peuvent en fait nous en éloigner: une étude mathématique sophistiquée doit alors être envisagée. En particulier, il faut que le problème soit

mathématiquement bien posé (existence et unicité de la solution),
numériquement bien posé (continuité suffisamment bonne par rapport aux conditions initiales)
bien conditionné = non raide (temps de calcul raisonnable)

Voici quelques exemples apparemment anodins tirés de la littérature.

Problème mal posé mathématiquement
Que pensez-vous du problème suivant? $$ \begin{cases} y'(t)=2\sqrt{|y(t)|}\\ y(0)=0 \end{cases} $$
Problème mal posé numériquement
Considérons le problème: $$ \begin{cases} y'(t)=5y(t)-5, &t\in[0;50]\\ y(0)=1 \end{cases} $$ et imaginons une petite perturbation de la valeur initiale: $y(0)=1+\varepsilon$. Estimez $y_\varepsilon(50)-y(50)$ pour $\varepsilon=10^{-17}$ (de l’ordre de l’epsilon de Python).
Problème mal conditionné
Cette fois, le problème est le suivant: $$ \begin{cases} y'(t)=-150y(t)-30\\ y(0)=\frac{1}{5}. \end{cases} $$ Vérifier qu’il est bien posé mathématiquement et numériquement. Si on l'approche à l’aide de la méthode d’Euler explicite on obtient $$ u_n=\frac{1}{5}+\left(1-\frac{150}{N}\right)^n\left(u_0-\frac{1}{5}\right). $$ Qu'obtenons-nous pour $N=50$ et $u_0=\frac{1}{5}+\varepsilon$?
Problèmes raides
Il existe des problèmes qui résistent à toutes les méthodes usuelles. On trouve par exemple ces phénomènes en chimie où deux réactions ont des échelles de temps très différentes. Considérons par exemple $$ \begin{cases} y'(t)=-1000\big(y(t)-e^{-t}\big)-e^{-t}, & t\in[0;1]\\ y(0)=0. \end{cases} $$

1 Propriétés du problème de Cauchy
2 Problème mathématiquement bien posé
- 2.1 Exemple de problème mathématiquement mal posé
- 2.2 Exercice
3 Problème numériquement bien posé
4 Problème bien conditionné
- 4.1 Exemple de problème mal conditionné (= raide = stiff)
- 4.2 Exemple de problème raide et approximation numérique
5 Propriétés d'un schéma d'approximation: consistance, stabilité, convergence
6 Étude de la zéro-stabilité pour les schémas à un pas
- 6.1 Exemple: Zéro-stabilité de la méthode d'Euler explicite
7 Étude de la A-stabilité pour les schémas à un pas
8 Tableau recapitulatif
9 Exemple : A-stabilité des schémas classiques à un pas

Propriétés du problème de Cauchy¶

Problème mathématiquement bien posé
Un problème de Cauchy est un problème (mathématiquement) bien posé si sa solution existe et est unique.
Problème numériquement bien posé
Un problème de Cauchy est numériquement bien posé si sa solution n’est pas trop perturbée par une erreur initiale ou des perturbations du second membre.
Problème bien conditionné
Un problème est bien conditionné si les méthodes numériques usuelles peuvent donner sa solution en un nombre raisonnable d’opérations. Sinon on parle de problème raide.

Cf. page 235-236 Jean-Pierre DEMAILLY, Analyse Numérique et équations différentielles

Problème mathématiquement bien posé¶

On dit qu'un problème de Cauchy est mathématiquement bien posé s'il existe une et une seule solution.

Si ce n'est pas le cas, on dit qu'il est mathématiquement mal posé.

Exemple de problème mathématiquement mal posé¶

Par exemple, on cherche une fonction $y\colon t\in\mathbb{R}^+\mapsto y(t)\in\mathbb{R}$ qui satisfait $$ \begin{cases} y'(t) = \sqrt[3]{y(t)}, &\forall t>0,\\ y(0) = 0. \end{cases} $$ Ce problème de Cauchy ne satisfait pas les hypothèse du théorème de Cauchy-Lipschitz.

En effet, $\varphi(t,y)=\sqrt[3]{y}$. Elle n'est pas lipschitizienne par rapport à la deuxième variable en $y=0$ car $\partial_y \varphi(t,y)=\frac{1}{3\sqrt[3]{y^2}} \xrightarrow[y\to0]{}+\infty$.

Il est simple de vérifier que, pour tout $t\ge0$, les trois fonctions

$y_1(t)=0$ pour tout $t$ et
$y_{2,3}(t)=\pm\sqrt{\frac{8}{27}t^3}=\pm\frac{2\sqrt{2}}{3\sqrt{3}}\sqrt{t^3}=\pm\frac{2\sqrt{6}}{9}\sqrt{t^3}$

sont 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.

# On demande à sympy de calculer les solutions exactes : il ne trouve pas la solution constante = 0
import sympy as sym
sym.init_printing()

t=sym.Symbol('t')
y=sym.Function('y')

edo=sym.Eq( sym.diff(y(t),t) , y(t)**(1/sym.S(3))  )
display(edo)

solgenLIST=sym.dsolve(edo,y(t))
display(solgenLIST)

for solgen in solgenLIST:
    t0,y0=0,0
    consts= sym.solve( [solgen.rhs.subs(t,t0)-y0 ], dict=True)[0]
    solpar=solgen.subs(consts).simplify()
    display(solpar)

Exercice¶

Montrer que le problème de Cauchy $$ \begin{cases} y'(t)=2\sqrt{|y(t)|}, & t>0\\ y(0)=0, \end{cases} $$ admet une infinité de solutions de classe $\mathcal{C}^1(\mathbb{R})$.

Parmi ces solutions, quelle solution approche-t-on si on utilise la méthode d'Euler explicite? et la méthode d'Euler implicite?

Que se passe-t-il si la donnée initiale est $y(0)=y_0>0$?

Correction
Sur $\mathbb{R}^+$ la fonction $\varphi(t,y)=2\sqrt{y}$ n'est pas lipschitzienne par rapport à $y$ au voisinage de $(t,0)$ (car $\partial_y\varphi=1/\sqrt{y}\to\infty$ lorsque $y\to0$), donc le théorème d'existence et unicité locale n'est pas valable au voisinage de $(0,0)$. Même raisonnement sur $\mathbb{R}^-$.

L'EDO est à variables séparables, on peut donc calculer explicitement toutes les solutions du problème de Cauchy:

elle admet une seule solution constante, la fonction $y(t)=0$ pour tout $t\in\mathbb{R}^+$,
et des solutions de la forme $y(t)=(t+c)^{2}$ pour tout $t\ge c$.

En imposant la CI on trouve que, pour tout $b\in\mathbb{R}^+$, les fonctions $$ y_b(t)= \begin{cases} 0 , &\text{si }0\le t\le b, \\ (t-b)^{2} , &\text{si }t\ge b, \end{cases} $$ sont de classe $\mathcal{C}^1(\mathbb{R}^+)$ et sont solution du problème de Cauchy donné.

Traçons quelques courbes:

%matplotlib inline
from matplotlib.pylab import *

xx =linspace(0,2,101)
yyr=[0 for x in xx]

f = lambda x,b : 0 if (x<=b) else (x-b)**2

yy1 =[f(x,1)    for x in xx]
yy75=[f(x,0.75) for x in xx]
yy0 =[f(x,0)    for x in xx]

figure(figsize=(10,7))
plot(xx,yyr,'r-',xx,yy1,'b:',xx,yy75,'g*',xx,yy0,'yd');
legend([r'$y(t)=0$ $\forall$ $t$', '$y_{b=1}(t)$', '$y_{b=0.75}(t)$', '$y_{b=0}(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)=0$ pour tout $t\in\mathbb{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-sûr $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\mathbb{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.

Problème numériquement bien posé¶

Une fois calculée la solution numérique $\{u_n\}_{n=1}^{N}$ d'un problème de Cauchy mathématiquement bien posé, il est légitime de chercher à savoir dans quelle mesure l'erreur $|y(t_n)-u_n|$ est petite. Cela dépend bien sur du schéma choisi mais aussi du problème en question.

Étant donné que les erreurs d’arrondi amènent toujours à résoudre un problème perturbé, il est important de savoir si la solution d'un problème perturbé est proche de la solution du problème non perturbé.

On dit qu'un problème de Cauchy est numériquement bien posé si la solution d'un problème faiblement perturbé (second membre ou condition initiale) possède une solution proche de celle du problème original.

Exemple de problème numériquement mal posé¶

On se donne $\varphi(t,y)=3t-3y$ et $y(0)=\alpha$ (un nombre quelconque). On cherche une fonction $y\colon t\in\mathbb{R}\mapsto y(t)\in\mathbb{R}$ qui satisfait $$ \begin{cases} y'(t) = 3y(t)-3t, &\forall t\in\mathbb{R},\\ y(0) = \alpha. \end{cases} $$ Sa solution, définie sur $\mathbb{R}$, est donnée par $$ y(t)=\left(\alpha-\frac{1}{3}\right)e^{3t}+t+\frac{1}{3}. $$ Calculons $y$ en $t=10$:

si $\alpha=\frac{1}{3}$ alors $y(10)=\frac{31}{3}$,
si $\alpha=0.333333$ alors $y(10)=(0.333333-1/3)e^{30}+10+1/3=-e^{30}/3000000+31/3\approx31/3+10^7/3$

Vérifions-le:

import sympy as sym
sym.init_printing()

t=sym.Symbol('t')
y=sym.Function('y')

edo=sym.Eq( sym.diff(y(t),t) , 3*y(t)-3*t  )
#display(edo)

solgen=sym.dsolve(edo,y(t))
#display(solgen)

alpha=sym.Symbol('alpha')
t0,y0=0,alpha
consts= sym.solve( [solgen.rhs.subs(t,t0)-y0 ], dict=True)[0]
solpar=solgen.subs(consts).simplify()
display(solpar)

sol1=solpar.subs(t,10).subs(alpha,sym.Rational(1,3))
sol2=solpar.subs(t,10).subs(alpha,0.333333)
display(sol1)
display(sol2)
print('Difference:')
display((sol1.rhs-sol2.rhs).evalf())

Difference:

Si nous cherchons à résoudre le problème de Cauchy jusqu'à $t=10$ avec $\alpha=1/3$, nous obtenons $y(10)=31/3$.
Par contre, si nous faisons le calcul avec l'approximation $\alpha=0.333333$ au lieu de $1/3$, nous avons $y(10)=31/3-e^{30}/3000000$ 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 calcule 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é

Exemple de problème numériquement bien posé¶

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 *

y0 = 1

y = lambda t,eps : (y0+eps)*exp(-t)

tt  = linspace(0,5,101)
yye = [y(t,0.1) for t in tt]
yy0 = [y(t,0) for t in tt]
#xkcd()
figure(figsize=(10,7))
plot(tt,yye,tt,yy0);
legend(["$y_0=1.1$","$y_0=1$"])
axis([0,5,0,1.1]);

Exemple de problème numériquement mal posé¶

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 *

y0 = 0

y = lambda t,eps : (y0+eps)*exp(t)

tt  = linspace(0,5,101)
yye = [y(t,0.01) for t in tt]
yy0 = [y(t,0) for t in tt]
#xkcd()
figure(figsize=(10,7)) 
plot(tt,yye,lw='2')
plot(tt,yy0,lw='5');
legend(["$y_0=0.01$","$y_0=0$"])
axis([0,5,0,1.1]);

Exemple de problème numériquement bien posé¶

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 *

y0 = -1

y = lambda t,eps : (y0+eps)*exp(1-cos(t))

tt=linspace(0,10,101)
yye = [y(t,0.1) for t in tt]
yy0 = [y(t,0) for t in tt]
#xkcd()
figure(figsize=(10,7)) 
plot(tt,yye,tt,yy0);
legend(["$y_0=-0.9$","$y_0=-1$"]);

Exercice¶

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}$$

Calculer l'unique solution exacte à ce problème de Cauchy en fonction de $a$ pour $t>0$. Remarquer qu'avec le changement de variable $u(t)=\frac{y(t)}{t}$ on obtient une EDO linéaire d'ordre 1

Tracer la solution exacte sur l'intervalle $[1;5]$ pour $a=1$, $a=1-\frac{1}{1000}$ et $a=1-\frac{1}{100}$.

Pour $a=1$, sur un autre graphique tracer la solution exacte et les solutions obtenues par la méthode d'Euler explicite avec $N=20$, $N=100$, $N=300$, $N=5000$. Comment expliquez-vous ce comportement du schéma?

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

$A(t)=\int \frac{b(t)}{a(t)}\mathrm{d}t=\int -2/t\mathrm{d}t= -2\ln(t)=\ln(t^{-2})$,
$B(t)=\int \frac{g(t)}{a(t)}e^{A(t)}\mathrm{d}t=\int -5t^{-4} e^{A(t)}\mathrm{d}t=\int -5t^{-6} \mathrm{d}t = t^{-5}$.

On conclut que $u(t)=ce^{-A(t)}+B(t)e^{-A(t)}=ct^2+t^{-3}$ et $y(t)=tu(t)=ct^3+t^{-2}$.

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}. $$

Vérifions ce calcul avec sympy:

import sympy as sym
sym.init_printing()

t=sym.Symbol('t')
y=sym.Function('y')

edo=sym.Eq( sym.diff(y(t),t) ,  3*y(t)/t-5/t**3  )
display(edo)

solgen=sym.dsolve(edo,y(t))
display(solgen)

a=sym.Symbol('a')
t0,y0=1,a
consts= sym.solve( [solgen.rhs.subs(t,t0)-y0 ], dict=True)[0]
solpar=solgen.subs(consts)
display(solpar)

%matplotlib inline
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]
figure(figsize=(10,7))
for y0 in Y0:
    yy=[sol_exacte(t,y0) for t in tt]
    plot(tt,yy,label=('$y_0=$'+str(y0)))
legend()
grid()
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éveille" le terme dormant $t^{3}$.

Lorsque l'on essaye d'approcher la solution avec les méthodes d'Euler explicite et implicite, bien sûr lorsque $N\to+\infty$ on converge vers la solution exacte, cependant en pratique $N$ est fini et on voit qu'inévitablement nous allons approcher une autre solution que celle correspondant à $a=1$.

from scipy.optimize import fsolve

phi = lambda t,y : 3*y/t-5/t**3
    
t0     = 1
tfinal = 5
y0     = 1


def EE(phi,tt,y0):
    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 EI(phi,tt,y0):
    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

def EM(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = [y0]
    for i in range(len(tt)-1):
        utilde=uu[i]+h/2*phi(tt[i],uu[i])
        uu.append(uu[i]+h*phi(tt[i]+h/2,utilde))
    return uu

def CN(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = [y0]
    for i in range(len(tt)-1):
        uu.append(fsolve(lambda x : -x +uu[i]+h/2*(phi(tt[i],uu[i])+phi(tt[i+1],x)),uu[i]))
    return uu

def HE(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = [y0]
    for i in range(len(tt)-1):
        utilde=uu[i]+h*phi(tt[i],uu[i])
        uu.append(uu[i]+h/2*(phi(tt[i],uu[i])+phi(tt[i+1],utilde)))
    return uu

NN=[20,100,300,5000]

figure(1, figsize=(18, 12))

subplot(2,3,1)
for N in NN:
	tt=linspace(t0,tfinal,N)
	uu=EE(phi,tt,y0)
	plot(tt,uu,'.-')
plot(tt,[sol_exacte(t,y0) for t in tt],'c-',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(2,3,3)
for N in NN:
	tt=linspace(t0,tfinal,N)
	uu=EI(phi,tt,y0)
	plot(tt,uu,'.-')
plot(tt,[sol_exacte(t,y0) for t in tt],'c-',lw=2)
legend(['N=20','N=100','N=300','N=5000','Exacte'])
title('$y_0=1$ - Exacte vs EI')
axis([t0,tfinal,-c,c]);

subplot(2,3,4)
for N in NN:
	tt=linspace(t0,tfinal,N)
	uu=EM(phi,tt,y0)
	plot(tt,uu,'.-')
plot(tt,[sol_exacte(t,y0) for t in tt],'c-',lw=2)
legend(['N=20','N=100','N=300','N=5000','Exacte'])
title('$y_0=1$ - Exacte vs EM')
axis([t0,tfinal,-c,c]);

subplot(2,3,5)
for N in NN:
	tt=linspace(t0,tfinal,N)
	uu=HE(phi,tt,y0)
	plot(tt,uu,'.-')
plot(tt,[sol_exacte(t,y0) for t in tt],'c-',lw=2)
legend(['N=20','N=100','N=300','N=5000','Exacte'])
title('$y_0=1$ - Exacte vs HE')
axis([t0,tfinal,-c,c]);

subplot(2,3,6)
for N in NN:
	tt=linspace(t0,tfinal,N)
	uu=CN(phi,tt,y0)
	plot(tt,uu,'.-')
plot(tt,[sol_exacte(t,y0) for t in tt],'c-',lw=2)
legend(['N=20','N=100','N=300','N=5000','Exacte'])
title('$y_0=1$ - Exacte vs CN')
axis([t0,tfinal,-c,c]);

Problème bien conditionné¶

On dit qu'un problème de Cauchy est bien conditionné si les méthodes numériques usuelles peuvent donner sa solution en un nombre raisonnable d’opérations.

Dans le cas contraire on parle de problème raide.

Exemple de problème mal conditionné (= raide = stiff)¶

On se donne $\varphi(t,y)=-\beta y$ et $y(0)=1$. On cherche une fonction $y\colon t\in\mathbb{R}\mapsto y(t)\in\mathbb{R}$ qui satisfait $$ \begin{cases} y'(t) = -\beta y(t), &\forall t\in\mathbb{R},\\ y(0) = 1. \end{cases} $$ Ce problème est

mathématiquement bien posé pour tout $\beta\in\mathbb{R}$ : il existe une et une seule solution, elle définie sur $\mathbb{R}$ et est donnée par $y(t)=e^{-\beta t}$;
numériquement bien posé pour tout $\beta>0$ : l'erreur sur la solution sera au pire de l'ordre de l'erreur sur la condition initiale;
mal conditionné ou stiff: si nous cherchons à résoudre le problème de Cauchy lorsque $\beta$ est très grand, il faut prendre un pas $h$ très petit et ce quel qu'il soit le schéma choisi car la solution devient de plus en plus raide (stiff).

Voici ci-dessous un exemple avec des valeurs de $\beta$ de $1$ à $100$:

%reset -f
%matplotlib inline
from matplotlib.pylab import *

exacte = lambda t,b : exp(-b*t) 

t0 = 0
tfinal = 1
tt=linspace(t0,tfinal,501)

figure(1, figsize=(15, 3))
for i,b in enumerate([1,10,50,100]):
    subplot(1,4,i+1)
    yy=[exacte(t,b) for t in tt]
    plot(tt,yy,label=('$b=$'+str(b)));
    axis([0,1,0,1])
    legend()
    grid()

Exemple de problème raide et approximation numérique¶

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}$$

Solution exacte

Donner l'unique solution exacte à ce problème de Cauchy et tracer le graphe $t\mapsto y(t)$ pour $t\ge0$ en fonction de $y_0$.

Soient $b=g=10$. Tracer avec python les solutions exactes sur l'intervalle $[0;1]$ pour $y_0=\frac{g}{b}$ et $y(0)=\frac{g}{b}+10^{-8}$ pour $b=g=10$.

Sur un graphe à coté tracer les mêmes courbes avec $b=g=40$. Que peut-on dire lorsque $b$ devient de plus en plus grand?

Méthode d'Euler explicite

Vérifier que la méthode génère une suite arithmético-géométrique et calculer analytiquement pour quelles valeurs de $h$ la suite ainsi construite converge et pour quelles valeurs la convergence est monotone.

Soient $b=g=40$ et $y_0=\frac{g}{b}$. Tracer les solutions approchées obtenues avec $N=19$, $N=23$, $N=45$, $N=100$.

Répéter pour $y_0=\frac{g}{b}+10^{-8}$ et expliquer chaque résultat.

Méthode d'Euler implicite

Répéter le même exercice pour la méthode d'Euler implicite. Noter que, étant $\varphi$ linéaire, la méthode peut être rendu explicite par un calcul élémentaire.

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\mathbb{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.

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

La méthode d'Euler explicite donne la suite définie par récurrence $u_{n+1}=Au_n+B$ et $u_0=y_0$ avec $A=1-bh$ et $B=gh$. Pour calculer $u_n$ sans passer par la récurrence, on introduit une suite auxiliaire définie par $v_n=u_n-\alpha$ ainsi $u_n=v_n+\alpha$ pour tout $n\in\mathbb{N}$. On a $v_{n+1}+\alpha=A(v_n+\alpha)+B=Av_n+A\alpha+B$. En posant $\alpha=\frac{B}{1-A}\ (A\neq1)$, la suite $(v_n)_{n\in \mathbb{N}}$ est géométrique de raison $A$ donc $v_{n}=v_0A^n$. En remplaçant $v_n$ par $u_n-\alpha$ et $v_0$ par $u_0-\alpha$ on en déduit que $u_n=A^n(u_0-\alpha)+\alpha,\quad \alpha=\frac{B}{1-A}$ soit encore $$ u_n-\frac{g}{b}=(1-hb)^n\left( y_0 -\frac{g}{b}\right) $$
Comparons la solution exacte et la solution approchée: \begin{align*} y_n&=\frac{g}{b}+\left( y_0 -\frac{g}{b}\right)e^{-bt_n}\\ u_n&=\frac{g}{b}+\left( y_0 -\frac{g}{b}\right)(1-hb)^n \end{align*} Pour que $u_n\simeq y_n$ il faut que $(1-hb)^n\simeq e^{-bt_n} \in]0,1]$, donc il faut que $0\le(1-hb)<1$, soit encore $h<\frac{1}{b}$. Comme $h=\frac{T-t-0}{N}=\frac{1}{N}$, cela implique $N\ge b$. Dans cet exemple, pour obtenir une bonne convergence il est nécessaire de prendre un pas $h$ de pus en plus petit que $b$ est grand, c'est-à-dire un coût en calcul plus élevé. On parle de problème raide ou mal conditionné.

%reset -f
%matplotlib inline
from matplotlib.pylab import *

phi = lambda t,y,b,g : -b*y+g

def EE(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(figsize=(18, 7))
subplot(1,2,1)
y0=Y0[0]
for N in NN:
    tt=linspace(t0,tfinal,N+1)
    yy=EE(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=EE(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

La méthode d'Euler implicite donne la suite définie par récurrence $u_{n+1}=Au_n+B$ et $u_0=y_0$ avec $A=\frac{1}{1+bh}$ et $B=\frac{gh}{1+bh}$. On introduit une suite auxiliaire définie par $v_n=u_n-\alpha$ ainsi $u_n=v_n+\alpha$ pour tout $n\in\mathbb{N}$. On a $v_{n+1}+\alpha=A(v_n+\alpha)+B=Av_n+A\alpha+B$. En posant $\alpha=\frac{B}{1-A}\ (A\neq1)$, la suite $(v_n)_{n\in \mathbb{N}}$ est géométrique de raison $A$ donc $v_{n}=v_0A^n$. En remplaçant $v_n$ par $u_n-\alpha$ et $v_0$ par $u_0-\alpha$ on en déduit que $u_n=A^n(u_0-\alpha)+\alpha,\quad \alpha=\frac{B}{1-A}$ soit encore $$ u_n-\frac{g}{b}=\frac{1}{(1+hb)^n}\left( y_0 -\frac{g}{b}\right) $$
Comparons la solution exacte et la solution approchée: \begin{align*} y_n&=\frac{g}{b}+\left( y_0 -\frac{g}{b}\right)e^{-bt_n}\\ u_n&=\frac{g}{b}+\left( y_0 -\frac{g}{b}\right)\frac{1}{(1+hb)^n} \end{align*} Pour que $u_n\simeq y_n$ il faut que $\frac{1}{(1+hb)^n}\simeq e^{-bt_n} \in]0,1]$, donc il faut que $0\le\frac{1}{1+hb}<1$ c'est-à-dire pour n'importe quelle valeur de $h$. Dans cet exemple, pour obtenir la convergence il n'est pas nécessaire de limiter le pas $h$. Cependant, pour une edo quelconque, à chaque pas de temps on doit résoudre une équation non-linéaire.

%reset -f
%matplotlib inline
from matplotlib.pylab import *
from scipy.optimize import fsolve

phi = lambda t,y,b,g : -b*y+g

def EI(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(figsize=(18,7))
subplot(1,2,1)
y0=Y0[0]
for N in NN:
    tt=linspace(t0,tfinal,N+1)
    yy=EI(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=EI(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()

Propriétés d'un schéma d'approximation: consistance, stabilité, convergence¶

Le théorème d’équivalence de Lax-Ritchmyer affirme qu'un schéma consistant est convergent ssi il est zéro-stable:

Consistance + zéro-stabilité = convergence

Rappelons ci-dessous la définition de schéma convergent et de schéma consistant, avant d'introduire la (les) notion(s) de stabilité d'un schéma.

Convergence¶

Une méthode numérique est convergente si $$ |y_n-u_n|\le C(h)\xrightarrow[h\to0]{}0 \qquad\forall n=0,\dots,N $$ Si $C(h) = \mathcal{O}(h^p)$ pour $p > 0$, on dit que la convergence de la méthode est d'ordre $p$.

Remarque:

$N\to+\infty$ lorsque $h\to0$ pour $T$ fixé (notion liée à la zéro-stabilité)
$N\to+\infty$ lorsque $T\to+\infty$ pour $h$ fixé (notion liée à la A-stabilité).

Soit $u_{n+1}^*$ la solution numérique au temps $t_{n+1}$ qu'on obtiendrait en insérant de la solution exacte dans le schéma. Pour vérifier qu'une méthode converge, on écrit l'erreur ainsi \begin{equation}\label{quart7.9} e_n \equiv y_n - u_n = (y_n - u_n^*) + (u_n^* - u_n). \end{equation}

Le terme $y_n - u_n^*$ représente l'erreur engendrée par une seule itération du schéma.

Le terme $u_{n+1}^* - u_{n+1}$ représente la propagation de $t_{n}$ à $t_{n+1}$ de l'erreur accumulée au temps précédent $t_{n}$.

Si les deux termes $(y_n - u_n^*)$ et $(u_n^* - u_n)$ tendent vers zéro quand $h \to 0$ alors la méthode converge.

Consistance¶

On appelle erreur de troncature locale la quantité $$ \tau_{n+1}(h)\equiv\frac{y_{n+1}-u_{n+1}^*}{h} $$ et erreur de troncature $$ \tau(h)=\max_{n=0,\dots,N}|\tau_n(h)|. $$ Si $\lim_{h\to0} \tau (h) = 0$ on dit que la méthode est consistante. On dit qu'elle est consistante d'ordre $p$ si $\tau (h) = \mathcal{O}(h^p)$ pour un certain $p\ge1$.

Stabilité : zéro-stabilité et A-stabilité¶

De manière général, un schéma numérique est dit stable s’il permet de contrôler la solution quand on perturbe les données. Il existe de nombreuses notions de stabilité.

Considérons le problème de Cauchy:

trouver une fonction $y \colon I\subset \mathbb{R} \to \mathbb{R}$ définie sur un intervalle $I$ telle que $$ \begin{cases} y'(t) = \varphi(t,y(t)), &\forall t \in I=]t_0,T[,\\ y(t_0) = y_0, \end{cases} $$

avec $y_0$ une valeur donnée et supposons que l'on ait montré l'existence et l'unicité d'une solution $y$ pour $t\in I$ (donc il est mathématiquement bien posé) et qu'il soit numériquement bien posé.

Deux questions naturelles se posent:

que se passe-t-il lorsqu'on fixe le temps final $T$ et on fait tendre le pas $h$ vers $0$?
que se passe-t-il lorsqu'on fixe le pas $h>0$ mais on fait tendre $T$ vers l'infini?

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é:

Zéro-stabilité:
la propriété de zéro-stabilité assure que la méthode numérique est peu sensible aux petites perturbations des données (encore faut-il que le problème soit numériquement bien posé). L'exigence d'avoir une méthode numérique stable provient avant tout de la nécessité de contrôler les (inévitables) erreurs introduites par l'arithmétique finie des ordinateurs. En effet, si la méthode numérique n'était pas zéro-stable, les erreurs d'arrondi commises sur $y_0$ et sur le calcul de $\varphi(t_n, y_n)$ rendraient la solution calculée inutilisable.
A-stabilité:
la zéro-stabilité s'intéresse à la résolution du problème de Cauchy sur des intervalles bornés. Dans ce cadre, le nombre $N_h$ de sous-intervalles ne tend vers l'infini que quand $h$ tend vers zéro. Il existe cependant de nombreuses situations dans lesquelles le problème de Cauchy doit être intégré sur des intervalles en temps très grands ou même infini. Dans ce cas, même pour $h$ fixé, $N_h$ tend vers l'infini. On s'intéresse donc à des méthodes capables d'approcher la solution pour des intervalles en temps arbitrairement grands, même pour des pas de temps $h$ "assez grands" (si on peut choisir $h>0$ quelconque, on dira que la méthode est inconditionnellement A-stable, sinon on aura une condition de stabilité sur $h$ qui détérminera le temps de calcul). Considérons un problème de Cauchy dont la solution exacte vérifie la propriété $y(t)\xrightarrow[t\to+\infty]{}0$. Soit $h>0$ fixé et considérons la limite $T\to+\infty$ (ainsi $N_h\to+\infty$). On dit que la méthode est A-stable si $u_n^{(h)}\xrightarrow[n\to+\infty]{}0$.

Étude de la zéro-stabilité pour les schémas à un pas¶

La zéro-stabilité garantit que, sur un intervalle borné, des petites perturbations des données entraînent des perturbations bornées de la solution numérique quand $h\to0$.

Plus précisément, une méthode numérique pour approcher le problème de Cauchy sur l'intervalle fini $I=[t_0,T]$ è zéro-stable si $$ \exists h_0>0,\ \exists C > 0,\ \exists \varepsilon_0> 0 \text{ tel que } \forall h \in ]0,h_0],\ \forall \varepsilon \in ]0,\varepsilon_0],\ si |\varrho_n| \le \varepsilon,\ 0 \le n \le N, $$ alors $$ |z_n - u_n | \le C\varepsilon,\ 0 \le n \le N $$ où

$C$ est une costante qui ne dépend pas de $h$,
$u_n$ est une solution obtenue par une schéma,
$z_n$ est la solution obtenue par le même schéma sur le problème perturbé,
$\varrho_n$ est la perturbation au $n$-ième pas,
$\varepsilon_0$ est la perturbation maximale.

Naturellement, $\varepsilon_0$ doit être assez petit pour que le problème perturbé ait encore une unique solution sur l’intervalle d’intégration donné.

Pour une méthode consistante à un pas, on peut prouver que la zéro-stabilité est une conséquence du fait que $\varphi$ est lipschitzienne par rapport à sa deuxième variable. Dans ce cas, la constante $C$ qui apparaît dépend de $\exp((T − t_0 )L)$, où $L$ est la constante de Lipschitz. Cependant, ceci n’est pas toujours vrai pour les autres familles de méthodes.

Exemple: Zéro-stabilité de la méthode d'Euler explicite¶

Soit $u_n$ l'approximation obtenue avec la méthode d'Euler explicite du problème donné: $$ \begin{cases} u_0=y_0,\\ u_{n+1}=u_n+h\varphi(t_n,u_n),&n=0,\dots,N-1 \end{cases} $$ et $z_n$ l'approximation obtenue avec la méthode d'Euler explicite du problème perturbé: $$ \begin{cases} z_0=y_0+\varrho_0,\\ z_{n+1}=z_n+h\left(\varphi(t_n,z_n)+\varrho_n\right),&n=0,\dots,N-1 \end{cases} $$ On a \begin{align} |u_n - z_n| &\le |u_{n-1}-z_{n-1}| + h|\varphi(t_{n-1},u_{n-1})-\varphi(t_{n-1},z_{n-1})|+h|\varrho_n| \\ &\le (1+hL)|u_{n-1}-z_{n-1}| + h\varepsilon \\ &\le (1+hL)^2|u_{n-2}-z_{n-2}| + (1+(1+hL))h\varepsilon \\ &\dots \\ &\le (1+hL)^n|\varrho_0|+ \left( 1+(1+hL)+\dots+(1+hL)^{n-1} \right)h\varepsilon \\ &=(1+hL)^n|\varrho_0|+ \left(\sum_{i=0}^{n-1}(1+hL)^i\right)h\varepsilon \\ &=(1+hL)^n|\varrho_0|+ \frac{(1+hL)^n-1}{L}\varepsilon \\ &\le (e^{hL})^n|\varrho_0|+\frac{(e^{hL})^n-1}{hL}h\varepsilon \qquad\text{car }(1+x)\le e^x \\ &= e^{nhL}|\varrho_0| +\frac{(e^{nhL})-1}{L}\varepsilon \\ &\le \left(e^{nhL}+\frac{(e^{nhL})-1}{L}\right)\varepsilon \\ &= \frac{(1+L)e^{L(T-t_0)}-1}{L}\varepsilon \qquad\text{car } t_n-t_0=nh \end{align}

Étude de la A-stabilité pour les schémas à un pas¶

Dans la section précédente, on a considéré la résolution du problème de Cauchy sur des intervalles bornés. Dans ce cadre, le nombre $N$ de sous-intervalles ne tend vers l’infini que quand $h$ tend vers zéro. Il existe cependant de nombreuses situations dans lesquelles le problème de Cauchy doit être intégré sur des intervalles en temps très grands ou même infini. Dans ce cas, même pour $h$ fixé, $N$ tend vers l’infini, et un résultat comme la zéro-stabilité n’a plus de sens puisque le membre de droite contient une quantité non bornée ($T-t_0$). On s’intéresse donc à des méthodes capables d’approcher la solution pour des intervalles en temps arbitrairement grands, même pour des pas de temps $h$ "assez grands".

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\mathbb{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. $$

A-stabilité du schéma d'Euler explicite¶

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$ à $$ \boxed{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.
De plus, cette convergence est monotone ssi $1-\beta h>0$, i.e. ssi $h<\frac{1}{\beta}$.

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}}. $$

A-stabilité du schéma d'Euler implicite¶

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$.

A-stabilité du schéma d'Euler modifié¶

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 +\frac{1}{2}(\beta h)^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]
#xkcd()
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 $$ \boxed{ 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é.

Notons qu'on a la même condition de A-stabilité que le schéma d'Euler explicite. Cependant, $q(x)>0$ pour tout $x$, donc cette convergence est monotone (tandis qu'avec le schéma d'Euler explicite la convergence n'est monotone que si $h<\frac{1}{\beta}$).

A-stabilité du schéma de Cranck-Nicolson¶

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\mathbb{R}_+^*$, donc $|q(x)|<1$ pour tout $x\in\mathbb{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)$');

A-stabilité du schéma de Heun¶

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 $$ \boxed{ h <\frac{2}{\beta} } $$ et cette convergence est monotone.
Cette condition de stabilité limite le pas $h$ d'avance en $t$ lorsqu'on utilise le schéma de Heun.

A-stabilité du schéma de Simpson implicite¶

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]
#xkcd()
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\mathbb{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$.

A-stabilité du schéma de Simpson explicite¶

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]
#xkcd()
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$.

Tableau recapitulatif¶

		Type	Ordre	Conditions pour A-stabilité	A-stable et suite positive
EE	$u_{n+1}=u_n+h\varphi(t_n,u_n)$	Explicite	1	$0<h<\frac{2}{\beta}$	$0<h<\frac{1}{\beta}$
EI	$u_{n+1}=u_n+h\varphi(t_{n+1},u_{n+1})$	Implicite	1	$\forall h>0$	$\forall h>0$
EM	$u_{n+1}=u_n+h\varphi\left(t_{n}+\frac{h}{2},u_n+\frac{h}{2}\varphi(t_n,u_n)\right)$	Explicite	2	$0<h<\frac{2}{\beta}$	$0<h<\frac{2}{\beta}$
CN	$u_{n+1}=u_n+\frac{h}{2}\left(\varphi(t_n,u_n)+\varphi(t_{n+1},u_{n+1})\right)$	Implicite	2	$\forall h>0$	$0<h<\frac{2}{\beta}$
H	$u_{n+1}=u_n+\frac{h}{2}\left(\varphi(t_n,u_n)+\varphi(t_{n+1},u_n+h\varphi(t_n,u_n))\right)$	Explicite	2	$0<h<\frac{2}{\beta}$	$0<h<\frac{2}{\beta}$

Exemple : A-stabilité des schémas classiques à un pas¶

On considère le problème de Cauchy $$ \begin{cases} y'(t)=-y(t),\\ y(0)=1, \end{cases} $$ sur l'intervalle $[0;12]$.

Solution exacte¶

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$ $\mathbb{R}$.

Solution approchée par la méthode d'Euler explicite¶

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

si $0<h<1$ alors la solution numérique est stable et convergente,
si $h=1$ alors la solution numérique est stationnaire $u_n=0$ pour tout $n\in\mathbb{N}^*$,
si $1<h<2$ alors la solution numérique oscille mais est encore convergente,
si $h=2$ alors la solution numérique oscille, plus précisément on a $u_{2n}=1$ et $u_{2n+1}=-1$ pour tout $n\in\mathbb{N}^*$,
si $h>2$ alors la solution numérique oscille et diverge. En effet, la méthode est A-stable si et seulement si $|1 - h| < 1$.

Remarque: la suite obtenue est une suite géométrique de raison $q=1-h$. On sait qu'une telle suite

diverge si $|q|>1$ ou $q=-1$,
est stationnaire si $q=1$,
converge vers $0$ si $|q|<1$.
converge de façon monotone vers $0$ si $0<q<1$.

Voyons graphiquement ce que cela donne avec différentes valeurs de $h>0$:

si $h=4$ alors $t_n=4n$ et $u_{n}=\left(-4\right)^{n}$ tandis que $y(t_{n})=e^{-4n}$,
si $h=2$ alors $t_n=2n$ et $u_{n}=\left(-1\right)^{n}$ tandis que $y(t_{n})=e^{-2n}$,
si $h=1$ alors $t_n=n$ et $u_{n}=0$ tandis que $y(t_{n})=e^{-n}$,
si $h=\frac{1}{2}$ alors $t_n=\frac{n}{2}$ et $u_{n}=\left(\dfrac{1}{2}\right)^{n}$ tandis que $y(t_{n})=e^{-n/2}$,
si $h=\frac{1}{4}$ alors $t_n=\frac{n}{4}$ et $u_{n}=\left(\dfrac{3}{4}\right)^{n}$ tandis que $y(t_{n})=e^{-n/4}$. Ci-dessous sont tracées sur l'intervalle $[0;10]$, les courbes représentatives de la solution exacte et de la solution calculée par la méthode d'Euler explicite. En faisant varier le pas $h$ nous pouvons constater que si $h>1$ l'erreur commise entre la solution exacte et la solution calculée est amplifiée d'un pas à l'autre.
(Remarque: la première itérée a la même pente quelque soit le pas $h$ (se rappeler de la construction géométrique de la méthode d'Euler)).

%matplotlib inline
from matplotlib.pylab import *

# SCHEMA
def EE(phi,tt,y0):
    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(EE(phi,tt[-1],y0))

# 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();

Solution approchée par la méthode d'Euler implicite¶

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=\frac{1}{1+h}\in]0;1[$.

%matplotlib inline
from matplotlib.pylab import *
from scipy.optimize import fsolve

# SCHEMA
def EI(phi,tt,y0):
    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(EI(phi,tt[-1],y0))

# 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();