M62 Étude de schémas numériques pour l’approximation de problèmes de Cauchy - TP 3 - Implémentation de schémas

TP 3 - Implémentation de schémas

Auteur·rice

Gloria FACCANONI

Date de publication

8 avril 2024

Code 
import sys #only needed to determine Python version number
print('Python version ' + sys.version)
Python version 3.10.12 (main, Nov 20 2023, 15:14:05) [GCC 11.4.0]

Objectif de ce TP : compléter le notebook en ajoutant l’implémentation des schémas indiqués et en vérifiant l’impléméntation sur l’exemple donné.

1 Implémentation des schémas

Remarques pour l’implémentation

On écrit les schémas numériques : + les \(N+1\) nœuds d’intégration \([t_0,t_1,\dots,t_{N}]\) sont contenus dans le vecteur tt + les \(N+1\) valeurs \([u_0,u_1,\dots,u_{N}]\) pour chaque méthode sont contenues dans le vecteur uu.

Comme len(tt) \(=N+1\) et range(1,N) produit 1,2,3,N-1 : - pour un schéma à un pas \(u_{n+1}=F(u_n)\) on initialise \(u_0\) et on calcule \(u_{n+1}\) 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 \(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) - etc.

Code 
%reset -f
%matplotlib inline
%autosave 300

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

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(t_n,u_n)& n=0,1,2,\dots N-1 \end{cases} \)

Code 
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(t_0,u_0),\\ u_{n+1}=u_n+\dfrac{h}{2}\Bigl(3\varphi(t_n,u_n)-\varphi(t_{n-1},u_{n-1})\Bigr)& n=1,2,3,4,5,\dots N-1 \end{cases} \)

Code 
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(t_0,u_0),\\ u_{2}=u_1+\dfrac{h}{2}\Bigl(3\varphi(t_1,u_1)-\varphi(t_{0},u_{0})\Bigr),\\ u_{n+1}=u_n+\dfrac{h}{12}\Bigl(23\varphi(t_n,u_n)-16\varphi(t_{n-1},u_{n-1})+5\varphi(t_{n-2},u_{n-2})\Bigr)& n=2,3,4,5,\dots N-1 \end{cases} \)

Code 
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(t_0,u_0),\\ u_{2}=u_1+\dfrac{h}{2}\Bigl(3\varphi(t_1,u_1)-\varphi(t_{0},u_{0})\Bigr),\\ u_{3}=u_2+\dfrac{h}{12}\Bigl(23\varphi(t_2,u_2)-16\varphi(t_{1},u_{1})+5\varphi(t_{0},u_{0})\Bigr),\\ u_{n+1}=u_n+\dfrac{h}{24}\Bigl(55\varphi(t_n,u_n)-59\varphi(t_{n-1},u_{n-1})+37\varphi(t_{n-2},u_{n-2})-9\varphi(t_{n-3},u_{n-3})\Bigr)& n=3,4,5,\dots N-1 \end{cases} \)

Code 
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(t_0,u_0),\\ u_{2}=u_1+\dfrac{h}{2}\Bigl(3\varphi(t_1,u_1)-\varphi(t_{0},u_{0})\Bigr),\\ u_{3}=u_2+\dfrac{h}{12}\Bigl(23\varphi(t_2,u_2)-16\varphi(t_{1},u_{1})+5\varphi(t_{0},u_{0})\Bigr),\\ u_{4}=u_3+\dfrac{h}{24}\Bigl(55\varphi(t_3,u_3)-59\varphi(t_{2},u_{2})+37\varphi(t_{1},u_{1})-9\varphi(t_{0},u_{0})\Bigr),\\ u_{n+1}=u_n+\dfrac{h}{720}\Bigl(1901\varphi(t_n,u_n)-2774\varphi(t_{n-1},u_{n-1})+2616\varphi(t_{n-2},u_{n-2})-1274\varphi(t_{n-3},u_{n-3})+251\varphi(t_{n-4},u_{n-4})\Bigr)& n=4,5,\dots N-1 \end{cases} \)

Code 
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(t_0,u_0),\\ u_{n+1}=u_{n-1}+2h\varphi(t_{n},u_{n})& n=1,2,3,4,5,\dots N-1 \end{cases} \)

Code 
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(t_0,u_0),\\ u_{2}=u_{0}+2h\varphi(t_{1},u_{1}),\\ u_{n+1}=u_{n-1}+\frac{h}{3}\Bigl(7\varphi(t_{n},u_{n})-2\varphi(t_{n-1},u_{n-1})+\varphi(t_{n-2},u_{n-2})\Bigr)& n=2,3,4,5,\dots N-1 \end{cases} \)

Code 
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(t_0,u_0),\\ u_{2}=u_{0}+2h\varphi(t_{1},u_{1}),\\ u_{3}=u_{1}+\frac{h}{3}\Bigl(7\varphi(t_{2},u_{2})-2\varphi(t_{1},u_{1})+\varphi(t_{0},u_{0})\Bigr),\\ u_{n+1}=u_{n-1}+\frac{h}{3}\Bigl(8\varphi(t_{n},u_{n})-5\varphi(t_{n-1},u_{n-1})+4\varphi(t_{n-2},u_{n-2})-\varphi(t_{n-3},u_{n-3})\Bigr)& n=3,4,5,\dots N-1 \end{cases} \)

Code 
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(t_n,u_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,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]+h*phi(tt[n]+h/2,uu[n]+k1*h/2)
    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} \)

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

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

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

Code 
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

\( \begin{cases} u_0=y_0,\\ u_{n+1}=u_n+\frac{h}{2}\Bigl(\varphi(t_n,u_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(t_n,u_n)+\varphi(t_{n+1},x))\)

Code 
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(t_1,u_1)+\varphi(t_{0},u_{0})\Bigr),\\ u_{n+1}=u_n+\frac{h}{12}\Bigl(5\varphi(t_{n+1},u_{n+1})+8\varphi(t_n,u_n)-\varphi(t_{n-1},u_{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(t_n,u_n)-\varphi(t_{n-1},u_{n-1})\Bigr). \)

Code 
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(t_1,u_1)+\varphi(t_{0},u_{0})\Bigr),\\ u_{2}=u_1+\frac{h}{12}\Bigl(5\varphi(t_{2},u_{2})+8\varphi(t_1,u_1)-\varphi(t_{0},u_{0})\Bigr),\\ u_{n+1}=u_n+\frac{h}{24}\Bigl(9\varphi(t_{n+1},u_{n+1})+19\varphi(t_n,u_n)-5\varphi(t_{n-1},u_{n-1})+\varphi(t_{n-2},u_{n-2})\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(t_n,u_n)-5\varphi(t_{n-1},u_{n-1})+\varphi(t_{n-2},u_{n-2})\Bigr). \)

Code 
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(t_1,u_1)+\varphi(t_{0},u_{0})\Bigr),\\ u_{2}=u_1+\frac{h}{12}\Bigl(5\varphi(t_{2},u_{2})+8\varphi(t_1,u_1)-\varphi(t_{0},u_{0})\Bigr),\\ u_{3}=u_2+\frac{h}{24}\Bigl(9\varphi(t_{3},u_{3})+19\varphi(t_2,u_2)-5\varphi(t_{1},u_{1})+\varphi(t_{0},u_{0})\Bigr),\\ u_{n+1}=u_n+\frac{h}{720}\Bigl(251\varphi(t_{n+1},u_{n+1})+646\varphi(t_n,u_n)-264\varphi(t_{n-1},u_{n-1})+106\varphi(t_{n-2},u_{n-2})-19\varphi(t_{n-3},u_{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(t_n,u_n)-264\varphi(t_{n-1},u_{n-1})+106\varphi(t_{n-2},u_{n-2})-19\varphi(t_{n-3},u_{n-3})\Bigr). \)

Code 
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(t_1,u_1)+\varphi(t_{0},u_{0})\Bigr),\\ u_{2}=u_1+\frac{h}{12}\Bigl(5\varphi(t_{2},u_{2})+8\varphi(t_1,u_1)-\varphi(t_{0},u_{0})\Bigr),\\ u_{3}=u_2+\frac{h}{24}\Bigl(9\varphi(t_{3},u_{3})+19\varphi(t_2,u_2)-5\varphi(t_{1},u_{1})+\varphi(t_{0},u_{0})\Bigr),\\ u_{4}=u_3+\frac{h}{720}\Bigl(251\varphi(t_{4},u_{4})+646\varphi(t_3,u_3)-264\varphi(t_{2},u_{2})+106\varphi(t_{1},u_{1})-19\varphi(t_{0},u_{0})\Bigr),\\ u_{n+1}=u_n+\frac{h}{1440}\Bigl(475\varphi(t_{n+1},u_{+1})+1427\varphi(t_n,u_n)-798\varphi(t_{n-1},u_{n-1})+482\varphi(t_{n-2},u_{n-2})-173\varphi(t_{n-3},u_{n-3})+27\varphi(t_{n-4},u_{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(t_n,u_n)-798\varphi(t_{n-1},u_{n-1})+482\varphi(t_{n-2},u_{n-2})-173\varphi(t_{n-3},u_{n-3})+27\varphi(t_{n-4},u_{n-4})\Bigr). \)

Code 
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(t_1,u_1),\\ u_{n+1}=\frac{4}{3}u_n-\frac{1}{3}u_{n-1}+\frac{2}{3}h\varphi(t_{n+1},u_{n+1})& n=1,2,\dots N-1 \end{cases} \)

Code 
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(t_1,u_1),\\ u_{2}=\frac{4}{3}u_1-\frac{1}{3}u_{0}+\frac{2}{3}h\varphi(t_{2},u_{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} \)

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

Code 
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.5 Schémas predicteur-correcteur

1.5.1 Schéma de Heun

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

Code 
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(t_0,u_0),\\ \tilde u=u_n+h\varphi(t_n,u_n),\\ u_{n+1}=u_n+\frac{h}{12}\Bigl(5\varphi(t_{n+1},\tilde u)+8\varphi(t_n,u_n)-\varphi(t_{n-1},u_{n-1})\Bigr)& n=1,2,\dots N-1 \end{cases} \)

Code 
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):
        pred = uu[n] + h*phi(tt[n],uu[n])
        uu[n+1] = uu[n]+h*(5*phi(tt[n+1],pred)+8*phi(tt[n],uu[n])-phi(tt[n-1],uu[n-1]))/12
    return uu

1.5.3 Schéma AM-3 AB-2

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

Code 
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] )
        pred = uu[n] + (3*k1-k2)*h/2
        uu[n+1] = uu[n]+h*(5*phi(tt[n+1],pred)+8*phi(tt[n],uu[n])-phi(tt[n-1],uu[n-1]))/12
    return uu

2 Comparaison sur un exemple

Considérons le problème de Cauchy

trouver la fonction \(y \colon I\subset \mathbb{R} \to \mathbb{R}\) définie sur l’intervalle \(I=[0,1]\) telle que \( \begin{cases} y'(t) = \sin(t)+y(t), &\forall t \in I=[0,1],\\ y(0) = 0 \end{cases} \)

  1. Calculer la solution exacte en utilisant le module sympy.
  2. Calculer la solution approchée obtenue avec la méthode d’Euler explicite avec \(h=1/N\) et \(N=8\) (pour bien visualiser les erreurs);
  3. même exercice pour les méthodes multipas explicites;
  4. même exercice pour les méthodes multipas implicites;
  5. même exercice pour les méthodes predictor-corrector.
  6. Pour chaque méthode, afficher solution exacte vs solution approchée ainsi que le maximum de la valeur absolue de l’erreur.

Correction 1
Calculons la solutions exacte en utilisant le module sympy. Attention, dans cette partie le sinus est la fonction de sympy et non pas celle de numpy.

Code 
t0 = 0
y0 = 0

t  = sym.Symbol('t')
y  = sym.Function('y')
edo = sym.Eq( sym.diff(y(t),t) , sym.sin(t)+y(t) )
display(edo)

# solgen = sym.dsolve(edo,y(t))
# display(solgen)

# cond_init = sym.Eq( y0, solgen.rhs.subs(t,t0))
# display(cond_init)

# consts = sym.solve( cond_init , dict=True)[0]
# display(consts)

# solpar = solgen.subs(consts).simplify()
# display(solpar)

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(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:

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

Correction de 2 à 5

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(t, y)=\sin(t)+y\) dépendant des variables \(t\) et \(y\). Attention, le sinus maintenant est celui du module numpy et non pas celui du module sympy.

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.
On a \(N+1\) points espacé de \(h=\frac{t_N-t_0}{N}\).

Code 
t0     = 0
tfinal = 1
y0 = 0

phi = lambda t,y : np.sin(t) + y

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

Code 
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_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_heun   = heun(phi,tt,y0)
uu_AM2AB1 = AM2AB1(phi,tt,y0)
uu_AM3AB2 = AM3AB2(phi,tt,y0)
Code 
fig, ax = plt.subplots(nrows=8, ncols=3, figsize=(18,28))
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_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_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_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_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.

Code 
yy = sol_exacte(tt)

schemas = ['EE','AB2','AB3','AB4','AB5','N2','N3','N4', 'EM', 'RK1_M','RK4', 
           'EI', 'CN', 'AM2', 'AM3', 'AM4', 'AM5', 'BDF2', 'BDF3', 
           '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=8, ncols=3, figsize=(18, 28), 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 3

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

Code 
yy = sol_exacte(tt)

schemas = [EE, AB2, AB3, AB4, AB5, N2, N3, N4, EM, RK1_M, RK4, 
           EI, CN, AM2, AM3, AM4, AM5, BDF2, BDF3, 
           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=8, ncols=3, figsize=(18, 28), 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()

Retour au sommet