TP 4 - Étude de la convergence

Auteur·rice

Gloria Faccanoni

Date de publication

11 mars 2026

Code

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
plt.rcParams.update({'font.size': 12})
from IPython.display import display, Markdown

Exercice (à compléter) : Estimation de l’Ordre de Convergence des Schémas Numériques

🎯 Objectif : compléter l’implémentation des schémas numériques indiqués et estimer l’ordre de convergence de chaque schéma sur le problème de Cauchy

\( \begin{cases} y'(t) = y(t), &\forall t \in I = [0,1], \\ y(0) = 1. \end{cases} \)

dont la solution exacte est \(y(t) = e^{t}\).

🛠️ Résultats attendus :

Méthodes \(\text{AB}_q\) et \(\text{N}_q\) (explicites à \(q\) pas) : convergence à l’ordre \(q\).
Méthodes \(\text{AM}_q\) et \(\text{MS}_q\) (implicites à \(q\) pas) : convergence à l’ordre \(q+1\) (donc sub-optimale lorsque \(q\) est pair). Noter que \(\text{MS}_2 \equiv \text{MS}_3\).
Méthodes \(\text{BDF}_q\) (implicites à \(q\) pas) : convergence à l’ordre \(q\) (donc sub-optimale quelle que soit la parité de \(q\)). Rappelle : elles sont zéro-stables pour \(q \le 6\) et instables pour \(q \ge 7\).

De plus, on vérifiera que les autres méthodes ont les ordres de convergence suivants :

Méthodes predictor-corrector : convergence à l’ordre \(\min(\omega_E + 1, \omega_I)\).
- \(\omega_E\) : ordre du schéma predictor (explicite),
- \(\omega_I\) : ordre du schéma corrector (implicite).
Méthodes de Runge-Kutta à \(s\) étages : convergence à l’ordre …
- \(\le s\) si elles sont explicites,
- \(\le 2s\) si elles sont implicites.

🔎 Remarque :

Puisque la fonction \(\varphi(t, y) = y\) est linéaire, toute méthode implicite peut être rendue explicite par un calcul élémentaire.
Cela consiste à expliciter directement, pour chaque schéma, l’expression de \(u_{n+1}\). Cependant, nous pouvons utiliser le module SciPy pour résoudre ces équations sans modifier l’implémentation des schémas. Attention : l’utilisation de fsolve peut réduire l’ordre de convergence global au même ordre que celui de fsolve lui-même. Dans notre exemple ce n’est pas le cas car la fonction \(\varphi\) est linéaire, mais il est important de le garder à l’esprit.
Initialisation des schémas multi-step : les schémas multi-step nécessitent l’initialisation de plusieurs pas pour pouvoir démarrer. Plutôt que d’utiliser un schéma d’ordre inférieur pour l’initialisation (ce qui dégraderait l’ordre de précision global), nous utiliserons la solution exacte.
Comme d’habitude, on notera \(\varphi_j = \varphi(t_j, u_j)\).
Selon la barrière de Dahlquist, les MPL à \(q\) pas sont au mieux
- d’ordre \(q\) si explicite,
- d’ordre \(q+1\) si implicite et \(q\) impair,
- d’ordre \(q+2\) si implicite et \(q\) pair. C’est pour cela que les méthodes \(\text{AM}_q\) et \(\text{MS}_q\) sont sub-optimales lorsque \(q\) est pair, tandis que les méthodes \(\text{BDF}_q\) sont sub-optimales pour tout \(q\).

Rappels de la démarche pour le calcul de l’orde de convergence.

On écrit les schémas numériques :
- les nœuds d’intégration \([t_0,t_1,\dots,t_{N}]\) sont contenus dans le vecteur tt (qui change en fonction de h)
- les valeurs \([u_0,u_1,\dots,u_{N}]\) pour chaque méthode sont contenues dans le vecteur uu.
Pour chaque schéma, on calcule la solution approchée avec différentes valeurs de \(h_k=1/N_k\). On sauvegarde les valeurs de \(h_k\) dans le vecteur H.
Pour chaque valeur de \(h_k\), on calcule le maximum de la valeur absolue de l’erreur et on sauvegarde toutes ces erreurs dans le vecteur err_schema de sort que err_schema[k] contient \(e_k=\max_{i=0,\dots,N_k}|y(t_i)-u_{i}|\).
Pour afficher l’ordre de convergence on utilise une échelle logarithmique, i.e. on représente \(\ln(h)\) sur l’axe des abscisses et \(\ln(\text{err})\) sur l’axe des ordonnées.
En effet, si \(\text{err}=Ch^p\) alors \(\ln(\text{err})=\ln(C)+p\ln(h)\).
En échelle logarithmique, \(p\) représente donc la pente de la ligne droite \(\ln(\text{err})\).

Nota Bene Le choix de la norme doit être cohérent avec la norme continue que l’on cherche à approcher.

Pour la norme infinie continue \(\|e\|_{L^\infty}=\sup_{t\in[0,1]}|e(t)|\), la norme discrète naturelle est \(\max_i |e_i|\), qui ne dépend pas du nombre de points du maillage.

En revanche, pour approcher la norme \(L^2\) continue \(\|e\|_{L^2}=\left(\int_0^1 e(t)^2\,dt\right)^{1/2}\), la norme euclidienne discrète \(\|e\|_2=\sqrt{\sum_i e_i^2}\) doit être multipliée par \(\sqrt{h}\) (où \(h\) est le pas de temps), car l’intégrale est approchée par une somme de Riemann : \(\int_0^1 e(t)^2 dt \approx h \sum_i e_i^2\). On obtient donc \(\|e\|_{L^2} \approx \sqrt{h}\,\|e\|_2\). Sans ce facteur \(\sqrt{h}\), la norme dépend artificiellement du nombre de points et l’ordre de convergence mesuré est incorrect.

Schémas explicites

Schéma EE

Schéma d’Adam-Bashforth à 1 pas = schéma d’Euler progressif

\( \begin{cases} u_0=y_0,\\ u_{n+1}=u_n+h \varphi_n & n=0,1,2,\dots N-1 \end{cases} \)

Code

def EE(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initialisation
    uu[0] = sol_exacte(tt[0])
    # récurrence
    for n in range(len(tt)-1):
        k1 = phi(tt[n],uu[n])
        uu[n+1] = uu[n]+h*k1
    return uu

Schéma AB-2

Schéma d’Adam-Bashforth à 2 pas

\( \begin{cases} u_0=y_0,\\ u_{1}=y_1,\\ u_{n+1}=u_n+\frac{h}{2}\Bigl(3 \varphi_n - \varphi_{n-1} \Bigr)& n=1,2,3,4,5,\dots N-1 \end{cases} \)

Code

def AB2(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initialisation
    q = 2
    uu[:q] = sol_exacte(tt[:q]) # on initialise uu[0] et uu[1]
    # récurrence
    for n in range(q-1,len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        k2 = phi( tt[n-1], uu[n-1] )
        uu[n+1] = uu[n] + (3*k1-k2)*h/2
    return uu

Schéma AB-3

Schéma d’Adam-Bashforth à 3 pas

\( \begin{cases} u_0=y_0,\\ u_{1}=y_1,\\ u_{2}=y_2,\\ u_{n+1}=u_n+\frac{h}{12}\Bigl(23 \varphi_n -16 \varphi_{n-1} +5 \varphi_{n-2} \Bigr)& n=2,3,4,5,\dots N-1 \end{cases} \)

Code

def AB3(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initialisation
    q = 3
    uu[:q] = sol_exacte(tt[:q])
    # récurrence
    for n in range(q-1,len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        k2 = phi( tt[n-1], uu[n-1] )
        k3 = phi( tt[n-2], uu[n-2] )
        uu[n+1] = uu[n] + (23*k1-16*k2+5*k3)*h/12
    return uu

Schéma AB-4

Schéma d’Adam-Bashforth à 4 pas

\( \begin{cases} u_0=y_0,\\ u_{1}=y_1,\\ u_{2}=y_2,\\ u_{3}=y_3,\\ u_{n+1}=u_n+\frac{h}{24}\Bigl(55 \varphi_n -59 \varphi_{n-1} +37 \varphi_{n-2} -9 \varphi_{n-3} \Bigr)& n=3,4,5,\dots N-1 \end{cases} \)

Code

def AB4(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initialisation
    q = 4
    uu[:q] = sol_exacte(tt[:q])
    # récurrence
    for n in range(q-1,len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        k2 = phi( tt[n-1], uu[n-1] )
        k3 = phi( tt[n-2], uu[n-2] )
        k4 = phi( tt[n-3], uu[n-3] )
        uu[n+1] = uu[n] + (55*k1-59*k2+37*k3-9*k4)*h/24
    return uu

Schéma AB-5

Schéma d’Adam-Bashforth à 5 pas

\( \begin{cases} u_0=y_0,\\ u_{1}=y_1,\\ u_{2}=y_2,\\ u_{3}=y_3,\\ u_{4}=y_4,\\ u_{n+1}=u_n+\frac{h}{720}\Bigl(1901 \varphi_n -2774 \varphi_{n-1} +2616 \varphi_{n-2} -1274 \varphi_{n-3} +251 \varphi_{n-4} \Bigr)& n=4,5,\dots N-1 \end{cases} \)

Code

def AB5(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initialisation
    q = 5
    uu[:q] = sol_exacte(tt[:q])
    # récurrence
    for n in range(q-1,len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        k2 = phi( tt[n-1], uu[n-1] )
        k3 = phi( tt[n-2], uu[n-2] )
        k4 = phi( tt[n-3], uu[n-3] )
        k5 = phi( tt[n-4], uu[n-4] )
        uu[n+1] = uu[n] + (1901*k1-2774*k2+2616*k3-1274*k4+251*k5)*h/720
    return uu

Schéma N-2

Schéma de Nylström à 2 pas

\( \begin{cases} u_0=y_0,\\ u_{1}=y_1,\\ u_{n+1}=u_{n-1}+2h \varphi_{n} & n=1,2,3,4,5,\dots N-1 \end{cases} \)

Code

def N2(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initialisation
    q = 2
    uu[:q] = sol_exacte(tt[:q])
    # récurrence
    for n in range(q-1,len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        uu[n+1] = uu[n-1] + 2*h*k1
    return uu

Schéma N-3

Schéma de Nylström à 3 pas

\( \begin{cases} u_0=y_0,\\ u_{1}=y_1,\\ u_{2}=y_2,\\ u_{n+1}=u_{n-1}+\frac{h}{3}\Bigl(7 \varphi_{n} -2 \varphi_{n-1} + \varphi_{n-2} \Bigr)& n=2,3,4,5,\dots N-1 \end{cases} \)

Code

def N3(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initialisation
    q = 3
    uu[:q] = sol_exacte(tt[:q])
    # récurrence
    for n in range(q-1,len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        k2 = phi( tt[n-1], uu[n-1] )
        k3 = phi( tt[n-2], uu[n-2] )
        uu[n+1] = uu[n-1] + (7*k1-2*k2+k3)*h/3
    return uu

Schéma N-4

Schéma de Nylström à 4 pas

\( \begin{cases} u_0=y_0,\\ u_{1}=y_1,\\ u_{2}=y_2,\\ u_{3}=y_3,\\ u_{n+1}=u_{n-1}+\frac{h}{3}\Bigl(8 \varphi_{n} -5 \varphi_{n-1} +4 \varphi_{n-2} - \varphi_{n-3} \Bigr)& n=3,4,5,\dots N-1 \end{cases} \)

Code

def N4(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initialisation
    q = 4
    uu[:q] = sol_exacte(tt[:q])
    # récurrence
    for n in range(q-1,len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        k2 = phi( tt[n-1], uu[n-1] )
        k3 = phi( tt[n-2], uu[n-2] )
        k4 = phi( tt[n-3], uu[n-3] )
        uu[n+1] = uu[n-1] + (8*k1-5*k2+4*k3-k4)*h/3
    return uu

Schéma EM

Schéma d’Euler modifié (RK explicite à 2 pas)

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

Code

def EM(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initialisation
    uu[0] = sol_exacte(tt[0])
    # récurrence
    for n in range(len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        uu[n+1] = uu[n]+h*phi(tt[n]+h/2,uu[n]+k1*h/2)
    return uu

Schéma RK4-1

Schéma de Runge-Kutta RK4-1

\( \begin{cases} u_0 = y_0 \\ K_1 = \varphi\left(t_n,u_n\right)\\ K_2 = \varphi\left(t_n+\frac{h}{2},u_n+\frac{h}{2} K_1)\right)\\ K_3 = \varphi\left(t_n+\frac{h}{2},u_n+\frac{h}{2}K_2\right)\\ K_4 = \varphi\left(t_{n+1},u_n+h K_3\right)\\ u_{n+1} = u_n + \frac{h}{6}\left(K_1+2K_2+2K_3+K_4\right) & n=0,1,\dots N-1 \end{cases} \)

Code

def RK4(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initialisation
    uu[0] = sol_exacte(tt[0])
    # récurrence
    for n in range(len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        k2 = phi( tt[n]+h/2 , uu[n]+h*k1/2 )
        k3 = phi( tt[n]+h/2 , uu[n]+h*k2/2 )
        k4 = phi( tt[n+1]    , uu[n]+h*k3 )
        uu[n+1] = uu[n] + (k1+2*k2+2*k3+k4)*h/6
    return uu

Schéma RK6_5

Schéma de Runge-Kutta RK6_5 (Butcher à 6 étages d’ordre 5)

\(\begin{array}{c|cccccc} 0 & 0 & 0 & 0 & 0 &0&0 \\ \frac{1}{4} & \frac{1}{4} & 0&0&0&0&0\\ \frac{1}{4} & \frac{1}{8} &\frac{1}{8}&0&0&0&0\\ \frac{1}{2} & 0 &-\frac{1}{2}&1&0&0&0\\ \frac{3}{4} & \frac{3}{16} &0&0&\frac{9}{16}&0&0\\ 1 & \frac{-3}{7} &\frac{2}{7}&\frac{12}{7}&\frac{-12}{7}&\frac{8}{7}&0\\ \hline & \frac{7}{90} & 0&\frac{32}{90} & \frac{12}{90} & \frac{32}{90} & \frac{7}{90} \end{array}\)

Code

def RK6_5(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initialisation
    uu[0] = sol_exacte(tt[0])
    # récurrence
    for n in range(len(tt)-1):
        k1 = phi( tt[n]      , uu[n] )
        k2 = phi( tt[n]+h/4  , uu[n]+h*k1/4 )
        k3 = phi( tt[n]+h/4  , uu[n]+h*(k1+k2)/8 )
        k4 = phi( tt[n]+h/2  , uu[n]+h*(-k2+2*k3)/2 )
        k5 = phi( tt[n]+h*3/4, uu[n]+h*(3*k1+9*k4)/16 )
        k6 = phi( tt[n+1]    , uu[n]+h*(-3*k1+2*k2+12*k3-12*k4+8*k5)/7 )
        uu[n+1] = uu[n] + (7*k1+32*k3+12*k4+32*k5+7*k6)*h/90
    return uu

Schéma RK7_6

Schéma de Runge-Kutta RK7_6 (Butcher à 7 étages d’ordre 6)

\( \begin{array}{c|ccccccc} 0 & 0 & 0 & 0 & 0 &0&0&0 \\ \frac{1}{2} & \frac{1}{2} & 0&0&0&0&0&0\\ \frac{2}{3} & \frac{2}{9} &\frac{4}{9}&0&0&0&0&0\\ \frac{1}{3} & \frac{7}{36} &\frac{2}{9}&-\frac{1}{12}&0&0&0&0\\ \frac{5}{6} & -\frac{35}{144} &-\frac{55}{36}&\frac{35}{48}&\frac{15}{8}&0&0&0\\ \frac{1}{6} & -\frac{1}{360} &-\frac{11}{36}&-\frac{1}{8}&\frac{1}{2}&\frac{1}{10}&0&0\\ 1 & \frac{-41}{260} &\frac{22}{13}&\frac{43}{156}&-\frac{118}{39}&\frac{32}{195}&\frac{80}{39}&0\\ \hline & \frac{13}{200} & 0&\frac{11}{40} & \frac{11}{40} & \frac{4}{25} & \frac{4}{25} & \frac{13}{200} \end{array}\)

Code

def RK7_6(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initialisation
    uu[0] = sol_exacte(tt[0])
    # récurrence
    for n in range(len(tt)-1):
        k1 = phi( tt[n]      , uu[n] )
        k2 = phi( tt[n]+h/2  , uu[n]+h*(        k1                                                     )/2 )
        k3 = phi( tt[n]+h*2/3, uu[n]+h*(      2*k1     +4*k2                                           )/9 )
        k4 = phi( tt[n]+h/3  , uu[n]+h*(      7*k1     +8*k2      -3*k3                                )/36 )
        k5 = phi( tt[n]+h*5/6, uu[n]+h*(    -35*k1   -220*k2    +105*k3    +270*k4                     )/144 )
        k6 = phi( tt[n]+h/6  , uu[n]+h*(       -k1   -110*k2     -45*k3    +180*k4     +36*k5          )/360 )
        k7 = phi( tt[n+1]    , uu[n]+h*(-41/260*k1 +22/13*k2 +43/156*k3 -118/39*k4 +32/195*k5 +80/39*k6)     )
        uu[n+1] = uu[n] + (13/200*k1+11/40*k3+11/40*k4+4/25*k5+4/25*k6+13/200*k7)*h
    return uu

Schémas implicites

Schéma EI

Schéma d’Euler régressif

\( \begin{cases} u_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} \)

avec \(u_{n+1}\) zéro de la fonction

\( x\mapsto -x+u_n+h\varphi(t_{n+1},x). \)

Code

def EI(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initialisation
    uu[0] = sol_exacte(tt[0])
    # récurrence
    for n in range(len(tt)-1):
        temp = fsolve(lambda x: -x+uu[n]+h*phi(tt[n+1],x), uu[n])
        uu[n+1] = temp[0]
    return uu

Schéma CN

Schéma de Crank-Nicolson

\( \begin{cases} u_0=y_0,\\ u_{n+1}=u_n+\frac{h}{2}\Bigl( \varphi_n +\varphi(t_{n+1},u_{n+1})\Bigr)& n=0,1,2,\dots N-1 \end{cases} \)