CM3 - Construction de schémas 👣 multipas 👣 linéaires

Auteur·rice

Gloria FACCANONI

Date de publication

12 avril 2024

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

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

L’ensemble de $N+1$ valeurs $\{t_0, t_1=t_0+h,\dots , t_{N}=T \}$ représente les points de la discrétisation.
L’ensemble de $N+1$ valeurs $\{y_0, y_1,\dots , y_{N} \}$ représente la solution exacte discrète.
L’ensemble de $N+1$ valeurs $\{u_0 = y_0, u_1,\dots , u_{N} \}$ représente la solution numérique obtenue en construisant une suite récurrente.

Les schémas que l’on va construire permettent de calculer (explicitement ou implicitement) $u_{n+1}$ à partir de $u_n, u_{n-1}, ..., u_{n-p}$ et il est donc possible de calculer successivement $u_1$, $u_2$,…, en partant de $u_0$ par une formule de récurrence de la forme

$\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-p}), &\forall n=p,p+1,\dots,N-1. \end{cases}$

Plus précisement, nous considérerons des schémas à $q=p+1$ étapes (=pas) linéaires qui s’écrivent sous la forme générale

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

où les $\{a_j\}$ et $\{b_j\}$ sont des coefficients donnés et $p\ge0$ un entier.

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_j$, $j\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).

Exemples Les schémas d’Euler explicite, d’Euler implicite et de Crank-Nicolson sont des schémas à $q=1$ pas ($p=0$) linéaires

$\begin{cases} u_0=y_0,\\ u_{n+1} = a_0u_{n}+hb_0\varphi(t_{n},u_{n})+hb_{-1}\varphi(t_{n+1},u_{n+1}), &\forall n=0,1,\dots,N-1. \end{cases}$

avec
$ \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(t_n,u_n) & \text{Explicite} & 1 & 1 & 0 \\ \text{EI} & u_{n+1}=u_n+h\varphi(t_{n+1},u_{n+1}) & \text{Implicite} & 1 & 0 & 1 \\ \text{CN} & u_{n+1}=u_n+\frac{h}{2}\left(\varphi(t_n,u_n)+\varphi(t_{n+1},u_{n+1})\right) & \text{Implicite} & 1 & \frac{1}{2} & \frac{1}{2} \\ \hline \end{array} $

Remarques pour l’implémentation

range(5) produit 0,1,2,3,4
range(1,5) produit 1,2,3,4
$\{t_0, t_1,\dots , t_{N} \}$ $N+1$ points de discrétisation
$\{u_0, u_1,\dots , u_{N} \}$ $N+1$ valeurs à calculer.
len(tt) $=N+1$

Pour un schéma à un pas, $u_{n+1}=F(u_n)$ donc on initialise $u_0$ et on calcule $u_{n+1}$ pour $n$ de $0$ jusqu’à $N-1$, autrement dit n in range(N) soit encore n in range(len(tt)-1)

Pour un schéma à deux pas, $u_{n+1}=F(u_n,u_{n-1})$ donc on initialise $u_0$ et $u_1$ et on calcule $u_{n+1}$ pour $n$ de $1$ jusqu’à $N-1$, autrement dit n in range(1,N) soit encore n in range(1,len(tt)-1)

Pour un schéma à trois pas, $u_{n+1}=F(u_n,u_{n-1},u_{n-2})$ donc on initialise $u_0$, $u_1$ et $u_2$ et on calcule $u_{n+1}$ pour $n$ de $2$ jusqu’à $N-1$, autrement dit n in range(2,N) soit encore n in range(2,len(tt)-1)

…

Pour un schéma à $q$ pas, $u_{n+1}=F(u_n,u_{n-1},u_{n-2},\dots,u_{n-q+1})$ donc on initialise $u_0$, $u_1$ et $u_{q-1}$ et on calcule $u_{n+1}$ pour $n$ de $q-1$ jusqu’à $N-1$, autrement dit n in range(q-1,N) soit encore n in range(q-1,len(tt)-1)

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.

Code

from IPython.display import display, Math, Markdown, Latex


import sympy as symb
# symb.init_printing()
symb.init_printing(use_latex='mathjax', use_unicode=True)


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(r'\phi_{n+1}')
phi_n  =symb.Symbol(r'\phi_{n}')
phi_nm1=symb.Symbol(r'\phi_{n-1}')
phi_nm2=symb.Symbol(r'\phi_{n-2}')
phi_nm3=symb.Symbol(r'\phi_{n-3}')
phi_nm4=symb.Symbol(r'\phi_{n-4}')
phi_nm5=symb.Symbol(r'\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}')

1 Schémas de Adam

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. $ Différents schémas peuvent être construits en utilisant des formules de quadrature pour approcher le membre de droite. Cette solution approchée sera obtenue en construisant une suite récurrente comme suit:

$\begin{cases} u_0=y_0,\\ u_{n+1}=u_n+\displaystyle\int_{t_n}^{t_{n+1}} \tilde f(t) \mathrm{d}t \quad \text{où }\tilde f(t) \approx \varphi(t,y(t)) \end{cases} $

Les schémas d’Adam choisissent $\tilde f(t)$ comme 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}]$.
Différents schémas peuvent être construits en choisissant différents points d’interpolation.

Ces schémas peuvent généralement être exprimés sous la forme suivante :

$\begin{cases} u_0=y_0,\\ u_1\approx y_1,\\ \vdots\\ u_{p}\approx y_{p},\\ \displaystyle u_{n+1}=u_n+h\sum_{j=0}^{p}b_j\varphi(t_j,u_{n-j})+hb_{-1}\varphi(t_{n+1},u_{n+1}) \qquad n=p,p+1,\dots,N-1 \end{cases} $

Ils se divisent en deux familles:

les schémas d’Adam-Bashforth à $q$ pas (AB-$q$) sont
- explicites,
- approchent l’intégrale $\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t$ par l’intégrale $\int_{t_n}^{t_{n+1}} p(t) \mathrm{d}t$ où $p$ est le polynôme interpolant $\varphi$ en $ \{t_n,t_{n-1},\dots,t_{n-q+1}\} \text{ avec } q\ge1; $
les schémas d’Adam-Moulton à $q$ pas (AM-$q$) sont
- implicites,
- approchent l’intégrale $\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t$ par l’intégrale $\int_{t_n}^{t_{n+1}} p(t) \mathrm{d}t$ où $p$ est le polynôme interpolant $\varphi$ en $ \{t_{n+1},t_n,t_{n-1},\dots,t_{n-q+1}\} \text{ avec } q\ge0. $

Remarques

La notation AM-$k$ peut varier selon les auteurs. Par exemple, dans les livres de Quarteroni et al., le symbole $k$ indique l’ordre de la méthode et non le nombre de pas. Les deux notations coincident pour les méthodes AB mais diffèrent pour les méthodes AM.
Noter que, pour calculer successivement $u_{q}$, $u_{q+1}$, …, il est nécessaire d’initialiser d’abord $u_0, u_1, \dots, u_{q-1}$. 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 de l’ordre. En particulier, une méthode explicitement A-stable est d’ordre au plus 2 (seconde barrière de Dahlquist).

1.1 AB-1

On a

Points d’interpolation: $\{(t_n,\varphi(t_n,y_n))\}$
Polynôme: $p(t)=\varphi(t_n,y_n)$
Approximation: $y_{n+1}-y_n \approx \int_{t_n}^{t_{n+1}} p(t) \mathrm{d}t =h\varphi(t_n,y_n)$

et on obtient le schéma

$ \text{($\text{AB}_1$)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ u_{n+1}=u_{n}+h \varphi(t_{n},u_{n})& n=0,1,\dots N-1 \end{cases}$

La méthode $\text{AB}_1$ coïncide avec la méthode d’Euler progressive.

Code

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))

$\displaystyle u_{n+1} = \phi_{n} h + u_{n}$

1.2 AB-2

On a

Points d’interpolation: $\{(t_n,\varphi(t_n,y_n));(t_{n-1},\varphi(t_{n-1},y_{n-1}))\}$
Polynôme: $p(t)=\varphi(t_n,y_n)\frac{t-t_{n-1}}{t_n-t_{n-1}}+\varphi(t_{n-1},y_{n-1})\frac{t-t_n}{t_{n-1}-t_n} =\varphi(t_n,y_n)\frac{t-t_{n-1}}{h}+\varphi(t_{n-1},y_{n-1})\frac{t-t_n}{-h}$
Approximation: $y_{n+1}-y_n \approx \int_{t_n}^{t_{n+1}} p(t) \mathrm{d}t =\frac{h}{2}\left(3\varphi(t_n,y(t_n))-\varphi(t_{n-1},y_{n-1})\right)$

et 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(t_{0},u_{0})\approx y(t_1)\\ u_{n+1}=u_{n}+ \frac{h}{2} \big(3\varphi(t_{n},u_{n})-\varphi(t_{n-1},u_{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{AB}_1$.

Code

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))

$\displaystyle u_{n+1} = \frac{h \left(- \phi_{n-1} + 3 \phi_{n}\right)}{2} + u_{n}$

1.3 AB-3

On a

Points d’interpolation: $\{(t_n,\varphi(t_n,y_n));(t_{n-1},\varphi(t_{n-1},y_{n-1}));(t_{n-2},\varphi(t_{n-2},y_{n-2}))\}$
Polynôme:
$\begin{aligned} p(t)&=\varphi(t_n,y_n)\frac{(t-t_{n-1})(t-t_{n-2})}{(t_n-t_{n-1})(t_n-t_{n-2})} +\varphi(t_{n-1},y(t_{n-1}))\frac{(t-t_{n})(t-t_{n-2})}{(t_{n-1}-t_{n})(t_{n-1}n-t_{n-2})} +\varphi(t_{n-2},y(t_{n-2}))\frac{(t-t_{n})(t-t_{n-1})}{(t_{n-2}-t_{n})(t_{n-2}n-t_{n-1})} \\ &=\frac{\varphi(t_{n-2},y_{n-2})}{2h^2}(t-t_{n-1})(t-t_{n}) -\frac{\varphi(t_{n-1},y_{n-1})}{h^2}(t-t_{n-2})(t-t_{n}) +\frac{\varphi(t_{n},y_{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}} p(t) \mathrm{d}t =\frac{h}{12}\left(23\varphi(t_n,y(t_n))-16\varphi(t_{n-1},y(t_{n-1}))+5\varphi(t_{n-2},y(t_{n-2}))\right)$

et 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(t_{0},u_{0})\approx y(t_1)\\ u_2 = u_1 + \dfrac{h}{2} \big(3\varphi(t_1,u_1)-\varphi(t_0,u_0)\big)\approx y(t_2)\\ u_{n+1} = u_{n} + \dfrac{h}{12}\big(23\varphi(t_n,u_n)-16\varphi(t_{n-1},u_{n-1})+5\varphi(t_{n-2},u_{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$.

Code

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))

$\displaystyle u_{n+1} = \frac{h \left(- 16 \phi_{n-1} + 5 \phi_{n-2} + 23 \phi_{n}\right)}{12} + u_{n}$

1.4 AB-4

On a

Points d’interpolation: $\{(t_n,\varphi(t_n,y_n));(t_{n-1},\varphi(t_{n-1},y_{n-1}));(t_{n-2},\varphi(t_{n-2},y_{n-2}));(t_{n-3},\varphi(t_{n-3},y_{n-3}))\}$
Polynôme: $p(t)=\displaystyle\sum_{i=n-3}^{n}\left(\varphi(t_i,y_i)\prod_{\substack{j=n-2,...,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}} p(t) \mathrm{d}t =\frac{h}{24}\left(55\varphi(t_n,y(t_n))-59\varphi(t_{n-1},y(t_{n-1}))+37\varphi(t_{n-2},y(t_{n-2})-9\varphi(t_{n-3},y(t_{n-3}))\right)$

et 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(t_{0},u_{0})\approx y(t_1)\\ u_2 = u_1 + \dfrac{h}{2} \big(3\varphi(t_1,u_1)-\varphi(t_0,u_0)\big)\approx y(t_2)\\ u_{3} = u_{2} + \dfrac{h}{12}\big(23\varphi(t_2,u_2)-16\varphi(t_{1},u_{1})+5\varphi(t_{0},u_{0})\big)\approx y(t_3) \\ u_{n+1} = u_{n} + \dfrac{h}{24}\big(55\varphi(t_n,u_n)-59\varphi(t_{n-1},u_{n-1})+37\varphi(t_{n-2},u_{n-2})-9\varphi(t_{n-3},u_{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.

Code

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))

$\displaystyle u_{n+1} = \frac{h \left(- 59 \phi_{n-1} + 37 \phi_{n-2} - 9 \phi_{n-3} + 55 \phi_{n}\right)}{24} + u_{n}$

Pour la construction des autres schémas $\text{AB}_5$ et $\text{AB}_6$ on verra plus bas dans le récapitulatif.

1.5 AM-0

On a

Points d’interpolation: $\{(t_{n+1},\varphi(t_{n+1},y_{n+1}))\}$
Polynôme: $p(t)=\varphi(t_{n+1},y(t_{n+1}))$
Approximation: $y_{n+1}-y_n \approx \int_{t_n}^{t_{n+1}} p(t) \mathrm{d}t =h\varphi(t_{n+1},y(t_{n+1}))$

et 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(t_{n+1},u_{n+1})& n=0,1,\dots N-1 \end{cases}$

La méthode $\text{AM}_0$ coïncide avec la méthode d’Euler regressive.

Code

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))

$\displaystyle u_{n+1} = \phi_{n+1} h + u_{n}$

1.6 AM-1

On a

Points d’interpolation: $\{(t_{n+1},\varphi(t_{n+1},y_{n+1})),(t_{n},\varphi(t_{n},y_{n}))\}$
Polynôme: $p(t)=\varphi(t_{n+1},y_{n+1})\frac{t-t_n}{t_{n+1}-t_n}+\varphi(t_{n},y_{n})\frac{t-t_{n+1}}{t_n-t_{n+1}}$
Approximation: $y_{n+1}-y_n \approx \int_{t_n}^{t_{n+1}} p(t) \mathrm{d}t =\frac{h}{2}\left(\varphi(t_{n+1},y_{n+1})+\varphi(t_n,y_n)\right)$

et 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(t_{n},u_{n})+\varphi(t_{n+1},u_{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.

Code

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))

$\displaystyle u_{n+1} = \frac{h \left(\phi_{n+1} + \phi_{n}\right)}{2} + u_{n}$

1.7 AM-2

On a

Points d’interpolation: $\{(t_{n+1},\varphi(t_{n+1},y_{n+1})),(t_{n},\varphi(t_{n},y_{n})),(t_{n-1},\varphi(t_{n-1},y_{n-1}))\}$
Polynôme: $p(t)=\varphi(t_{n+1},y_{n+1})\frac{(t-t_n)(t-t_{n-1})}{(t_{n+1}-t_n)(t_{n+1}-t_{n-1})} +\varphi(t_{n},y_{n})\frac{(t-t_{n+1})(t-t_{n-1})}{(t_n-t_{n+1})(t_n-t_{n-1})}+\varphi(t_{n-1},y_{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}} p(t) \mathrm{d}t =\frac{h}{12}\left(5\varphi(t_{n+1},y_{n+1})+8\varphi(t_n,y_n)-\varphi(t_{n-1},y_{n-1}))\right)$

et 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(t_{1},u_{1})+\varphi(t_{0},u_{0})\big)\\ u_{n+1}=u_{n}+\dfrac{h}{12}\big(5\varphi(t_{n+1},u_{n+1})+8\varphi(t_n,u_n)-\varphi(t_{n-1},u_{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$.

Code

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))

$\displaystyle u_{n+1} = \frac{h \left(5 \phi_{n+1} - \phi_{n-1} + 8 \phi_{n}\right)}{12} + u_{n}$

Pour la construction des autres schémas $\text{AM}_3$, $\text{AM}_4$, $\text{AM}_5$ et $\text{AM}_6$ on verra plus bas dans le récapitulatif.

1.8 Calcul systématique des coefficients des méthodes AB et AM

Code

display(Markdown(f"#### Adam-Bashforth"))
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())
    display(Markdown(f"**AB-{q}: explicite, à {q} pas**"))
    display(symb.Eq(u_np1, u_n + AB))
    #display(Markdown("---"))


display(Markdown(f"#### Adam-Moulton"))
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())
    display(Markdown(f"**AM-{q}: implicite, à {q} pas**"))
    display(symb.Eq(u_np1, u_n + AM))
    #(Markdown("---"))

1.8.1 Adam-Bashforth

AB-1: explicite, à 1 pas

$\displaystyle u_{n+1} = \phi_{n} h + u_{n}$

AB-2: explicite, à 2 pas

$\displaystyle u_{n+1} = \frac{h \left(- \phi_{n-1} + 3 \phi_{n}\right)}{2} + u_{n}$

AB-3: explicite, à 3 pas

$\displaystyle u_{n+1} = \frac{h \left(- 16 \phi_{n-1} + 5 \phi_{n-2} + 23 \phi_{n}\right)}{12} + u_{n}$

AB-4: explicite, à 4 pas

$\displaystyle u_{n+1} = \frac{h \left(- 59 \phi_{n-1} + 37 \phi_{n-2} - 9 \phi_{n-3} + 55 \phi_{n}\right)}{24} + u_{n}$

AB-5: explicite, à 5 pas

$\displaystyle u_{n+1} = \frac{h \left(- 2774 \phi_{n-1} + 2616 \phi_{n-2} - 1274 \phi_{n-3} + 251 \phi_{n-4} + 1901 \phi_{n}\right)}{720} + u_{n}$

AB-6: explicite, à 6 pas

$\displaystyle u_{n+1} = \frac{h \left(- 7923 \phi_{n-1} + 9982 \phi_{n-2} - 7298 \phi_{n-3} + 2877 \phi_{n-4} - 475 \phi_{n-5} + 4277 \phi_{n}\right)}{1440} + u_{n}$

1.8.2 Adam-Moulton

AM-0: implicite, à 0 pas

$\displaystyle u_{n+1} = \phi_{n+1} h + u_{n}$

AM-1: implicite, à 1 pas

$\displaystyle u_{n+1} = \frac{h \left(\phi_{n+1} + \phi_{n}\right)}{2} + u_{n}$

AM-2: implicite, à 2 pas

$\displaystyle u_{n+1} = \frac{h \left(5 \phi_{n+1} - \phi_{n-1} + 8 \phi_{n}\right)}{12} + u_{n}$

AM-3: implicite, à 3 pas

$\displaystyle u_{n+1} = \frac{h \left(9 \phi_{n+1} - 5 \phi_{n-1} + \phi_{n-2} + 19 \phi_{n}\right)}{24} + u_{n}$

AM-4: implicite, à 4 pas

$\displaystyle u_{n+1} = \frac{h \left(251 \phi_{n+1} - 264 \phi_{n-1} + 106 \phi_{n-2} - 19 \phi_{n-3} + 646 \phi_{n}\right)}{720} + u_{n}$

AM-5: implicite, à 5 pas

$\displaystyle u_{n+1} = \frac{h \left(475 \phi_{n+1} - 798 \phi_{n-1} + 482 \phi_{n-2} - 173 \phi_{n-3} + 27 \phi_{n-4} + 1427 \phi_{n}\right)}{1440} + u_{n}$

AM-6: implicite, à 6 pas

$\displaystyle u_{n+1} = \frac{h \left(19087 \phi_{n+1} - 46461 \phi_{n-1} + 37504 \phi_{n-2} - 20211 \phi_{n-3} + 6312 \phi_{n-4} - 863 \phi_{n-5} + 65112 \phi_{n}\right)}{60480} + u_{n}$

2 Schémas de Nyström et de Milne-Simpson

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:

$\begin{cases} u_0=y_0,\\ u_1=?,\\ u_{n+1}=u_{n-1}+\displaystyle\int_{t_{n-1}}^{t_{n+1}} \tilde f(t) \mathrm{d}t \quad \text{où }\tilde f(t)\approx\varphi(t,y(t)) \end{cases} $

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

$\begin{cases} u_0=y_0,\\ u_1\approx y_1,\\ \vdots\\ u_{p}\approx y_{p},\\ \displaystyle u_{n+1}=u_{n-1}+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} $

Ils se divisent en deux familles:

les schémas de Nyström à $q$ pas (N-$q$) sont
- explicites,
- approchent l’intégrale $\int_{t_{n-1}}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t$ par l’intégrale $\int_{t_{n-1}}^{t_{n+1}} p(t) \mathrm{d}t$ où $p$ est le polynôme interpolant $\varphi$ en $ \{t_n,t_{n-1},\dots,t_{n-p}\} \text{ avec } p\ge0. $
les schémas de Milne-Simpson à $q$ pas (MS-$q$) sont
- implicites,
- approchent l’intégrale $\int_{t_{n-1}}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t$ par l’intégrale $\int_{t_{n-1}}^{t_{n+1}} p(t) \mathrm{d}t$ où $p$ est le polynôme interpolant $\varphi$ en $ \{t_{n+1},t_n,t_{n-1},\dots,t_{n-p}\} \text{ avec } p\ge-1. $

2.1 N-1

On a

Points d’interpolation: $\{(t_n,\varphi(t_n,y_n))\}$
Polynôme: $p(t)=\varphi(t_n,y(t_n))$
Approximation: $y_{n+1}-y_{n-1} \approx \int_{t_{n-1}}^{t_{n+1}} p(t) \mathrm{d}t =2h\varphi(t_n,y(t_n))$

et on obtient le schéma

$ \text{(N$_1$)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ u_1=u_0+h\varphi(t_{0},u_{0})\approx y(t_1)\\ u_{n+1}=u_{n-1}+2h \varphi(t_{n},u_{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.

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

Code

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))

$\displaystyle u_{n+1} = 2 \phi_{n} h + u_{n-1}$

2.2 N-2

On a

Points d’interpolation: $\{(t_n,\varphi(t_n,y_n)),(t_{n-1},\varphi(t_{n-1},y_{n-1}))\}$
Polynôme: $p(t)=\varphi(t_n,y(t_n))\frac{t-t_{n-1}}{t_n-t_{n-1}}+\varphi(t_{n-1},y_{n-1})\frac{t-t_n}{t_{n-1}-t_n}$
Approximation: $y_{n+1}-y_{n-1} \approx \int_{t_{n-1}}^{t_{n+1}} p(t) \mathrm{d}t =2h\varphi(t_n,y(t_n))$

et on obtient le schéma

$ \text{(N$_2$)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ u_1=u_0+h\varphi(t_{0},u_{0})\approx y(t_1)\\ u_{n+1}=u_{n-1}+2h \varphi(t_{n},u_{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.

On retrouve la méthode N$_1$.

Code

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))

$\displaystyle u_{n+1} = 2 \phi_{n} h + u_{n-1}$

2.3 N-3

On a

Points d’interpolation: $\{(t_n,\varphi(t_n,y_n));(t_{n-1},\varphi(t_{n-1},y_{n-1}));(t_{n-2},\varphi(t_{n-2},y_{n-2}))\}$
Polynôme: $\begin{aligned}[t] p(t)&=\varphi(t_n,y_n)\frac{(t-t_{n-1})(t-t_{n-2})}{(t_n-t_{n-1})(t_n-t_{n-2})} +\varphi(t_{n-1},y_{n-1})\frac{(t-t_{n})(t-t_{n-2})}{(t_{n-1}-t_{n})(t_{n-1}n-t_{n-2})} +\varphi(t_{n-2},y_{n-2})\frac{(t-t_{n})(t-t_{n-1})}{(t_{n-2}-t_{n})(t_{n-2}n-t_{n-1})} \\ &=\frac{\varphi(t_{n-2},y_{n-2})}{2h^2}(t-t_{n-1})(t-t_{n}) -\frac{\varphi(t_{n-1},y_{n-1})}{h^2}(t-t_{n-2})(t-t_{n}) +\frac{\varphi(t_{n},y_{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}} p(t) \mathrm{d}t =\frac{h}{3}\left(7\varphi(t_n,y_n)-2\varphi(t_{n-1},y_{n-1})+\varphi(t_{n-2},y_{n-2})\right)$

et on obtient le schéma

$ \text{(N$_3$)}\quad \begin{cases} u_0 = y(t_0) = y_0,\\ u_1 = u_0 + h\varphi(t_{0},u_{0})\approx y(t_1)\\ u_2 = u_{0} + 2h\varphi(t_{1},u_{1})\approx y(t_2)\\ u_{n+1} = u_{n-1} + \dfrac{h}{3}\big(7\varphi(t_n,u_n)-2\varphi(t_{n-1},u_{n-1})+\varphi(t_{n-2},u_{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 d’Euler explicite et $u_{2}$ est une approximation de $y(t_2)$ obtenue en utilisant la méthode N$_2$.

Code

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))

$\displaystyle u_{n+1} = \frac{h \left(- 2 \phi_{n-1} + \phi_{n-2} + 7 \phi_{n}\right)}{3} + u_{n-1}$

Pour la construction des autres schémas N$_4$, N$_5$ et N$_6$ on verra plus bas dans le récapitulatif.

2.4 MS-0

On a

Points d’interpolation: $\{(t_{n+1},\varphi(t_{n+1},y_{n+1}))\}$
Polynôme: $p(t)=\varphi(t_{n+1},y(t_{n+1}))$
Approximation: $y_{n+1}-y_{n-1} \approx \int_{t_{n-1}}^{t_{n+1}} p(t) \mathrm{d}t =2h\varphi(t_{n+1},y(t_{n+1}))$

et on obtient le schéma

$ \text{(MS$_0$)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ u_1=u_0+h\varphi(t_{1},u_{1})\approx y(t_1)\\ u_{n+1}=u_{n-1}+2h \varphi(t_{n+1},u_{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 il est à 2 pas car il dépend de $u_{n-1}$.

Code

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))

$\displaystyle u_{n+1} = 2 \phi_{n+1} h + u_{n-1}$

2.5 MS-1

On a

Points d’interpolation: $\{(t_{n+1},\varphi(t_{n+1},y_{n+1})),(t_{n},\varphi(t_{n},y_{n}))\}$
Polynôme: $p(t)=\varphi(t_{n+1},y_{n+1})\frac{t-t_n}{t_{n+1}-t_n} +\varphi(t_{n},y_{n})\frac{t-t_{n+1}}{t_n-t_{n+1}}$
Approximation: $y_{n+1}-y_{n-1} \approx \int_{t_{n-1}}^{t_{n+1}} p(t) \mathrm{d}t =2h\varphi(t_n,y_n)$

et on obtient le schéma

$ \text{(MS$_1$)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ u_1=u_0+h\varphi(t_{1},u_{1})\approx y(t_1)\\ u_{n+1}=u_{n-1}+2h \varphi(t_{n},u_{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.

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

Code

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))

$\displaystyle u_{n+1} = 2 \phi_{n} h + u_{n-1}$

2.6 MS-2

On a

Points d’interpolation: $\{(t_{n+1},\varphi(t_{n+1},y_{n+1})),(t_{n},\varphi(t_{n},y_{n})),(t_{n-1},\varphi(t_{n-1},y_{n-1}))\}$
Polynôme: $p(t)=\varphi(t_{n+1},y_{n+1})\frac{(t-t_n)(t-t_{n-1})}{(t_{n+1}-t_n)(t_{n+1}-t_{n-1})}+ \varphi(t_{n},y_{n})\frac{(t-t_{n+1})(t-t_{n-1})}{(t_n-t_{n+1})(t_n-t_{n-1})}+ \varphi(t_{n-1},y_{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-1} \approx \int_{t_{n-1}}^{t_{n+1}} p(t) \mathrm{d}t =\frac{h}{3}\left(\varphi(t_{n+1},y_{n+1})+4\varphi(t_n,y_n)+\varphi(t_{n-1},y_{n-1}))\right)$

et 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(t_{1},u_{1})+\varphi(t_{0},u_{0})\big)\\ u_{n+1}=u_{n-1}+\dfrac{h}{3}\big( \varphi(t_{n+1},u_{n+1})+4\varphi(t_n,u_n)+\varphi(t_{n-1},u_{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).

Code

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))

$\displaystyle u_{n+1} = \frac{h \left(\phi_{n+1} + \phi_{n-1} + 4 \phi_{n}\right)}{3} + u_{n-1}$

Pour la construction des autres schémas MS$_3$, MS$_4$, MS$_5$ et MS$_6$ on verra plus bas dans le récapitulatif.

2.7 Calcul systématique des coefficients des méthodes N et MS

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("---"))

2.7.1 Nyström

N-1: explicite, à 1 pas

$\displaystyle u_{n+1} = 2 \phi_{n} h + u_{n-1}$

N-2: explicite, à 2 pas

$\displaystyle u_{n+1} = 2 \phi_{n} h + u_{n-1}$

N-3: explicite, à 3 pas

$\displaystyle u_{n+1} = \frac{h \left(- 2 \phi_{n-1} + \phi_{n-2} + 7 \phi_{n}\right)}{3} + u_{n-1}$

N-4: explicite, à 4 pas

$\displaystyle u_{n+1} = \frac{h \left(- 5 \phi_{n-1} + 4 \phi_{n-2} - \phi_{n-3} + 8 \phi_{n}\right)}{3} + u_{n-1}$

N-5: explicite, à 5 pas

$\displaystyle u_{n+1} = \frac{h \left(- 266 \phi_{n-1} + 294 \phi_{n-2} - 146 \phi_{n-3} + 29 \phi_{n-4} + 269 \phi_{n}\right)}{90} + u_{n-1}$

N-6: explicite, à 6 pas

$\displaystyle u_{n+1} = \frac{h \left(- 406 \phi_{n-1} + 574 \phi_{n-2} - 426 \phi_{n-3} + 169 \phi_{n-4} - 28 \phi_{n-5} + 297 \phi_{n}\right)}{90} + u_{n-1}$

2.7.2 Milne-Simpson

MS-0: implicite, à 0 pas

$\displaystyle u_{n+1} = 2 \phi_{n+1} h + u_{n-1}$

MS-1: implicite, à 1 pas

$\displaystyle u_{n+1} = 2 \phi_{n} h + u_{n-1}$

MS-2: implicite, à 2 pas

$\displaystyle u_{n+1} = \frac{h \left(\phi_{n+1} + \phi_{n-1} + 4 \phi_{n}\right)}{3} + u_{n-1}$

MS-3: implicite, à 3 pas

$\displaystyle u_{n+1} = \frac{h \left(\phi_{n+1} + \phi_{n-1} + 4 \phi_{n}\right)}{3} + u_{n-1}$

MS-4: implicite, à 4 pas

$\displaystyle u_{n+1} = \frac{h \left(29 \phi_{n+1} + 24 \phi_{n-1} + 4 \phi_{n-2} - \phi_{n-3} + 124 \phi_{n}\right)}{90} + u_{n-1}$

MS-5: implicite, à 5 pas

$\displaystyle u_{n+1} = \frac{h \left(28 \phi_{n+1} + 14 \phi_{n-1} + 14 \phi_{n-2} - 6 \phi_{n-3} + \phi_{n-4} + 129 \phi_{n}\right)}{90} + u_{n-1}$

MS-6: implicite, à 6 pas

$\displaystyle u_{n+1} = \frac{h \left(1139 \phi_{n+1} + 33 \phi_{n-1} + 1328 \phi_{n-2} - 807 \phi_{n-3} + 264 \phi_{n-4} - 37 \phi_{n-5} + 5640 \phi_{n}\right)}{3780} + u_{n-1}$

3 Schémas BDF (Backward Differentiation Formulae)

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:

$ \begin{cases} u_0=y_0,\\ u_{n+1} \text{ solution de l'équation }\varphi(t_{n+1},u_{n+1}) = p'(t_{n+1}) \quad \text{où $p(t)$ est un polynôme interpolant }y(t). \end{cases} $

Toutes ces méthodes peuvent s’écrire sous la forme

$\begin{cases} u_0=y_0,\\ u_1\approx y_1,\\ \vdots\\ u_{q-1}\approx y_{q-1},\\ \displaystyle u_{n+1}=\sum_{j=n-q+1}^{n}a_j u_j+hb_{-1}\varphi(t_{n+1},u_{n+1}). \end{cases} $

On verra que les méthodes BDF-$q$ à $q$ pas sont d’ordre $q$ et sont zéro-stables ssi $q\le6$.

3.1 BDF-1

On a

Points d’interpolation: $\{(t_{n+1},y_{n+1}),(t_{n},y_{n})\}$
Polynôme: $p(t)=\frac{y_{n+1}-y_{n}}{h}(t-t_{n})+y_{n}$
Dérivée: $p'(t)=\frac{y_{n+1}-y_{n}}{h}$
Approximation: $\varphi(t_{n+1},y_{n+1}) \approx p'(t_{n+1}) = \frac{y_{n+1}-y_{n}}{h}$

et on obtient le schéma $\begin{cases} u_0=y(t_0)=y_0,\\ u_{n+1} \text{ solution de l'équation }\varphi(t_{n+1},u_{n+1}) = \frac{u_{n+1}-u_{n}}{h} & n=1,2,\dots N-1 \end{cases}$ c’est-à-dire le schéma

$ \text{($\text{BDF}_1$)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ u_{n+1}=u_n+h\varphi(t_{n+1},u_{n+1}) & n=0,1,\dots N-1. \end{cases}$

La méthode $\text{BDF}_1$ coïncide avec la méthode d’Euler implicite.

Code

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))

$\displaystyle p{\left(t \right)} = \frac{h u_{n} + t u_{n+1} - t u_{n} - t_{n} u_{n+1} + t_{n} u_{n}}{h}$

$\displaystyle \left. \frac{d}{d t} p{\left(t \right)} \right|_{\substack{ t=h + t_{n} }} = \frac{u_{n+1} - u_{n}}{h}$

$\displaystyle u_{n+1} = \phi_{n+1} h + u_{n}$

3.2 BDF-2

On a

Points d’interpolation: $\{(t_{n+1},y_{n+1}),(t_{n},y_{n}),(t_{n-1},y_{n-1})\}$
Polynôme:
$\begin{aligned} p(t)&=y_{n+1}\frac{(t-t_n)(t-t_{n-1})}{(t_{n+1}-t_n)(t_{n+1}-t_{n-1})}+y_{n}\frac{(t-t_{n+1})(t-t_{n-1})}{(t_{n}-t_{n+1})(t_{n}-t_{n-1})}+y_{n-1}\frac{(t-t_{n+1})(t-t_{n})}{(t_{n-1}-t_{n+1})(t_{n-1}-t_{n})} \\ &=\frac{(t-t_{n})(t-t_{n-1})}{2h^2}y_{n+1}+\frac{(t-t_{n+1})(t-t_{n-1})}{-h^2}y_{n}+\frac{(t-t_{n+1})(t-t_{n})}{-2h^2}y_{n-1} \end{aligned}$
Dérivée: $p'(t)=\frac{(t-t_{n})+(t-t_{n-1})}{2h^2}y_{n+1}+\frac{(t-t_{n+1})+(t-t_{n-1})}{-h^2}y_{n}+\frac{(t-t_{n+1})+(t-t_{n})}{-2h^2}y_{n-1}$
Approximation: $\varphi(t_{n+1},y_{n+1}) \approx p'(t_{n+1}) = \frac{3}{2h}y_{n+1}-\frac{2}{h}y_{n}+\frac{1}{2h}y_{n-1}$

et on obtient le schéma $\begin{cases} u_0=y(t_0)=y_0,\\ u_{n+1} \text{ solution de l'équation }\varphi(t_{n+1},u_{n+1}) = \frac{3}{2h}u_{n+1}-\frac{2}{h}u_{n}+\frac{1}{2h}u_{n-1} & n=2,3,\dots N-1 \end{cases}$ c’est-à-dire le schéma

$ \text{($\text{BDF}_2$)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ u_1=u_0+h\varphi(t_{1},u_{1}) \\ 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}) & n=1,2,3,\dots N-1. \end{cases}$

Code

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))

$\displaystyle p{\left(t \right)} = \frac{h^{2} u_{n} + \frac{h \left(t u_{n+1} - t u_{n-1} - t_{n} u_{n+1} + t_{n} u_{n-1}\right)}{2} + \frac{t^{2} u_{n+1}}{2} + \frac{t^{2} u_{n-1}}{2} - t^{2} u_{n} - t t_{n} u_{n+1} - t t_{n} u_{n-1} + 2 t t_{n} u_{n} + \frac{t_{n}^{2} u_{n+1}}{2} + \frac{t_{n}^{2} u_{n-1}}{2} - t_{n}^{2} u_{n}}{h^{2}}$

$\displaystyle \left. \frac{d}{d t} p{\left(t \right)} \right|_{\substack{ t=h + t_{n} }} = \frac{\frac{h \left(u_{n+1} - u_{n-1}\right)}{2} - t_{n} u_{n+1} - t_{n} u_{n-1} + 2 t_{n} u_{n} + u_{n+1} \left(h + t_{n}\right) + u_{n-1} \left(h + t_{n}\right) - 2 u_{n} \left(h + t_{n}\right)}{h^{2}}$

$\displaystyle u_{n+1} = \frac{2 \phi_{n+1} h}{3} - \frac{u_{n-1}}{3} + \frac{4 u_{n}}{3}$

3.3 Calcul systématique des coefficients des méthodes BDF

Code

display(Markdown(f"#### BDF"))
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(Markdown(f"**BDF-{q}: implicite, à {q} pas**"))
    display(symb.Eq(u_np1, sol[0]))
    # display(Markdown("---"))

3.3.1 BDF

BDF-1: implicite, à 1 pas

$\displaystyle u_{n+1} = \phi_{n+1} h + u_{n}$

BDF-2: implicite, à 2 pas

$\displaystyle u_{n+1} = \frac{2 \phi_{n+1} h}{3} - \frac{u_{n-1}}{3} + \frac{4 u_{n}}{3}$

BDF-3: implicite, à 3 pas

$\displaystyle u_{n+1} = \frac{6 \phi_{n+1} h}{11} - \frac{9 u_{n-1}}{11} + \frac{2 u_{n-2}}{11} + \frac{18 u_{n}}{11}$

BDF-4: implicite, à 4 pas

$\displaystyle u_{n+1} = \frac{12 \phi_{n+1} h}{25} - \frac{36 u_{n-1}}{25} + \frac{16 u_{n-2}}{25} - \frac{3 u_{n-3}}{25} + \frac{48 u_{n}}{25}$

BDF-5: implicite, à 5 pas

$\displaystyle u_{n+1} = \frac{60 \phi_{n+1} h}{137} - \frac{300 u_{n-1}}{137} + \frac{200 u_{n-2}}{137} - \frac{75 u_{n-3}}{137} + \frac{12 u_{n-4}}{137} + \frac{300 u_{n}}{137}$

BDF-6: implicite, à 6 pas

$\displaystyle u_{n+1} = \frac{20 \phi_{n+1} h}{49} - \frac{150 u_{n-1}}{49} + \frac{400 u_{n-2}}{147} - \frac{75 u_{n-3}}{49} + \frac{24 u_{n-4}}{49} - \frac{10 u_{n-5}}{147} + \frac{120 u_{n}}{49}$

4 Schémas (multi-pas) de type predictor-corrector

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.

On considère une méthode implicite
$ u_{n+1}=\displaystyle\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},{\color{blue}u_{n+1}}) \text{ pour $n=p,p+1,\dots,N-1$.} $

On modifie la méthode implicite comme suit:

étape de prédiction : on calcule $\tilde u_{n+1}$ une approximation de $u_{n+1}$ par une méthode explicite, typiquement une autre méthode multipas linéaire,
étape de correction : on “corrige” la méthode implicite en remplacant $ (t_{n+1},{u_{n+1}}) $ par $ (t_{n+1} , {u_{n+1}}) $.

Cela donne une récurrence du type

$\begin{cases} {\color{blue}\tilde u_{n+1}}=\displaystyle\sum_{j=0}^{\tilde p} \tilde a_ju_{n-j} + h\sum_{j=0}^{\tilde p} \tilde b_j\varphi(t_{n-j},u_{n-j}) \\ u_{n+1}=\displaystyle\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},{\color{blue}\tilde u_{n+1}}) \end{cases}$

4.1 Remarques

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

$\omega_{[C]}$ l’ordre du corrector
$\omega_{[P]}$ l’ordre du predictor

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.

4.2 Exemple : schéma de Heun

Considérons les deux schémas PC suivants:

predictor: méthode d’Euler explicite (ou $\text{AB}_1$) $\tilde u_{n+1}=u_n+h\varphi(t_n,u_n)$
corrector: méthode de Crank-Nicolson (ou $\text{AM}_1$) $u_{n+1}=u_n+\frac{h}{2}\left(\varphi(t_n,u_n)+\varphi(t_{n+1},u_{n+1})\right)$

On retrouve la méthode de Heun

$\text{(H)}\quad \begin{cases} u_0 = y_0 \\ \tilde u_{n+1} = u_n + h\varphi(t_n,u_n) \\ u_{n+1} = u_n + \frac{h}{2} \left(\varphi(t_{n},u_{n})+\varphi(t_{n+1},\tilde u_{n+1})\right) & n=1,2,\dots N-1 \end{cases}$

4.3 Exemple : AB2-AM2

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

predictor: méthode $\text{AB}_2$ $u_{n+1}=u_{n}+ \frac{h}{2} \left(3\varphi(t_{n},u_{n})-\varphi(t_{n-1},u_{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(t_n,u_n)-\varphi(t_{n-1},u_{n-1})\right)$

On obtient la méthode $\text{AB}_2-\text{AM}_2$:

$\text{($\text{AB}_2$-$\text{AM}_2$)}\quad \begin{cases} u_0 = y_0 \\ u_1 = u_0 +h\varphi(t_0,u_0),\\ \tilde u_{n+1} = u_n +\frac{3}{2}h\left( \varphi(t_n,u_n)-\varphi(t_{n-1},u_{n-1}) \right) \\ u_{n+1} = u_{n}+\frac{h}{12}\left(5\varphi(t_{n+1},\tilde u_{n+1})+8\varphi(t_n,u_n)-\varphi(t_{n-1},u_{n-1})\right) & n=2,3,\dots N-1 \end{cases}$

4.4 Exemple : AB2-BDF2

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

predictor: méthode $\text{AB}_2$ $u_{n+1}=u_{n}+ \frac{h}{2} \left(3\varphi(t_{n},u_{n})-\varphi(t_{n-1},u_{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})$

On obtient la méthode $\text{AB}_2$-$\text{BDF}_2$:

$\text{($\text{AB}_2$-$\text{BDF}_2$)}\quad \begin{cases} u_0 = y_0 \\ u_1 = u_0 +h\varphi(t_0,u_0),\\ \tilde u_{n+1} = u_n +\frac{3}{2}h\left( \varphi(t_n,u_n)-\varphi(t_{n-1},u_{n-1}) \right)\\ 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}) & n=1,2,3,\dots N-1. \end{cases}$

Retour au sommet