TP 3 - Implémentation de schémas

import sys #only needed to determine Python version number
print('Python version ' + sys.version)

Python version 3.12.3 (main, Feb  4 2025, 14:48:35) [GCC 13.3.0]

1 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 : # !!!!!!

Remarques 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)

%reset -f
%matplotlib inline
%autosave 300

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">12</span>})

from scipy.optimize import fsolve

import sympy as sym
sym.init_printing()

Autosaving every 300 seconds

1.1 Schémas explicites

1.1.1 Schéma d’Euler progressif = de 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.zeros_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

1.1.2 Schéma de 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.zeros_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

1.1.3 Schéma de 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.zeros_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

1.1.4 Schéma de 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.zeros_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

1.1.5 Schéma de 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.zeros_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

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

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

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

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

def EM(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.zeros_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

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

def RK4(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.zeros_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

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

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

1.4 Schémas implicites

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

def EI(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.zeros_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

1.4.2 Schéma de Crank-Nicolson (AM-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.zeros_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

1.4.3 Schéma de AM-2

\( \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.zeros_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

1.4.4 Schéma de AM-3

\( \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.zeros_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

1.4.5 Schéma de AM-4

\( \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.zeros_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

1.4.6 Schéma de AM-5

\( \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.zeros_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

1.4.7 Schéma BDF2

\( \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.zeros_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   

1.4.8 Schéma BDF3

\( \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.zeros_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

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

def RK1_M(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = np.zeros_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

1.4.10 Schéma DIRK-2 (RK 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.zeros_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

1.4.11 Schéma Randau-IIa-2 (RK 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.zeros_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

1.5 Schémas predicteur-correcteur

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

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

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

1.5.4 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')

1.5.5 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_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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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|)={</span>err<span class="sc">: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', 'DIRK2', 'Randau',
           'heun', 'AM2AB1', 'AM3AB2']

uu = { s : <span class="bu">eval</span>(s)(phi,tt,y0) <span class="cf">for</span> s <span class="kw">in</span> schemas }
err= { s : np.linalg.norm(uu[s]<span class="op">-</span>yy, np.inf) <span class="cf">for</span> s <span class="kw">in</span> 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'{</span>key<span class="sc">}\nmax(|err|)={</span>err[key]<span class="sc">:g}' )
    idx += 1

fig.tight_layout() 

Méthode 2 bis

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, DIRK2, Randau,
           heun, AM2AB1, AM3AB2]

uu = { s.<span class="va">__name__</span> : s(phi,tt,y0) <span class="cf">for</span> s <span class="kw">in</span> schemas }
err= { s.<span class="va">__name__</span> : np.linalg.norm(uu[s.<span class="va">__name__</span>]<span class="op">-</span>yy, np.inf) <span class="cf">for</span> s <span class="kw">in</span> 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'{</span>key<span class="sc">}\nmax(|err|)={</span>err[key]<span class="sc">:g}' )
    idx += 1

fig.tight_layout()