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

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

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