Schémas “classiques” à un pas

Auteur·rice

Gloria Faccanoni

Date de publication

14 février 2026

Code

import matplotlib.pyplot as plt
import numpy as np
# plt.rcdefaults()
plt.rcParams.update({'font.size': 14})
from scipy.optimize import fsolve

Considérons le problème de Cauchy

$\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’existence et l’unicité d’une solution $y$ sur $I$ aient été établies.

Soit $N\ge 1$ un entier et posons le pas de discrétisation $h=\frac{T-t_0}{N}>0$. On définit alors les $N+1$ nœuds $t_n=t_0+nh$, pour $n=0,1,\dots , N$ qui divisent l’intervalle $I$ en $N$ sous-intervalles $I_n=[t_n,t_{n+1}]$ de même longueur $h$.

Pour chaque nœud $t_n$, on cherche une valeur $u_n$ qui approchera 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 \}$ constitue les points de la discrétisation.
L’ensemble de $N+1$ valeurs $\{y_0, y_1,\dots , y_{N} \}$ constitue la solution exacte discrète.
L’ensemble de $N+1$ valeurs $\{u_0 = y_0, u_1,\dots , u_{N} \}$ constitue la solution numérique obtenue via une suite récurrente.

Les schémas que l’on va construire permettent de calculer (explicitement ou implicitement) $u_{n+1}$ à partir de $u_n, u_{n-1}, ..., u_{n-\kappa}$. Il est ainsi possible de calculer successivement $u_{\kappa+1}, u_{\kappa+2}, \dots, u_N$ à partir des valeurs initiales $u_0, \dots, u_\kappa$ via une formule de récurrence de la forme

$ \begin{cases} u_0 = y_0,\\ u_1 \approx y_1,\\ \vdots\\ u_\kappa \approx y_\kappa,\\ u_{n+1}=\Phi(u_{n+1},u_n, u_{n-1}, \dots, u_{n-\kappa}),&\text{pour } n=\kappa,\kappa+1,\dots,N-1. \end{cases} $

Méthodes explicites vs implicites

Une méthode est dite explicite si la valeur $u_{n+1}$ peut être calculée directement à partir des valeurs précédentes $u_k$, $k\le n$ (ou d’une partie d’entre elles).
Une méthode est dite implicite s’il faut résoudre une équation impliquant $u_{n+1}$ et la fonction $\varphi$ pour obtenir la valeur $u_{n+1}$.

Méthodes à un pas vs multi-pas

Une méthode est dite à un pas si $u_{n+1}$ ne dépend que de $u_n$ (et éventuellement de $u_{n+1}$ lui-même).
Sinon, le schéma est une méthode multi-pas, ou à pas multiples, si $u_{n+1}$ dépend de plusieurs valeurs précédentes.

Construction de schémas à un pas

Si nous intégrons l’EDO $y'(t)=\varphi(t,y(t))$ entre $t_n$ et $t_{n+1}$ nous obtenons

$ \int_{t_n}^{t_{n+1}}y'(t)\mathrm{d}t=\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t $

c’est-à-dire

$ y_{n+1}-y_n=\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t. $

On peut construire différents schémas selon la formule d’approximation utilisée pour approcher le membre de droite : $\begin{cases} u_0=y_0,\\ u_{n+1}=u_n+\displaystyle \text{approximation de }\left(\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t\right) \end{cases} $

Pour construire cette approximation, plusieurs stratégies sont possibles. La plus simple est d’utiliser une formule de quadrature numérique. Pour cela, on remplace la fonction $\varphi(t,y(t))$ par un polynôme interpolant $\varphi(t,y(t))$ en certains points de l’intervalle $[t_n,t_{n+1}]$ (on verra plus tard qu’on peut même construire le polynôme en utilisant des points en dehors de cet intervalle). On obtient ainsi le schéma général suivant :

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

Rappels sur les polynômes interpolants et les formules de quadrature

Schéma d’Euler explicite

Si on remplace une fonction $f$ par une constante égale à la valeur de $f$ dans la borne gauche de l’intervalle $[a;b]$
(polynôme qui interpole $f$ au point $(a,f(a))$ et donc de degré $0$), on a

$ \begin{aligned} \tilde f(x)&=f(a)\\ \int_a^b f(x)\mathrm{d}x &\approx \int_a^b \tilde f(x)\mathrm{d}x=(b-a)f(a). \end{aligned} $

Cette formule est dite formule de quadrature du rectangle à gauche.

En utilisant cette formule pour approcher la fonction $t\mapsto \varphi(t,y(t))$ on a

$ \int_{t_n}^{t_{n+1}} \varphi(t,y(t)) \mathrm{d}t \approx h \varphi(t_n,y(t_n)) $

et on reconnait le schéma d’Euler progressif (= Euler explicite = Euler en avant)

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

Il s’agit d’un schéma à 1 pas explicite, car il permet d’expliciter $u_{n+1}$ directement en fonction de $u_n$.

Code

def EE(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = np.empty_like(tt)
    N = len(tt)-1
    uu[0] = y0
    for n in range(N):
        uu[n+1] = uu[n]+h*phi(tt[n],uu[n])
    return uu

Schéma d’Euler implicite

Si on remplace une fonction $f$ par une constante égale à la valeur de $f$ en la borne droite de l’intervalle $[a;b]$ (polynôme qui interpole $f$ au point $(b,f(b))$ et donc de degré $0$), on a

$ \begin{aligned} \tilde f(x)&=f(b)\\ \int_a^b f(x)\mathrm{d}x &\approx \int_a^b \tilde f(x)\mathrm{d}x=(b-a)f(b). \end{aligned} $

Cette formule est dite formule de quadrature du rectangle à droite.

En utilisant cette formule pour approcher la fonction $t\mapsto \varphi(t,y(t))$ on a

$ \int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t\approx (t_{n+1}-t_n) \varphi(t_{n+1},y(t_{n+1})) = h \varphi(t_{n+1},y(t_{n+1})) $

et on obtient le schéma d’Euler rétrograde (=Euler implicite = Euler en arrière)

$ \text{(EI)}\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,2,\dots N-1 \end{cases} $

Il s’agit d’un schéma - à 1 pas (car il permet de calculer $u_{n+1}$ en fonction de $u_n$) - implicite car, pour calculer $u_{n+1}$, il faudra utiliser un schéma pour le calcul du zéro d’une fonction quelconque. En effet, $u_{n+1}$ est solution de l’équation $x=u_n+h\varphi(t_{n+1},x)$, c’est-à-dire un zéro de la fonction (en général non linéaire) $x\mapsto -x+u_n+h\varphi(t_{n+1},x).$

Code

def EI(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = np.empty_like(tt)
    N = len(tt)-1
    uu[0] = y0
    for n in range(N):
        fct_implicite = lambda x: -x+uu[n]+h*phi(tt[n+1],x)
        tmp = fsolve(fct_implicite, uu[n])
        uu[n+1] = tmp[0] # on récupère la valeur scalaire
    return uu

Schémas d’Euler modifié (EM) ou de type Runge-Kutta (RK1_M)

Si on remplace une fonction $f$ par une constante égale à la valeur de $f$ au milieu de l’intervalle $[a;b]$ (polynôme qui interpole $f$ au point $\left(\tfrac{a+b}{2},f\left(\tfrac{a+b}{2}\right)\right)$ et donc de degré $0$), on a

$ \begin{aligned} \tilde f(x)&=f\left(\tfrac{a+b}{2}\right)\\ \int_a^b f(x)\mathrm{d}x &\approx \int_a^b \tilde f(x)\mathrm{d}x=(b-a)f\left(\tfrac{a+b}{2}\right). \end{aligned} $

Cette formule est dite formule de quadrature du rectangle ou du point milieu.

En utilisant cette formule pour approcher la fonction $t\mapsto \varphi(t,y(t))$ on a

$ \int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t\approx h \varphi\left(t_n+\frac{h}{2},y\left(t_n+\frac{h}{2}\right)\right) $

et on obtient

$\begin{cases} u_0=y(t_0)=y_0,\\ u_{n+1}=u_n+h \varphi\left(t_n+\frac{h}{2},u_{n+1/2}\right)& n=0,1,2,\dots N-1 \end{cases}$

où $u_{n+1/2}$ est une approximation de $y(t_n+h/2)$.

Cependant, $u_{n+1/2}$ est une valeur inconnue, on doit alors en calculer une approximation $\tilde u_{n+1/2}$. Pour cela nous pouvons utiliser deux stratégies différentes :

Utiliser une prédiction d’Euler progressive sur un demi-pas: $\tilde u_{n+1/2}=u_n+(h/2) \varphi(t_{n},u_{n})$.

On peut ainsi remplacer $\varphi(t_{n}+h/2,u_{n+1/2})$ par $\varphi(t_{n}+h/2,\tilde u_{n+1/2})$.

Nous avons construit ainsi un nouveau schéma appelé schéma d’Euler modifié qui s’écrit comme suit et qui est à 1 pas explicite car il permet d’expliciter $u_{n+1}$ en fonction de $u_n$. Cette technique d’utiliser un schéma pour prédire une valeur de la solution en un point qui n’est pas un nœud de la discrétisation est très utilisée et constitue l’une des bases des méthodes de type Runge-Kutta. Ici on a un schéma de type Runge-Kutta à un étage.

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

Utiliser une moyenne: $\tilde u_{n+1/2}=\dfrac{u_n+u_{n+1}}{2}$.

On peut ainsi remplacer $\varphi(t_{n}+h/2,u_{n+1/2})$ par $\varphi(t_{n}+h/2,\tilde u_{n+1/2})$.

Nous avons construit ainsi un nouveau schéma qu’on appelera schéma RK1_M qui s’écrit comme suit et qui est à 1 pas implicite car il ne permet pas d’expliciter $u_{n+1}$ en fonction de $u_n$. Pour calculer $u_{n+1}$ il faudra utiliser un schéma pour le calcul du zéro d’une fonction quelconque. En effet, $u_{n+1}$ est solution de l’équation $x=u_n+h\varphi(t_{n}+\frac{h}{2},\frac{u_n+x}{2})$, c’est-à-dire un zéro de la fonction (en général non linéaire) $x\mapsto -x+u_n+h\varphi\left(t_{n}+\frac{h}{2},\frac{u_n+x}{2}\right).$ À nouveau, il s’agit d’un schéma de type Runge-Kutta à un étage.

$\text{(RK1\_M)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ u_{n+1}=u_n+h \varphi\left(\dfrac{t_n+t_{n+1}}{2},\dfrac{u_n+u_{n+1}}{2}\right)& n=0,1,2,\dots N-1 \end{cases} $

Code

def EM(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = np.empty_like(tt)
    N = len(tt)-1
    uu[0] = y0
    for n in range(N):
        u_tilde = uu[n]+0.5*h*phi(tt[n],uu[n])
        uu[n+1] = uu[n]+h*phi(tt[n]+h/2,u_tilde)
    return uu

def RK1_M(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = np.empty_like(tt)
    N = len(tt)-1
    uu[0] = y0
    for n in range(N):
        fct_implicite = lambda x: -x+uu[n]+h*phi( (tt[n]+tt[n+1])/2,(uu[n]+x)/2 )
        tmp = fsolve(fct_implicite, uu[n])
        uu[n+1] = tmp[0]
    return uu

Schéma du trapèze (ou de Crank-Nicolson)

John Crank (1916-2006) Phyllis Nicolson (1917-1968)

Si on remplace une fonction $f$ par le segment qui relie $(a,f(a))$ à $(b,f(b))$ (polynôme qui interpole $f$ en les points $(a,f(a))$ et $(b,f(b))$ et donc de degré $1$), on a

$ \begin{aligned} \tilde f(x)&=\dfrac{f(b)-f(a)}{b-a}(x-a)+f(a)\\ \int_a^b f(x)\mathrm{d}x &\approx \int_a^b \tilde f(x)\mathrm{d}x=\frac{b-a}{2}\left(f(a)+f(b)\right). \end{aligned} $

Cette formule est dite formule de quadrature du trapèze.

En utilisant cette formule pour approcher la fonction $t\mapsto \varphi(t,y(t))$ on a

\( {t_n}^{t</em>{n+1}} (t,y(t))t((t_{n},y(t_{n}))+(t_{n+1},y(t_{n+1})))

\) et on obtient le schéma du trapèze ou de Crank-Nicolson.

$ \text{(CN)}\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,2,\dots N-1 \end{cases} $

En fait, ce schéma fait la moyenne des schémas d’Euler progressif et rétrograde.

Il s’agit à nouveau d’un schéma à 1 pas implicite, car il ne permet pas d’expliciter directement $u_{n+1}$ en fonction de $u_n$ lorsque la fonction $\varphi$ n’est pas triviale. Pour calculer $u_{n+1}$ il faudra utiliser un schéma pour le calcul du zéro d’une fonction quelconque. En effet, $u_{n+1}$ est solution de l’équation $x=u_n+\frac{h}{2}\left(\varphi(t_n,u_n)+\varphi(t_{n+1},x)\right)$, c’est-à-dire un zéro de la fonction (en général non linéaire) $x\mapsto -x+u_n+\frac{h}{2}\Big(\varphi(t_n,u_n)+\varphi(t_{n+1},x)\Big).$

Code

def CN(phi,tt,y0): h = tt[1]-tt[0] uu = np.empty_like(tt) N = len(tt)-1 uu[0] = y0 for n in range(N): fct_implicite = lambda x: -x+uu[n]+0.5*h*( phi(tt[n],uu[n])+phi(tt[n+1],x) ) tmp = fsolve(fct_implicite, uu[n]) uu[n+1] = tmp[0] return uu

Schéma de Heun

Karl Heun (1859-1929)

Pour éviter le calcul implicite de $u_{n+1}$ dans le schéma du trapèze, nous pouvons utiliser une prédiction d’Euler progressive et remplacer le $u_{n+1}$ dans le terme $\varphi(t_{n+1},u_{n+1})$ par $\tilde u_{n+1}=u_n+h \varphi(t_{n},u_{n})$. Nous avons construit ainsi un nouveau schéma appelé schéma de Heun. Plus précisément, la méthode de Heun s’écrit

$ \text{(H)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ \tilde u_{n+1}=u_n+h \varphi(t_{n},u_{n}),\\ u_{n+1}=u_n+\dfrac{h}{2} \Big(\varphi(t_{n},u_{n}) + \varphi(t_{n+1},\tilde u_{n+1})\Big) & n=0,1,2,\dots N-1 \end{cases} $

On a ainsi construit un schéma à 1 pas explicite. Cette technique de prédiction-correction est très utilisée pour transformer un schéma implicite en schéma explicite.

Code

def heun(phi,tt,y0): h = tt[1]-tt[0] uu = np.empty_like(tt) N = len(tt)-1 uu[0] = y0 for n in range(N): k0 = phi( tt[n], uu[n] ) k1 = phi( tt[n+1], uu[n] + h*k0 ) uu[n+1] = uu[n] + 0.5*h*(k0+k1) return uu

Remarques

À première vue, il semble que les schémas d’Euler progressif et de Heun soient préférable aux schémas d’Euler rétrograde et de Crank-Nicolson puisque ces derniers ne sont pas explicites. Cependant, on verra que les méthodes d’Euler implicite et de Crank-Nicolson sont inconditionnellement A-stable (autrement dit, on peut choisir n’importe quel pas $h>0$, contrairement aux méthodes explicites où le pas doit être suffisamment petit pour garantir la stabilité). C’est aussi le cas de nombreuses autres méthodes implicites. Cette propriété rend les méthodes implicites attractives, bien qu’elles soient plus coûteuses que les méthodes explicites.

Pour la mise en application d’un schéma, il faut aussi prendre en compte l’influence des erreurs d’arrondi. En effet, afin de minimiser l’erreur globale théorique, on pourrait être tenté d’appliquer une méthode avec un pas très petit, par exemple de l’ordre de $10^{-16}$, mais ce faisant, outre que le temps de calcul deviendrait irréaliste, très rapidement les erreurs d’arrondi feraient diverger la solution approchée. En pratique, il faut prendre $h$ assez petit pour que la méthode converge, mais pas trop petit non plus pour que les erreurs d’arrondi ne donnent pas lieu à des résultats incohérents et pour que les calculs puissent être effectués en un temps raisonnable.

Annexes

Schémas basés sur la formule de quadrature de Simpson

Thomas Simpson (1710-1761)

Si on remplace une fonction $f$ par la parabole qui passe par $(a,f(a))$, $(b,f(b))$ et $\left(\tfrac{a+b}{2},f\left(\tfrac{a+b}{2}\right)\right)$ (polynôme qui interpole $f$ en les points $(a,f(a))$, $(b,f(b))$ et $\left(\tfrac{a+b}{2},f\left(\tfrac{a+b}{2}\right)\right)$ et donc de degré $2$), dans la base de Lagrange on a

$ \begin{aligned} \tilde f(x)&=f(a)\frac{(x-b)\left(x-\tfrac{a+b}{2}\right)}{(a-b)\left(a-\tfrac{a+b}{2}\right)}\\ & +f\left(\tfrac{a+b}{2}\right)\frac{(x-a)(x-b)}{\left(\tfrac{a+b}{2}-a\right)\left(\tfrac{a+b}{2}-b\right)} \\ & +f(b)\frac{(x-a)\left(x-\tfrac{a+b}{2}\right)}{(b-a)\left(b-\tfrac{a+b}{2}\right)}\\ \int_a^b f(x)\mathrm{d}x &\approx \int_a^b \tilde f(x)\mathrm{d}x = \frac{b-a}{6} \left( f(a)+4f\left(\tfrac{a+b}{2}\right)+ f(b)\right). \end{aligned} $

Cette formule est dite formule de Simpson.

Simpson’s Rule

Schéma de Simpson implicite

En utilisant cette formule pour approcher la fonction $t\mapsto \varphi(t,y(t))$ on a

$ \int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t\approx \frac{h}{6} \left( \varphi(t_n,y(t_n))+4\varphi\left(t_n+\frac{h}{2},y\left(t_n+\frac{h}{2}\right)\right)+ \varphi(t_{n+1},y(t_{n+1})) \right) $

et on obtient

$\begin{cases} u_0=y(t_0)=y_0,\\ u_{n+1}=u_n+\frac{h}{6} \left(\varphi(t_n,u_n)+4\varphi\left(t_n+\frac{h}{2},u_{n+1/2}\right)+\varphi(t_{n+1},u_{n+1})\right)& n=0,1,2,\dots N-1 \end{cases}$

où $u_{n+1/2}$ est une approximation de $y(t_n+h/2)$.

Puisque $u_{n+1/2}$ est une valeur inconnue, on doit en calculer une approximation $\tilde u_{n+1/2}$. Pour cela nous pouvons utiliser soit la moyenne $\tilde u_{n+1/2}=\dfrac{u_n+u_{n+1}}{2}$ (comme dans le schéma RK1_M), soit une prédiction d’Euler progressive sur un demi-pas $\tilde u_{n+1/2}=u_n+(h/2) \varphi(t_{n},u_{n})$ (comme dans le schéma EM).

Si par exemple on utilise la moyenne, on obtient le schéma suivant

$ \text{(SI-moy)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ u_{n+1}=u_n+\frac{h}{6} \left(\varphi(t_n,u_n)+4\varphi\left(t_n+\frac{h}{2},\dfrac{u_n+u_{n+1}}{2}\right)+\varphi(t_{n+1},u_{n+1})\right)& n=0,1,2,\dots N-1 \end{cases} $

Il s’agit d’un schéma à 1 pas implicite.
Si par exemple on utilise la prédiction d’Euler progressive sur un demi-pas, on obtient le schéma suivant, appelé schéma de Simpson implicite, qui s’écrit

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

Il s’agit d’un schéma à 1 pas implicite.

Schéma de Simpson explicite

Pour éviter le calcul implicite de $u_{n+1}$ dans le schéma de SI, nous pouvons utiliser une prédiction d’Euler progressive et remplacer le $u_{n+1}$ dans le terme $\varphi(t_{n+1},u_{n+1})$ par $\check u_{n+1}=u_n+h \varphi(t_{n},u_{n})$. Nous avons construit ainsi un nouveau schéma qu’on appellera schéma de Simpson explicite RK et qui s’écrit

$ \text{(SE-RK)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ \tilde u_{n+1/2}=u_n+\frac{h}{2}\varphi(t_{n},u_{n}),\\ \check u_{n+1}=u_n+h \varphi(t_{n},u_{n}),\\ u_{n+1}=u_n+\frac{h}{6} \left(\varphi(t_n,u_n)+4\varphi\left(t_n+\frac{h}{2},\tilde u_{n+1/2}\right)+\varphi(t_{n+1},\check u_{n+1})\right)& n=0,1,2,\dots N-1 \end{cases} $

Il s’agit à nouveau d’un schéma à 1 pas explicite de type Runge-Kutta à deux étages.

Schéma de Simpson explicite-variante

Notons qu’on aurait pu remplacer le $u_{n+1}$ dans le terme $\varphi(t_{n+1},u_{n+1})$ par une approximation utilisant $\tilde u_{n+1/2}$ comme une prédiction d’Euler progressive à partir de $t_n+h/2$, ce qui donne $\hat u_{n+1}=\tilde u_{n+1/2}+\frac{h}{2}\varphi(t_{n}+h/2,\tilde u_{n+1/2})$ :

$ \text{(SEV-RK)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ \tilde u_{n+1/2}=u_n+\frac{h}{2}\varphi(t_{n},u_{n}),\\ \hat u_{n+1}=\tilde u_{n+1/2}+\frac{h}{2}\varphi(t_{n}+\frac{h}{2},\tilde u_{n+1/2}),\\ u_{n+1}=u_n+\frac{h}{6} \left(\varphi(t_n,u_n)+4\varphi\left(t_n+\frac{h}{2},\tilde u_{n+1/2}\right)+\varphi(t_{n+1},\hat u_{n+1})\right)& n=0,1,2,\dots N-1 \end{cases} $

Schémas de Simpson multipas

Au lieu de travailler sur $[t_n,t_{n+1}]$ on peut répéter la même construction sur $[t_{n-1},t_{n+1}]$ (et donc intégrer entre $t_{n-1}$ et $t_{n+1}$). Il faut alors initialiser le schéma avec $u_0=y_0$ et $u_1$ calculé par exemple avec le schéma d’Euler explicite.

$ \text{(SMI à deux pas)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ u_{1}=u_0+h\varphi(t_{0},u_{0}),\\ 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,3,\dots N-1 \end{cases} $

On peut ensuite éviter le calcul implicite de $u_{n+1}$ dans le schéma de SMI en utilisant une prédiction d’Euler progressive (dans l’esprit du schéma de Heun) et remplacer le $u_{n+1}$ dans le terme $\varphi(t_{n+1},u_{n+1})$ par $\tilde u_{n+1}=u_n+h \varphi(t_{n},u_{n})$.:

$ \text{(SME à deux pas)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ u_{1}=u_0+h\varphi(t_{0},u_{0}),\\ \tilde u_{n+1}=u_{n}+h\varphi(t_{n},u_{n}),\\ 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},\tilde u_{n+1})\right)& n=1,2,3,\dots N-1 \end{cases} $

Annexe : automatisation de la construction de schémas

Pour vérifier nos calculs d’interpolation puis intégration, nous pouvons utiliser le package de calcul formel SymPy.

On notera

t_n : $t_n$

t_nm1 : $t_{n-1} = t_n-h$

t_np1 : $t_{n+1} = t_n+h$

t_np12 : $t_{n+1/2} = t_n+\frac{h}{2}$

phi_n : $\varphi_n = \varphi(t_n,y(t_n))$

phi_nm1 : $\varphi_{n-1} = \varphi(t_{n-1},y(t_{n-1}))$

phi_np1 : $\varphi_{n+1} = \varphi(t_{n+1},y(t_{n+1}))$

phi_np12 : $\varphi_{n+1/2} = \varphi(t_{n+1/2},y(t_{n+1/2}))$

Code

import sympy as sym from IPython.display import display, Math, Markdown sym.init_printing() phi_n = sym.Symbol(r'\varphi_{n}') phi_np12 = sym.Symbol(r'\varphi_{n+\frac{1}{2}}') phi_np1 = sym.Symbol(r'\varphi_{n+1}') phi_nm1 = sym.Symbol(r'\varphi_{n-1}') sym.var('h,t,t_n') t_np1 = t_n+h t_np12 = t_n+h/2 t_nm1 = t_n-h u_np1=sym.Symbol('u_{n+1}') u_n =sym.Symbol('u_{n}') u_nm1=sym.Symbol('u_{n-1}') display(Markdown("- **EE**:")) p = sym.interpolate([(t_n,phi_n)], t) # display(Math(r't\mapsto \tilde f(t)='+sym.latex(p))) EE = sym.integrate(p,(t,t_n,t_np1)).simplify() # display(Math(r'\int_{t_n}^{t_{n+1}}\tilde f(t) \mathrm{d}t='+sym.latex(EE))) display(sym.Eq(u_np1, u_n + EE)) display(Markdown("- **EI**:")) p = sym.interpolate([(t_np1,phi_np1)], t) # display(Math(r't\mapsto \tilde f(t)='+sym.latex(p))) EI = sym.integrate(p,(t,t_n,t_np1)).simplify() # display(Math(r'\int_{t_n}^{t_{n+1}}\tilde f(t) \mathrm{d}t='+sym.latex(EI))) display(sym.Eq(u_np1, u_n + EI)) display(Markdown("- **EM-init**:")) p = sym.interpolate([(t_np12,phi_np12)], t) # display(Math(r't\mapsto \tilde f(t)='+sym.latex(p))) PM = sym.integrate(p,(t,t_n,t_np1)).simplify() # display(Math(r'\int_{t_n}^{t_{n+1}}\tilde f(t) \mathrm{d}t='+sym.latex(PM))) display(sym.Eq(u_np1, u_n + PM)) display(Markdown("- **CN**:")) p = sym.interpolate([(t_np1,phi_np1),(t_n,phi_n)], t) # display(Math(r't\mapsto \tilde f(t)='+sym.latex(p))) CN = sym.integrate(p,(t,t_n,t_np1)).simplify() # display(Math(r'\int_{t_n}^{t_{n+1}}\tilde f(t) \mathrm{d}t='+sym.latex(CN))) display(sym.Eq(u_np1, u_n + CN)) display(Markdown("- **SI-init**:")) p = sym.interpolate([(t_np1,phi_np1),(t_np12,phi_np12),(t_n,phi_n)], t) # display(Math(r't\mapsto \tilde f(t)='+sym.latex(p))) SI = sym.integrate(p,(t,t_n,t_np1)).simplify() # display(Math(r'\int_{t_n}^{t_{n+1}}\tilde f(t) \mathrm{d}t='+sym.latex(SI))) display(sym.Eq(u_np1, u_n + SI)) display(Markdown("- **SMI à deux pas**:")) p = sym.interpolate([(t_np1,phi_np1),(t_n,phi_n),(t_nm1,phi_nm1)], t) # display(Math(r't\mapsto \tilde f(t)='+sym.latex(p))) SMI = sym.integrate(p,(t,t_nm1,t_np1)).simplify() # display(Math(r'\int_{t_{n-1}}^{t_{n+1}}\tilde f(t) \mathrm{d}t='+sym.latex(SMI))) display(sym.Eq(u_np1, u_nm1 + SMI))

EE:

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

EI:

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

EM-init:

$\displaystyle u_{n+1} = \varphi_{n+\frac{1}{2}} h + u_{n}$

CN:

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

SI-init:

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

SMI à deux pas:

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