TP 3 - Implémentation de schémas

Exercice : Implémentation et Comparaison de Schémas

🎯 Objectif : implémenter les méthodes indiquées et les tester progressivement sur le problème de Cauchy :

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

📌 Travail à réaliser :

1. Implémenter les schémas numériques indiqués, classés par catégories : schémas explicites, schémas implicites, schémas predictor-corrector. Comme d’habitude, on notera \(\varphi_j = \varphi(t_j, u_j)\).

2. Calculer la solution exacte du problème de Cauchy en utilisant le module sympy.

3. Pour chaque schéma, calculer la solution approchée avec un pas \(h = \frac{1}{N}\) et \(N = 8\) (une grille grossière permet de mieux visualiser les erreurs).

4. Analyse et comparaison des résultats :

   • Pour chaque schéma, tracer la solution exacte et la solution approchée.
     
   • Calculer et afficher l’erreur maximale définie par : \(e = \max_{i=0,\dots,N} |y(t_i) - u_i|\)

🚀 Démarche recommandée :

• Calculer la solution exacte et tester les deux schémas déjà implémentés (Étape 0).
• Implémenter un seul schéma (Étape 1).
  
• Tester immédiatement ce schéma sur le problème de Cauchy donné (Étape 2).
  
• Répéter cette démarche pour chaque méthode demandée (commencer par une méthode explicite, puis une méthode implicite, puis une méthode predictor-corrector et recommencer).

Les sections où vous devez ajouter du code sont signalées par les commentaires commençant par : # !!!!!!

AvertissementRemarques sur les boucles pour l’implémentation des suites définies par récurrence :

Rappels :

1. En python,
   • L’instruction range(5) produit 0,1,2,3,4
     
   • L’instruction range(1,5) produit 1,2,3,4
2. Selon notre notation,
   • la suite \(\{t_0, t_1,\dots , t_{N} \}\) contient les \(N+1\) points de discrétisation
     
   • la suite \(\{u_0, u_1,\dots , u_{N} \}\) contient les \(N+1\) valeurs à calculer.

\(\leadsto\) len(tt) \(=N+1\).

Itérations pour calculer \(u_N\) :

• pour un schéma à un pas, \(u_{n+1}=F(u_n)\) ; on initialise \(u_0\) et on calcule \(F(u_{n})\) pour \(n\) de \(0\) jusqu’à \(N-1\), autrement dit n in range(N) soit encore n in range(len(tt)-1) ;
• pour un schéma à deux pas, \(u_{n+1}=F(u_n,u_{n-1})\) ; on initialise \(u_0\) et \(u_1\) et on calcule \(F(u_{n},u_{n-1})\) pour \(n\) de \(1\) jusqu’à \(N-1\), autrement dit n in range(1,N) soit encore n in range(1,len(tt)-1);
• pour un schéma à trois pas, \(u_{n+1}=F(u_n,u_{n-1},u_{n-2})\) ; on initialise \(u_0\), \(u_1\) et \(u_2\) et on calcule \(u_{n+1}\) pour \(n\) de \(2\) jusqu’à \(N-1\), autrement dit n in range(2,N) soit encore n in range(2,len(tt)-1)
• \(\dots\)
• pour un schéma à \(q=p+1\) pas, \(u_{n+1}=F(u_n,u_{n-1},u_{n-2},\dots,u_{n-p})\) ; on initialise \(u_0\), \(u_1\) et \(u_{p}\) et on calcule \(u_{n+1}\) pour \(n\) de \(p\) jusqu’à \(N-1\), autrement dit n in range(p,N) soit encore n in range(p,len(tt)-1)

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 12})

from scipy.optimize import fsolve

import sympy as sym
sym.init_printing()

Autosaving every 300 seconds

Schémas explicites

Schéma AB-1

Schéma d’Euler progressif = d’Adam-Bashforth à 1 pas \( \begin{cases} u_0=y_0,\\ u_{n+1}=u_n+h\varphi_{n}& n=0,1,2,\dots N-1 \end{cases} \)

def EE(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial condition
    uu[0] = y0
    # time-stepping loop
    for n in range(len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        uu[n+1] = uu[n]+k1*h
    return uu

Schéma AB-2

Schéma d’Adam-Bashforth à 2 pas \( \begin{cases} u_0=y_0,\\ u_{1}=u_0+h\varphi_{0},\\ u_{n+1}=u_n+\dfrac{h}{2}\Bigl(3\varphi_{n}-\varphi_{n-1}\Bigr)& n=1,2,3,4,5,\dots N-1 \end{cases} \)

def AB2(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial conditions
    uu[0] = y0
    uu[1] = uu[0]+h*phi(tt[0],uu[0])
    # time-stepping loop
    for n in range(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}=u_0+h\varphi_{0},\\ u_{2}=u_1+\dfrac{h}{2}\Bigl(3\varphi_{1}-\varphi_{0}\Bigr),\\ u_{n+1}=u_n+\dfrac{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} \)

def AB3(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial conditions
    uu[0] = y0
    uu[1] = uu[0]+h*phi(tt[0],uu[0])
    uu[2] = uu[1]+h*(3*phi(tt[1],uu[1])-phi(tt[0],uu[0]))/2
    # time-stepping loop
    for n in range(2,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}=u_0+h\varphi_{0},\\ u_{2}=u_1+\dfrac{h}{2}\Bigl(3\varphi_{1}-\varphi_{0}\Bigr),\\ u_{3}=u_2+\dfrac{h}{12}\Bigl(23\varphi_{2}-16\varphi_{1}+5\varphi_{0}\Bigr),\\ u_{n+1}=u_n+\dfrac{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} \)

def AB4(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial conditions
    uu[0] = y0
    uu[1] = uu[0]+h*phi(tt[0],uu[0])
    uu[2] = uu[1]+h*(3*phi(tt[1],uu[1])-phi(tt[0],uu[0]))/2
    uu[3] = uu[2]+h*(23*phi(tt[2],uu[2])-16*phi(tt[1],uu[1])+5*phi(tt[0],uu[0]))/12
    # time-stepping loop
    for n in range(3,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}=u_0+h\varphi_{0},\\ u_{2}=u_1+\dfrac{h}{2}\Bigl(3\varphi_{1}-\varphi_{0}\Bigr),\\ u_{3}=u_2+\dfrac{h}{12}\Bigl(23\varphi_{2}-16\varphi_{1}+5\varphi_{0}\Bigr),\\ u_{4}=u_3+\dfrac{h}{24}\Bigl(55\varphi_{3}-59\varphi_{2}+37\varphi_{1}-9\varphi_{0}\Bigr),\\ u_{n+1}=u_n+\dfrac{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} \)

def AB5(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial conditions
    uu[0] = y0
    uu[1] = uu[0]+h*phi(tt[0],uu[0])
    uu[2] = uu[1]+h*(3*phi(tt[1],uu[1])-phi(tt[0],uu[0]))/2
    uu[3] = uu[2]+h*(23*phi(tt[2],uu[2])-16*phi(tt[1],uu[1])+5*phi(tt[0],uu[0]))/12
    uu[4] = uu[3]+h*(55*phi(tt[3],uu[3])-59*phi(tt[2],uu[2])+37*phi(tt[1],uu[1])-9*phi(tt[0],uu[0]))/24
    # time-stepping loop
    for n in range(4,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}=u_0+h\varphi_{0},\\ u_{n+1}=u_{n-1}+2h\varphi_{n}& n=1,2,3,4,5,\dots N-1 \end{cases} \)

def N2(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial conditions
    uu[0] = y0
    uu[1] = uu[0]+h*phi(tt[0],uu[0])
    # time-stepping loop
    for n in range(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}=u_0+h\varphi_{0},\\ u_{2}=u_{0}+2h\varphi_{1},\\ 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} \)

def N3(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial conditions
    uu[0] = y0
    uu[1] = uu[0]+h*phi(tt[0],uu[0])
    uu[2] = uu[0]+2*h*phi(tt[1],uu[1])
    # time-stepping loop
    for n in range(2,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}=u_0+h\varphi_{0},\\ u_{2}=u_{0}+2h\varphi_{1},\\ u_{3}=u_{1}+\frac{h}{3}\Bigl(7\varphi_{2}-2\varphi_{1}+\varphi_{0}\Bigr),\\ 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} \)

def N4(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial conditions
    uu[0] = y0
    uu[1] = uu[0]+h*phi(tt[0],uu[0])
    uu[2] = uu[0]+2*h*phi(tt[1],uu[1])
    uu[3] = uu[1]+(7*h*phi(tt[2],uu[2])-2*h*phi(tt[1],uu[1])+h*phi(tt[0],uu[0]))/3
    # time-stepping loop
    for n in range(3,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é (qui est un schéma RK à 2 étages d’ordre 2) \( \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} \)

def EM(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial condition
    uu[0] = y0
    # time-stepping loop
    for n in range(len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        k2 = phi( tt[n]+h/2, uu[n]+h/2*k1 )
        uu[n+1] = uu[n]+h*k2
    return uu

Schéma RK-4.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} \)

def RK4(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial condition
    uu[0] = y0
    # time-stepping loop
    for n in range(len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        k2 = phi( tt[n]+h/2 , uu[n]+k1*h/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 RK-6.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}\)

def RK6_5(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initialisation
    uu[0] = y0
    # recurrence
    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 RK-7.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}\)

def RK7_6(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initialisation
    uu[0] = y0
    # recurrence
    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 (AM-0) \( \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)\)

def EI(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial condition
    uu[0] = y0
    # time-stepping loop
    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 = AM-1

Schéma de Crank-Nicolson (AM-1 = BDF-1)

\( \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} \)

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

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

def CN(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial condition
    uu[0] = y0
    # time-stepping loop
    for n in range(len(tt)-1):
        temp = fsolve(lambda x: -x+uu[n]+0.5*h*( phi(tt[n+1],x)+phi(tt[n],uu[n]) ), uu[n])
        uu[n+1] = temp[0]
    return uu

Schéma AM-2

Schéma d’Adam-Moulton à 2 pas

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

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

\( x\mapsto -x+u_n+\frac{h}{12}\Bigl(5\varphi(t_{n+1},x)+8\varphi_{n}-\varphi_{n-1}\Bigr). \)

def AM2(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial condition
    uu[0] = y0
    temp = fsolve(lambda x: -x+uu[0]+0.5*h*( phi(tt[1],x)+phi(tt[0],uu[0]) ), uu[0])
    uu[1] = temp[0]
    # time-stepping loop
    for n in range(1,len(tt)-1):
        temp = fsolve(lambda x: -x+uu[n]+h*( 5*phi(tt[n+1],x)+8*phi(tt[n],uu[n])-phi(tt[n-1],uu[n-1]) )/12, uu[n])
        uu[n+1] = temp[0]
    return uu

Schéma AM-3

Schéma d’Adam-Moulton à 3 pas

\( \begin{cases} u_0=y_0,\\ u_1=u_0+\frac{h}{2}\Bigl(\varphi_{1}+\varphi_{0}\Bigr),\\ u_{2}=u_1+\frac{h}{12}\Bigl(5\varphi_{2}+8\varphi_{1}-\varphi_{0}\Bigr),\\ u_{n+1}=u_n+\frac{h}{24}\Bigl(9\varphi(t_{n+1},u_{n+1})+19\varphi_{n}-5\varphi_{n-1}+\varphi_{n-2}\Bigr)& n=2,3,\dots N-1 \end{cases} \)

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

\( x\mapsto -x+u_n+\frac{h}{24}\Bigl(9\varphi(t_{n+1},x)+19\varphi_{n}-5\varphi_{n-1}+\varphi_{n-2}\Bigr). \)

def AM3(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial condition
    uu[0] = y0
    temp = fsolve(lambda x: -x+uu[0]+0.5*h*( phi(tt[1],x)+phi(tt[0],uu[0]) ), uu[0])
    uu[1] = temp[0]
    temp = fsolve(lambda x: -x+uu[1]+h*( 5*phi(tt[2],x)+8*phi(tt[1],uu[1])-phi(tt[0],uu[0]) )/12, uu[1])
    uu[2] = temp[0]
    # time-stepping loop
    for n in range(2,len(tt)-1):
        temp = fsolve(lambda x: -x+uu[n]+h*( 9*phi(tt[n+1],x)+19*phi(tt[n],uu[n])-5*phi(tt[n-1],uu[n-1])+phi(tt[n-2],uu[n-2]) )/24, uu[n])
        uu[n+1] = temp[0]
    return uu

Schéma AM-4

Schéma d’Adam-Moulton à 4 pas

\( \begin{cases} u_0=y_0,\\ u_1=u_0+\frac{h}{2}\Bigl(\varphi_{1}+\varphi_{0}\Bigr),\\ u_{2}=u_1+\frac{h}{12}\Bigl(5\varphi_{2}+8\varphi_{1}-\varphi_{0}\Bigr),\\ u_{3}=u_2+\frac{h}{24}\Bigl(9\varphi_{3}+19\varphi_{2}-5\varphi_{1}+\varphi_{0}\Bigr),\\ u_{n+1}=u_n+\frac{h}{720}\Bigl(251\varphi(t_{n+1},u_{n+1})+646\varphi_{n}-264\varphi_{n-1}+106\varphi_{n-2}-19\varphi_{n-3}\Bigr)& n=3,4,\dots N-1 \end{cases} \)

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

\( x\mapsto -x+u_n+\frac{h}{720}\Bigl(251\varphi(t_{n+1},x)+646\varphi_{n}-264\varphi_{n-1}+106\varphi_{n-2}-19\varphi_{n-3}\Bigr). \)

def AM4(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    
    # initial condition
    uu[0] = y0
    
    temp = fsolve(lambda x: -x+uu[0]+0.5*h*( phi(tt[1],x)+phi(tt[0],uu[0]) ), uu[0])
    uu[1] = temp[0]
    
    temp = fsolve(lambda x: -x+uu[1]+h*( 5*phi(tt[2],x)+8*phi(tt[1],uu[1])-phi(tt[0],uu[0]) )/12, uu[1])
    uu[2] = temp[0]
    
    temp = fsolve(lambda x: -x+uu[2]+h*( 9*phi(tt[3],x)+19*phi(tt[2],uu[2])-5*phi(tt[1],uu[1])+phi(tt[0],uu[0]) )/24, uu[2])
    uu[3] = temp[0]
    
    # time-stepping loop
    for n in range(3,len(tt)-1):
        temp = fsolve(lambda x: -x+uu[n]+h*( 251*phi(tt[n+1],x)+646*phi(tt[n],uu[n])-264*phi(tt[n-1],uu[n-1])+106*phi(tt[n-2],uu[n-2])-19*phi(tt[n-3],uu[n-3]) )/720, uu[n])
        uu[n+1] = temp[0]
    return uu

Schéma AM-5

Schéma d’Adam-Moulton à 5 pas

\( \begin{cases} u_0=y_0,\\ u_1=u_0+\frac{h}{2}\Bigl(\varphi_{1}+\varphi_{0}\Bigr),\\ u_{2}=u_1+\frac{h}{12}\Bigl(5\varphi_{2}+8\varphi_{1}-\varphi_{0}\Bigr),\\ u_{3}=u_2+\frac{h}{24}\Bigl(9\varphi_{3}+19\varphi_{2}-5\varphi_{1}+\varphi_{0}\Bigr),\\ u_{4}=u_3+\frac{h}{720}\Bigl(251\varphi_{4}+646\varphi_{3}-264\varphi_{2}+106\varphi_{1}-19\varphi_{0}\Bigr),\\ u_{n+1}=u_n+\frac{h}{1440}\Bigl(475\varphi(t_{n+1},u_{+1})+1427\varphi_{n}-798\varphi_{n-1}+482\varphi_{n-2}-173\varphi_{n-3}+27\varphi_{n-4}\Bigr),& n=4,5,\dots N-1 \end{cases} \)

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

\( x\mapsto -x+u_n+\frac{h}{1440}\Bigl(475\varphi(t_{n+1},x)+1427\varphi_{n}-798\varphi_{n-1}+482\varphi_{n-2}-173\varphi_{n-3}+27\varphi_{n-4}\Bigr). \)

def AM5(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    
    # initial condition
    uu[0] = y0
    
    temp = fsolve(lambda x: -x+uu[0]+0.5*h*( phi(tt[1],x)+phi(tt[0],uu[0]) ), uu[0])
    uu[1] = temp[0]
    
    temp = fsolve(lambda x: -x+uu[1]+h*( 5*phi(tt[2],x)+8*phi(tt[1],uu[1])-phi(tt[0],uu[0]) )/12, uu[1])
    uu[2] = temp[0]
    
    temp = fsolve(lambda x: -x+uu[2]+h*( 9*phi(tt[3],x)+19*phi(tt[2],uu[2])-5*phi(tt[1],uu[1])+phi(tt[0],uu[0]) )/24, uu[2])
    uu[3] = temp[0]

    temp = fsolve(lambda x: -x+uu[3]+h*( 251*phi(tt[4],x)+646*phi(tt[3],uu[3])-264*phi(tt[2],uu[2])+106*phi(tt[1],uu[1])-19*phi(tt[0],uu[0]) )/720, uu[3])
    uu[4] = temp[0]
    
    # time-stepping loop
    for n in range(4,len(tt)-1):
        temp = fsolve(lambda x: -x+uu[n]+h*( 475*phi(tt[n+1],x)+1427*phi(tt[n],uu[n])-798*phi(tt[n-1],uu[n-1])+482*phi(tt[n-2],uu[n-2])-173*phi(tt[n-3],uu[n-3])+27*phi(tt[n-4],uu[n-4]) )/1440, uu[n])
        uu[n+1] = temp[0]
    return uu

Schéma BDF-2

Schéma Backward Differentiation Formula à 2 pas

\( \begin{cases} u_0=y_0,\\ u_1=u_0+h\varphi_{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,\dots N-1 \end{cases} \)

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

\( x\mapsto -x+\frac{4}{3}u_n-\frac{1}{3}u_{n-1}+\frac{2}{3}h\varphi(t_{n+1},x). \)

def BDF2(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial condition
    uu[0] = y0
    temp = fsolve(lambda x: -x+uu[0]+h*phi(tt[1],x), uu[0])
    uu[1] = temp[0]
    # time-stepping loop
    for n in range(1,len(tt)-1):
        temp = fsolve(lambda x: -x+4/3*uu[n]-1/3*uu[n-1] + 2/3*h*phi(tt[n+1],x) , uu[n])
        uu[n+1] = temp[0]
    return uu   

Schéma BDF-3

Schéma Backward Differentiation Formula à 3 pas

\( \begin{cases} u_0=y_0,\\ u_1=u_0+h\varphi_{1},\\ u_{2}=\frac{4}{3}u_1-\frac{1}{3}u_{0}+\frac{2}{3}h\varphi_{2},\\ u_{n+1}=\frac{18}{11}u_n-\frac{9}{11}u_{n-1}+\frac{2}{11}u_{n-2}+\frac{6}{11}h\varphi(t_{n+1},u_{n+1})& n=2,3,\dots N-1 \end{cases} \)

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

\( x\mapsto -x+\frac{18}{11}u_n-\frac{9}{11}u_{n-1}+\frac{2}{11}u_{n-2}+\frac{6}{11}h\varphi(t_{n+1},x). \)

def BDF3(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial condition
    uu[0] = y0
    temp = fsolve(lambda x: -x+uu[0]+h*phi(tt[1],x), uu[0])
    uu[1] = temp[0]
    temp = fsolve(lambda x: -x+4/3*uu[1]-1/3*uu[0] + 2/3*h*phi(tt[2],x), uu[1])
    uu[2] = temp[0]
    # time-stepping loop
    for n in range(2,len(tt)-1):
        temp = fsolve(lambda x: -x+18/11*uu[n]-9/11*uu[n-1] + 2/11*uu[n-2]+6/11*h*phi(tt[n+1],x) , uu[n])
        uu[n+1] = temp[0]
    return uu

Schéma RK1_M

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

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

\( x\mapsto -x+u_n+h\varphi\left(\frac{t_n+t_{n+1}}{2},\frac{u_n+x}{2}\right). \)

Rappelons-nous que ce schéma est bien un schéma de Runge-Kutta:

\( \begin{array}{c|cc} \frac{1}{2} & \frac{1}{2}\\ \hline & 1 \end{array} \quad \implies \begin{cases} u_0 = y_0 \\ K_1 = \varphi\left(t_n+\frac{h}{2},u_n+\frac{h}{2}K_1\right)\\ u_{n+1} = u_n + h K_1 & n=0,1,\dots N-1 \end{cases} \)

avec \(K_1\) un zéro de la fonction

\( x\mapsto -x+\varphi\left(t_n+\frac{h}{2},u_n+\frac{h}{2}x\right). \)

def RK1_M(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial condition
    uu[0] = y0
    # time-stepping loop
    for n in range(len(tt)-1):
        uu[n+1] = fsolve(lambda x: -x+uu[n]+h*phi( (tt[n]+tt[n+1])/2,(uu[n]+x)/2 ), uu[n])[0]
    return uu

def RK1_M_bis(phi,tt,sol_exacte):
    h  = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial condition
    uu[0] = y0
    # time-stepping loop
    for n in range(len(tt)-1):
        k1 = fsolve( lambda x : -x + phi(tt[n]+h/2, uu[n]+h*x/2), phi(tt[n],uu[n]) )[0]    
        uu[n+1] = uu[n]+h*k1
    return uu

Schéma DIRK-2

Schéma de Runge-Kutta semi-implicite à 2 étages

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

avec \(K_1\) zéro de la fonction

\( x\mapsto -x+\varphi\left(t_n+\frac{h}{3},u_n+\frac{h}{3}x\right). \)

Le point de départ du fsolve n’est donc pas \(u_n\) mais \(\varphi_{n}\).

def DIRK2(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial condition
    uu[0] = y0
    # time-stepping loop
    for n in range(len(tt)-1):
        k1 = fsolve(lambda kappa : -kappa+phi(tt[n]+h/3,uu[n]+h/3*kappa), phi(tt[n],uu[n])  )[0]
        k2 = phi( tt[n+1], uu[n]+h*k1 )
        uu[n+1] = uu[n]+h/4*(3*k1+k2)
    return uu

Schéma Randau-IIa-2

Schéma de Runge-Kutta implicite à 2 étages

\(\begin{cases} u_0 = y_0 \\ K_1 = \varphi\left(t_n+\frac{h}{3},u_n+\frac{h}{12}(5K_1-K_2)\right)\\ K_2 = \varphi\left(t_{n+1},u_n+\frac{h}{4}(3K_1+K_2)\right)\\ u_{n+1} = u_n + \frac{h}{4}\left(3K_1+K_2\right) & n=0,1,\dots N-1 \end{cases} \)

avec \(K_1\) et \(K_2\) solutions des équations

\( \begin{cases} K_1 = \varphi\left(t_n+\frac{h}{3},u_n+\frac{h}{12}(5K_1-K_2)\right)\\ K_2 = \varphi\left(t_{n+1},u_n+\frac{h}{4}(3K_1+K_2)\right) \end{cases} \)

soit encore zéro de la fonction de \(\mathbb{R}^2\) dans \(\mathbb{R}^2\)

\( (x,y)\mapsto \begin{pmatrix} x-\varphi\left(t_n+\frac{h}{3},u_n+\frac{h}{12}(5x-y)\right)\\ y-\varphi\left(t_{n+1},u_n+\frac{h}{4}(3x+y)\right)\end{pmatrix}. \)

Le point de départ du fsolve n’est donc pas \((u_n,u_n)\) mais \((\varphi_{n},\varphi_{n})\).

def Randau(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial condition
    uu[0] = y0
    # time-stepping loop
    for n in range(len(tt)-1):
          systeme = lambda K: [  K[0] - phi(tt[n] + h / 3, uu[n] + h / 12 * (5 * K[0] - K[1])),
                               K[1] - phi(tt[n + 1], uu[n] + h / 4 * (3 * K[0] + K[1]))
                               ]
          K_start = phi(tt[n], uu[n])
          k1, k2 = fsolve(systeme, [K_start, K_start])
          uu[n+1] = uu[n]+h/4*(3*k1+k2)
    return uu

Schémas predicteur-correcteur

Schéma de Heun

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

def heun(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial condition
    uu[0] = y0
    # time-stepping loop
    for n in range(len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        k2 = phi( tt[n+1], uu[n] + h*k1  )
        uu[n+1] = uu[n] + (k1+k2)*h/2
    return uu

Schéma AM-2 AB-1

\( \begin{cases} u_0=y_0,\\ u_1=u_0+h\varphi_{0},\\ \tilde u=u_n+h\varphi_{n},\\ u_{n+1}=u_n+\frac{h}{12}\Bigl(5\varphi(t_{n+1},\tilde u)+8\varphi_{n}-\varphi_{n-1}\Bigr)& n=1,2,\dots N-1 \end{cases} \)

def AM2AB1(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial condition
    uu[0] = y0
    uu[1] = uu[0]+h*phi(tt[0],uu[0])
    # time-stepping loop
    for n in range(1,len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        k2 = phi( tt[n-1], uu[n-1] )
        pred = uu[n] + h*k1
        uu[n+1] = uu[n]+h*(5*phi(tt[n+1],pred)+8*k1-k2)/12
    return uu

Schéma AM-3 AB-2

\( \begin{cases} u_0=y_0,\\ u_1=u_0+h\varphi_{0},\\ u_2=u_0+\frac{h}{2}(3\varphi_{1}-\varphi_{0}),\\ \tilde u=u_n+\frac{h}{2}(3\varphi_{n}-\varphi_{n-1}),\\ u_{n+1}=u_n+\frac{h}{24}\Bigl(9\varphi(t_{n+1},\tilde u)+19\varphi_{n}-5\varphi_{n-1}+\varphi_{n-2}\Bigr)& n=2,3,\dots N-1 \end{cases} \)

def AM3AB2(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initial condition
    uu[0] = y0
    uu[1] = uu[0]+h*phi(tt[0],uu[0])
    uu[2] = uu[0]+0.5*h*(3*phi(tt[1],uu[1])-phi(tt[0],uu[0]))
    # time-stepping loop
    for n in range(2,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] )
        pred = uu[n] + (3*k1-k2)*h/2
        uu[n+1] = uu[n]+h*(9*phi(tt[n+1],pred)+19*k1-5*k2+k3)/24
    return uu

Solution exacte

Nous allons calculer la solution exacte en utilisant le module sympy.

t  = sym.Symbol('t')
y  = sym.Function('y')

edo = sym.Eq( sym.diff(y(t),t) , sym.sin(t)+y(t) ) # Attention : c'est le sinus de sympy
display(edo)

t0 = 0
y0 = 0
display(sym.Eq(y(t0),y0))

solpar = sym.dsolve(edo,y(t), ics={y(t0):y0})
display(solpar)

\(\displaystyle \frac{d}{d t} y{\left(t \right)} = y{\left(t \right)} + \sin{\left(t \right)}\)

\(\displaystyle y{\left(0 \right)} = 0\)

\(\displaystyle y{\left(t \right)} = \frac{e^{t}}{2} - \frac{\sin{\left(t \right)}}{2} - \frac{\cos{\left(t \right)}}{2}\)

On définit la solution exacte à utiliser pour estimer les erreurs :

sol_exacte = sym.lambdify(t,solpar.rhs,'numpy')

Tests

On initialise le problème de Cauchy et on définit l’équation différentielle : phi est une fonction python qui contient la fonction mathématique \(\varphi \colon [t_0;T]\times\mathbb{R}\to\mathbb{R}\) définie per \(\varphi(t, y)=\sin(t)+y\) dépendant des variables \(t\) et \(y\).

On introduit la discrétisation : les \(N+1\) nœuds d’intégration \([t_0,t_1,\dots,t_{N}]\) sont contenus dans le vecteur tt et sont espacés de \(h=\frac{t_N-t_0}{N}\).

t0     = 0
tfinal = 1
y0 = 0

phi = lambda t,y : np.sin(t) + y # Attention : c'est le sinus de numpy

N = 8
tt = np.linspace(t0, tfinal, N+1)

On calcule les solutions exacte et approchées et on les compare.

Nous pouvons utiliser deux méthodes différentes pour calculer et afficher les solutions.

Méthode 1

La première méthode est celle utilisée lors des deux premiers TP: - on crée autant de liste que de solutions approchée - on compare les graphes des solutions exacte et approchées - on affiche le maximum de l’erreur

yy = sol_exacte(tt)

uu_EE     = EE(phi,tt,y0)
uu_AB2    = AB2(phi,tt,y0)
uu_AB3    = AB3(phi,tt,y0)
uu_AB4    = AB4(phi,tt,y0)
uu_AB5    = AB5(phi,tt,y0)
uu_N2     = N2(phi,tt,y0)
uu_N3     = N3(phi,tt,y0)
uu_N4     = N4(phi,tt,y0)
uu_EM     = EM(phi,tt,y0)
uu_RK4    = RK4(phi,tt,y0)
uu_RK6_5  = RK6_5(phi,tt,y0)
uu_RK7_6  = RK7_6(phi,tt,y0)

uu_EI   = EI(phi,tt,y0)
uu_CN   = CN(phi,tt,y0)
uu_AM2  = AM2(phi,tt,y0)
uu_AM3  = AM3(phi,tt,y0)
uu_AM4  = AM4(phi,tt,y0)
uu_AM5  = AM5(phi,tt,y0)
uu_BDF2 = BDF2(phi,tt,y0)
uu_BDF3 = BDF3(phi,tt,y0)
uu_RK1_M  = RK1_M(phi,tt,y0)
uu_RK1_M_bis = RK1_M_bis(phi,tt,y0)
uu_DIRK2  = DIRK2(phi,tt,y0)
uu_Randau = Randau(phi,tt,y0)

uu_heun   = heun(phi,tt,y0)
uu_AM2AB1 = AM2AB1(phi,tt,y0)
uu_AM3AB2 = AM3AB2(phi,tt,y0)

fig, ax = plt.subplots(nrows=9, ncols=3, figsize=(18,32))
ax = ax.reshape(-1)


i = 0
ax[i].plot(tt,yy,'b-',tt,uu_EE,'r-D')
err = np.linalg.norm(uu_EE-yy, np.inf)
ax[i].set_title(f'AB1=EE\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_AB2,'r-D')
err = np.linalg.norm(uu_AB2-yy, np.inf)
ax[i].set_title(f'AB2\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_AB3,'r-D')
err = np.linalg.norm(uu_AB3-yy, np.inf)
ax[i].set_title(f'AB3\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_AB4,'r-D')
err = np.linalg.norm(uu_AB4-yy, np.inf)
ax[i].set_title(f'AB4\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_AB5,'r-D')
err = np.linalg.norm(uu_AB5-yy, np.inf)
ax[i].set_title(f'AB5\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_N2,'r-D')
err = np.linalg.norm(uu_N2-yy, np.inf)
ax[i].set_title(f'N2\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_N3,'r-D')
err = np.linalg.norm(uu_N3-yy, np.inf)
ax[i].set_title(f'N3\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_N4,'r-D')
err = np.linalg.norm(uu_N4-yy, np.inf)
ax[i].set_title(f'N4\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_EM,'r-D')
err = np.linalg.norm(uu_EM-yy, np.inf)
ax[i].set_title(f'EM\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_RK4,'r-D')
err = np.linalg.norm(uu_RK4-yy, np.inf)
ax[i].set_title(f'RK4\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_RK6_5,'r-D')
err = np.linalg.norm(uu_RK6_5-yy, np.inf)
ax[i].set_title(f'RK6_5\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_RK7_6,'r-D')
err = np.linalg.norm(uu_RK7_6-yy, np.inf)
ax[i].set_title(f'RK7_6\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_EI,'r-D')
err = np.linalg.norm(uu_EI-yy, np.inf)
ax[i].set_title(f'AM0=EI\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_CN,'r-D')
err = np.linalg.norm(uu_CN-yy, np.inf)
ax[i].set_title(f'AM1=CN\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_AM2,'r-D')
err = np.linalg.norm(uu_AM2-yy, np.inf)
ax[i].set_title(f'AM2\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_AM3,'r-D')
err = np.linalg.norm(uu_AM3-yy, np.inf)
ax[i].set_title(f'AM3\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_AM4,'r-D')
err = np.linalg.norm(uu_AM4-yy, np.inf)
ax[i].set_title(f'AM4\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_AM5,'r-D')
err = np.linalg.norm(uu_AM5-yy, np.inf)
ax[i].set_title(f'AM5\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_BDF2,'r-D')
err = np.linalg.norm(uu_BDF2-yy, np.inf)
ax[i].set_title(f'BDF2\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_BDF3,'r-D')
err = np.linalg.norm(uu_BDF3-yy, np.inf)
ax[i].set_title(f'BDF3\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_RK1_M,'r-D')
err = np.linalg.norm(uu_RK1_M-yy, np.inf)
ax[i].set_title(f'RK1_M\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_RK1_M_bis,'r-D')
err = np.linalg.norm(uu_RK1_M_bis-yy, np.inf)
ax[i].set_title(f'RK1_M_bis\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_DIRK2,'r-D')
err = np.linalg.norm(uu_DIRK2-yy, np.inf)
ax[i].set_title(f'DIRK2\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_Randau,'r-D')
err = np.linalg.norm(uu_Randau-yy, np.inf)
ax[i].set_title(f'Randau\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_heun,'r-D')
err = np.linalg.norm(uu_heun-yy, np.inf)
ax[i].set_title(f'Heun\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_AM2AB1,'r-D')
err = np.linalg.norm(uu_AM2AB1-yy, np.inf)
ax[i].set_title(f'AM2AB1\nmax(|err|)={err:g}')

i += 1
ax[i].plot(tt,yy,'b-',tt,uu_AM3AB2,'r-D')
err = np.linalg.norm(uu_AM3AB2-yy, np.inf)
ax[i].set_title(f'AM3AB2\nmax(|err|)={err:g}')


fig.tight_layout() ;

Méthode 2

La deuxième méthode fait appelle à la notion de dictionnaire, elle est compacte mais peut-être plus difficile à comprendre. On crée

• une liste avec les noms des schémas (sous forme de chaînes de caractères),
• un dictionnaire avec comme clé les noms des schémas et comme valeur l’array solution approchée
• un dictionnaire avec comme clé les noms des schémas et comme valeur le maximum de la valeur absolue de l’erreur.

Pour appeler les fonctions, on utilise eval qui permet d’évaluer une chaîne de caractères comme une commande python.

yy = sol_exacte(tt)

schemas = ['EE','AB2','AB3','AB4','AB5','N2','N3','N4', 'EM', 'RK4', 'RK6_5', 'RK7_6',
           'EI', 'CN', 'AM2', 'AM3', 'AM4', 'AM5', 'BDF2', 'BDF3', 'RK1_M', 'RK1_M_bis', 'DIRK2', 'Randau',
           'heun', 'AM2AB1', 'AM3AB2']

uu = { s : eval(s)(phi,tt,y0) for s in schemas }
err= { s : np.linalg.norm(uu[s]-yy, np.inf) for s in schemas }

fig, ax = plt.subplots(nrows=9, ncols=3, figsize=(18, 32), constrained_layout=False)
ax = ax.reshape(-1)

idx = 0
for key in uu:
    ax[idx].plot(tt,yy,'b-',tt,uu[key],'r-D')    
    ax[idx].set_title( f'{key}\nmax(|err|)={err[key]:g}' )
    idx += 1

fig.tight_layout() 

Méthode 2bis

La troisième méthode ressemble à la précédente sans utiliser eval(). On crée

• une liste avec les noms des schémas (il s’agit d’une liste de fonctions cette fois-ci),
• un dictionnaire avec comme clé les noms des schémas et comme valeur l’array solution approchée
• un dictionnaire avec comme clé les noms des schémas et comme valeur le maximum de la valeur absolue de l’erreur.

Pour récupérer le nom des fonctions, on utilise .__name__.

yy = sol_exacte(tt)


schemas = [EE, AB2, AB3, AB4, AB5, N2, N3, N4, EM, RK4, RK6_5, RK7_6,
           EI, CN, AM2, AM3, AM4, AM5, BDF2, BDF3, RK1_M, RK1_M_bis, DIRK2, Randau,
           heun, AM2AB1, AM3AB2]

uu = { s.__name__ : s(phi,tt,y0) for s in schemas }
err= { s.__name__ : np.linalg.norm(uu[s.__name__]-yy, np.inf) for s in schemas }

fig, ax = plt.subplots(nrows=9, ncols=3, figsize=(18, 32), constrained_layout=False)
ax = ax.reshape(-1)

idx = 0
for key in uu:
    ax[idx].plot(tt,yy,'b-',tt,uu[key],'r-D')    
    ax[idx].set_title( f'{key}\nmax(|err|)={err[key]:g}' )
    idx += 1

fig.tight_layout()