from IPython.display import display, Latex
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)
M62_CM5 : Convergence et A-stabilité des schémas multipas linéaires¶
$\newcommand{\R}{\mathbb{R}}$ $\newcommand{\N}{\mathbb{N}}$ $\newcommand{\dt}{\mathrm{d}t}$ $\newcommand{\dx}{\mathrm{d}x}$ $\newcommand{\dtau}{\mathrm{d}\tau}$ $\newcommand{\CC}[1]{\mathscr{C}}$
Convergence d'un schéma multipas linéaire¶
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=[t_0,T]$ 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$.
Les données initiales $u_0, u_1, \dots, u_p$ doivent être fournies.
Mis à part $u_0$, qui est égale à $y_0$ donné, les autres valeurs, $u_1, \dots, u_p$ peuvent être obtenues à l’aide de méthodes suffisament précises.
Si $b_{-1}=0$ la méthode est explicite, sinon implicite.
Les méthodes AB$_q$, AM$_q$, N$_q$, MS$_q$ et BFD$_q$ sont toutes des méthodes numériques à $q=p+1$ étapes (=pas) linéaires. En revanche, les méthodes EM et Heun ne rentrent pas dans ce cadre: elles sont des méthodes de la famille RK.
Exemples :
- Euler explicite
$$ u_{n+1}=u_n+h\varphi(t_n,u_n) \quad\implies\quad p=0,\ a_0=1,\ b_{-1}=0,\ b_0=1; $$ - Euler implicite
$$ u_{n+1}=u_n+h\varphi(t_{n+1},u_{n+1}) \quad\implies\quad p=0,\ a_0=1,\ b_{-1}=1,\ b_0=0; $$ - Crank-Nicolson $$ u_{n+1}=u_n+\frac{h}{2}\left(\varphi(t_n,u_n)+\varphi(t_{n+1},u_{n+1})\right) \quad\implies\quad p=0,\ a_0=1,\ b_{-1}=\frac{1}{2},\ b_0=\frac{1}{2}. $$
Polynômes caractéristiques¶
Réécrivons le schéma multistep comme suit: $$ u_{n+1} -\sum_{j=0}^p a_ju_{n-j}= h\left(\sum_{j=0}^p b_j\varphi(t_{n-j},u_{n-j})+b_{-1}\varphi(t_{n+1},u_{n+1})\right). $$ Formellement,
- pour construire le premier polynômes caractéristiques on ne regarde que le LHS et les $u_{n-j}$ deviennent des $r^{p-j}$,
- pour construire le second polynômes caractéristiques on ne regarde que le RHS divisé par $h$ et les $\varphi(t_{n-j},u_{n-j})$ deviennent des $r^{p-j}$.Définition: le polynôme $$ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j} $$ est appelé premier polynôme caractéristique associé à la méthode numérique.Définition: le polynôme $$ \sigma(r)=b_{-1}r^{p+1}+\sum_{j=0}^p b_jr^{p-j} $$ est appelé second polynôme caractéristique associé à la méthode numérique.
Exemples :
- Euler explicite: $\varrho(r)=r-1$ et $\sigma(r)=1$
- Euler implicite: $\varrho(r)=r-1$ et $\sigma(r)=r$
- Crank-Nicolson: $\varrho(r)=r-1$ et $\sigma(r)=\frac{1}{2}r+\frac{1}{2}$
- AB$_q$ ou AM$_q$: $\varrho(r)=r^{q+1}-r^{q}$
- N$_q$ ou MS$_q$: $\varrho(r)=r^{q+1}-r^{q-1}$
- BDF$_2$: $\varrho(r)=r^2-\frac{4}{3}r+\frac{1}{3}$ et $\sigma(r)=\frac{2}{3}$
Consistance¶
L’erreur de troncature locale de la méthode multi-pas s'écrit $$ \tau_n(h)=\frac{1}{h}\left(y_{n+1}-\sum_{j=0}^p a_jy_{n-j}-h\sum_{j=0}^p b_j\varphi(t_{n-j},y_{n-j})-hb_{-1}\varphi(t_{n+1},y_{n+1})\right). $$ La méthode est dite consistante si $\tau(h)=\max_n |\tau_n(h)|$ tend vers zéro quand $h$ tend vers zéro. Par un développement de Taylor assez fastidieux, on peut montrer que cette condition est équivalente à la condition suivante:
Exemple: schéma AB$_2$
- Schéma: $u_{n+1}=u_{n}+ \frac{h}{2} \big(3\varphi_{n}-\varphi_{n-1}\big)$
- Coefficients: $q=2$, $p=1$, $a_0=1$, $a_1=0$, $b_0=\frac{3}{2}$, $b_1=-\frac{1}{2}$, $b_{-1}=0$
- Consistance: $\begin{cases}\sum_{j=0}^p a_j=1+0=1\\-\sum_{j=0}^p ja_j+\sum_{j=-1}^p b_j=-(0a_0+1a_1)+(b_{-1}+b_0+b_1)=1\end{cases}$
- Ordre $\omega$:
- pour $i=2$ on a $\sum_{j=0}^p (-j)^{i}a_j+i\sum_{j=-1}^p (-j)^{i-1}b_j=\sum_{j=0}^p j^2 a_j+2\sum_{j=-1}^p j b_j=(0^2a_0+1^2a_1)+2(-1b_{-1}+0b_0+1b_1)=1$ donc $\omega\ge2$,
- pour $i=3$ on a $\sum_{j=0}^p (-j)^{i}a_j+i\sum_{j=-1}^p (-j)^{i-1}b_j=\sum_{j=0}^p (-j)^3 a_j+3\sum_{j=-1}^p j^2 b_j=(0^3a_0-1^3a_1)+3((-1)^2b_{-1}+0^2b_0+1^2b_1)=3$ donc $\omega<3$
- donc l'ordre $\omega$ est $2$
Zéro-stabilité¶
En particulier:
- les méthodes à un pas (EE,EI,CN) vérifient $\varrho(r)=r-1$: le polynôme $\varrho$ n'a qu'une racine, c'est $r=1$. Toutes ces méthodes sont donc zéro-stables. On peut donc conclure que pour les méthodes à un pas, "consistance $=$ convergence";
- les méthodes de Adam à $q$ pas vérifient $\varrho(r)=r^{p+1}-r^p=r^p(r-1)$: il n'y a que deux racines, la racine $1$ de multiplicité $1$ et la racine $0$ de multiplicité $p$. Toutes ces méthodes sont donc zéro-stables;
- les méthodes de Nystrom et Milne-Simpson à $q$ pas vérifient $\varrho(r)=r^{p+1}-r^{p-1}=r^{p-1}(r-1)(r+1)$: il n'y a que trois racines, les racines $\pm1$ de multiplicité $1$ et la racine $0$ de multiplicité $p-1$. Toutes ces méthodes sont donc zéro-stables;
- enfin, on peut montrer que les méthodes BDF sont zéro-stables ssi $p\le5$ (voir les calculs
sympyci-dessous).
La zero stabilité implique une restriction significative sur l’ordre atteignable pour une méthode multipas.
Exemples :
- Soit une méthode explicite à un pas, alors $q=1$, $p=0$, $b_{-1}=0$ et son ordre de consistance est au plus $\omega=p+1=1$.
(Exemple: EE qui a ordre 1) - Soit une méthode implicite à un pas, alors $q=1$, $p=0$, $b_{-1}\neq0$ et son ordre de consistance est au plus $\omega=p+2=2$.
(Exemple: EI qui a ordre 1 et CN qui a ordre 2) - Soit une méthode explicite à deux pas, alors $q=2$, $p=1$, $b_{-1}=0$ et son ordre de consistance est au plus $\omega=p+1=2$.
(Exemple: AB$_2$ qui a ordre 2) - Soit une méthode implicite à deux pas, alors $q=2$, $p=1$, $b_{-1}\neq0$ et son ordre de consistance est au plus $\omega=p+3=4$.
(Exemple: AM$_2$ qui a ordre 3 et MS$_2$ qui a ordre 4)
%reset -f
from matplotlib.pylab import *
from fractions import Fraction
from IPython.display import display, Math
import sympy as sym
sym.init_printing()
r=sym.Symbol('r')
sc = { "EE_AB1" : { "aa":[1], "bb":[ 1], "bm1":0 },
"EI_AM0_BDF1" : { "aa":[1], "bb":[ 0], "bm1":1 },
"CN_AM1" : { "aa":[1], "bb":[Fraction(1,2)], "bm1":1/2 },
"AB2": { "aa":[1,0], "bb":[Fraction(3,2),Fraction(-1,2)], "bm1":0 },
"AB3": { "aa":[1,0,0], "bb":[Fraction(23,12),Fraction(-16,12), Fraction(5,12)], "bm1":0 },
"AB4": { "aa":[1,0,0,0], "bb":[Fraction(55,24),Fraction(-59,24), Fraction(37,24), Fraction(-9,24)], "bm1":0 },
"AB5": { "aa":[1,0,0,0,0], "bb":[Fraction(1901,720), Fraction(-2774,720), Fraction(2616,720), Fraction(-1274,720), Fraction(251,720)], "bm1":0 },
"AM2": { "aa":[1,0], "bb":[Fraction(8,12),Fraction(-1,12)], "bm1":Fraction(5,12) },
"AM3": { "aa":[1,0,0], "bb":[Fraction(19,24),Fraction(-5,24),Fraction(1,24)], "bm1":Fraction(9,24) },
"AM4": { "aa":[1,0,0,0], "bb":[Fraction(646,720),Fraction(-264,720),Fraction(106,720),Fraction(-19,720)], "bm1":Fraction(251,720) },
"BDF2": { "aa":[Fraction(4,3),Fraction(-1,3)], "bb":[0,0], "bm1":Fraction(2,3) },
"BDF3": { "aa":[Fraction(18,11),Fraction(-9,11),Fraction(2,11)], "bb":[0,0,0], "bm1":Fraction(6,11) },
"BDF4": { "aa":[Fraction(48,25),Fraction(-36,25),Fraction(16,25),Fraction(-3,25)], "bb":[0,0,0,0], "bm1":Fraction(12,25) },
"BDF5": { "aa":[Fraction(300,137),Fraction(-300,137),Fraction(200,137),Fraction(-75,137),Fraction(12,137)], "bb":[0,0,0,0,0], "bm1":Fraction(60,137) },
"BDF6": { "aa":[Fraction(120,49),Fraction(-150,49),Fraction(400,147),Fraction(-75,49),Fraction(24,49),Fraction(-10,147)], "bb":[0,0,0,0,0,0], "bm1":Fraction(20,49) },
"MS2": { "aa":[0,1], "bb":[Fraction(4,3),Fraction(1,3)], "bm1":Fraction(1,3) }
}
for schema in sc.keys():
print("\n",schema)
q=len(sc[schema]["aa"])
p=q-1
print(f"\tq = {q}")
aa=sc[schema]["aa"]
bb=sc[schema]["bb"]
bm1=sc[schema]["bm1"]
if bm1==0:
print("\tExplicite")
else:
print("\tImplicite")
print("\tConsistance: ",end="")
ordre=[]
ordre.append( sum(aa)==1 and sum([-j*aa[j] for j in range(len(aa))])+bm1+sum(bb)==1 )
print(ordre[-1])
i=1
while ordre[-1]:
i+=1
ordre.append( sum([(-j)**i*aa[j]+i*(-j)**(i-1)*bb[j] for j in range(len(aa))])+i*bm1==1)
print(f"\tOrdre {len(ordre)-1}")
rho=array([1])
rho=append(rho,-array(aa))
display(Math(r'\varrho(r)='+sym.latex(sym.expand(poly1d(rho)(r)))))
# print(f"\tZéro-stabilité: |r_j|= {abs(roots(rho))}")
# print("\tRacines: ",end="")
# display(sym.solve(sym.expand(poly1d(rho)(r))))
display( Math(r' \{|r_j|\} = '+sym.latex( [abs(rac).simplify().evalf() for rac in sym.solve(sym.expand(poly1d(rho)(r)))] )))
A-stabilité¶
Dans la section précédente, on a considéré la résolution du problème de Cauchy sur des intervalles bornés.
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 $\lambda\in\mathbb{C}$ un nombre complexe tel que $\Re(\lambda)=-\beta<0$ et considérons le problème de Cauchy $$\begin{cases} y'(t)=\lambda y(t), &\text{pour }t>0,\\ y(0)=1 \end{cases}$$ Sa solution est $y(t)=e^{\lambda t}$ et comme $\Re(\lambda)<0$ alors $$\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 (absolument stable).
Région de stabilité absolue¶
On appelle région de stabilité absolue $A$ d’une méthode numérique l’ensemble des nombres complexes $z = h\lambda$ pour lesquels la méthode est absolument stable: $$ A=\Big\{z = h\lambda \in \mathbb{C} \ \Big| \ |u_n|\xrightarrow[n\to+\infty]{}0 \Big\} $$ Une méthode numérique est inconditionellement A-stable si $\mathbb{C}^-\subset A$ où $\mathbb{C}^-=\{z\in\mathbb{C}|\Re(z)<0\}$.
Une méthode à $q=p+1$ pas appliquée à ce problème modèle conduit à léquation récurrente à $$ u_{n+1}=\sum_{j=0}^p a_j u_{n-j} + h\lambda\sum_{j=0}^p b_ju_{n-j} + h\lambda b_{-1}u_{n+1} $$ dont le polynôme caracthéristique associé $$ p(r) = \varrho(r) - (h\lambda) \sigma(r) $$ est appelé polynôme de A-stabilité associé à la méthode numérique.
Régions de A-stabilité de diverses méthodes d'Adams-Bashforth (à gauche) et d'Adams-Moulton (à droite).
Source: A. Quarteroni, F. Saleri, P. Gervasio Calcul Scientifique: Cours, exercices corrigés et illustrations en Matlab et Octave.
On remarque que les régions de stabilité absolue des méthodes d'Adams-Bashforth et d'Adams-Moulton sont bornées et diminuent quand l'ordre augmente.
Régions de A-stabilité de diverses méthodes BDF.
Source: A. Quarteroni, F. Saleri, P. Gervasio Calcul Scientifique: Cours, exercices corrigés et illustrations en Matlab et Octave.
On remarque que les régions de stabilité absolue des méthodes BDF, pour $q\le6$, sont non bornées. En particulier, elles contiennent toujours les réels négatifs.
Tracer les régions de A-stabilité.
Posons $z=h\lambda$.
Il est en général difficile de déterminer la région de stabilité absolue car on doit décider pour quels $z\in\mathbb{C}$ les racines du polynôme de stabilité vérifient la condition de racines stricte (i.e. $|r| < 1$).
Il est plus aisé de déterminer la frontière de cette région car au moins une des racines de $p$ à la frontière est de module $1$.
On a donc $|r| = 1$ sur $\partial A$.
La frontière de $A$ est alors un sous ensemble des points $z\in\mathbb{C}$ pour lesquels $r=e^{i\vartheta}$, $\vartheta\in\mathbb{R}$.
On remplace alors $r =e^{i\vartheta}$ dans $p$ et comme $e^{i\vartheta}$ est racine, $p(e^{i\vartheta}) = 0$ et on résoud l’équation. On obtient $z = z(\vartheta)$ ce qui donne une courbe dans $\mathbb{C}$ en fonction de l’angle $\vartheta$, ce qui décrit donc une courbe polaire.
Par conséquent, pour obtenir une représentation approchée de $\partial A$, il suffit d’évaluer le second membre pour diverses valeurs de $r$ sur le cercle unité, diverses valeurs de $\vartheta\in\mathbb{R}$ et $r=e^{i\vartheta}$.
%reset -f
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
sp.init_printing()
def polyCar(schema):
sp.var('r z')
p,aa,bb,bm1=schema["p"],schema["aa"],schema["bb"],schema["bm1"]
rho=r**(p+1)-sum(aa[j]*r**(p-j)for j in range(p+1))
sigma=bm1*r**(p+1)+sum(bb[j]*r**(p-j)for j in range(p+1))
display(f'rho(r)={rho}')
display(f'sigma(r)={sigma}')
return
def polyCarIMPL(schema):
sp.var('r z')
p,aa,bb,bm1=schema["p"],schema["aa"],schema["bb"],schema["bm1"]
rho=r**(p+1)-sum(aa[j]*r**(p-j)for j in range(p))
sigma=bm1*r**(p+1)+sum(bb[j]*r**(p-j)for j in range(p))
display(f'rho(r)={rho}')
display(f'sigma(r)={sigma}')
return
# Returns > 1 if argument is not in region of absolute stability
def stabilityFunction(z,schema):
pol=np.array(schema["aa"])+z*np.array(schema["bb"])
pol=np.append([1-z*schema["bm1"]],-pol)
return max(abs(np.roots(pol)))
def Astabilite(schema):
p,aa,bb,bm1=schema["p"],schema["aa"],schema["bb"],schema["bm1"]
# cc liste de points sur le circle unitaire
rr = np.linspace(0,2*np.pi,201)
cc = np.cos(rr)+1j*np.sin(rr)
## ff fonction vectorisée qui donne zz
ff = lambda cc,p,aa,bb,bm1 : (cc**(p+1)-np.polyval(aa,cc))/(bm1*cc**(p+1)+np.polyval(bb,cc))
zz = ff(cc,p,aa,bb,bm1)
plt.figure(figsize=(18,7))
plt.subplot(1,2,1)
plt.plot(zz.real,zz.imag,'r')
plt.grid(True)
# plt.axis('equal')
plt.xlabel(r'$\Re(h\beta)$')
plt.ylabel(r'$\Im(h\beta)$');
plt.title(schema["Nom"])
Lx,Ly = plt.xlim(),plt.ylim()
plt.subplot(1,2,2)
x = np.linspace(Lx[0],Lx[1], 300)
y = np.linspace(Ly[0],Ly[1], 300)
[X,Y] = np.meshgrid(x,y)
Z = np.zeros(X.shape)
for m in range(X.shape[0]):
for n in range(X.shape[1]):
Z[m,n] = stabilityFunction(X[m,n] + 1j * Y[m,n], schema)
plt.contour(X, Y, Z, [1], colors='k')
plt.contourf(X, Y, Z, [0,1], colors=[[1, 0.5, 0.8]])
return None
EE = { "Nom":"EE=AB1", "p":0, "aa":[1], "bb":[1], "bm1":0 }
polyCar(EE)
Astabilite(EE)
AB2 = { "Nom":"AB2", "p":1, "aa":[1,0], "bb":[3/2,-1/2], "bm1":0 }
polyCar(AB2)
Astabilite(AB2)
AB3 = { "Nom":"AB3", "p":2, "aa":[1,0,0], "bb":[23/12,-16/12, 5/12], "bm1":0 }
polyCar(AB3)
Astabilite(AB3)
AB4 = { "Nom":"AB4", "p":3, "aa":[1,0,0,0], "bb":[55/24,-59/24, 37/24, -9/24], "bm1":0 }
polyCar(AB4)
Astabilite(AB4)
AB5 = { "Nom":"AB5", "p":4, "aa":[1,0,0,0,0], "bb":[1901/720, -2774/720, 2616/720, -1274/720, 251/720], "bm1":0 }
polyCar(AB5)
Astabilite(AB5)
AM1 = { "Nom":"AM1=CN", "p":1, "aa":[1], "bb":[1/2], "bm1":1/2 }
polyCarIMPL(AM1)
#Astabilite(AM1)
AM2 = { "Nom":"AM2", "p":1, "aa":[1,0], "bb":[8/12,-1/12], "bm1":5/12 }
polyCar(AM2)
Astabilite(AM2)
AM3 = { "Nom":"AM3", "p":2, "aa":[1,0,0], "bb":[19/24,-5/24,1/24], "bm1":9/24 }
polyCar(AM3)
Astabilite(AM3)
AM4 = { "Nom":"AM4", "p":3, "aa":[1,0,0,0], "bb":[646/720,-264/720,106/720,-19/720], "bm1":251/720 }
polyCar(AM4)
Astabilite(AM4)
MS2 = { "Nom":"MS2", "p":1, "aa":[0,1], "bb":[4/3,1/3], "bm1":1/3 }
polyCar(MS2)
#Astabilite(MS2)
BDF1 = { "Nom":"BDF1=EI", "p":1, "aa":[1], "bb":[0], "bm1":1 }
polyCarIMPL(BDF1)
Astabilite(BDF1)
BDF2 = { "Nom":"BDF2", "p":1, "aa":[4/3,-1/3], "bb":[0], "bm1":2/3 }
polyCarIMPL(BDF2)
Astabilite(BDF2)
BDF3 = { "Nom":"BDF3", "p":2, "aa":[18/11,-9/11,2/11], "bb":[0], "bm1":6/11 }
# polyCarIMPL(BDF3)
Astabilite(BDF3)
BDF4 = { "Nom":"BDF4", "p":3, "aa":[48/25,-36/25,16/25,-3/25], "bb":[0], "bm1":12/25 }
# polyCarIMPL(BDF4)
Astabilite(BDF4)