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)
Considérons le problème de Cauchy
Pour $h>0$ soit $t_n\equiv t_0+nh$ avec $n=0,1,2,\dots,N$ une suite de $N+1$ nœuds de $I$ induisant une discrétisation de $I$ en $N$ sous-intervalles $I_n=[t_n;t_{n+1}]$ chacun de longueur $h=\frac{T-t_0}{N}>0$ (appelé le pas de discrétisation).
Pour chaque nœud $t_n$, on cherche la valeur inconnue $u_n$ qui approche la valeur exacte $y_n\equiv y(t_n)$.
Méthodes explicites et méthodes implicites
Une méthode est dite explicite si la valeur $u_{n+1}$ peut être calculée directement à l'aide des valeurs précédentes $u_k$, $k\le n$ (ou d'une partie d'entre elles). Pour un schéma linéaire cela équivaut à $b_{-1}=0$.
Une méthode est dite implicite si $u_{n+1}$ n'est défini que par une relation implicite faisant intervenir la fonction $\varphi$. Pour un schéma linéaire cela équivaut à $b_{-1}\neq0$.
Méthodes à un pas et méthodes multi-pas
Une méthode numérique pour l'approximation du problème de Cauchy est dite à un pas si pour tout $n\in\mathbb{N}$, $u_{n+1}$ ne dépend que de $u_n$ et éventuellement de lui-même.
Autrement, on dit que le schéma est une méthode multi-pas (ou à pas multiples).
Les schémas d'Euler explicite, d'Euler implicite et de Crank-Nicolson sont des schémas à un pas ($p=0$) linéaires
$$ u_{n+1} = a_0u_{n}+hb_0\varphi(t_{n},u_{n})+hb_{-1}\varphi(t_{n+1},u_{n+1}),$$avec
| Type | $a_0$ | $b_0$ | $b_{-1}$ | ||
|---|---|---|---|---|---|
| EE | $u_{n+1}=u_n+h\varphi(t_n,u_n)$ | Explicite | 1 | 1 | 0 |
| EI | $u_{n+1}=u_n+h\varphi(t_{n+1},u_{n+1})$ | Implicite | 1 | 0 | 1 |
| CN | $u_{n+1}=u_n+\frac{h}{2}\left(\varphi(t_n,u_n)+\varphi(t_{n+1},u_{n+1})\right)$ | Implicite | 1 | $$\tfrac{1}{2}$$ | $$\tfrac{1}{2}$$ |
Dans la suite nous allons écrire explicitement ces schémas.
Pour vérifier nos calculs d'interpolation puis intégration, nous allons utiliser le package de calcul formel SymPy.
import sympy as symb
symb.init_printing()
symb.var('h,t,t_n')
t_np1 = t_n + h
t_nm1 = t_n - h
t_nm2 = t_n - 2 * h
t_nm3 = t_n - 3 * h
t_nm4 = t_n - 4 * h
t_nm5 = t_n - 5 * h
#
y_np1=symb.Symbol('y_{n+1}')
y_n =symb.Symbol('y_{n}')
y_nm1=symb.Symbol('y_{n-1}')
y_nm2=symb.Symbol('y_{n-2}')
y_nm3=symb.Symbol('y_{n-3}')
y_nm4=symb.Symbol('y_{n-4}')
y_nm5=symb.Symbol('y_{n-5}')
y_np1=symb.Symbol('y_{n+1}')
# Pour Adam et N-MS
phi_np1=symb.Symbol('\phi_{n+1}')
phi_n =symb.Symbol('\phi_{n}')
phi_nm1=symb.Symbol('\phi_{n-1}')
phi_nm2=symb.Symbol('\phi_{n-2}')
phi_nm3=symb.Symbol('\phi_{n-3}')
phi_nm4=symb.Symbol('\phi_{n-4}')
phi_nm5=symb.Symbol('\phi_{n-5}')
# Pour BDF
p = symb.Function('p')
u_np1=symb.Symbol('u_{n+1}')
u_n =symb.Symbol('u_{n}')
u_nm1=symb.Symbol('u_{n-1}')
u_nm2=symb.Symbol('u_{n-2}')
u_nm3=symb.Symbol('u_{n-3}')
u_nm4=symb.Symbol('u_{n-4}')
u_nm5=symb.Symbol('u_{n-5}')
Si nous intégrons l'EDO $y'(t)=\varphi(t,y(t))$ entre $t_n$ et $t_{n+1}$ nous obtenons $$ y_{n+1}-y_n=\int_{t_n}^{t_{n+1}} \varphi(t,y(t))dt. $$ On peut construire différents schémas selon la formule de quadrature utilisée pour approcher le membre de droite. Cette solution approchée sera obtenue en construisant une suite récurrente comme suit:
Les schémas d'Adam approchent l'intégrale $\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t$ par l'intégrale d'un polynôme $p$ interpolant $\varphi$ en des points donnés qui peuvent être à l'extérieur de l'intervalle d'intégration $[t_{n};t_{n+1}]$.
On peut construire différents schémas selon les points d'interpolation choisis.
Toutes ces méthodes peuvent s'écrire sous la forme
Ils se divisent en deux familles:
Remarques
On a
</div> La méthode AB$_1$ coïncide avec la méthode d'Euler progressive.
pol = symb.interpolate([(t_n, phi_n)], t)
#display(pol)
integrale = symb.integrate(pol, (t, t_n, t_np1)).simplify()
display(symb.Eq(u_np1, u_n + integrale))
On a
pol = symb.interpolate([(t_n, phi_n), (t_nm1, phi_nm1)], t)
# display(pol)
integrale = symb.integrate(pol, (t, t_n, t_np1)).simplify()
display(symb.Eq(u_np1, u_n + integrale))
On a
pol = symb.interpolate([(t_n, phi_n), (t_nm1, phi_nm1), (t_nm2, phi_nm2)], t)
# display(pol)
integrale = symb.integrate(pol, (t, t_n, t_np1)).simplify()
display(symb.Eq(u_np1, u_n + integrale))
On a
pol = symb.interpolate([(t_n, phi_n), (t_nm1, phi_nm1), (t_nm2, phi_nm2),
(t_nm3, phi_nm3)], t)
integrale = symb.integrate(pol, (t, t_n, t_np1)).simplify()
display(symb.Eq(u_np1, u_n + integrale))
Pour la construction des autres schémas AB$_5$ et AB$_6$ on verra plus bas dans le récapitulatif.
On a
</div> La méthode AM$_0$ coïncide avec la méthode d'Euler regressive.
pol = symb.interpolate([(t_np1, phi_np1)], t)
integrale = symb.integrate(pol, (t, t_n, t_np1)).simplify()
display(symb.Eq(u_np1, u_n + integrale))
On a
</div> La méthode AM$_1$ coïncide avec la méthode de Crank-Nicolson.
pol = symb.interpolate([(t_np1, phi_np1), (t_n, phi_n)], t)
integrale = symb.integrate(pol, (t, t_n, t_np1)).simplify()
display(symb.Eq(u_np1, u_n + integrale))
On a
pol = symb.interpolate([(t_np1, phi_np1), (t_n, phi_n), (t_nm1, phi_nm1)], t)
integrale = symb.integrate(pol, (t, t_n, t_np1)).simplify()
display(symb.Eq(u_np1, u_n + integrale))
Pour la construction des autres schémas AM$_3$, AM$_4$, AM$_5$ et AM$_6$ on verra plus bas dans le récapitulatif.
Points = [(t_n, phi_n), (t_nm1, phi_nm1), (t_nm2, phi_nm2), (t_nm3, phi_nm3),
(t_nm4, phi_nm4), (t_nm5, phi_nm5)]
for q in range(1, len(Points) + 1):
pol = symb.interpolate(Points[:q], t)
AB = (symb.integrate(pol, (t, t_n, t_np1)).simplify())
print(f"AB-{q}: explicite, à {q} pas, d'ordre {q}")
display(symb.Eq(u_np1, u_n + AB))
print("\n")
print("\n")
Points = [(t_np1, phi_np1), (t_n, phi_n), (t_nm1, phi_nm1), (t_nm2, phi_nm2),
(t_nm3, phi_nm3), (t_nm4, phi_nm4), (t_nm5, phi_nm5)]
for q in range(len(Points)):
pol = symb.interpolate(Points[0:q + 1], t)
AM = (symb.integrate(pol, (t, t_n, t_np1)).simplify())
print(f"AM-{q}: implicite, à {q} pas, d'ordre {q+1}")
display(symb.Eq(u_np1, u_n + AM))
print("\n")
Les méthodes d'Adam peuvent être facilement généralisées en intégrant l'EDO $y'(t)=\varphi(t,y(t))$ entre $t_{n-r}$ et $t_{n+1}$ avec $r\ge1$.
Avec $r=1$, si nous intégrons l'EDO $y'(t)=\varphi(t,y(t))$ entre $t_{n-1}$ et $t_{n+1}$ nous obtenons $$ y_{n+1}-y_{n-1}=\int_{t_{n-1}}^{t_{n+1}} \varphi(t,y(t))dt. $$ On peut construire différents schémas selon la formule de quadrature utilisée pour approcher le membre de droite.
Cette solution approchée sera obtenue en construisant une suite récurrente comme suit:
Ces schémas approchent l'intégrale $\int_{t_{n-1}}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t$ par l'intégrale d'un polynôme $p$ interpolant $\varphi$ en des points donnés qui peuvent être à l'extérieur de l'intervalle d'intégration $[t_{n-1};t_{n+1}]$.
On peut construire différents schémas selon les points d'interpolation choisis.
Toutes ces méthodes peuvent s'écrire à nouveau sous la forme
Ils se divisent en deux familles:
On a
La méthode N$_1$ coïncide avec la méthode du point milieu (appelée aussi Saute-mouton ou Leapfrog).
Remarque: on l'appelle N$_1$ mais il est à 2 pas car il dépend de $u_{n-1}$.
pol = symb.interpolate([(t_n, phi_n)], t).simplify()
integrale = symb.integrate(pol, (t, t_nm1, t_np1)).simplify()
display(symb.Eq(u_np1, u_nm1 + integrale))
On a
où $u_{1}$ est une approximation de $y(t_1)$ obtenue en utilisant une prédiction d'Euler explicite.
</div> On retrouve la méthode N$_1$.
pol = symb.interpolate([(t_n, phi_n), (t_nm1, phi_nm1)], t).simplify()
integrale = symb.integrate(pol, (t, t_nm1, t_np1)).simplify()
display(symb.Eq(u_np1, u_nm1 + integrale))
On a
pol = symb.interpolate([(t_n, phi_n), (t_nm1, phi_nm1), (t_nm2, phi_nm2)],
t).simplify()
integrale = symb.integrate(pol, (t, t_nm1, t_np1)).simplify()
display(symb.Eq(u_np1, u_nm1 + integrale))
Pour la construction des autres schémas N$_4$, N$_5$ et N$_6$ on verra plus bas dans le récapitulatif.
On a
Remarque: on l'appelle MS$_0$ mais il est à 2 pas car il dépend de $u_{n-1}$.
pol = symb.interpolate([(t_np1, phi_np1)], t).simplify()
integrale = symb.integrate(pol, (t, t_nm1, t_np1)).simplify()
display(symb.Eq(u_np1, u_nm1 + integrale))
On a
où $u_{1}$ est une approximation de $y(t_1)$ obtenue en utilisant une prédiction d'Euler implicite.
</div> On retrouve la méthode N$_1$: le schéma est devenu explicite.
Remarque: on l'appelle MS$_1$ mais il est à 2 pas car il dépend de $u_{n-1}$.
pol = symb.interpolate([(t_np1, phi_np1), (t_n, phi_n)], t).simplify()
integrale = symb.integrate(pol, (t, t_nm1, t_np1)).simplify()
display(symb.Eq(u_np1, u_nm1 + integrale))
On a
pol = symb.interpolate([(t_np1, phi_np1), (t_n, phi_n), (t_nm1, phi_nm1)],
t).simplify()
integrale = symb.integrate(pol, (t, t_nm1, t_np1)).simplify()
display(symb.Eq(u_np1, u_nm1 + integrale))
Pour la construction des autres schémas MS$_3$, MS$_4$, MS$_5$ et MS$_6$ on verra plus bas dans le récapitulatif.
Points = [(t_n, phi_n), (t_nm1, phi_nm1), (t_nm2, phi_nm2), (t_nm3, phi_nm3),
(t_nm4, phi_nm4), (t_nm5, phi_nm5)]
for q in range(1, len(Points) + 1):
pol = symb.interpolate(Points[:q], t)
N = (symb.integrate(pol, (t, t_nm1, t_np1)).simplify())
print(f"N-{q}: explicite, à {q} pas")
display(symb.Eq(u_np1, u_nm1 + N))
print("\n")
print("\n")
Points = [(t_np1, phi_np1), (t_n, phi_n), (t_nm1, phi_nm1), (t_nm2, phi_nm2),
(t_nm3, phi_nm3), (t_nm4, phi_nm4), (t_nm5, phi_nm5)]
for q in range(len(Points)):
pol = symb.interpolate(Points[:q + 1], t)
MS = (symb.integrate(pol, (t, t_nm1, t_np1)).simplify())
print(f"MS-{q}: implicite, à {q} pas")
display(symb.Eq(u_np1, u_nm1 + MS))
print("\n")
Soit $p$ le polynôme qui interpole la fonction inconnue $t\mapsto y(t)$ en les points $$ \{t_{n+1},t_n,\dots, t_{n-q+1}\}\text{ avec } q\ge 1. $$
Au lieu d'évaluer l'EDO $\varphi(t,y) = y'(t)$ en $t_{n+1}$, nous évaluons l'EDO approchée $\varphi(t,y) = p'(t)$ en $t_{n+1}$ ce qui donne la relation $$ \varphi(t_{n+1},y_{n+1}) \simeq p'(t_{n+1}). $$
On peut construire différents schémas implicites selon les points d'interpolation utilisés pour approcher la fonction inconnue $t\mapsto y(t)$. Cette solution approchée sera obtenue en construisant une suite récurrente comme suit:
Toutes ces méthodes peuvent s'écrire sous la forme
On verra que les méthodes BDF-$q$ à $q$ pas sont d'ordre $q$ et sont zéro-stables ssi $q\le6$.
On a
</div> La méthode BDF$_1$ coïncide avec la méthode d'Euler implicite.
pol = symb.interpolate([(t_np1, u_np1), (t_n, u_n)], t).factor()
# display(pol)
display(symb.Eq(p(t),pol))
dp = symb.diff(pol, t).subs(t, t_np1)
display(symb.Eq(symb.diff(p(t),t).subs(t, t_np1),dp))
rhs=symb.solve(dp - phi_np1, u_np1)[0]
display(symb.Eq(u_np1,rhs))
On a
pol = symb.interpolate([(t_np1, u_np1), (t_n, u_n), (t_nm1, u_nm1)],
t).simplify()
display(symb.Eq(p(t),pol))
dp = symb.diff(pol, t).subs(t, t_np1)
display(symb.Eq(symb.diff(p(t),t).subs(t, t_np1),dp))
rhs=symb.solve(dp - phi_np1, u_np1)[0]
display(symb.Eq(u_np1,rhs))
Points = [(t_np1, u_np1), (t_n, u_n), (t_nm1, u_nm1), (t_nm2, u_nm2),
(t_nm3, u_nm3), (t_nm4, u_nm4), (t_nm5, u_nm5)]
for q in range(1, len(Points)):
pol = symb.interpolate(Points[0:q + 1], t)
dp = symb.diff(pol, t).subs(t, t_np1)
sol = symb.solve(dp - phi_np1, u_np1)
#display(sol)
print(f"BDF-{q}: implicite, à {q} pas")
display(symb.Eq(u_np1, sol[0]))
print("\n")
Lorsqu'on utilise une méthode implicite, pour calculer $u_{n+1}$ on doit résoudre une équation (en générale non-linéaire), par exemple avec une méthode de point fixe.
Une approche différente qui permet de s'affranchir de cette étape est donnée par les méthodes predictor-corrector. Une méthode predictor-corrector est une méthode qui permet de calculer $u_{n+1}$ de façon explicite à partir d'une méthode implicite.
Ordre Parmi les méthodes linéaires à $q$ pas, les méthodes implicites sont, à nombre de pas égal, plus précises et en général plus stables (on verra ces propriétés lors du prochaine cours). Cependant, elles sont aussi plus coûteuses car le caractère implicite doit être résolu par une méthode itérative, de type point fixe par exemple. D’où l’idée, pour diminuer le coût, d’utiliser une méthode explicite pour calculer un bon prédicteur de la solution et de ne faire qu’une itération (ou quelques itérations) de la méthode implicite utilisée alors comme correcteur explicite.
Ces méthodes héritent de l’ordre de précision du correcteur. (La contribution à l’erreur du prédicteur est d’un ordre inférieur à celle du correcteur). Pour avoir une méthode d’ordre $\omega$ on pourra prendre un prédicteur d’ordre $\omega-1$ et un correcteur d’ordre $\omega$ (typiquement deux méthodes de Adam ayant le même nombre de pas).
Schématiquement, si on note
alors l'ordre du predictor-corrector est $\omega_{[PC]}=\min\{\omega_{[C]},\omega_{[P]}+1\}$
Zéro-stabilité On démontre que la zéro-stabilité du prédicteur n’influe pas sur la zéro-stabilité de la méthode prédicteur-correcteur. On peut même utiliser un prédicteur instable.
A-stabilité Elles sont soumises à une condition de A-stabilité qui est typiquement celle du prédicteur. Elles ne sont donc pas adaptées à la résolution des problèmes de Cauchy sur des intervalles non bornés ou des problèmes stiff.
Considérons les deux schémas PC suivants:
On retrouve la méthode de Heun
Des méthodes de type predictor-corrector sont souvent construites en utilisant une prédiction d'Adam-Bashforth suivie d'une correction d'Adam-Moulton.
Par exemple, si on considère les deux étapes suivantes
On obtient la méthode AB$_2$-AM$_2$:
Des méthodes de type predictor-corrector sont souvent construites en utilisant une prédiction d'Adam-Bashforth suivie d'une correction BDF.
Par exemple, si on considère les deux étapes suivantes
On obtient la méthode AB$_2$-BDF$_2$: