None TP3-Multisteps-impl
In [21]:
from IPython.display import display, Latex
from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read()) 
Out[21]:
In [22]:
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é.

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.
In [23]:
%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} $$
In [24]:
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} $$
In [25]:
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} $$
In [26]:
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} $$
In [27]:
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} $$
In [28]:
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} $$
In [29]:
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} $$
In [30]:
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} $$
In [31]:
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} $$
In [32]:
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} $$
In [33]:
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)$$

In [34]:
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))$$

In [35]:
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} $$
In [36]:
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} $$
In [37]:
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} $$
In [38]:
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} $$
In [39]:
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} $$
In [40]:
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} $$
In [41]:
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} $$
In [42]:
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} $$
In [43]:
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} $$
In [44]:
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}$$

  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:

In [45]:
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)
$\displaystyle \frac{d}{d t} y{\left(t \right)} = y{\left(t \right)} + \sin{\left(t \right)}$
$\displaystyle y{\left(t \right)} = C_{1} e^{t} - \frac{\sin{\left(t \right)}}{2} - \frac{\cos{\left(t \right)}}{2}$
$\displaystyle \left\{ C_{1} : \frac{1}{2}\right\}$
$\displaystyle y{\left(t \right)} = \frac{e^{t}}{2} - \frac{\sqrt{2} \sin{\left(t + \frac{\pi}{4} \right)}}{2}$

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

In [46]:
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}$.

In [72]:
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
In [71]:
# 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.
In [73]:
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() 
In [ ]: