In [1]:
from IPython.display import display, Latex
from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read()) 
Out[1]:
In [2]:
import sys #only needed to determine Python version number
print('Python version ' + sys.version)
Python version 3.8.5 (default, Jan 27 2021, 15:41:15) 
[GCC 9.3.0]

M62_CM4 : Construction de schémas multipas

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 qu'on va construire permettent de calculer (explicitement ou implicitement) $u_{n+1}$ à partir de $u_n, u_{n-1}, ..., u_{n-k}$ 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} u_0=y_0,\\ \vdots\\ u_\kappa\approx y_\kappa,\\ u_{n+1}=\Phi(u_{n+1},u_n, u_{n-1}, \dots, u_{n-k}),&\forall n=\kappa,\kappa+1,\dots,N-1. \end{cases}$$ Plus précisement, **nous considérerons des schémas à $k=p+1$ étapes (=pas) linéaires** qui s'écrivent sous la forme générale $$ 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_k\}$ et $\{b_k\}$ 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_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 2 $\frac{1}{2}$ $\frac{1}{2}$

Remarques pour l'implémentation

  1. range(5) produit 0,1,2,3,4
    range(1,5) produit 1,2,3,4

  2. $\{t_0, t_1,\dots , t_{N} \}$ $N+1$ points de discrétisation
    $\{u_0, u_1,\dots , u_{N} \}$ $N+1$ valeurs à calculer.

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

In [3]:
import sympy as symb
symb.init_printing()

symb.var('h,t,t_n')
# Pour Adam et N-MS
symb.var('phi_np1,phi_n,phi_nm1,phi_nm2,phi_nm3,phi_nm4,phi_nm5')
# Pour BDF
symb.var('  y_np1,  y_n,  y_nm1,  y_nm2,  y_nm3,  y_nm4,  y_nm5')

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

Schémas d'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. $$ On peut construire différentes 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_{n+1}=u_n+\displaystyle\int_{t_n}^{t_{n+1}} \tilde f(t) \mathrm{d}t \quad \text{où $\tilde f(t)$ est un polynôme interpolant }\varphi(t,y(t)) \end{cases} $$

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érentes schémas selon les points d'interpolation choisis.

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}=u_n+h\sum_{j=n-q+1}^{n}b_j\varphi(t_j,u_j)+hb_{-1}\varphi(t_{n+1},u_{n+1}) \qquad n=q-1,q,\dots,N-1 \end{cases} $$

Ils se divisent en deux familles:

  • les schémas d'Adam-Bashforth à $q$ pas (AB-$q$) sont
    • explicites,
    • d'ordre $q$,
    • 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,
    • d'ordre $q+1$,
    • 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

  1. La notation AM-$k$ n'est pas toujours la même selon les auteurs. Par exemple, dans les livres de Quarteroni et al., le $k$ indique l'ordre de la méthode et non le nombre de pas. Les deux notations coincident pour les methodes AB mais pas pour les méthodes AM.
  2. Notons que, pour calculer successivement $u_{q}$, $u_{q+1}$, ..., il faut d'abord initialiser $u_0,u_1,\dots,u_{q-1}$. Cette initialisation doit être faite par des approximations adéquates qui ne dégradent pas l'ordre de la méthode.
  3. On verra plus tard que
    • la condition de zéro-stabilité est toujours satisfaite,
    • la condition de A-stabilité est de plus en plus contraignante quand l’ordre augmente. Plus particulièrement, une méthode explicite A-stable a au plus ordre 2 (seconde barrière de Dahlquist).

AB-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 \approx \int_{t_n}^{t_{n+1}} p(t) \mathrm{d}t =h\varphi(t_n,y(t_n))$ et on obtient le schéma
    $$ \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 AB$_1$ coïncide avec la méthode d'Euler progressive.
In [4]:
p=symb.interpolate([(t_n,phi_n)], t)
display(p)
symb.integrate(p,(t,t_n,t_np1)).simplify()
$\displaystyle \phi_{n}$
Out[4]:
$\displaystyle h \phi_{n}$

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(t_{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(t_{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(t_{n-1}))\right)$ et on obtient le schéma
    $$ \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}+ \dfrac{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 AB$_1$.
In [5]:
p=symb.interpolate([(t_n,phi_n),(t_nm1,phi_nm1)], t)
# display(p)
symb.integrate(p,(t,t_n,t_np1)).simplify()
Out[5]:
$\displaystyle \frac{h \left(3 \phi_{n} - \phi_{nm1}\right)}{2}$

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{align} 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{align}
  • 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{(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 AB$_1$ et $u_{2}$ est une approximation de $y(t_2)$ obtenue en utilisant la méthode AB$_2$.
In [6]:
p=symb.interpolate([(t_n,phi_n),(t_nm1,phi_nm1),(t_nm2,phi_nm2)], t)
# display(p)
symb.integrate(p,(t,t_n,t_np1)).simplify()
Out[6]:
$\displaystyle \frac{h \left(23 \phi_{n} - 16 \phi_{nm1} + 5 \phi_{nm2}\right)}{12}$

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{(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 AB$_1$, $u_{2}$ est une approximation de $y(t_2)$ obtenue en utilisant la méthode AB$_2$ etc.
In [7]:
p=symb.interpolate([(t_n,phi_n),(t_nm1,phi_nm1),(t_nm2,phi_nm2),(t_nm3,phi_nm3)], t)
symb.integrate(p,(t,t_n,t_np1)).simplify()
Out[7]:
$\displaystyle \frac{h \left(55 \phi_{n} - 59 \phi_{nm1} + 37 \phi_{nm2} - 9 \phi_{nm3}\right)}{24}$

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

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{(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 AM$_0$ coïncide avec la méthode d'Euler regressive.
In [8]:
p=symb.interpolate([(t_np1,phi_np1)], t)
symb.integrate(p,(t,t_n,t_np1)).simplify()
Out[8]:
$\displaystyle h \phi_{np1}$

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{(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 AM$_1$ coïncide avec la méthode de Crank-Nicolson.
In [9]:
p=symb.interpolate([(t_np1,phi_np1),(t_n,phi_n)], t)
symb.integrate(p,(t,t_n,t_np1)).simplify()
Out[9]:
$\displaystyle \frac{h \left(\phi_{n} + \phi_{np1}\right)}{2}$

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{(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 AM$_1$.
In [10]:
p=symb.interpolate([(t_np1,phi_np1),(t_n,phi_n),(t_nm1,phi_nm1)], t)
symb.integrate(p,(t,t_n,t_np1)).simplify()
Out[10]:
$\displaystyle \frac{h \left(8 \phi_{n} - \phi_{nm1} + 5 \phi_{np1}\right)}{12}$

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

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

In [11]:
import sympy as symb
symb.init_printing()

symb.var('h,t,t_n')
symb.var('phi_np1,phi_n,phi_nm1,phi_nm2,phi_nm3,phi_nm4,phi_nm5')
symb.var('  u_np1,  u_n,  u_nm1,  u_nm2,  u_nm3,  u_nm4,  u_nm5')

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

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):
    p=symb.interpolate(Points[:q], t)
    AB=(symb.integrate(p,(t,t_n,t_np1)).simplify())
    print("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)):
    p=symb.interpolate(Points[0:q+1], t)
    AM=(symb.integrate(p,(t,t_n,t_np1)).simplify())
    print("AM-",q, ": implicite, à",q,"pas, d'ordre", q+1)
    display(symb.Eq(u_np1,u_n+AM))
    print("\n")
AB- 1 : explicite, à 1 pas, d'ordre 1
$\displaystyle u_{np1} = h \phi_{n} + u_{n}$

AB- 2 : explicite, à 2 pas, d'ordre 2
$\displaystyle u_{np1} = \frac{h \left(3 \phi_{n} - \phi_{nm1}\right)}{2} + u_{n}$

AB- 3 : explicite, à 3 pas, d'ordre 3
$\displaystyle u_{np1} = \frac{h \left(23 \phi_{n} - 16 \phi_{nm1} + 5 \phi_{nm2}\right)}{12} + u_{n}$

AB- 4 : explicite, à 4 pas, d'ordre 4
$\displaystyle u_{np1} = \frac{h \left(55 \phi_{n} - 59 \phi_{nm1} + 37 \phi_{nm2} - 9 \phi_{nm3}\right)}{24} + u_{n}$

AB- 5 : explicite, à 5 pas, d'ordre 5
$\displaystyle u_{np1} = \frac{h \left(1901 \phi_{n} - 2774 \phi_{nm1} + 2616 \phi_{nm2} - 1274 \phi_{nm3} + 251 \phi_{nm4}\right)}{720} + u_{n}$

AB- 6 : explicite, à 6 pas, d'ordre 6
$\displaystyle u_{np1} = \frac{h \left(4277 \phi_{n} - 7923 \phi_{nm1} + 9982 \phi_{nm2} - 7298 \phi_{nm3} + 2877 \phi_{nm4} - 475 \phi_{nm5}\right)}{1440} + u_{n}$



AM- 0 : implicite, à 0 pas, d'ordre 1
$\displaystyle u_{np1} = h \phi_{np1} + u_{n}$

AM- 1 : implicite, à 1 pas, d'ordre 2
$\displaystyle u_{np1} = \frac{h \left(\phi_{n} + \phi_{np1}\right)}{2} + u_{n}$

AM- 2 : implicite, à 2 pas, d'ordre 3
$\displaystyle u_{np1} = \frac{h \left(8 \phi_{n} - \phi_{nm1} + 5 \phi_{np1}\right)}{12} + u_{n}$

AM- 3 : implicite, à 3 pas, d'ordre 4
$\displaystyle u_{np1} = \frac{h \left(19 \phi_{n} - 5 \phi_{nm1} + \phi_{nm2} + 9 \phi_{np1}\right)}{24} + u_{n}$

AM- 4 : implicite, à 4 pas, d'ordre 5
$\displaystyle u_{np1} = \frac{h \left(646 \phi_{n} - 264 \phi_{nm1} + 106 \phi_{nm2} - 19 \phi_{nm3} + 251 \phi_{np1}\right)}{720} + u_{n}$

AM- 5 : implicite, à 5 pas, d'ordre 6
$\displaystyle u_{np1} = \frac{h \left(1427 \phi_{n} - 798 \phi_{nm1} + 482 \phi_{nm2} - 173 \phi_{nm3} + 27 \phi_{nm4} + 475 \phi_{np1}\right)}{1440} + u_{n}$

AM- 6 : implicite, à 6 pas, d'ordre 7
$\displaystyle u_{np1} = \frac{h \left(65112 \phi_{n} - 46461 \phi_{nm1} + 37504 \phi_{nm2} - 20211 \phi_{nm3} + 6312 \phi_{nm4} - 863 \phi_{nm5} + 19087 \phi_{np1}\right)}{60480} + u_{n}$

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érentes 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+\displaystyle\int_{t_{n-1}}^{t_{n+1}} \tilde f(t) \mathrm{d}t \quad \text{où $\tilde f(t)$ est un polynôme interpolant }\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érentes 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_{q-1}\approx y_{q-1},\\ \displaystyle u_{n+1}=u_{n-1}+h\sum_{j=n-q+1}^{n}b_j\varphi(t_j,u_j)+hb_{-1}\varphi(t_{n+1},u_{n+1}) \qquad n=q-1,q,\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-q+1}\} \text{ avec } q\ge1. $$
  • 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-q+1}\} \text{ avec } q\ge0. $$

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

In [12]:
p=symb.interpolate([(t_n,phi_n)], t).simplify()
symb.integrate(p,(t,t_nm1,t_np1)).simplify()
Out[12]:
$\displaystyle 2 h \phi_{n}$

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$.
In [13]:
p=symb.interpolate([(t_n,phi_n),(t_nm1,phi_nm1)], t).simplify()
symb.integrate(p,(t,t_nm1,t_np1)).simplify()
Out[13]:
$\displaystyle 2 h \phi_{n}$

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$.
In [14]:
p=symb.interpolate([(t_n,phi_n),(t_nm1,phi_nm1),(t_nm2,phi_nm2)], t).simplify()
symb.integrate(p,(t,t_nm1,t_np1)).simplify()
Out[14]:
$\displaystyle \frac{h \left(7 \phi_{n} - 2 \phi_{nm1} + \phi_{nm2}\right)}{3}$

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

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

In [15]:
p=symb.interpolate([(t_np1,phi_np1)], t).simplify()
symb.integrate(p,(t,t_nm1,t_np1)).simplify()
Out[15]:
$\displaystyle 2 h \phi_{np1}$

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

In [16]:
p=symb.interpolate([(t_np1,phi_np1),(t_n,phi_n)], t).simplify()
symb.integrate(p,(t,t_nm1,t_np1)).simplify()
Out[16]:
$\displaystyle 2 h \phi_{n}$

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 AM$_1$ (=Crank-Nicolson).
In [17]:
p=symb.interpolate([(t_np1,phi_np1),(t_n,phi_n),(t_nm1,phi_nm1)], t).simplify()
symb.integrate(p,(t,t_nm1,t_np1)).simplify()
Out[17]:
$\displaystyle \frac{h \left(4 \phi_{n} + \phi_{nm1} + \phi_{np1}\right)}{3}$

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

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

In [18]:
import sympy as symb
symb.init_printing()
symb.var('h,t,t_n')
symb.var('phi_np1,phi_n,phi_nm1,phi_nm2,phi_nm3,phi_nm4,phi_nm5')
symb.var('  u_np1,  u_n,  u_nm1,  u_nm2,  u_nm3,  u_nm4,  u_nm5')
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

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):
    p=symb.interpolate(Points[:q], t)
    N=(symb.integrate(p,(t,t_nm1,t_np1)).simplify())
    print("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)):
    p=symb.interpolate(Points[:q+1], t)
    MS=(symb.integrate(p,(t,t_nm1,t_np1)).simplify())
    print("MS-",q, ": implicite, à",q,"pas")
    display(symb.Eq(u_np1,u_nm1+MS))
    print("\n")
N- 1 : explicite, à 1 pas
$\displaystyle u_{np1} = 2 h \phi_{n} + u_{nm1}$

N- 2 : explicite, à 2 pas
$\displaystyle u_{np1} = 2 h \phi_{n} + u_{nm1}$

N- 3 : explicite, à 3 pas
$\displaystyle u_{np1} = \frac{h \left(7 \phi_{n} - 2 \phi_{nm1} + \phi_{nm2}\right)}{3} + u_{nm1}$

N- 4 : explicite, à 4 pas
$\displaystyle u_{np1} = \frac{h \left(8 \phi_{n} - 5 \phi_{nm1} + 4 \phi_{nm2} - \phi_{nm3}\right)}{3} + u_{nm1}$

N- 5 : explicite, à 5 pas
$\displaystyle u_{np1} = \frac{h \left(269 \phi_{n} - 266 \phi_{nm1} + 294 \phi_{nm2} - 146 \phi_{nm3} + 29 \phi_{nm4}\right)}{90} + u_{nm1}$

N- 6 : explicite, à 6 pas
$\displaystyle u_{np1} = \frac{h \left(297 \phi_{n} - 406 \phi_{nm1} + 574 \phi_{nm2} - 426 \phi_{nm3} + 169 \phi_{nm4} - 28 \phi_{nm5}\right)}{90} + u_{nm1}$



MS- 0 : implicite, à 0 pas
$\displaystyle u_{np1} = 2 h \phi_{np1} + u_{nm1}$

MS- 1 : implicite, à 1 pas
$\displaystyle u_{np1} = 2 h \phi_{n} + u_{nm1}$

MS- 2 : implicite, à 2 pas
$\displaystyle u_{np1} = \frac{h \left(4 \phi_{n} + \phi_{nm1} + \phi_{np1}\right)}{3} + u_{nm1}$

MS- 3 : implicite, à 3 pas
$\displaystyle u_{np1} = \frac{h \left(4 \phi_{n} + \phi_{nm1} + \phi_{np1}\right)}{3} + u_{nm1}$

MS- 4 : implicite, à 4 pas
$\displaystyle u_{np1} = \frac{h \left(124 \phi_{n} + 24 \phi_{nm1} + 4 \phi_{nm2} - \phi_{nm3} + 29 \phi_{np1}\right)}{90} + u_{nm1}$

MS- 5 : implicite, à 5 pas
$\displaystyle u_{np1} = \frac{h \left(129 \phi_{n} + 14 \phi_{nm1} + 14 \phi_{nm2} - 6 \phi_{nm3} + \phi_{nm4} + 28 \phi_{np1}\right)}{90} + u_{nm1}$

MS- 6 : implicite, à 6 pas
$\displaystyle u_{np1} = \frac{h \left(5640 \phi_{n} + 33 \phi_{nm1} + 1328 \phi_{nm2} - 807 \phi_{nm3} + 264 \phi_{nm4} - 37 \phi_{nm5} + 1139 \phi_{np1}\right)}{3780} + u_{nm1}$

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

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{(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 BDF$_1$ coïncide avec la méthode d'Euler implicite.
In [19]:
p=symb.interpolate([(t_np1,y_np1),(t_n,y_n)], t).factor()
print("p(t)")
display(p)
dp=symb.diff(p,t).subs(t,t_np1)
print("p'(t)")
display(dp)
print("y_np1=")
symb.solve(dp-phi_np1,y_np1)[0]
p(t)
$\displaystyle - \frac{- h y_{n} + t y_{n} - t y_{np1} - t_{n} y_{n} + t_{n} y_{np1}}{h}$
p'(t)
$\displaystyle - \frac{y_{n} - y_{np1}}{h}$
y_np1=
Out[19]:
$\displaystyle h \phi_{np1} + y_{n}$

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}[t] 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{(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}$$
In [20]:
p=symb.interpolate([(t_np1,y_np1),(t_n,y_n),(t_nm1,y_nm1)], t).simplify()
print("p(t)")
display(p)
dp=symb.diff(p,t).subs(t,t_np1)
print("p'(t)")
display(dp)
print("y_np1=")
symb.solve(dp-phi_np1,y_np1)[0]
p(t)
$\displaystyle \frac{h^{2} y_{n} + \frac{h \left(- t y_{nm1} + t y_{np1} + t_{n} y_{nm1} - t_{n} y_{np1}\right)}{2} - t^{2} y_{n} + \frac{t^{2} y_{nm1}}{2} + \frac{t^{2} y_{np1}}{2} + 2 t t_{n} y_{n} - t t_{n} y_{nm1} - t t_{n} y_{np1} - t_{n}^{2} y_{n} + \frac{t_{n}^{2} y_{nm1}}{2} + \frac{t_{n}^{2} y_{np1}}{2}}{h^{2}}$
p'(t)
$\displaystyle \frac{\frac{h \left(- y_{nm1} + y_{np1}\right)}{2} + 2 t_{n} y_{n} - t_{n} y_{nm1} - t_{n} y_{np1} - 2 y_{n} \left(h + t_{n}\right) + y_{nm1} \left(h + t_{n}\right) + y_{np1} \left(h + t_{n}\right)}{h^{2}}$
y_np1=
Out[20]:
$\displaystyle \frac{2 h \phi_{np1}}{3} + \frac{4 y_{n}}{3} - \frac{y_{nm1}}{3}$

Calcul systématique des coefficients des méthodes BDF

In [21]:
import sympy as symb
symb.init_printing()
symb.var('h,t,t_n')
symb.var('phi_np1')
symb.var('  u_np1,  u_n,  u_nm1,  u_nm2,  u_nm3,  u_nm4,  u_nm5')
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

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)):
    p=symb.interpolate(Points[0:q+1], t)
    dp=symb.diff(p,t).subs(t,t_np1)
    sol=symb.solve(dp-phi_np1,u_np1)
    #display(sol)
    print("BDF-",q, ": implicite, à",q,"pas")
    display(symb.Eq(u_np1,sol[0]))
    print("\n")
BDF- 1 : implicite, à 1 pas
$\displaystyle u_{np1} = h \phi_{np1} + u_{n}$

BDF- 2 : implicite, à 2 pas
$\displaystyle u_{np1} = \frac{2 h \phi_{np1}}{3} + \frac{4 u_{n}}{3} - \frac{u_{nm1}}{3}$

BDF- 3 : implicite, à 3 pas
$\displaystyle u_{np1} = \frac{6 h \phi_{np1}}{11} + \frac{18 u_{n}}{11} - \frac{9 u_{nm1}}{11} + \frac{2 u_{nm2}}{11}$

BDF- 4 : implicite, à 4 pas
$\displaystyle u_{np1} = \frac{12 h \phi_{np1}}{25} + \frac{48 u_{n}}{25} - \frac{36 u_{nm1}}{25} + \frac{16 u_{nm2}}{25} - \frac{3 u_{nm3}}{25}$

BDF- 5 : implicite, à 5 pas
$\displaystyle u_{np1} = \frac{60 h \phi_{np1}}{137} + \frac{300 u_{n}}{137} - \frac{300 u_{nm1}}{137} + \frac{200 u_{nm2}}{137} - \frac{75 u_{nm3}}{137} + \frac{12 u_{nm4}}{137}$

BDF- 6 : implicite, à 6 pas
$\displaystyle u_{np1} = \frac{20 h \phi_{np1}}{49} + \frac{120 u_{n}}{49} - \frac{150 u_{nm1}}{49} + \frac{400 u_{nm2}}{147} - \frac{75 u_{nm3}}{49} + \frac{24 u_{nm4}}{49} - \frac{10 u_{nm5}}{147}$

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 comme suit.

On considère une méthode implicite $u_{n+1}=u_n+hG(u_n,u_{n-1}\dots,u_{n-q+1})+hb_{-1}\varphi(t_{n+1},u_{n+1})$:
  • étape de prédiction: on calcule $\tilde u_{n+1}$ une approximation de $u_{n+1}$ par une méthode explicite;
  • étape de correction: on "corrige" la méthode implicite en approchant $\varphi(t_{n+1},u_{n+1})$ par $\varphi(t_{n+1},\tilde u_{n+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.

Exemple : schéma de Heun

Considérons les deux schémas PC suivants:

  • predictor: méthode d'Euler explicite (ou AB$_1$) $\tilde u_{n+1}=u_n+h\varphi(t_n,u_n)$
  • corrector: méthode de Crank-Nicolson (ou 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) & n=1,2,\dots N-1\\ 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}$$

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 AB$_2$ $\tilde 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 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 AB$_2$-AM$_2$:
    $$ \text{(AB$_2$-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)& n=2,3,\dots N-1\\ 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}$$