Considérons le problème de Cauchy dont on suppose l’existence et l’unicité de la solution \(y\) sur un intervalle \(I=[t_0,T]\) :
\(\begin{cases}
y'(t) = \varphi(t,y(t)), &\forall t \in I,\\
y(t_0) = y_0,
\end{cases}\)
où \(\varphi : I \times \mathbb{R} \to \mathbb{R}\) est une fonction donnée et \(y_0\) est une donnée initiale également donnée.
Discrétisation de l’intervalle \(I\) : posons \(h=\frac{T-t_0}{N}>0\) (appelé le pas de discrétisation) et définissons \(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 l’intervalle \(I\) en \(N\) sous-intervalles \(I_n=[t_n;t_{n+1}]\) chacun de longueur \(h\).
- L’ensemble de \(N+1\) valeurs \(\{t_0, t_1=t_0+h,\dots , t_{N}=T \}\) représente les points (ou nœud) de la discrétisation.
Discrétisation de la solution : pour chaque nœud \(t_n\), on cherche la valeur inconnue \(u_n\) qui approche la valeur exacte \(y(t_n)\).
- L’ensemble de \(N+1\) valeurs \(\{y_0, y_1,\dots , y_{N} \}\) avec \(y_n\equiv y(t_n)\) représente la solution exacte discrète (en général inconnue).
- L’ensemble de \(N+1\) valeurs \(\{u_0 = y_0, u_1,\dots , u_{N} \}\) représente la solution numérique. Les schémas que l’on va construire permettent de calculer (explicitement ou implicitement) \(u_{n+1}\) à partir de \(u_n, u_{n-1}, \dots\) et il est donc possible de calculer successivement \(u_1\), \(u_2\),…, en partant de \(u_0\) par une formule de récurrence du type \(\begin{cases}
\left. \begin{array}{l}
u_0=y_0, \text{ connu}\\
u_1 \approx y_1,\\
\vdots\\
u_\kappa\approx y_\kappa,
\end{array}\right] \text{ initialisation}\\
u_{n+1}=\Phi(u_{n+1},u_n, u_{n-1}, \dots, u_{n-\kappa}), &\forall n=\kappa,\kappa+1,\dots,N-1.
\end{cases}\)
MĂ©thodes explicites VS 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_j\), \(j\le n\) (ou d’une partie d’entre elles).
- Une méthode est dite implicite si \(u_{n+1}\) n’est défini que par une relation implicite faisant intervenir la fonction \(\varphi\) évaluée en \(t_{n+1}\) et \(u_{n+1}\).
Méthodes à un pas VS méthodes multi-pas et initialisation(s)
- 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 (donc \(\kappa=0\)).
- Autrement, on dit que le schéma est une méthode multi-pas (ou à pas multiples ou multistep). Noter que, pour calculer successivement \(u_{\kappa+1}\), \(u_{\kappa+2}\), …, il est nécessaire d’initialiser d’abord \(u_0, u_1, \dots, u_{\kappa}\). Cette initialisation doit être effectuée par des approximations appropriées qui ne dégradent pas l’ordre de la méthode.
Définition d’un schéma multi-pas linéaire à \(q\) pas
Soit \(p\ge0\) un entier. Un schéma multi-pas linéaire à \(q=p+1\) pas (= étapes) pour l’approximation du problème de Cauchy consiste à définir une suite de valeurs \(\{u_n\}_{n=0}^N\) par la relation de récurrence
\(\begin{cases}
\left. \begin{array}{l}
u_0=y_0, \text{ connu}\\
u_1 \approx y_1,\\
\vdots\\
u_p\approx y_p,
\end{array}\right] \text{ initialisation}\\
\displaystyle
u_{n+1} = \sum_{j=0}^p a_ju_{n-j} + h\sum_{j=0}^p b_j\varphi(t_{n-j},u_{n-j}) + hb_{-1}\varphi(t_{n+1},u_{n+1}), \qquad
n=p,p+1,\dots,N-1
\end{cases}\)
où les \(\{a_j\}_{j=0}^p\) et \(\{b_j\}_{j=-1}^p\) sont \(2p+1\) coefficients donnés déterminant le schéma.
Le terme linéaire traduit le fait que le schéma est construit comme une combinaison linéaire des valeurs \(u_{n-j}\) et des valeurs de la fonction \(\varphi\) évaluée en ces même points.
Notation : la fonction \(\varphi\) est évaluée seulement en \(t_{n-j}\) pour \(j=-1,0,1,\dots,p\), qui sont des points de la grille de discrétisation, ainsi sans ambiguïté, on peut écrire \(\varphi_{n-j}\) au lieu de \(\varphi(t_{n-j},y_{n-j})\).
Notons qu’on peut réécrire la formule précédente en regroupant les termes dans d’autres manières équivalentes en privilégiant les termes en \(u_{n-j}\) pour \(j\in\{0,1,\dots,p\}\) et en \(u_{n+1}\) pour bien visualiser si le schéma est explicite ou implicite ou en privilégiant la distinction entre les termes en \(u_{n-j}\) et ceux en \(\varphi_{n-j}\). Par exemple, on peut écrire
\(
u_{n+1}
= \sum_{j=0}^p \Big(a_j u_{n-j} + h b_j\varphi_{n-j} \Big) + hb_{-1}\varphi_{n+1}
\)
ou encore
\(
\left(u_{n+1}-\sum_{j=0}^p a_ju_{n-j}\right) = h\sum_{j=-1}^p b_j\varphi_{n-j}.
\)
- Méthodes explicites et méthodes implicites : si \(b_{-1}=0\), le schéma est explicite; sinon, il est dit implicite.
- Méthodes à un pas et méthodes multi-pas : si \(p=0\), le schéma est à un pas (\(q=p+1=1\)); sinon, il est multi-pas.
Cas particuliers :
- Schémas de type Adams : \(a_0=1\) et \(a_j=0\) pour \(j\neq0\), i.e. de la forme \(\left(u_{n+1}-u_{n}\right) = h\sum_{j=-1}^p b_j\varphi_{n-j}.\)
- Schéma de Nyström et de Milne-Simpson : \(a_1=1\) et \(a_j=0\) pour \(j\neq1\), i.e. de la forme \(\left(u_{n+1}-u_{n-1}\right) = h\sum_{j=-1}^p b_j\varphi_{n-j}.\)
- Schémas de type BDF (Backward Differentiation Formula) : \(b_{-1}=1\) et \(b_j=0\) pour \(j\neq-1\), i.e. de la forme \(\left(u_{n+1}-\sum_{j=0}^p a_ju_{n-j}\right) = h b_{-1}\varphi_{n+1}.\)
Exemples
Les schémas d’Euler explicite (EE), d’Euler implicite (EI) et de Crank-Nicolson (CN) sont des schémas linéaires à \(q=1\) pas (\(p=0\)) car de la forme :
\(\begin{cases}
u_0=y_0,\\
u_{n+1} = a_0u_{n}+hb_0\varphi_{n}+hb_{-1}\varphi_{n+1}, &\forall n=0,1,\dots,N-1.
\end{cases}\)
En effet, on a les relations suivantes :
\(
\begin{array}{|l|l|c|c|c|c|}
\hline
& \text{RĂ©currence} & \text{Type} & a_0 & b_0 & b_{-1} \\
\hline
\text{EE} & u_{n+1}=u_n+h\varphi_n & \text{Explicite} & 1 & 1 & 0 \\
\text{EI} & u_{n+1}=u_n+h\varphi_{n+1} & \text{Implicite} & 1 & 0 & 1 \\
\text{CN} & u_{n+1}=u_n+\frac{h}{2}\left(\varphi_n+\varphi_{n+1}\right) & \text{Implicite} & 1 & \frac{1}{2} & \frac{1}{2} \\
\hline
\end{array}
\)
NB :
- le schéma d’Euler modifié (EM) et le schéma (RK1_M) ne sont pas des schémas multi-pas linéaires car ils nécéssitent l’évaluation de \(\varphi\) en un point qui n’appartient pas à la grille de discrétisation, à savoir en \(t_n+\frac{h}{2}\). On verra que ces types de schémas sont des schémas de la famille Runge-Kutta (on les verras au dernier chapitre) ;
- le schéma de Heun non plus n’est pas un schéma multi-pas linéaire car il évalue \(\varphi\) non pas en \((t_{n+1},u_{n+1})\) mais en \((t_{n+1},\tilde u_{n+1})\) où \(\tilde u_{n+1}\approx u_{n+1}\). On verra que ce type de schéma peut être considéré soit comme un schéma de type predicteur-corrector (dernière section de ce chapitre) ou encore comme un schéma de la famille Runge-Kutta.
Dans la suite de ce chapitre nous allons nous intéresser à la construction de schémas linéaires à \(q\) pas bien connus en littérature et concrètement utilisés pour l’approximation numérique de problèmes de Cauchy.
En particulier, nous allons voir 2 stratégies pour construire ces schémas :
- intégration sur un intervalle de temps d’une approximation polynomiale de la fonction \(\varphi\) :
- méthodes de Adams : intégration entre \([t_n,t_{n+1}]\) d’une approximation de la fonction \(\varphi\) par un polynôme de degré \(p\) (méthodes de type Adams-Bashforth et Adams-Moulton) ; dans ces schémas \(a_j=0\) pour \(j\neq0\);
- méthodes de Nyström et de Milne-Simpson : intégration entre \([t_{n-1},t_{n+1}]\) d’une approximation de la fonction \(\varphi\) par un polynôme de degré \(p\) ; dans ces schémas \(a_j=0\) pour \(j\neq1\);
- dérivation d’une approximation polynomiale de la fonction \(y\) puis évaluation en \(t_{n+1}\) (méthodes de type BDF); dans ces schémas \(b_j=0\) pour \(j\neq -1\).
On verra ensuite une approche dite “predictor-corrector” qui combine deux schémas multi-pas linéaires (l’un implicite, l’autre explicite) permettant la construction de schémas d’ordre plus élevé tout en gardant le caractère explicite.
Après avoir introduit les notions de convergence, consistance et stabilité(s) (dans les prochains chapitres), nous pourrons enfin étudier l’ordre de convergence et la stabilité de ces schémas. Dans les TD, nous étudierons des schémas multi-pas linéaires à \(q\) pas plus généraux, sans nous intéresser à leur construction, mais seulement à leur analyse de convergence et stabilité.
Sympy pour la construction de schémas multi-pas linéaires
Dans la suite, nous allons construire explicitement plusieurs schémas multi-pas linéaires à \(q\) pas pour l’approximation du problème de Cauchy. Cette construction est basée sur l’interpolation de fonction \(\varphi\) par un polynôme et son intégration sur un intervalle donné. Pour vérifier nos calculs d’interpolation puis intégration, nous allons utiliser le package de calcul formel SymPy. Nous allons donc l’importer et définir les variables symboliques nécessaires.
Code
# Affichage des sorties en Latex + Markdown
from IPython.display import display, Math, Markdown, Latex
# Calcul symbolique
import sympy as symb
symb.init_printing(use_latex='mathjax', use_unicode=True)
h = symb.Symbol('h') # pas de discrétisation
t = symb.Symbol('t') # temps continu utilisé pour les polynômes d'interpolation
f = symb.Symbol(r'\tilde f(t)') # polynĂ´me d'interpolation
t_n = symb.Symbol('t_n') # temps discret
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
t_nm6 = t_n - 6 * h
t_nm7 = t_n - 7 * 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_nm6=symb.Symbol('y_{n-6}')
y_nm7=symb.Symbol('y_{n-7}')
y_np1=symb.Symbol('y_{n+1}')
# Pour Adams et N-MS
phi_np1=symb.Symbol(r'\varphi_{n+1}')
phi_n =symb.Symbol(r'\varphi_{n}')
phi_nm1=symb.Symbol(r'\varphi_{n-1}')
phi_nm2=symb.Symbol(r'\varphi_{n-2}')
phi_nm3=symb.Symbol(r'\varphi_{n-3}')
phi_nm4=symb.Symbol(r'\varphi_{n-4}')
phi_nm5=symb.Symbol(r'\varphi_{n-5}')
phi_nm6=symb.Symbol(r'\varphi_{n-6}')
phi_nm7=symb.Symbol(r'\varphi_{n-7}')
# 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}')
u_nm6=symb.Symbol('u_{n-6}')
u_nm7=symb.Symbol('u_{n-7}')
Schémas de Adam : intégration sur \([t_{n},t_{n+1}]\)
Idée : intégrer l’EDO \(y'(t)=\varphi(t,y(t))\) entre \(t_n\) et \(t_{n+1}\) puis approcher l’intégrale par une formule de quadrature
\(
y_{n+1}-y_n=\int_{t_n}^{t_{n+1}} \varphi(t,y(t))dt \approx \int_{t_n}^{t_{n+1}} \tilde f(t) \mathrm{d}t
\quad \text{oĂą }\tilde f(t) \approx \varphi(t,y(t)).
\)
Les schémas d’Adams choisissent \(\tilde f(t)\) comme un polynôme interpolant \(\varphi\) en des points donnés qui peuvent être à l’extérieur de l’intervalle d’intégration \([t_{n};t_{n+1}].\)
Schémas de Adams à \(q\) pas
Soit \(p\ge0\) un entier. Un schéma de Adams à \(q=p+1\) pas (= étapes) est un cas particulier de schéma multi-pas linéaire
\(\begin{cases}
\left. \begin{array}{l}
u_0=y_0, \text{ connu}\\
u_1 \approx y_1,\\
\vdots\\
u_p\approx y_p,
\end{array}\right] \text{ initialisation}\\
\displaystyle
u_{n+1}=u_n+h\sum_{j=-1}^{p}b_j\varphi_{n-j} \qquad n=p,p+1,\dots,N-1
\end{cases}
\)
Autrement dit, un schéma de Adams est un schéma multi-pas linéaire à \(q\) pas où
- les coefficients \(a_j\) sont tous nuls pour \(j\ge1\),
- les \(\{b_j\}_{j=-1}^p\) sont \(p+1\) coefficients donnés déterminant le schéma.
Sous-classes de schémas de Adams
Ils se divisent en deux familles :
- schémas explicites dits d’Adams-Bashforth (AB-\(q\)) où \(b_{-1}=0\);
- schémas implicites dits d’Adams-Moulton (AM-\(q\)) où \(b_{-1}\ne0\) (avec un abus de notation car AM-\(0\) est toujours un schéma à un pas).
Caveat : la notation AM-\(k\) peut varier selon les auteurs. Par exemple, dans les livres de Quarteroni et al., le terme \(k\) indique l’ordre de la méthode et non le nombre de pas. Les deux notations coïncident pour les méthodes AB, mais diffèrent pour les méthodes AM.
Remarques
Polynômes d’interpolation :
- dans les schémas explicites (AB-\(q\)), \(\tilde f\) est le polynôme \(p\in\mathbb{R}^p[t]\) interpolant \(\varphi\) en les \(p+1\) points \(\{t_n,t_{n-1},\dots,t_{n-q+1}\}\) avec \(q=p+1\ge1\);
- dans les schémas implicites (AM-\(q\)), \(\tilde f\) est le polynôme \(p\in\mathbb{R}^{p+1}[t]\) interpolant \(\varphi\) en les \(p+2\) points \(\{t_{n+1},t_n,t_{n-1},\dots,t_{n-q+1}\}\) avec \(q=p+1\ge0\).
Noter que, pour calculer successivement \(u_{p+1}\), \(u_{p+2}\), …, il est nécessaire d’initialiser d’abord \(u_0, u_1, \dots, u_{p}\). Cette initialisation doit être effectuée par des approximations appropriées qui ne dégradent pas l’ordre de la méthode.
Nous verrons ultérieurement que :
- la condition de zéro-stabilité est toujours satisfaite ;
- la condition de A-stabilité devient de plus en plus contraignante avec l’augmentation du nombre de pas. En particulier, une méthode explicitement A-stable est d’ordre au plus 2 (seconde barrière de Dahlquist).
Construction des schémas de Adams-Bashford (AB-q)
Construction du schéma de Adams-Bashforth à 1 pas (AB-1)
\(q=1\) pas donc \(p=q-1=0\)
Points d’interpolation : (un seul point) \($ \{ (t_n,\varphi_n) \} \)
Polynôme d’interpolation : (degré 0) \( \tilde f(t)=\varphi_n \)
Approximation : (coĂŻncide avec la formule du rectangle Ă gauche) \(
y_{n+1}-y_n \approx \int_{t_n}^{t_{n+1}} \tilde f(t) \mathrm{d}t =h\varphi_n
\)
Code
pol = symb.interpolate([(t_n, phi_n)], t)
display(symb.Eq(f, pol))
integrale = symb.integrate(pol, (t, t_n, t_np1)).simplify()
display(symb.Eq(u_np1, u_n + integrale))
\(\displaystyle \tilde f(t) = \varphi_{n}\)
\(\displaystyle u_{n+1} = \varphi_{n} h + u_{n}\)
Ainsi, on obtient le schéma :
\(
\text{($\text{AB}_1$)}\quad
\begin{cases}
u_0 = y_0,\\
u_{n+1} = u_{n}+h \varphi_n, & n=0,1,\dots,N-1.
\end{cases}
\)
La méthode \(\text{AB}_1\) coïncide donc avec la méthode d’Euler explicite.
Construction du schéma de Adam-Bashforth à 2 pas (AB-2)
\(q=2\) pas donc \(p=q-1=1\)
Points d’interpolation : (deux points) \( \{(t_n,\varphi_n);(t_{n-1},\varphi_{n-1})\} \)
Polynôme d’interpolation : (degré 1) \( \tilde f(t)=\varphi_n\frac{t-t_{n-1}}{t_n-t_{n-1}}+\varphi_{n-1}\frac{t-t_n}{t_{n-1}-t_n}=\varphi_n\frac{t-t_{n-1}}{h}+\varphi_{n-1}\frac{t-t_n}{-h} \)
Approximation : \(y_{n+1}-y_n \approx \int_{t_n}^{t_{n+1}} \tilde f(t) \mathrm{d}t =\frac{h}{2}\left(3\varphi_n-\varphi_{n-1}\right)\)
Code
pol = symb.interpolate([(t_n, phi_n),
(t_nm1, phi_nm1)],
t)
display(symb.Eq(f, pol.factor(t)))
integrale = symb.integrate(pol, (t, t_n, t_np1)).simplify()
display(symb.Eq(u_np1, u_n + integrale))
\(\displaystyle \tilde f(t) = - \frac{- \varphi_{n-1} t_{n} - \varphi_{n} h + \varphi_{n} t_{n} + t \left(\varphi_{n-1} - \varphi_{n}\right)}{h}\)
\(\displaystyle u_{n+1} = \frac{h \left(- \varphi_{n-1} + 3 \varphi_{n}\right)}{2} + u_{n}\)
Ainsi, on obtient le schéma :
\(
\text{($\text{AB}_2$)}\quad
\begin{cases}
u_0=y(t_0)=y_0,\\
u_1=u_0+h\varphi_0\approx y(t_1),\\
u_{n+1}=u_{n}+ \frac{h}{2} \big(3\varphi_n-\varphi_{n-1}\big), & n=1,2,\dots,N-1.
\end{cases}
\)
où \(u_{1}\) est une approximation de \(y(t_1)\) obtenue en utilisant ici une prédiction \(\text{AB}_1\).
Construction du schéma de Adams-Bashforth à 3 pas (AB-3)
Points d’interpolation : (trois points) \(\{(t_n,\varphi_n);(t_{n-1},\varphi_{n-1});(t_{n-2},\varphi_{n-2})\}\)
Polynôme d’interpolation : (degré 2) \(
\begin{aligned}
\tilde f(t)&=\varphi_n\frac{(t-t_{n-1})(t-t_{n-2})}{(t_n-t_{n-1})(t_n-t_{n-2})}
+\varphi_{n-1}\frac{(t-t_{n})(t-t_{n-2})}{(t_{n-1}-t_{n})(t_{n-1}-t_{n-2})}
+\varphi_{n-2}\frac{(t-t_{n})(t-t_{n-1})}{(t_{n-2}-t_{n})(t_{n-2}-t_{n-1})} \\
&=\frac{\varphi_{n-2}}{2h^2}(t-t_{n-1})(t-t_{n})
-\frac{\varphi_{n-1}}{h^2}(t-t_{n-2})(t-t_{n})
+\frac{\varphi_n}{2h^2}(t-t_{n-2})(t-t_{n-1})
\end{aligned}
\)
Approximation : \(y_{n+1}-y_n \approx \int_{t_n}^{t_{n+1}} \tilde f(t) \mathrm{d}t =\frac{h}{12}\left(23\varphi_n -16\varphi_{n-1} +5\varphi_{n-2}\right)\)
Code
pol = symb.interpolate([(t_n, phi_n),
(t_nm1, phi_nm1),
(t_nm2, phi_nm2)],
t)
display(symb.Eq(f, pol.factor(t)))
integrale = symb.integrate(pol, (t, t_n, t_np1)).simplify()
display(symb.Eq(u_np1, u_n + integrale))
\(\displaystyle \tilde f(t) = - \frac{- 4 \varphi_{n-1} h t_{n} + 2 \varphi_{n-1} t_{n}^{2} + \varphi_{n-2} h t_{n} - \varphi_{n-2} t_{n}^{2} - 2 \varphi_{n} h^{2} + 3 \varphi_{n} h t_{n} - \varphi_{n} t_{n}^{2} + t^{2} \cdot \left(2 \varphi_{n-1} - \varphi_{n-2} - \varphi_{n}\right) + t \left(4 \varphi_{n-1} h - 4 \varphi_{n-1} t_{n} - \varphi_{n-2} h + 2 \varphi_{n-2} t_{n} - 3 \varphi_{n} h + 2 \varphi_{n} t_{n}\right)}{2 h^{2}}\)
\(\displaystyle u_{n+1} = \frac{h \left(- 16 \varphi_{n-1} + 5 \varphi_{n-2} + 23 \varphi_{n}\right)}{12} + u_{n}\)
On obtient le schéma
\(
\text{($\text{AB}_3$)}\quad
\begin{cases}
u_0 = y(t_0) = y_0,\\
u_1 = u_0 + h\varphi_0 \approx y(t_1),\\
u_2 = u_1 + \dfrac{h}{2} \big(3\varphi_1 - \varphi_0 \big) \approx y(t_2),\\
u_{n+1} = u_{n} + \dfrac{h}{12} \big(23\varphi_n -16\varphi_{n-1} +5\varphi_{n-2} \big), & n=2,3,\dots,N-1.
\end{cases}
\)
où \(u_{1}\) est une approximation de \(y(t_1)\) obtenue en utilisant une prédiction \(\text{AB}_1\) et \(u_{2}\) est une approximation de \(y(t_2)\) obtenue en utilisant la méthode \(\text{AB}_2\).
Construction du schéma de Adams-Bashforth à 4 pas (AB-4)
Points d’interpolation : (quatre points) \(\{(t_n,\varphi_n);(t_{n-1},\varphi_{n-1});(t_{n-2},\varphi_{n-2});(t_{n-3},\varphi_{n-3})\}\)
Polynôme d’interpolation : (degré 3) \(
\tilde f(t)=\sum_{i=n-3}^{n} \left(\varphi_i \prod_{\substack{j=n-2,\dots,n\\j\neq i}} \frac{t-t_j}{t_i-t_j} \right)
\)
Approximation :
\(y_{n+1}-y_n \approx \int_{t_n}^{t_{n+1}} \tilde f(t) \mathrm{d}t =\frac{h}{24}\left(55\varphi_n - 59\varphi_{n-1} + 37\varphi_{n-2} - 9\varphi_{n-3}\right)\)
Code
pol = symb.interpolate([(t_n, phi_n),
(t_nm1, phi_nm1),
(t_nm2, phi_nm2),
(t_nm3, phi_nm3)],
t)
# display(symb.Eq(f, pol.factor(t)))
integrale = symb.integrate(pol, (t, t_n, t_np1)).simplify()
display(symb.Eq(u_np1, u_n + integrale))
\(\displaystyle u_{n+1} = \frac{h \left(- 59 \varphi_{n-1} + 37 \varphi_{n-2} - 9 \varphi_{n-3} + 55 \varphi_{n}\right)}{24} + u_{n}\)
On obtient le schéma
\(
\text{($\text{AB}_4$)}\quad
\begin{cases}
u_0 = y(t_0) = y_0,\\
u_1 = u_0 + h\varphi_0 \approx y(t_1),\\
u_2 = u_1 + \dfrac{h}{2} \big(3\varphi_1 - \varphi_0 \big) \approx y(t_2),\\
u_3 = u_2 + \dfrac{h}{12} \big(23\varphi_2 -16\varphi_1 +5\varphi_0 \big) \approx y(t_3),\\
u_{n+1} = u_n + \dfrac{h}{24} \big(55\varphi_n - 59\varphi_{n-1} + 37\varphi_{n-2} - 9\varphi_{n-3} \big), & n=3,4,\dots,N-1.
\end{cases}
\)
où \(u_{1}\) est une approximation de \(y(t_1)\) obtenue en utilisant une prédiction \(\text{AB}_1\), \(u_{2}\) est une approximation de \(y(t_2)\) obtenue en utilisant la méthode \(\text{AB}_2\) etc.
Construction des schémas de Adams-Moulton (AM-q)
Construction du schéma de Adams-Moulton à 0 pas (AM-0) [abus de notation car AM-0 est toujours un schéma à un pas]
Points d’interpolation : (un seul point) \(\{(t_{n+1},\varphi_{n+1})\}\)
Polynôme d’interpolation : (degré 0) \(\tilde f(t) = \varphi_{n+1}\)
Approximation : (coĂŻncide avec la formule du rectangle Ă droite) \(y_{n+1} - y_n \approx \int_{t_n}^{t_{n+1}} \tilde f(t) \mathrm{d}t = h\varphi_{n+1}\)
Code
pol = symb.interpolate([(t_np1, phi_np1)], t)
display(symb.Eq(f, pol.factor(t)))
integrale = symb.integrate(pol, (t, t_n, t_np1)).simplify()
display(symb.Eq(u_np1, u_n + integrale))
\(\displaystyle \tilde f(t) = \varphi_{n+1}\)
\(\displaystyle u_{n+1} = \varphi_{n+1} h + u_{n}\)
Ainsi, on obtient le schéma
\(
\text{($\text{AM}_0$)}\quad
\begin{cases}
u_0 = y(t_0) = y_0,\\
u_{n+1} = u_n + h \varphi_{n+1}, & n=0,1,\dots,N-1.
\end{cases}
\)
La méthode \(\text{AM}_0\) coïncide avec la méthode d’Euler implicite.
Construction du schéma de Adam-Multon à 1 pas (AM-1)
Points d’interpolation : (deux points) \(\{(t_{n+1},\varphi_{n+1}),(t_n,\varphi_n)\}\)
Polynôme d’interpolation : (degré 1)
\(\tilde f(t) = \varphi_{n+1} \frac{t - t_n}{h} + \varphi_n \frac{t - t_{n+1}}{-h}\)
Approximation : (coïncide avec la formule du trapèze)
\(y_{n+1} - y_n \approx \int_{t_n}^{t_{n+1}} \tilde f(t) \mathrm{d}t = \frac{h}{2} (\varphi_{n+1} + \varphi_n)\)
Code
pol = symb.interpolate([(t_np1, phi_np1),
(t_n, phi_n)],
t)
display(symb.Eq(f, pol.factor(t)))
integrale = symb.integrate(pol, (t, t_n, t_np1)).simplify()
display(symb.Eq(u_np1, u_n + integrale))
\(\displaystyle \tilde f(t) = \frac{- \varphi_{n+1} t_{n} + \varphi_{n} h + \varphi_{n} t_{n} + t \left(\varphi_{n+1} - \varphi_{n}\right)}{h}\)
\(\displaystyle u_{n+1} = \frac{h \left(\varphi_{n+1} + \varphi_{n}\right)}{2} + u_{n}\)
Ainsi, on obtient le schéma
\(
\text{($\text{AM}_1$)}\quad
\begin{cases}
u_0 = y(t_0) = y_0,\\
u_{n+1} = u_n + \dfrac{h}{2} \big(\varphi_n + \varphi_{n+1}\big), & n=0,1,\dots,N-1.
\end{cases}
\)
La méthode \(\text{AM}_1\) coïncide avec la méthode de Crank-Nicolson.
Construction du schéma de Adam-Multon à 2 pas (AM-2)
Points d’interpolation : (trois points) \(\{(t_{n+1},\varphi_{n+1}),(t_n,\varphi_n),(t_{n-1},\varphi_{n-1})\}\)
Polynôme d’interpolation : (degré 2) \(\tilde f(t) = \varphi_{n+1} \frac{(t - t_n)(t - t_{n-1})}{(t_{n+1} - t_n)(t_{n+1} - t_{n-1})}
+ \varphi_n \frac{(t - t_{n+1})(t - t_{n-1})}{(t_n - t_{n+1})(t_n - t_{n-1})}
+ \varphi_{n-1} \frac{(t - t_{n+1})(t - t_n)}{(t_{n-1} - t_{n+1})(t_{n-1} - t_n)}\)
Approximation :
\(y_{n+1} - y_n \approx \int_{t_n}^{t_{n+1}} \tilde f(t) \mathrm{d}t = \frac{h}{12} \left(5\varphi_{n+1} + 8\varphi_n - \varphi_{n-1}\right)\)
Code
pol = symb.interpolate([(t_np1, phi_np1),
(t_n, phi_n),
(t_nm1, phi_nm1)],
t)
display(symb.Eq(f, pol.factor(t)))
integrale = symb.integrate(pol, (t, t_n, t_np1)).simplify()
display(symb.Eq(u_np1, u_n + integrale))
\(\displaystyle \tilde f(t) = \frac{- \varphi_{n+1} h t_{n} + \varphi_{n+1} t_{n}^{2} + \varphi_{n-1} h t_{n} + \varphi_{n-1} t_{n}^{2} + 2 \varphi_{n} h^{2} - 2 \varphi_{n} t_{n}^{2} + t^{2} \left(\varphi_{n+1} + \varphi_{n-1} - 2 \varphi_{n}\right) + t \left(\varphi_{n+1} h - 2 \varphi_{n+1} t_{n} - \varphi_{n-1} h - 2 \varphi_{n-1} t_{n} + 4 \varphi_{n} t_{n}\right)}{2 h^{2}}\)
\(\displaystyle u_{n+1} = \frac{h \left(5 \varphi_{n+1} - \varphi_{n-1} + 8 \varphi_{n}\right)}{12} + u_{n}\)
On obtient le schéma
\(
\text{($\text{AM}_2$)}\quad
\begin{cases}
u_0 = y(t_0) = y_0,\\
u_1 = u_0 + \dfrac{h}{2} \big(\varphi_1 + \varphi_0\big),\\
u_{n+1} = u_n + \dfrac{h}{12} \big(5\varphi_{n+1} + 8\varphi_n - \varphi_{n-1}\big), & n=1,2,\dots,N-1.
\end{cases}
\)
où \(u_{1}\) est une approximation de \(y(t_1)\) obtenue en utilisant une prédiction \(\text{AM}_1\).
Automatisation de la construction des schémas de Adam
Nous allons maintenant automatiser le calcul des coefficients des schémas de Adam en utilisant la bibliothèque SymPy.
Code
display(Markdown(f"#### Adam-Bashforth (explictes)"))
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), (t_nm6, phi_nm6), (t_nm7, phi_nm7)]
for q in range(1, len(Points) + 1):
pol = symb.interpolate(Points[:q], t)
AB = (symb.integrate(pol, (t, t_n, t_np1)).simplify())
display(Markdown("- "+f"**AB-{q}: explicite, Ă {q} pas**"))
display(symb.Eq(u_np1, u_n + AB))
#display(Markdown("---"))
display(Markdown(f"#### Adam-Moulton (implicites)"))
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), (t_nm6, phi_nm6), (t_nm7, phi_nm7)]
for q in range(len(Points)):
pol = symb.interpolate(Points[0:q + 1], t)
AM = (symb.integrate(pol, (t, t_n, t_np1)).simplify())
display(Markdown("- "+f"**AM-{q}: implicite, Ă {q} pas**"))
display(symb.Eq(u_np1, u_n + AM))
#(Markdown("---"))
Adam-Bashforth (explictes)
\(\displaystyle u_{n+1} = \varphi_{n} h + u_{n}\)
\(\displaystyle u_{n+1} = \frac{h \left(- \varphi_{n-1} + 3 \varphi_{n}\right)}{2} + u_{n}\)
\(\displaystyle u_{n+1} = \frac{h \left(- 16 \varphi_{n-1} + 5 \varphi_{n-2} + 23 \varphi_{n}\right)}{12} + u_{n}\)
\(\displaystyle u_{n+1} = \frac{h \left(- 59 \varphi_{n-1} + 37 \varphi_{n-2} - 9 \varphi_{n-3} + 55 \varphi_{n}\right)}{24} + u_{n}\)
\(\displaystyle u_{n+1} = \frac{h \left(- 2774 \varphi_{n-1} + 2616 \varphi_{n-2} - 1274 \varphi_{n-3} + 251 \varphi_{n-4} + 1901 \varphi_{n}\right)}{720} + u_{n}\)
\(\displaystyle u_{n+1} = \frac{h \left(- 7923 \varphi_{n-1} + 9982 \varphi_{n-2} - 7298 \varphi_{n-3} + 2877 \varphi_{n-4} - 475 \varphi_{n-5} + 4277 \varphi_{n}\right)}{1440} + u_{n}\)
\(\displaystyle u_{n+1} = \frac{h \left(- 447288 \varphi_{n-1} + 705549 \varphi_{n-2} - 688256 \varphi_{n-3} + 407139 \varphi_{n-4} - 134472 \varphi_{n-5} + 19087 \varphi_{n-6} + 198721 \varphi_{n}\right)}{60480} + u_{n}\)
\(\displaystyle u_{n+1} = \frac{h \left(- 1152169 \varphi_{n-1} + 2183877 \varphi_{n-2} - 2664477 \varphi_{n-3} + 2102243 \varphi_{n-4} - 1041723 \varphi_{n-5} + 295767 \varphi_{n-6} - 36799 \varphi_{n-7} + 434241 \varphi_{n}\right)}{120960} + u_{n}\)
Adam-Moulton (implicites)
\(\displaystyle u_{n+1} = \varphi_{n+1} h + u_{n}\)
\(\displaystyle u_{n+1} = \frac{h \left(\varphi_{n+1} + \varphi_{n}\right)}{2} + u_{n}\)
\(\displaystyle u_{n+1} = \frac{h \left(5 \varphi_{n+1} - \varphi_{n-1} + 8 \varphi_{n}\right)}{12} + u_{n}\)
\(\displaystyle u_{n+1} = \frac{h \left(9 \varphi_{n+1} - 5 \varphi_{n-1} + \varphi_{n-2} + 19 \varphi_{n}\right)}{24} + u_{n}\)
\(\displaystyle u_{n+1} = \frac{h \left(251 \varphi_{n+1} - 264 \varphi_{n-1} + 106 \varphi_{n-2} - 19 \varphi_{n-3} + 646 \varphi_{n}\right)}{720} + u_{n}\)
\(\displaystyle u_{n+1} = \frac{h \left(475 \varphi_{n+1} - 798 \varphi_{n-1} + 482 \varphi_{n-2} - 173 \varphi_{n-3} + 27 \varphi_{n-4} + 1427 \varphi_{n}\right)}{1440} + u_{n}\)
\(\displaystyle u_{n+1} = \frac{h \left(19087 \varphi_{n+1} - 46461 \varphi_{n-1} + 37504 \varphi_{n-2} - 20211 \varphi_{n-3} + 6312 \varphi_{n-4} - 863 \varphi_{n-5} + 65112 \varphi_{n}\right)}{60480} + u_{n}\)
\(\displaystyle u_{n+1} = \frac{h \left(36799 \varphi_{n+1} - 121797 \varphi_{n-1} + 123133 \varphi_{n-2} - 88547 \varphi_{n-3} + 41499 \varphi_{n-4} - 11351 \varphi_{n-5} + 1375 \varphi_{n-6} + 139849 \varphi_{n}\right)}{120960} + u_{n}\)
\(\displaystyle u_{n+1} = \frac{h \left(1070017 \varphi_{n+1} - 4604594 \varphi_{n-1} + 5595358 \varphi_{n-2} - 5033120 \varphi_{n-3} + 3146338 \varphi_{n-4} - 1291214 \varphi_{n-5} + 312874 \varphi_{n-6} - 33953 \varphi_{n-7} + 4467094 \varphi_{n}\right)}{3628800} + u_{n}\)
Schémas de Nyström et de Milne-Simpson : intégration sur \([t_{n-1},t_{n+1}]\)
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\).
Idée : choisir \(r=1\); intégrer l’EDO \(y'(t)=\varphi(t,y(t))\) entre \(t_{n-1}\) et \(t_{n+1}\) puis approcher l’intégrale par une formule de quadrature
\(
y_{n+1}-y_{n-1}=\int_{t_{n-1}}^{t_{n+1}} \varphi(t,y(t))dt \approx \int_{t_{n-1}}^{t_{n+1}} \tilde f(t) \mathrm{d}t
\quad \text{oĂą }\tilde f(t) \approx \varphi(t,y(t)).
\)
À nouveau, \(\tilde f(t)\) sera un polynôme 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}]\).
Schémas de Nyström et de Milne-Simpson à \(q\ge2\) pas
Soit \(p\ge1\) un entier. Un schéma de Nystrom ou de Milne-Simpson à \(q=p+1\ge2\) pas (= étapes) est un cas particulier de schéma multi-pas linéaire
\(\begin{cases}
\left. \begin{array}{l}
u_0=y_0, \text{ connu}\\
u_1 \approx y_1,\text{ tj présent}\\
\vdots\\
u_p\approx y_p,
\end{array}\right] \text{ initialisation}\\
\displaystyle
u_{n+1}=u_{n-1}+h\sum_{j=-1}^{p}b_j\varphi_{n-j} \qquad n=p,p+1,\dots,N-1
\end{cases}
\)
Autrement dit, un schéma multi-pas linéaire à \(q\) pas où
- les coefficients \(a_j\) sont tous nuls pour \(j\neq1\),
- les \(\{b_j\}_{j=-1}^p\) sont \(p+1\) coefficients donnés déterminant le schéma,
- \(u_1\) doit toujours être calculé par une méthode de prédiction.
Sous-classes de schémas
- Schémas explicites dits de Nyström (N-\(q\)) où \(b_{-1}=0\);
- schémas implicites dits de Milne-Simpson (MS-\(q\)) où \(b_{-1}\ne0\).
Remarques
- Polynômes d’interpolation :
- dans les schémas explicites (N-\(q\)), \(\tilde f\) est le polynôme \(p\in\mathbb{R}^p[t]\) interpolant \(\varphi\) en les \(p+1\) points \(\{t_n,t_{n-1},\dots,t_{n-q+1}\}\) avec \(q=p+1\ge1\);
- dans les schémas implicites (MS-\(q\)), \(\tilde f\) est le polynôme \(p\in\mathbb{R}^{p+1}[t]\) interpolant \(\varphi\) en les \(p+2\) points \(\{t_{n+1},t_n,t_{n-1},\dots,t_{n-q+1}\}\) avec \(q=p+1\ge0\).
- Quel que soit le schéma, il est nécessaire d’initialiser \(u_0\) et \(u_1\). Cette initialisation doit être effectuée par des approximations appropriées qui ne dégradent pas l’ordre de la méthode.
Construction des schémas de Nyström (N-q)
Construction du schéma de Nyström à 1 pas (N-1) [abus de notation, car N-1 est en réalité un schéma à 2 pas]
Points d’interpolation: (un seul point) \(\{(t_n,\varphi_n)\}\)
Polynôme: (degré 0) \(
\tilde f(t) = \varphi_n
\)
Approximation: \(
y_{n+1}-y_{n-1} \approx \int_{t_{n-1}}^{t_{n+1}} \tilde f(t) \, \mathrm{d}t = 2h \varphi_n
\)
Code
pol = symb.interpolate([(t_n, phi_n)], t).simplify()
display(symb.Eq(f, pol))
integrale = symb.integrate(pol, (t, t_nm1, t_np1)).simplify()
display(symb.Eq(u_np1, u_nm1 + integrale))
\(\displaystyle \tilde f(t) = \varphi_{n}\)
\(\displaystyle u_{n+1} = 2 \varphi_{n} h + u_{n-1}\)
On obtient le schéma
\(
\text{(N$_1$)}\quad
\begin{cases}
u_0 = y(t_0) = y_0, \\
u_1 = u_0 + h \varphi_0 \approx y(t_1) \\
u_{n+1} = u_{n-1} + 2h \varphi_n & n = 1, 2, \dots N-1
\end{cases}
\)
où \(u_{1}\) est une approximation de \(y(t_1)\) obtenue en utilisant une prédiction d’Euler explicite.
Remarque : on l’appelle N-1 mais c’est un abus de notation car elle est à 2 pas (elle dépend de \(u_{n-1}\) directement et non pas à travers \(\varphi\)).
Construction du schéma de Nyström à 2 pas (N-2)
Points d’interpolation : (deux points) \(\{(t_n,\varphi_n);(t_{n-1},\varphi_{n-1})\}\)
Polynôme : (degré 1) \(
\tilde f(t) = \varphi_n \frac{t - t_{n-1}}{h} + \varphi_{n-1} \frac{t - t_n}{-h}
\)
Approximation : (coïncide avec la formule d’approximation trouvée précédemment) \(y_{n+1}-y_{n-1} \approx \int_{t_{n-1}}^{t_{n+1}} \tilde f(t) \mathrm{d}t = 2h\varphi_n\)
Code
pol = symb.interpolate([(t_n, phi_n),
(t_nm1, phi_nm1)],
t)
display(symb.Eq(f, pol.factor(t)))
integrale = symb.integrate(pol, (t, t_nm1, t_np1)).simplify()
display(symb.Eq(u_np1, u_nm1 + integrale))
\(\displaystyle \tilde f(t) = - \frac{- \varphi_{n-1} t_{n} - \varphi_{n} h + \varphi_{n} t_{n} + t \left(\varphi_{n-1} - \varphi_{n}\right)}{h}\)
\(\displaystyle u_{n+1} = 2 \varphi_{n} h + u_{n-1}\)
On obtient le schéma
\(
\text{(N$_2$)}\equiv\text{(N$_1$)}\quad
\begin{cases}
u_0=y(t_0)=y_0,\\
u_1 = u_0 + h \varphi_0 \approx y(t_1) \\
u_{n+1}=u_{n-1}+2h \varphi_{n}& n=1,2,\dots N-1
\end{cases}\)
où \(u_{1}\) est une approximation de \(y(t_1)\) obtenue en utilisant une prédiction d’Euler explicite.
Cette méthode coïncide avec la méthode du point milieu (aussi appelée Saute-mouton ou Leapfrog).
Construction du schéma de Nyström à 3 pas (N-3)
Points d’interpolation : (trois points) \(\{(t_n,\varphi_n);(t_{n-1},\varphi_{n-1});(t_{n-2},\varphi_{n-2})\}\)
Polynôme : (degré 2) \(\begin{aligned}
\tilde f(t) &= \varphi_n \frac{(t - t_{n-1})(t - t_{n-2})}{(t_n - t_{n-1})(t_n - t_{n-2})}
{}+ \varphi_{n-1} \frac{(t - t_n)(t - t_{n-2})}{(t_{n-1} - t_n)(t_{n-1} - t_{n-2})}
{}+ \varphi_{n-2} \frac{(t - t_n)(t - t_{n-1})}{(t_{n-2} - t_n)(t_{n-2} - t_{n-1})} \\
&= \frac{\varphi_{n-2}}{2h^2}(t - t_{n-1})(t - t_n) - \frac{\varphi_{n-1}}{h^2}(t - t_{n-2})(t - t_n) + \frac{\varphi_n}{2h^2}(t - t_{n-2})(t - t_{n-1})
\end{aligned}\)
Approximation : \(
y_{n+1} - y_{n-1} \approx \int_{t_{n-1}}^{t_{n+1}} \tilde f(t) \, \mathrm{d}t = \frac{h}{3}\left(7 \varphi_n - 2 \varphi_{n-1} + \varphi_{n-2}\right)
\)
Code
pol = symb.interpolate([(t_n, phi_n), (t_nm1, phi_nm1), (t_nm2, phi_nm2)],t).simplify()
display(symb.Eq(f, pol))
integrale = symb.integrate(pol, (t, t_nm1, t_np1)).simplify()
display(symb.Eq(u_np1, u_nm1 + integrale))
\(\displaystyle \tilde f(t) = \frac{- \varphi_{n-1} t^{2} + 2 \varphi_{n-1} t t_{n} - \varphi_{n-1} t_{n}^{2} + \frac{\varphi_{n-2} t^{2}}{2} - \varphi_{n-2} t t_{n} + \frac{\varphi_{n-2} t_{n}^{2}}{2} + \varphi_{n} h^{2} + \frac{\varphi_{n} t^{2}}{2} - \varphi_{n} t t_{n} + \frac{\varphi_{n} t_{n}^{2}}{2} + \frac{h \left(- 4 \varphi_{n-1} t + 4 \varphi_{n-1} t_{n} + \varphi_{n-2} t - \varphi_{n-2} t_{n} + 3 \varphi_{n} t - 3 \varphi_{n} t_{n}\right)}{2}}{h^{2}}\)
\(\displaystyle u_{n+1} = \frac{h \left(- 2 \varphi_{n-1} + \varphi_{n-2} + 7 \varphi_{n}\right)}{3} + u_{n-1}\)
On obtient le schéma
\(
\text{(N$_3$)}\quad
\begin{cases}
u_0 = y(t_0) = y_0, \\
u_1 = u_0 + h \varphi_0 \approx y(t_1) \\
u_2 = u_0 + 2h \varphi_1 \approx y(t_2) \\
u_{n+1} = u_{n-1} + \frac{h}{3}\left(7 \varphi_n - 2 \varphi_{n-1} + \varphi_{n-2}\right) & n = 2, 3, \dots N-1
\end{cases}
\)
où \(u_{1}\) est une approximation de \(y(t_1)\) obtenue en utilisant une prédiction d’Euler explicite et \(u_{2}\) est une approximation de \(y(t_2)\) obtenue en utilisant la méthode N-2.
Construction des schémas de Milne-Simpson (MS-q)
Construction du schéma de Milne-Simpson à 0 pas (MS-0) [abus de notation car MS-0 est toujours un schéma à 2 pas]
Points d’interpolation : (un seul point) \(
\{(t_{n+1}, \varphi_{n+1})\}
\)
Polynôme : (degré 0) \(
\tilde f(t) = \varphi_{n+1}
\)
Approximation : (coĂŻncide avec la formule du rectangle Ă droite) \(
y_{n+1} - y_{n-1} \approx \int_{t_{n-1}}^{t_{n+1}} \tilde f(t) \, \mathrm{d}t = 2h \varphi_{n+1}
\)
Code
pol = symb.interpolate([(t_np1, phi_np1)], t).simplify()
display(symb.Eq(f, pol))
integrale = symb.integrate(pol, (t, t_nm1, t_np1)).simplify()
display(symb.Eq(u_np1, u_nm1 + integrale))
\(\displaystyle \tilde f(t) = \varphi_{n+1}\)
\(\displaystyle u_{n+1} = 2 \varphi_{n+1} h + u_{n-1}\)
On obtient le schéma
\(
\text{(MS$_0$)}\quad
\begin{cases}
u_0 = y(t_0) = y_0,\\
u_1 = u_0 + h \varphi_1 \approx y(t_1)\\
u_{n+1} = u_{n-1} + 2h \varphi_{n+1} & n = 1, 2, \dots, N-1
\end{cases}
\)
où \(u_{1}\) est une approximation de \(y(t_1)\) obtenue en utilisant une prédiction d’Euler implicite.
Remarque : on l’appelle MS-0 mais c’est une schéma à 2 pas, car il dépend de \(u_{n-1}\).
Construction du schéma de Milne-Simpson à 1 pas (MS-1) [abus de notation car MS-1 est toujours un schéma à 2 pas]
Points d’interpolation : (deux points) \(
\{(t_{n+1}, \varphi_{n+1}), (t_n, \varphi_n)\}
\)
Polynôme : (degré 1) \(
\tilde f(t) = \varphi_{n+1} \frac{t - t_n}{h} + \varphi_n \frac{t - t_{n+1}}{-h}
\)
Approximation : (coĂŻncide avec les formules N-1 et N-2) \(
y_{n+1} - y_{n-1} \approx \int_{t_{n-1}}^{t_{n+1}} \tilde f(t) \, \mathrm{d}t = 2h\varphi(t_n,y_n)
\)
Code
pol = symb.interpolate([(t_np1, phi_np1),
(t_n, phi_n)],
t).simplify()
display(symb.Eq(f, pol.factor(t)))
integrale = symb.integrate(pol, (t, t_nm1, t_np1)).simplify()
display(symb.Eq(u_np1, u_nm1 + integrale))
\(\displaystyle \tilde f(t) = \frac{- \varphi_{n+1} t_{n} + \varphi_{n} h + \varphi_{n} t_{n} + t \left(\varphi_{n+1} - \varphi_{n}\right)}{h}\)
\(\displaystyle u_{n+1} = 2 \varphi_{n} h + u_{n-1}\)
On obtient le schéma
\(
\text{(MS$_1$)} \equiv \text{(N$_1$)} \equiv \text{(N$_2$)}\quad
\begin{cases}
u_0 = y(t_0) = y_0,\\
u_1 = u_0+h\varphi_1\approx y(t_1)\\
u_{n+1} = u_{n-1} + 2h \varphi_{n}& n=1,2,\dots N-1
\end{cases}\)
où \(u_{1}\) est une approximation de \(y(t_1)\) obtenue en utilisant une prédiction d’Euler implicite.
Nota bene : le schéma est en realité explicite.
Remarque : on l’appelle MS-1 mais c’est un schéma à 2 pas, car il dépend de \(u_{n-1}\).
Construction du schéma de Milne-Simpson à 2 pas (MS-2)
Points d’interpolation : (trois points) \(
\{(t_{n+1}, \varphi_{n+1}), (t_{n}, \varphi_{n}), (t_{n-1}, \varphi_{n-1})\}
\)
Polynôme : (degré 2) \(
\tilde f(t) = \varphi_{n+1} \frac{(t - t_n)(t - t_{n-1})}{(t_{n+1} - t_n)(t_{n+1} - t_{n-1})}
+ \varphi_n \frac{(t - t_{n+1})(t - t_{n-1})}{(t_n - t_{n+1})(t_n - t_{n-1})}
+ \varphi_{n-1} \frac{(t - t_{n+1})(t - t_n)}{(t_{n-1} - t_{n+1})(t_{n+1} - t_n)}
\)
Approximation : (coĂŻncide avec la formule de Simpson sur un intervalle de largeur \(2h\)) \(
y_{n+1} - y_{n-1} \approx \int_{t_{n-1}}^{t_{n+1}} \tilde f(t) \, \mathrm{d}t = \frac{h}{3} \left( \varphi_{n+1} + 4 \varphi_n + \varphi_{n-1} \right)
\)
Code
pol = symb.interpolate([(t_np1, phi_np1),
(t_n, phi_n),
(t_nm1, phi_nm1)],
t).simplify()
display(symb.Eq(f, pol.factor(t)))
integrale = symb.integrate(pol, (t, t_nm1, t_np1)).simplify()
display(symb.Eq(u_np1, u_nm1 + integrale))
\(\displaystyle \tilde f(t) = \frac{- \varphi_{n+1} h t_{n} + \varphi_{n+1} t_{n}^{2} + \varphi_{n-1} h t_{n} + \varphi_{n-1} t_{n}^{2} + 2 \varphi_{n} h^{2} - 2 \varphi_{n} t_{n}^{2} + t^{2} \left(\varphi_{n+1} + \varphi_{n-1} - 2 \varphi_{n}\right) + t \left(\varphi_{n+1} h - 2 \varphi_{n+1} t_{n} - \varphi_{n-1} h - 2 \varphi_{n-1} t_{n} + 4 \varphi_{n} t_{n}\right)}{2 h^{2}}\)
\(\displaystyle u_{n+1} = \frac{h \left(\varphi_{n+1} + \varphi_{n-1} + 4 \varphi_{n}\right)}{3} + u_{n-1}\)
On obtient le schéma
\(
\text{(MS$_2$)}\quad
\begin{cases}
u_0=y(t_0)=y_0,\\
u_1=u_{0}+ \dfrac{h}{2} \big(\varphi_{1}+\varphi_{0})\big)\\
u_{n+1}=u_{n-1}+\dfrac{h}{3}\big( \varphi_{n+1}+4\varphi_n+\varphi_{n-1}\big) & n=1,2,\dots N-1
\end{cases}\)
où \(u_{1}\) est une approximation de \(y(t_1)\) obtenue en utilisant une prédiction \(\text{AM}_1\) (= Crank-Nicolson).
Automatisation de la construction des schémas de Nyström et de Milne-Simpson
Code
display(Markdown(f"#### Nyström"))
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())
display(Markdown("- "+f"**N-{q}: explicite, Ă {q} pas**"))
display(symb.Eq(u_np1, u_nm1 + N))
#display(Markdown("---"))
display(Markdown(f"#### Milne-Simpson"))
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())
display(Markdown("- "+f"**MS-{q}: implicite, Ă {q} pas**"))
display(symb.Eq(u_np1, u_nm1 + MS))
##display(Markdown("---"))
\(\displaystyle u_{n+1} = 2 \varphi_{n} h + u_{n-1}\)
\(\displaystyle u_{n+1} = 2 \varphi_{n} h + u_{n-1}\)
\(\displaystyle u_{n+1} = \frac{h \left(- 2 \varphi_{n-1} + \varphi_{n-2} + 7 \varphi_{n}\right)}{3} + u_{n-1}\)
\(\displaystyle u_{n+1} = \frac{h \left(- 5 \varphi_{n-1} + 4 \varphi_{n-2} - \varphi_{n-3} + 8 \varphi_{n}\right)}{3} + u_{n-1}\)
\(\displaystyle u_{n+1} = \frac{h \left(- 266 \varphi_{n-1} + 294 \varphi_{n-2} - 146 \varphi_{n-3} + 29 \varphi_{n-4} + 269 \varphi_{n}\right)}{90} + u_{n-1}\)
\(\displaystyle u_{n+1} = \frac{h \left(- 406 \varphi_{n-1} + 574 \varphi_{n-2} - 426 \varphi_{n-3} + 169 \varphi_{n-4} - 28 \varphi_{n-5} + 297 \varphi_{n}\right)}{90} + u_{n-1}\)
\(\displaystyle u_{n+1} = 2 \varphi_{n+1} h + u_{n-1}\)
\(\displaystyle u_{n+1} = 2 \varphi_{n} h + u_{n-1}\)
\(\displaystyle u_{n+1} = \frac{h \left(\varphi_{n+1} + \varphi_{n-1} + 4 \varphi_{n}\right)}{3} + u_{n-1}\)
\(\displaystyle u_{n+1} = \frac{h \left(\varphi_{n+1} + \varphi_{n-1} + 4 \varphi_{n}\right)}{3} + u_{n-1}\)
\(\displaystyle u_{n+1} = \frac{h \left(29 \varphi_{n+1} + 24 \varphi_{n-1} + 4 \varphi_{n-2} - \varphi_{n-3} + 124 \varphi_{n}\right)}{90} + u_{n-1}\)
\(\displaystyle u_{n+1} = \frac{h \left(28 \varphi_{n+1} + 14 \varphi_{n-1} + 14 \varphi_{n-2} - 6 \varphi_{n-3} + \varphi_{n-4} + 129 \varphi_{n}\right)}{90} + u_{n-1}\)
\(\displaystyle u_{n+1} = \frac{h \left(1139 \varphi_{n+1} + 33 \varphi_{n-1} + 1328 \varphi_{n-2} - 807 \varphi_{n-3} + 264 \varphi_{n-4} - 37 \varphi_{n-5} + 5640 \varphi_{n}\right)}{3780} + u_{n-1}\)
Combiner des schéma multi-pas avec une approche predictor-corrector
En général, les méthodes multi-pas implicites sont, à nombre de pas équivalent, plus précises et plus stables que les méthodes explicites (nous aborderons ces propriétés dans les prochains cours). Cependant, elles sont aussi plus coûteuses, car le calcul de \(u_{n+1}\) nécessite la résolution d’une équation. En général, cette équation est non linéaire et on ne sait pas la résoudre explicitement en toute généralité, on doit alors en chercher une solution approchée, par exemple à l’aide d’une méthode itérative comme la méthode de point fixe (nous utilisons habituellement la méthode fsolve de scipy).
Une approche différente, qui permet de s’affranchir de cette étape coûteuse, est donnée par les méthodes predictor-corrector (PC). L’idée des méthodes PC est de diminuer le coût de chaque itération en remplaçant la résolution de l’équation non-linéaire par une prédiction obtenue en utilisant une méthode explicite, puis de n’effectuer qu’une itération (ou quelques itérations) de la méthode implicite, qui est alors utilisée comme correcteur explicite.
En somme, les méthodes predictor-corrector permettent de calculer \(u_{n+1}\) de façon explicite à partir d’une méthode implicite, réduisant ainsi le besoin de résoudre des équations non-linéaires à chaque étape.
Principe de la méthode predictor-corrector (PC)
On considère une méthode multi-pas implicite (de coefficients \(\{a_j\}\) et \(\{b_j\}\))
\(
u_{n+1} = \displaystyle\sum_{j=0}^p \Big(a_ju_{n-j} + h b_j\varphi_{n-j}\Big) + hb_{-1}\varphi_{n+1} \quad\text{ pour $n=p,p+1,\dots,N-1$.}
\)
On modifie la méthode implicite comme suit :
- étape de prédiction : on calcule \({\color{blue}{\tilde u_{n+1}}}\) une approximation de \(u_{n+1}\) obtenue par une méthode explicite, typiquement une autre méthode multi-pas linéaire (de coefficients \(\{\tilde a_j\}\) et \(\{\tilde b_j\}\));
- étape de correction : on “corrige” la méthode implicite en remplaçant \(\varphi_{n+1}\) par \(\varphi(t_{n+1},{\color{blue}{\tilde u_{n+1}}})\).
Cela donne une récurrence du type
\(\begin{cases}
{\color{blue}{\tilde u_{n+1}}} = \displaystyle\sum_{j=0}^{\tilde p} \Big(\tilde a_ju_{n-j} + h \tilde b_j\varphi_{n-j}\Big) \\
u_{n+1}=\displaystyle\sum_{j=0}^p \Big(a_ju_{n-j} + h b_j\varphi_{n-j} \Big) + hb_{-1}\varphi(t_{n+1},{\color{blue}{\tilde u_{n+1}}})
\end{cases}\)
Bien-sûr, il reste à définir les valeurs d’initialisation. Comme le schéma final est un schéma explicite, on utlisera des initialisations obtenues avec des méthodes explicites.
Exemple : schéma de Heun
Considérons les deux schémas multi-pas à 1 pas suivants :
- predictor: méthode d’Euler explicite (ou \(\text{AB}_1\)) \(\tilde u_{n+1}=u_n+h\varphi_n\)
- corrector: méthode de Crank-Nicolson (ou \(\text{AM}_1\)) \(u_{n+1}=u_n+\frac{h}{2}\left(\varphi_n+\varphi(t_{n+1},u_{n+1})\right)\)
Le schéma PC associé à ces deux schémas n’est rien d’autre que la méthode de Heun
\(\text{(H)}\quad
\begin{cases}
u_0 = y_0 \\
\text{Pour }n=0,1,\dots N-1:\\
\quad \tilde u_{n+1} = u_n + h\varphi_n \\
\quad u_{n+1} = u_n + \frac{h}{2} \left(\varphi(t_{n+1},\tilde u_{n+1})+\varphi_{n}\right)
\end{cases}\)
Des méthodes de type predictor-corrector sont souvent construites en utilisant une prédiction d’Adams-Bashforth suivie d’une correction d’Adams-Moulton.
Exemple : schéma AB2-AM2
Considérons les deux schémas multi-pas à 2 pas suivants :
- predictor: méthode \(\text{AB}_2\) \(u_{n+1}=u_{n}+ \frac{h}{2} \left(3\varphi_{n}-\varphi_{n-1}\right)\)
- corrector: méthode \(\text{AM}_2\) \(u_{n+1}=u_{n}+\frac{h}{12}\left(5\varphi(t_{n+1},u_{n+1})+8\varphi_n-\varphi_{n-1}\right)\)
Le schéma PC associé à ces deux schémas s’écrit
\(\text{($\text{AB}_2$-$\text{AM}_2$)}\quad
\begin{cases}
u_0 = y_0 \\
u_1 = u_0 +h\varphi_0,\\
\text{Pour }n=1,2,\dots N-1:\\
\quad\tilde u_{n+1} = u_n +\frac{1}{2}h\left( 3\varphi_n-\varphi_{n-1} \right) \\
\quad u_{n+1} = u_{n}+\frac{h}{12}\left(5\varphi(t_{n+1},\tilde u_{n+1})+8\varphi_n-\varphi_{n-1}\right)
\end{cases}\)
Comme c’est un schéma à 2 pas, on a besoin d’initialiser \(u_1\). Puisqu’au final nous voulons un schéma explicite, on utiliser alors un schéma explicite pour initialiser cette quantité, à savoir le schéma \(\text{AB}_1\) qui n’est rien d’autre que le schéma d’Euler explicite.
Astuce de programmation : on affecte les Ă©valuations de \(\varphi_j\) pour les utiliser dans les deux Ă©tapes :
\(\text{($\text{AB}_2$-$\text{AM}_2$)}\quad
\begin{cases}
u_0 = y_0 \\
u_1 = u_0 +h\varphi_0,\\
\text{Pour }n=1,2,\dots N-1:\\
\quad k_0 = \varphi_n, \\
\quad k_1 = \varphi_{n-1}\\
\quad \tilde u_{n+1} = u_n +\frac{1}{2}h\left( 3k_0-k_1 \right) \\
\quad u_{n+1} = u_{n}+\frac{h}{12}\left(5\varphi(t_{n+1},\tilde u_{n+1})+8k_0-k_1\right)
\end{cases}\)
On peut combiner des schémas multi-pas de différents types, comme une prédiction d’Adam-Bashforth suivie d’une correction BDF.
Exemple : schéma AB2-BDF2
Considérons les deux schémas multi-pas à 2 pas suivants :
- predictor: méthode \(\text{AB}_2\) \(u_{n+1}=u_{n}+ \frac{h}{2} \left(3\varphi_{n}-\varphi_{n-1}\right)\)
- corrector: méthode \(\text{BDF}_2\) \(u_{n+1}=\frac{4}{3}u_n-\frac{1}{3}u_{n-1}+\frac{2}{3}h\varphi(t_{n+1},u_{n+1})\)
Le schéma PC associé à ces deux schémas s’écrit
\(\text{($\text{AB}_2$-$\text{BDF}_2$)}\quad
\begin{cases}
u_0 = y_0 \\
u_1 = u_0 +h\varphi_0,\\
\text{Pour }n=1,2,\dots N-1:\\
\quad \tilde u_{n+1} = u_n +\frac{h}{2}\left( 3\varphi_n-\varphi_{n-1} \right)\\
\quad u_{n+1}=\frac{4}{3}u_n-\frac{1}{3}u_{n-1}+\frac{2}{3}h\varphi(t_{n+1},\tilde u_{n+1})
\end{cases}\)
