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]:

M62_CM6-7 : Schémas multistep $\newcommand{\R}{\mathbb{R}}$ $\newcommand{\N}{\mathbb{N}}$ $\newcommand{\dt}{\mathrm{d}t}$ $\newcommand{\dx}{\mathrm{d}x}$ $\newcommand{\dtau}{\mathrm{d}\tau}$ $\newcommand{\CC}[1]{\mathscr{C}}$

Table of Contents

Construction de schémas

Considérons le problème de Cauchy

trouver une fonction $y \colon I\subset \R \to \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\stackrel{\text{déf.}}{=} t_0+nh$ avec $n=0,1,2,\dots,N_h$ une suite de $N_h+1$ nœuds de $I$ induisant une discrétisation de $I$ en $N_h$ sous-intervalles $I_n=[t_n;t_{n+1}]$ chacun de longueur $h>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_h+1$ valeurs $\{t_0, t_1=t_0+h,... , t_{N_h}=T \}$ représente les points de la discrétisation.
  • L'ensemble de $N_h+1$ valeurs $\{y_0, y_1,... , y_{N_h} \}$ représente la solution exacte.
  • L'ensemble de $N_h+1$ valeurs $\{u_0 = y_0, u_1,... , u_{N_h} \}$ représente la solution numérique. Cette solution approchée sera 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,\\ u_{n+1}=\Phi(u_{n+1},u_n, u_{n-1}, ..., u_{n-k}),&\forall n\in\N. \end{cases} Les schémas qu'on va construire permettent ainsi 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$,\dots, en partant de $u_0$.

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).
Une méthode est dite implicite si $u_{n+1}$ n'est défini que par une relation implicite faisant intervenir la fonction $\varphi$.

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

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 [2]:
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('  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}} \text{un polynôme d'interpolation de }\varphi(t,u) \dt. \end{cases} Les schémas d'Adam approchent l'intégrale $\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\dt$ 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. Ils se divisent en deux familles: les méthodes d'Adam-Bashforth qui sont explicites et les méthodes d'Adam-Moulton qui sont implicites:

  • schémas AB-$q$: les schémas d'Adam-Bashforth d'ordre $q$ approchent l'intégrale $\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\dt$ par l'intégrale $\int_{t_n}^{t_{n+1}} p(t) \dt$ où $p$ est le polynôme interpolant $\varphi$ en $$ \{t_n,t_{n-1},t_{n+1-q}\} \text{ avec } q\ge1; $$
  • schémas AM-$q$: les schémas d'Adam-Moulton d'ordre $q$ approchent l'intégrale $\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\dt$ par l'intégrale $\int_{t_n}^{t_{n+1}} p(t) \dt$ où $p$ est le polynôme interpolant $\varphi$ en $$ \{t_{n+1},t_n,t_{n-1},t_{n+1-q}\} \text{ avec } q\ge0. $$

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 car seul $u_0$ est donné.

Remarque: la condition de A-stabilité est de plus en plus contraignante quand l’ordre augmente.

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) \dt =h\varphi(t_n,y(t_n))$ et on obtient le schéma \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 [3]:
p=symb.interpolate([(t_n,phi_n)], t)
symb.integrate(p,(t,t_n,t_np1)).simplify()
Out[3]:
$$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) \dt =\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 \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} \left(3\varphi(t_{n},u_{n})-\varphi(t_{n-1},u_{n-1})\right)& 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 [4]:
p=symb.interpolate([(t_n,phi_n),(t_nm1,phi_nm1)], t)
symb.integrate(p,(t,t_n,t_np1)).simplify()
Out[4]:
$$\frac{h}{2} \left(3 \phi_{n} - \phi_{nm1}\right)$$

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-1},\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) \dt =\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 \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 + \frac{h}{2} \left(3\varphi(t_1,u_1)-\varphi(t_0,u_0)\right)\approx y(t_2)\\ u_{n+1} = u_{n} + \frac{h}{12}\left(23\varphi(t_n,u_n)-16\varphi(t_{n-1},u_{n-1}+5\varphi(t_{n-2},u_{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 AB$_1$ et $u_{2}$ est une approximation de $y(t_2)$ obtenue en utilisant la méthode AB$_2$.
In [5]:
p=symb.interpolate([(t_n,phi_n),(t_nm1,phi_nm1),(t_nm2,phi_nm2)], t)
symb.integrate(p,(t,t_n,t_np1)).simplify()
Out[5]:
$$\frac{h}{12} \left(23 \phi_{n} - 16 \phi_{nm1} + 5 \phi_{nm2}\right)$$

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-1},\varphi(t_{n-2},y_{n-2}))\}$
  • Polynôme: $p(t)=\displaystyle\sum_{i=n-2}^{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) \dt =\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 \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 + \frac{h}{2} \left(3\varphi(t_1,u_1)-\varphi(t_0,u_0)\right)\approx y(t_2)\\ u_{3} = u_{2} + \frac{h}{12}\left(23\varphi(t_2,u_2)-16\varphi(t_{1},u_{1}+5\varphi(t_{0},u_{0})\right)\approx y(t_3) \\ u_{n+1} = u_{n} + \frac{h}{24}\left(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})\right) & 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 [6]:
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[6]:
$$\frac{h}{24} \left(55 \phi_{n} - 59 \phi_{nm1} + 37 \phi_{nm2} - 9 \phi_{nm3}\right)$$

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) \dt =h\varphi(t_{n+1},y(t_{n+1}))$ et on obtient le schéma \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 [7]:
p=symb.interpolate([(t_np1,phi_np1)], t)
symb.integrate(p,(t,t_n,t_np1)).simplify()
Out[7]:
$$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-tn}{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) \dt =\frac{h}{2}\left(\varphi(t_{n+1},y_{n+1})+\varphi(t_n,y_n)\right)$ et on obtient le schéma \begin{cases} u_0=y(t_0)=y_0,\\ u_{n+1}=u_{n}+ \frac{h}{2} \left(\varphi(t_{n},u_{n})+\varphi(t_{n+1},u_{n+1})\right)& n=0,1,\dots N-1 \end{cases} La méthode AM$_1$ coïncide avec la méthode de Crank-Nicolson.
In [8]:
p=symb.interpolate([(t_np1,phi_np1),(t_n,phi_n)], t)
symb.integrate(p,(t,t_n,t_np1)).simplify()
Out[8]:
$$\frac{h}{2} \left(\phi_{n} + \phi_{np1}\right)$$

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-tn)(t-t{n-1})}{(t_{n+1}-tn)(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) \dt =\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 \begin{cases} u_0=y(t_0)=y_0,\\ u_1=u_{0}+ \frac{h}{2} \left(\varphi(t_{1},u_{1})+\varphi(t_{0},u_{0})\right)\\ 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) & 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 [9]:
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[9]:
$$\frac{h}{12} \left(8 \phi_{n} - \phi_{nm1} + 5 \phi_{np1}\right)$$

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

In [10]:
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('  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

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)):
    p=symb.interpolate(Points[0:q+1], t)
    AB=(symb.integrate(p,(t,t_n,t_np1)).simplify())
    print("AB-",q)
    display(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)
    display(AM)
    print("\n")

AB- 1
$$\frac{h}{2} \left(3 \phi_{n} - \phi_{nm1}\right)$$ 
AB- 2
$$\frac{h}{12} \left(23 \phi_{n} - 16 \phi_{nm1} + 5 \phi_{nm2}\right)$$ 
AB- 3
$$\frac{h}{24} \left(55 \phi_{n} - 59 \phi_{nm1} + 37 \phi_{nm2} - 9 \phi_{nm3}\right)$$ 
AB- 4
$$\frac{h}{720} \left(1901 \phi_{n} - 2774 \phi_{nm1} + 2616 \phi_{nm2} - 1274 \phi_{nm3} + 251 \phi_{nm4}\right)$$ 
AB- 5
$$\frac{h}{1440} \left(4277 \phi_{n} - 7923 \phi_{nm1} + 9982 \phi_{nm2} - 7298 \phi_{nm3} + 2877 \phi_{nm4} - 475 \phi_{nm5}\right)$$ 


AM- 0
$$h \phi_{np1}$$ 
AM- 1
$$\frac{h}{2} \left(\phi_{n} + \phi_{np1}\right)$$ 
AM- 2
$$\frac{h}{12} \left(8 \phi_{n} - \phi_{nm1} + 5 \phi_{np1}\right)$$ 
AM- 3
$$\frac{h}{24} \left(19 \phi_{n} - 5 \phi_{nm1} + \phi_{nm2} + 9 \phi_{np1}\right)$$ 
AM- 4
$$\frac{h}{720} \left(646 \phi_{n} - 264 \phi_{nm1} + 106 \phi_{nm2} - 19 \phi_{nm3} + 251 \phi_{np1}\right)$$ 
AM- 5
$$\frac{h}{1440} \left(1427 \phi_{n} - 798 \phi_{nm1} + 482 \phi_{nm2} - 173 \phi_{nm3} + 27 \phi_{nm4} + 475 \phi_{np1}\right)$$ 
AM- 6
$$\frac{h}{60480} \left(65112 \phi_{n} - 46461 \phi_{nm1} + 37504 \phi_{nm2} - 20211 \phi_{nm3} + 6312 \phi_{nm4} - 863 \phi_{nm5} + 19087 \phi_{np1}\right)$$ 




















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-1}+\displaystyle\int_{t_{n-1}}^{t_{n+1}} \text{un polynôme d'interpolation de }\varphi(t,u) \dt.
\end{cases}


On obtient les schémas de Nyström, qui sont explicites, et les schémas  de Milne-Simpson, qui sont implicites:


    

  • schémas N-$q$: les schémas de Nyström d'ordre $q$ approchent l'intégrale $\int_{t_{n-1}}^{t_{n+1}} \varphi(t,y(t))\dt$ par l'intégrale $\int_{t_{n-1}}^{t_{n+1}} p(t) \dt$ où $p$ est le polynôme interpolant $\varphi$ en 
$$
\{t_n,t_{n-1},t_{n+1-q}\} \text{ avec } q\ge1;
$$
    • 

  • schémas MS-$q$: les schémas de Milne-Simpson d'ordre $q$ approchent l'intégrale $\int_{t_{n-1}}^{t_{n+1}} \varphi(t,y(t))\dt$ par l'intégrale $\int_{t_{n-1}}^{t_{n+1}} p(t) \dt$ où $p$ est le polynôme interpolant $\varphi$  en 
$$
\{t_{n+1},t_n,t_{n-1},t_{n+1-q}\} \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) \dt =2h\varphi(t_n,y(t_n))$
et on obtient le schéma
\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). 
    • 














In [11]:



    


p=symb.interpolate([(t_n,phi_n)], t).simplify()
symb.integrate(p,(t,t_nm1,t_np1)).simplify()



    















    
Out[11]:







$$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) \dt =2h\varphi(t_n,y(t_n))$
et on obtient le schéma
\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 [12]:



    


p=symb.interpolate([(t_n,phi_n),(t_nm1,phi_nm1)], t).simplify()
symb.integrate(p,(t,t_nm1,t_np1)).simplify()



    















    
Out[12]:







$$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-1},\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) \dt =\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
\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} + \frac{h}{3}\left(7\varphi(t_n,u_n)-2\varphi(t_{n-1},u_{n-1}+\varphi(t_{n-2},u_{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$.  
    • 














In [13]:



    


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[13]:







$$\frac{h}{3} \left(7 \phi_{n} - 2 \phi_{nm1} + \phi_{nm2}\right)$$




















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) \dt =2h\varphi(t_{n+1},y(t_{n+1}))$
et on obtient le schéma
\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+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 explicite.
    • 














In [14]:



    


p=symb.interpolate([(t_np1,phi_np1)], t).simplify()
symb.integrate(p,(t,t_nm1,t_np1)).simplify()



    















    
Out[14]:







$$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-tn}{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) \dt =2h\varphi(t_n,y_n)$
et on obtient le schéma
\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.  
    • 














In [15]:



    


p=symb.interpolate([(t_np1,phi_np1),(t_n,phi_n)], t).simplify()
symb.integrate(p,(t,t_nm1,t_np1)).simplify()



    















    
Out[15]:







$$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-tn)(t-t{n-1})}{(t_{n+1}-tn)(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) \dt =\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
\begin{cases}
u_0=y(t_0)=y_0,\\
u_1=u_{0}+ \frac{h}{2} \left(\varphi(t_{1},u_{1})+\varphi(t_{0},u_{0})\right)\\
u_{n+1}=u_{n-1}+\frac{h}{3}\left( \varphi(t_{n+1},u_{n+1})+4\varphi(t_n,u_n)+\varphi(t_{n-1},u_{n-1}\right) & 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 [16]:



    


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[16]:







$$\frac{h}{3} \left(4 \phi_{n} + \phi_{nm1} + \phi_{np1}\right)$$




















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












In [17]:



    


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('  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

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)):
    p=symb.interpolate(Points[0:q+1], t)
    N=(symb.integrate(p,(t,t_nm1,t_np1)).simplify())
    print("N-",q)
    display(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[0:q+1], t)
    MS=(symb.integrate(p,(t,t_nm1,t_np1)).simplify())
    print("MS-",q)
    display(MS)
    print("\n")



    















    






N- 1










    







$$2 h \phi_{n}$$









    






N- 2










    







$$\frac{h}{3} \left(7 \phi_{n} - 2 \phi_{nm1} + \phi_{nm2}\right)$$









    






N- 3










    







$$\frac{h}{3} \left(8 \phi_{n} - 5 \phi_{nm1} + 4 \phi_{nm2} - \phi_{nm3}\right)$$









    






N- 4










    







$$\frac{h}{90} \left(269 \phi_{n} - 266 \phi_{nm1} + 294 \phi_{nm2} - 146 \phi_{nm3} + 29 \phi_{nm4}\right)$$









    






N- 5










    







$$\frac{h}{90} \left(297 \phi_{n} - 406 \phi_{nm1} + 574 \phi_{nm2} - 426 \phi_{nm3} + 169 \phi_{nm4} - 28 \phi_{nm5}\right)$$









    








MS- 0










    







$$2 h \phi_{np1}$$









    






MS- 1










    







$$2 h \phi_{n}$$









    






MS- 2










    







$$\frac{h}{3} \left(4 \phi_{n} + \phi_{nm1} + \phi_{np1}\right)$$









    






MS- 3










    







$$\frac{h}{3} \left(4 \phi_{n} + \phi_{nm1} + \phi_{np1}\right)$$









    






MS- 4










    







$$\frac{h}{90} \left(124 \phi_{n} + 24 \phi_{nm1} + 4 \phi_{nm2} - \phi_{nm3} + 29 \phi_{np1}\right)$$









    






MS- 5










    







$$\frac{h}{90} \left(129 \phi_{n} + 14 \phi_{nm1} + 14 \phi_{nm2} - 6 \phi_{nm3} + \phi_{nm4} + 28 \phi_{np1}\right)$$









    






MS- 6










    







$$\frac{h}{3780} \left(5640 \phi_{n} + 33 \phi_{nm1} + 1328 \phi_{nm2} - 807 \phi_{nm3} + 264 \phi_{nm4} - 37 \phi_{nm5} + 1139 \phi_{np1}\right)$$









    


























Schémas BDF (Backward Differentiation Formulae)
Si nous interpolons la fonction inconnue $t\mapsto y(t)$ par un polynôme $p$  en des points donnés 
$$
\{t_{n+1},t_n,\dots, t_{n+1-i}\}\text{ avec } q\ge 1,
$$ 
nous obtenons $y(t) \simeq p(t)$ et nous pouvons écrire
$$
\varphi(t,y) = y'(t) \simeq p'(t).
$$
En évaluant cette relation en $t_{n+1}$ nous obtenons 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}) \simeq p'(t_{n+1}).
\end{cases}
$$


Remarque: la condition de A-stabilité est de plus en plus contraignante quand l’ordre augmente.















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
\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 [18]:



    


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)










    







$$- \frac{1}{h} \left(- h y_{n} + t y_{n} - t y_{np1} - t_{n} y_{n} + t_{n} y_{np1}\right)$$









    






p'(t)










    







$$- \frac{1}{h} \left(y_{n} - y_{np1}\right)$$









    






y_np1=










    
Out[18]:







$$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
\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 [19]:



    


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)










    







$$\frac{1}{h^{2}} \left(h^{2} y_{n} + \frac{h}{2} \left(- t y_{nm1} + t y_{np1} + t_{n} y_{nm1} - t_{n} y_{np1}\right) - 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}\right)$$









    






p'(t)










    







$$\frac{1}{h^{2}} \left(\frac{h}{2} \left(- y_{nm1} + y_{np1}\right) + 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)\right)$$









    






y_np1=










    
Out[19]:







$$\frac{2 h}{3} \phi_{np1} + \frac{4 y_{n}}{3} - \frac{y_{nm1}}{3}$$




















Calcul systématique des coefficients des méthodes BDF












In [20]:



    


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)):
    print("BDF-",q)
    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("\n")



    















    






BDF- 1










    







$$\left [ h \phi_{np1} + u_{n}\right ]$$









    






BDF- 2










    







$$\left [ \frac{2 h}{3} \phi_{np1} + \frac{4 u_{n}}{3} - \frac{u_{nm1}}{3}\right ]$$









    






BDF- 3










    







$$\left [ \frac{6 h}{11} \phi_{np1} + \frac{18 u_{n}}{11} - \frac{9 u_{nm1}}{11} + \frac{2 u_{nm2}}{11}\right ]$$









    






BDF- 4










    







$$\left [ \frac{12 h}{25} \phi_{np1} + \frac{48 u_{n}}{25} - \frac{36 u_{nm1}}{25} + \frac{16 u_{nm2}}{25} - \frac{3 u_{nm3}}{25}\right ]$$









    






BDF- 5










    







$$\left [ \frac{60 h}{137} \phi_{np1} + \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}\right ]$$









    






BDF- 6










    







$$\left [ \frac{20 h}{49} \phi_{np1} + \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}\right ]$$









    


























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+1};u_n,...)$:


    

  • é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 $G(u_{n+1};u_n,...)$ par $G(\tilde u_{n+1})$.
    • 



Remarque

Ces méthodes héritent de l’ordre de précision du correcteur. Cependant, étant explicites, elles sont soumises à une condition de 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
La méthode de Heun
\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}
est construite par les deux étapes suivantes:


    

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



Exemple : AB-AM
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$:
\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,un)-\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} 
    • 
















Exercice : Interpolation, Quadrature et EDO
    

  • Soit $f$ une fonction de classe $\CC^1([-1,1])$. Écrire le polynôme $p\in\R_2[\tau]$ qui interpole $f$ aux points $-1$, $0$ et $1$.
    • 

  • Construire une méthode de quadrature comme suit:
$$
\int_{0}^{1}f(\tau)\dtau \approx \int_{0}^{1}p(\tau)\dtau.
$$
NB: on intègre sur $[0,1]$ mais on interpole en $-1$, $0$ et $1$.
    

    • 

  • À l'aide d'un changement de variable affine entre l'intervalle $[0,1]$ et l'intervalle $[a,b]$, en déduire une formule de quadrature pour l'intégrale
$$
\int_{a}^{b}f(x) \dx
$$
lorsque $f$ est une fonction de classe $\CC^1([2a-b,b])$.
    
Remarque: $[2a-b,b]=[a-(b-a),a+(b-a)]$
    

    • 

  • Considérons le problème de Cauchy:  trouver $y \colon [t_0,T]\subset \R \to \R$ tel que
  \begin{cases}
  y'(t) = \varphi(t,y(t)), &\forall t \in [t_0,T],\\
  y(t_0) = y_0,
  \end{cases}
  dont on suppose l'existence d'une unique solution $y$.
    

    On subdivise l'intervalle $[t_0;T]$ en $N$ intervalles $[t_{n};t_{n+1}]$ de largeur $h=\dfrac{T-t_0}{N}$ avec $t_n=t_0+nh$ pour $n=0,\dots,N$.
  Utiliser la formule obtenue au point 3 pour approcher l'intégrale
  $$
  \int_{t_n}^{t_{n+1}}\varphi(t,y(t))\dt.
  $$
  En déduire un schéma à deux pas implicite pour l'approximation de la solution du problème  de Cauchy.
    

    • 
















Correction
On cherche les coefficients $\alpha$, $\beta$ et $\gamma$ du polynôme $p(\tau)=\alpha+\beta \tau+\gamma \tau^2$ tels que
$$
\begin{cases}
p(-1)=f(-1),\\
p(0) =f(0),\\
p(1) =f(1),
\end{cases}
\qquad\text{c'est à dire}\qquad
\begin{cases}
\alpha-\beta+\gamma=f(-1),\\
\alpha=f(0),\\
\alpha+\beta+\gamma=f(1).
\end{cases}
$$
Donc $\alpha=f(0)$, $\beta=\frac{f(1)-f(-1)}{2}$ et $\gamma=\frac{f(1)-2f(0)+f(-1)}{2}$.













In [21]:



    


import sympy as symb
symb.init_printing()

symb.var('f_0,f_1,f_m1,')

p=symb.interpolate([(1,f_1),(0,f_0),(-1,f_m1)], t).factor(t)
p



    















    
Out[21]:







$$- \frac{1}{2} \left(- 2 f_{0} + t^{2} \left(2 f_{0} - f_{1} - f_{m1}\right) + t \left(- f_{1} + f_{m1}\right)\right)$$




















Correction
On en déduit la méthode de quadrature
$$
\int_{0}^{1}f(\tau)\dtau
\approx


\int_{0}^{1}p(\tau)\dtau
\alpha+\frac{\beta}{2}+\frac{\gamma}{3}
\frac{-f(-1)+8f(0)+5f(1)}{12}.
$$













In [22]:



    


q=(symb.integrate(p,(t,0,1)).simplify())
q



    















    
Out[22]:







$$\frac{2 f_{0}}{3} + \frac{5 f_{1}}{12} - \frac{f_{m1}}{12}$$




















Correction
Soit $x=m\tau+q$, alors
$$


\int_{a}^{b} f(x)\dx
m\int_{0}^{1}f(m\tau+q)\dtau
\quad\text{avec}\quad
\begin{cases}
a=q,\\
b=m+q,
\end{cases}
\quad\text{\ie}\quad
\begin{cases}
m=b-a,\\
q=a,
\end{cases}
$$
d'où le changement de variable $x=(b-a)\tau+a$. On en déduit la formule de quadrature
$$


\int_{a}^{b}!!!!f(x)\dx
(b-a)\int_{0}^{1}f\left((b-a)\tau+a\right)\dtau
\approx
(b-a)\frac{-f(2a-b)+8f(a)+5f(b)}{12}.
$$















Correction
On pose $a=t_n$ et $b=t_{n+1}$ d'où la formule de quadrature
$$
\int{t{n}}^{t_{n+1}}!!!!f(t)\dt
\approx


(t{n+1}-t{n}) \frac{-f(2t{n}-t{n+1})+8f(t{n})+5f(t{n+1})}{12}
h \frac{-f(t{n-1})+8f(t{n})+5f(t_{n+1})}{12}.
$$
En utilisant la formule de quadrature pour l'intégration de l'EDO $y'(t)=\varphi(t,y(t))$ entre $t_n$ et $t_{n+1}$ on obtient
$$
y(t_{n+1})=y(tn)+\int{tn}^{t{n+1}} \varphi(t,y(t))\dt
\approx
h \dfrac{-\varphi(t{n-1},y(t{n-1}))+8\varphi(t{n},y(t{n}))+5\varphi(t{n+1},y(t{n+1}))}{12}.
$$
Si on note $u_n$ une approximation de la solution $y$ au temps $t_n$, on obtient le schéma à deux pas implicite suivant:
$$
\begin{cases}
u_0=y(t_0)=y_0,\\
u_1\text{ à définir},\\
u_{n+1}=u_n+ h \dfrac{-\varphi(t_{n-1},u_{n-1})+8\varphi(t_{n},u_{n})+5\varphi(t_{n+1},u_{n+1})}{12} & n=1,2,\dots N-1
\end{cases}
$$
On peut utiliser une pédiction d'Euler explicite pour initialiser $u_1$:
$$
\begin{cases}
u_0=y(t_0)=y_0,\\
u_1=u_0+h\varphi(t_0,u_0),\\
u_{n+1}=u_n+ h \dfrac{-\varphi(t_{n-1},u_{n-1})+8\varphi(t_{n},u_{n})+5\varphi(t_{n+1},u_{n+1})}{12} & n=1,2,\dots N-1
\end{cases}
$$















Exercice : Interpolation, Quadrature et EDO
    

  1. Soit $h>0$ et $f\colon [a-h,a+h] \to\R$ une fonction de classe $\mathscr{C}^1([a-h,a+h])$. Écrire le polynôme $p\in\R_1[x]$ qui interpole $f$ aux points $a-h$ et $a$, i.e. l'équation de la droite $p\in\R_1[x]$ qui passe par les deux points $(a-h,f(a-h))$ et $(a,f(a))$.
    2. 

  2. Construire une méthode de quadrature comme suit:
$$
\int_{a}^{a+h}f(x)\dx \approx \int_{a}^{a+h}p(x)\dx.
$$
NB: on intègre sur $[a,a+h]$ mais on interpole en $a-h$ et $a$.
    3. 

  3. Considérons le problème de Cauchy:  trouver $y \colon [t_0,T]\subset \R \to \R$ tel que
 \begin{cases}
 y'(t) = \varphi(t,y(t)), &\forall t \in [t_0,T],\\
 y(t_0) = y_0,
 \end{cases}
 dont on suppose l'existence d'une unique solution $y$.
    
 On subdivise l'intervalle $[t_0;T]$ en $N$ intervalles $[t_{n};t_{n+1}]$ de largeur $h=\dfrac{T-t_0}{N}$ avec $t_n=t_0+nh$ pour $n=0,\dots,N$.
    
 Utiliser la formule obtenue au point 2 pour approcher l'intégrale
 $$
 \int_{t_n}^{t_{n+1}}\varphi(t,y(t))\dt.
 $$
 En déduire un schéma à deux pas explicite pour l'approximation de la solution du problème  de Cauchy.
    4. 
















Correction
$p(x)=\dfrac{f(a)-f(a-h)}{a-(a-h)}(x-a)+f(a)=\dfrac{f(a)-f(a-h)}{h}(x-a)+f(a)$.













In [23]:



    


import sympy as symb
symb.init_printing()

symb.var('a,h,f_a,f_amh')

p=symb.interpolate([(a,f_a),(a-h,f_amh)], t).factor(t)
p



    















    
Out[23]:







$$\frac{1}{h} \left(- a f_{a} + a f_{amh} + f_{a} h + t \left(f_{a} - f_{amh}\right)\right)$$




















Correction
On en déduit la méthode de quadrature
\begin{align*}
\int_{a}^{a+h}f(x)\dx
&\approx
\int_{a}^{a+h}p(x)\dx
\\&=
\dfrac{f(a)-f(a-h)}{h}\left[ \frac{(x-a)^2}{2} \right]_a^{a+h}+f(a)\left[ x \right]_a^{a+h}
\\&=
\dfrac{f(a)-f(a-h)}{2h}\left((a+h-a)^2-(a-a)^2 \right)+f(a)(a+h-a)
\\&=
\dfrac{f(a)-f(a-h)}{2h} h^2 + hf(a)
\\&=
h \dfrac{3f(a)-f(a-h)}{2}.
\end{align*}













In [24]:



    


q=(symb.integrate(p,(t,a,a+h)).simplify())
q



    















    
Out[24]:







$$\frac{h}{2} \left(3 f_{a} - f_{amh}\right)$$




















Correction
On pose $a=t_n$ et $a+h=t_{n+1}$ d'où la formule de quadrature
$$
\int{t{n}}^{t_{n+1}}!!!!f(t)\dt
\approx


(t{n+1}-t{n}) \frac{3f(tn)-f(2t{n}-t_{n+1})}{2}
h \frac{3f(tn)-f(t{n-1})}{2}.
$$
En utilisant la formule de quadrature pour l'intégration de l'EDO $y'(t)=\varphi(t,y(t))$ entre $t_n$ et $t_{n+1}$ on obtient
$$
y(t_{n+1})=y(tn)+\int{tn}^{t{n+1}} \varphi(t,y(t))\dt
\approx
h \dfrac{3\varphi(t{n},y(t{n}))-\varphi(t{n-1},y(t{n-1}))}{2}.
$$
Si on note $u_n$ une approximation de la solution $y$ au temps $t_n$, on obtient le schéma à deux pas implicite suivant:
$$
\begin{cases}
u_0=y(t_0)=y_0,\\
u_1\text{ à définir},\\
u_{n+1}=u_n+ h \dfrac{3\varphi(t_{n},u_{n})-\varphi(t_{n-1},u_{n-1})}{2} & n=1,2,\dots N-1
\end{cases}
$$
On peut utiliser une prédiction d'Euler explicite pour initialiser $u_1$:
$$
\begin{cases}
u_0=y(t_0)=y_0,\\
u_1=u_0+h\varphi(t_0,u_0),\\
u_{n+1}=u_n+ h \dfrac{3\varphi(t_{n},u_{n})-\varphi(t_{n-1},u_{n-1})}{2} & n=1,2,\dots N-1
\end{cases}
$$















Exercice: Construction d'un schéma
Considérons le problème de Cauchy


trouver une fonction $y \colon I\subset \R \to \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$.




On subdivise l'intervalle $I=[t_0;T]$, avec $T<+\infty$, en $N$ intervalles $[t_{n};t_{n+1}]$ de largeur $h=(T-t_0)/N$ avec $t_n=t_0+nh$ pour $n=0,\dots,N$.

Pour chaque nœud $t_n$, on note $y_n=y(t_n)$ et on cherche la valeur inconnue $u_n$ qui approche la valeur exacte $y_n$


Si nous intégrons l'EDO $y'(t)=\varphi(t,y(t))$ entre $t_{n-2}$ et $t_{n+1}$ nous obtenons
$$
y_{n+1}-y_{n-2}=\int_{t_{n-2}}^{t_{n+1}} \varphi(t,y(t))dt.
$$
Écrire le schéma obtenu en approchant l'intégrale $\int_{t_{n-2}}^{t_{n+1}} \varphi(t,y(t))\dt$ par l'intégrale $\int_{t_{n-2}}^{t_{n+1}} p(t) \dt$ où $p$ est le polynôme interpolant $\varphi$ en $t_{n}$ et $t_{n-1}$ (attention à bien initialiser la suite définie par récurrence pour qu'on puisse effectivement calculer tous les termes). Le schéma obtenu est-il explicite?















Correction
Le polynôme interpolant $\varphi$ en $t_{n}$ et $t_{n-1}$ a équation
$$
p(t) 
= \frac{\varphi(t_n,y_n)-\varphi(t_{n-1},y_{n-1})}{t_n-t_{n-1}}(t-t_n)+ \varphi(t_n,y_n).
$$
On intègre ce polynôme entre $t_{n-2}$ et $t_{n+1}$:
\begin{align*}
\int_{t_{n-2}}^{t_{n+1}} p(t) \dt
&= \frac{\varphi(t_n,y_n)-\varphi(t_{n-1},y_{n-1})}{h}\left[\frac{(t-t_n)^2}{2}\right]_{t_{n-2}}^{t_{n+1}}+ \varphi(t_n,y_n)\left[t\right]_{t_{n-2}}^{t_{n+1}}
\\
&=\frac{\varphi(t_n,y_n)-\varphi(t_{n-1},y_{n-1})}{2h}\left[ (t_{n+1}-t_n)^2-(t_{n-2}-t_n)^2 \right]+ \varphi(t_n,y_n)\left[t_{n+1}-t_{n-2}\right]
\\
&=\frac{3}{2}h\left(\varphi(t_n,y_n)+\varphi(t_{n-1},y_{n-1})\right).
\end{align*}
On obtient le schéma explicite
\begin{cases}
u_0=y_0,\\
u_1=u_0+h\varphi(t_0,u_0), \\
u_{n+1}=u_{n-2}+\frac{3}{2}h\left(\varphi(t_n,u_n)+\varphi(t_{n-1},u_{n-1})\right).
\end{cases}















Exercice : Déduction d'un schéma Predictor-Corrector
Écrire le schéma predictor-corrector basé sur AM-3 AB-2. Attention à bien initialiser la suite récurrente.















Correction
    

  • AB-1
$u_{n+1}=u_n+h\varphi(t_n,u_n)$
    • 

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

  • AB-3
$u_{n+1}=u_n+\frac{h}{12}\left(23\varphi(t_n,u_n)-16\varphi(t_{n-1},u_{n-1})+5\varphi(t_{n-2},u_{n-2})\right)$
    • 

  • AM-3
$u_{n+1}=u_n+\frac{h}{24}\left(9\varphi(t_{n+1},u_{n+1})+19\varphi(t_n,u_n)-5\varphi(t_{n-1},u_{n-1})+\varphi(t_{n-2},u_{n-2}\right)$
donc
\begin{cases}
u_0=y_0,\\
u_1=u_0+h\varphi(t_0,u_0)\\
u_{2}=u_1+\frac{h}{2}\left(3\varphi(t_1,u_1)-\varphi(t_0,u_0)\right)\\
\tilde u=u_n+\frac{h}{12}\left(23\varphi(t_n,u_n)-16\varphi(t_{n-1},u_{n-1})+5\varphi(t_{n-2},u_{n-2})\right)\\
u_{n+1}=u_n+\frac{h}{24}\left(9\varphi(t_{n+1},\tilde u)+19\varphi(t_n,u_n)-5\varphi(t_{n-1},u_{n-1})+\varphi(t_{n-2},u_{n-2}\right)
\end{cases} 
    • 
















Exercice : Ordre d'un schéma
Une méthode numérique à un pas a été utilisée pour résoudre une équation différentielle avec condition initiale. 
Les résultats obtenus par cette méthode en prenant des pas de temps $h = 0.1$, $h = 0.05$ et $h = 0.025$ sont donnés dans le tableau suivant 
(remarque: une valeur sur deux est affichée pour $h = 0.05$ et une valeur sur quatre est affichée pour $h = 0.025$)
$$\begin{array}{|c|c|c|c|}
\hline
t_i &y_i\text{ pour }h=0.1  &y_i\text{ pour }h=0.05 &y_i\text{ pour }h=0.025\\
\hline
1.0 &0.500000           &0.500000 &0.500000\\
1.1 &0.512084           &0.512242 &0.512280\\
1.2 &0.511698           &0.512101 &0.512196\\
1.3 &0.500927           &0.501559 &0.501704\\
1.4 &0.482686           &0.483447 &0.483619\\
1.5 &0.459861           &0.460633 &0.460804\\
\hline
\end{array}
$$
En calculant le rapport des erreurs et sachant que $y(1.5)=0.460857$, déterminer l’ordre de la méthode numérique utilisée.















Correction
La méthode utilisée est d’ordre $2$, en effet, pour $t=1.5$ nous avons
$$


\frac{E(h = 0. 05)}{E(h = 0. 025)}
\frac{|0.460633 - 0.460857|}{|0.460804 - 0.460857|}
= 4.226
\simeq
2^2. 
$$
Si on a accès à polyfit on peut tracer la courbe d'erreur en echèlle logaritmique et estimer la pente:













In [25]:



    


%reset -f
%matplotlib inline
%autosave 300

from matplotlib.pylab import *

H=[0.1,0.05,0.025]
err=[]
y=0.460857
err.append(abs(y-0.459861))
err.append(abs(y-0.460633))
err.append(abs(y-0.460804))

print ('Ordre\t %1.2f' %(polyfit(log(H),log(err), 1)[0]))

subplot(1,2,1)
loglog(H,err, 'r-o');
grid()
subplot(1,2,2)
plot(log(H),log(err), 'r-o');



    















    




















    






Autosaving every 300 seconds
Ordre	 2.12