from IPython.display import display, Latex
from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())

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

Python version 3.10.6 (main, Nov 14 2022, 16:10:14) [GCC 11.3.0]

M62 TP 3 - Implémentation de schémas¶

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
2 Comparaison sur un exemple

Implémentation des schémas¶

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.

%reset -f
%matplotlib inline
%autosave 300

from matplotlib.pylab import *
# rcdefaults()
rcParams.update({'font.size': 12})
from scipy.optimize import fsolve

Autosaving every 300 seconds

Schémas explicites¶

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

def EE(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    for i in range(len(tt)-1):
        k1 = phi( tt[i], uu[i] )
        uu.append(uu[i]+k1*h)
    return uu

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

def AB2(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    uu.append(uu[0]+h*phi(tt[0],uu[0]))
    for i in range(1,len(tt)-1):
        k1 = phi( tt[i], uu[i] )
        k2 = phi( tt[i-1], uu[i-1] )
        uu.append( uu[i] + (3*k1-k2)*h/2 )
    return uu

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

def AB3(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    uu.append(uu[0]+h*phi(tt[0],uu[0]))
    uu.append(uu[1]+h*(3*phi(tt[1],uu[1])-phi(tt[0],uu[0]))/2)
    for i in range(2,len(tt)-1):
        k1 = phi( tt[i], uu[i] )
        k2 = phi( tt[i-1], uu[i-1] )
        k3 = phi( tt[i-2], uu[i-2] )
        uu.append( uu[i] + (23*k1-16*k2+5*k3)*h/12 )
    return uu

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

def AB4(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    uu.append(uu[0]+h*phi(tt[0],uu[0]))
    uu.append(uu[1]+h*(3*phi(tt[1],uu[1])-phi(tt[0],uu[0]))/2)
    uu.append(uu[2]+h*(23*phi(tt[2],uu[2])-16*phi(tt[1],uu[1])+5*phi(tt[0],uu[0]))/12)
    for i in range(3,len(tt)-1):
        k1 = phi( tt[i], uu[i] )
        k2 = phi( tt[i-1], uu[i-1] )
        k3 = phi( tt[i-2], uu[i-2] )
        k4 = phi( tt[i-3], uu[i-3] )
        uu.append( uu[i] + (55*k1-59*k2+37*k3-9*k4)*h/24 )
    return uu

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

def AB5(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    uu.append(uu[0]+h*phi(tt[0],uu[0]))
    uu.append(uu[1]+h*(3*phi(tt[1],uu[1])-phi(tt[0],uu[0]))/2)
    uu.append(uu[2]+h*(23*phi(tt[2],uu[2])-16*phi(tt[1],uu[1])+5*phi(tt[0],uu[0]))/12)
    uu.append(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)
    for i in range(4,len(tt)-1):
        k1 = phi( tt[i], uu[i] )
        k2 = phi( tt[i-1], uu[i-1] )
        k3 = phi( tt[i-2], uu[i-2] )
        k4 = phi( tt[i-3], uu[i-3] )
        k5 = phi( tt[i-4], uu[i-4] )
        uu.append( uu[i] + (1901*k1-2774*k2+2616*k3-1274*k4+251*k5)*h/720 )
    return uu

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

def N2(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    uu.append(uu[0]+h*phi(tt[0],uu[0]))
    for i in range(1,len(tt)-1):
        k1 = phi( tt[i], uu[i] )
        uu.append( uu[i-1] + 2*h*k1 )
    return uu

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

def N3(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    uu.append(uu[0]+h*phi(tt[0],uu[0]))
    uu.append(uu[0]+2*h*phi(tt[1],uu[1]))
    for i in range(2,len(tt)-1):
        k1 = phi( tt[i], uu[i] )
        k2 = phi( tt[i-1], uu[i-1] )
        k3 = phi( tt[i-2], uu[i-2] )
        uu.append( uu[i-1] + (7*k1-2*k2+k3)*h/3 )
    return uu

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

def N4(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    uu.append(uu[0]+h*phi(tt[0],uu[0]))
    uu.append(uu[0]+2*h*phi(tt[1],uu[1]))
    uu.append(uu[1]+(7*h*phi(tt[2],uu[2])-2*h*phi(tt[1],uu[1])+h*phi(tt[0],uu[0]))/3)
    for i in range(3,len(tt)-1):
        k1 = phi( tt[i], uu[i] )
        k2 = phi( tt[i-1], uu[i-1] )
        k3 = phi( tt[i-2], uu[i-2] )
        k4 = phi( tt[i-3], uu[i-3] )
        uu.append( uu[i-1] + (8*k1-5*k2+4*k3-k4)*h/3 )
    return uu

Schéma d'Euler modifié¶

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

def EM(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    for i in range(len(tt)-1):
        k1 = phi( tt[i], uu[i] )
        uu.append( uu[i]+h*phi(tt[i]+h/2,uu[i]+k1*h/2) )
    return uu

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 = [y0]
    for i in range(len(tt)-1):
        k1 = phi( tt[i], uu[i] )
        k2 = phi( tt[i]+h/2 , uu[i]+k1*h/2 )
        k3 = phi( tt[i]+h/2 , uu[i]+h*k2/2 )
        k4 = phi( tt[i+1]   , uu[i]+h*k3 )
        uu.append( uu[i] + (k1+2*k2+2*k3+k4)*h/6 )
    return uu

Schémas implicites¶

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 = [y0]
    for i in range(len(tt)-1):
        temp = fsolve(lambda x: -x+uu[i]+h*phi(tt[i+1],x), uu[i])
        uu.append(temp[0])
    return uu

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

def CN(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    for i in range(len(tt)-1):
        temp = fsolve(lambda x: -x+uu[i]+0.5*h*( phi(tt[i+1],x)+phi(tt[i],uu[i]) ), uu[i])
        uu.append(temp[0])
    return uu

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

def AM2(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    temp = fsolve(lambda x: -x+uu[0]+0.5*h*( phi(tt[1],x)+phi(tt[0],uu[0]) ), uu[0])
    uu.append(temp[0])
    for i in range(1,len(tt)-1):
        temp = fsolve(lambda x: -x+uu[i]+h*( 5*phi(tt[i+1],x)+8*phi(tt[i],uu[i])-phi(tt[i-1],uu[i-1]) )/12, uu[i])
        uu.append(temp[0])
    return uu

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

def AM3(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    temp = fsolve(lambda x: -x+uu[0]+0.5*h*( phi(tt[1],x)+phi(tt[0],uu[0]) ), uu[0])
    uu.append(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.append(temp[0])
    for i in range(2,len(tt)-1):
        temp = fsolve(lambda x: -x+uu[i]+h*( 9*phi(tt[i+1],x)+19*phi(tt[i],uu[i])-5*phi(tt[i-1],uu[i-1])+phi(tt[i-2],uu[i-2]) )/24, uu[i])
        uu.append(temp[0])
    return uu

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

def AM4(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    temp = fsolve(lambda x: -x+uu[0]+0.5*h*( phi(tt[1],x)+phi(tt[0],uu[0]) ), uu[0])
    uu.append(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.append(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.append(temp[0])
    for i in range(3,len(tt)-1):
        temp = fsolve(lambda x: -x+uu[i]+h*( 251*phi(tt[i+1],x)+646*phi(tt[i],uu[i])-264*phi(tt[i-1],uu[i-1])+106*phi(tt[i-2],uu[i-2])-19*phi(tt[i-3],uu[i-3]) )/720, uu[i])
        uu.append(temp[0])
    return uu

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

def BDF2(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    temp = fsolve(lambda x: -x+uu[0]+h*phi(tt[1],x), uu[0])
    uu.append(temp[0])
    for i in range(1,len(tt)-1):
        temp = fsolve(lambda x: -x+4/3*uu[i]-1/3*uu[i-1] + 2/3*h*phi(tt[i+1],x) , uu[i])
        uu.append(temp[0])
    return uu

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

def BDF3(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    temp = fsolve(lambda x: -x+uu[0]+h*phi(tt[1],x), uu[0])
    uu.append(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.append(temp[0])
    for i in range(2,len(tt)-1):
        temp = fsolve(lambda x: -x+18/11*uu[i]-9/11*uu[i-1] + 2/11*uu[i-2]+6/11*h*phi(tt[i+1],x) , uu[i])
        uu.append(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} $$

def RK1_M(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = [y0]
    for i in range(len(tt)-1):
        uu.append( fsolve(lambda x: -x+uu[i]+h*phi( (tt[i]+tt[i+1])/2,(uu[i]+x)/2 ), uu[i])[0] )
    return uu

Schémas predicteur-correcteur¶

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

def heun(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    for i in range(len(tt)-1):
        k1 = phi( tt[i], uu[i] )
        k2 = phi( tt[i+1], uu[i] + h*k1  )
        uu.append( uu[i] + (k1+k2)*h/2 )
    return uu

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

def AM2AB1(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    uu.append(uu[0]+h*phi(tt[0],uu[0]))
    for i in range(1,len(tt)-1):
        pred = uu[i] + h*phi(tt[i],uu[i])
        uu.append(uu[i]+h*(5*phi(tt[i+1],pred)+8*phi(tt[i],uu[i])-phi(tt[i-1],uu[i-1]))/12)
    return uu

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

def AM3AB2(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    uu.append(uu[0]+h*phi(tt[0],uu[0]))
    uu.append(uu[0]+0.5*h*(3*phi(tt[1],uu[1])-phi(tt[0],uu[0])))
    for i in range(2,len(tt)-1):
        k1 = phi( tt[i], uu[i] )
        k2 = phi( tt[i-1], uu[i-1] )
        pred = uu[i] + (3*k1-k2)*h/2
        uu.append(uu[i]+h*(5*phi(tt[i+1],pred)+8*phi(tt[i],uu[i])-phi(tt[i-1],uu[i-1]))/12)
    return uu

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

Calculer la solution exacte en utilisant le module sympy.
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);
même exercice pour les méthodes multipas explicites;
même exercice pour les méthodes multipas implicites;
même exercice pour les méthodes predictor-corrector.
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:

import sympy as sym
sym.init_printing()

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)

t0=0
y0=0
consts = sym.solve( sym.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)[0]
display(consts)
solpar=solgen.subs(consts).simplify()
display(solpar)

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

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

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

t0     = 0
tfinal = 1
y0 = 0

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

N = 8
tt = 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

# Attention, les deux listes suivantes ne donnent pas les mêmes valeurs que la troisième

# yy_tmp1 = [sol_exacte(t) for t in tt]
# print(f'{yy_tmp1=}')

# yy_tmp2 = list(map(sol_exacte,tt))
# print(f'{yy_tmp2=}')

yy = sol_exacte(tt)
# print(f'{yy=}')


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_RK1_M  = RK1_M(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_BDF2 = BDF2(phi,tt,y0)
uu_BDF3 = BDF3(phi,tt,y0)

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


fig = figure(figsize=(28,12))

subplot(3,7,1)
plot(tt,yy,'b-',tt,uu_EE,'r-D')
err = norm(uu_EE-yy,inf)
title(f'AB1=EE\nmax(|err|)={err:g}')

subplot(3,7,2)
plot(tt,yy,'b-',tt,uu_AB2,'r-D')
err = norm(uu_AB2-yy,inf)
title(f'AB2\nmax(|err|)={err:g}')

subplot(3,7,3)
plot(tt,yy,'b-',tt,uu_AB3,'r-D')
err = norm(uu_AB3-yy,inf)
title(f'AB3\nmax(|err|)={err:g}')

subplot(3,7,4)
plot(tt,yy,'b-',tt,uu_AB4,'r-D')
err=norm(uu_AB4-yy,inf)
title(f'AB4\nmax(|err|)={err:g}')

subplot(3,7,5)
plot(tt,yy,'b-',tt,uu_AB5,'r-D')
err=norm(uu_AB5-yy,inf)
title(f'AB5\nmax(|err|)={err:g}')

subplot(3,7,6)
plot(tt,yy,'b-',tt,uu_N2,'r-D')
err=norm(uu_N2-yy,inf)
title(f'N2\nmax(|err|)={err:g}')

subplot(3,7,7)
plot(tt,yy,'b-',tt,uu_N3,'r-D')
err=norm(uu_N3-yy,inf)
title(f'N3\nmax(|err|)={err:g}')

subplot(3,7,8)
plot(tt,yy,'b-',tt,uu_N4,'r-D')
err=norm(uu_N4-yy,inf)
title(f'N4\nmax(|err|)={err:g}')

subplot(3,7,9)
plot(tt,yy,'b-',tt,uu_EM,'r-D')
err=norm(uu_EM-yy,inf)
title(f'EM\nmax(|err|)={err:g}')

subplot(3,7,10)
plot(tt,yy,'b-',tt,uu_RK1_M,'r-D')
err=norm(uu_RK1_M-yy,inf)
title(f'RK1_M\nmax(|err|)={err:g}')

subplot(3,7,11)
plot(tt,yy,'b-',tt,uu_RK4,'r-D')
err=norm(uu_RK4-yy,inf)
title(f'RK4\nmax(|err|)={err:g}')

subplot(3,7,12)
plot(tt,yy,'b-',tt,uu_EI,'r-D')
err=norm(uu_EI-yy,inf)
title(f'AM0=EI\nmax(|err|)={err:g}')

subplot(3,7,13)
plot(tt,yy,'b-',tt,uu_CN,'r-D')
err=norm(uu_CN-yy,inf)
title(f'AM1=CN\nmax(|err|)={err:g}')

subplot(3,7,14)
plot(tt,yy,'b-',tt,uu_AM2,'r-D')
err=norm(uu_AM2-yy,inf)
title(f'AM2\nmax(|err|)={err:g}')

subplot(3,7,15)
plot(tt,yy,'b-',tt,uu_AM3,'r-D')
err=norm(uu_AM3-yy,inf)
title(f'AM3\nmax(|err|)={err:g}')

subplot(3,7,16)
plot(tt,yy,'b-',tt,uu_AM4,'r-D')
err=norm(uu_AM4-yy,inf)
title(f'AM4\nmax(|err|)={err:g}')

subplot(3,7,17)
plot(tt,yy,'b-',tt,uu_BDF2,'r-D')
err=norm(uu_BDF2-yy,inf)
title(f'BDF2\nmax(|err|)={err:g}')

subplot(3,7,18)
plot(tt,yy,'b-',tt,uu_BDF3,'r-D')
err=norm(uu_BDF3-yy,inf)
title(f'BDF3\nmax(|err|)={err:g}')

subplot(3,7,19)
plot(tt,yy,'b-',tt,uu_heun,'r-D')
err=norm(uu_heun-yy,inf)
title(f'Heun\nmax(|err|)={err:g}')

subplot(3,7,20)
plot(tt,yy,'b-',tt,uu_AM2AB1,'r-D')
err=norm(uu_AM2AB1-yy,inf)
title(f'AM2AB1\nmax(|err|)={err:g}')

subplot(3,7,21)
plot(tt,yy,'b-',tt,uu_AM3AB2,'r-D')
err=norm(uu_AM3AB2-yy,inf)
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,
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.

yy = sol_exacte(tt)

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

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

fig, ax = subplots(nrows=3, ncols=7, figsize=(28, 12),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|)=\n{err[key]:g}' )
    idx += 1

fig.tight_layout()