None TP4-Multisteps-conv
In [1]:
from IPython.display import display, Latex
from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())
Out[1]:
In [2]:
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 4 - Étude de la convergence

In [3]:
%reset -f
%matplotlib inline
%autosave 300

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

Autosaving every 300 seconds

Exercice

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) = y(t), &\forall t \in I=[0,1],\\ y(0) = 1 \end{cases} $$ dont la solution est $y(t)=e^{t}$.

Estimer l'ordre de convergence des méthodes:

  1. EE, AB$_2$, AB$_3$, AB$_4$, AB$_5$, N$_2$, N$_3$, N$_4$, EM, RK4$_1$, RK6$_5$, RK7$_6$
  2. EI, CN, AM$_2$, AM$_3$, AM$_4$, AM$_5$, MA$_2$, BDF$_2$, BDF$_3$, RK1_M
  3. Heun, AM$_4$-AB$_1$, AM$_4$-AB$_2$, AM$_4$-AB$_3$

Remarque: puisque la fonction $\varphi(t,y)=y$ est linéaire, toute méthode implicite peut être rendue explicite par un calcul élémentaire en explicitant directement pour chaque schéma l'expression de $u_{n+1}$. Cependant, nous pouvons utiliser le le module SciPy sans modifier l'implémentation des schémas (mais on payera l'ordre de convergence de fsolve).

Attention: les schémas multistep ont besoin d'initialiser plusieurs pas de la suite définie pas récurrence pour pouvoir démarrer. Dans cette étude, au lieu d'utiliser un schéma d'ordre inférieur pour initialiser la suite, on utilisera la solution exacte (en effet, l'utilisation d'un schéma d'ordre inférieur dégrade l'ordre de précision).

Rappels de la démarche pour le calcul de l'orde de convergence. + On écrit les schémas numériques : + les nœuds d'intégration $[t_0,t_1,\dots,t_{N}]$ sont contenus dans le vecteur `tt` (qui change en fonction de `h`) + les valeurs $[u_0,u_1,\dots,u_{N}]$ pour chaque méthode sont contenues dans le vecteur `uu`. + Pour chaque schéma, on calcule la solution approchée avec différentes valeurs de $h_k=1/N_k$. On sauvegarde les valeurs de $h_k$ dans le vecteur `H`. + Pour chaque valeur de $h_k$, on calcule le maximum de la valeur absolue de l'erreur et on sauvegarde toutes ces erreurs dans le vecteur `err_schema` de sort que `err_schema[k]` contient $e_k=\max_{i=0,\dots,N_k}|y(t_i)-u_{i}|$. + Pour afficher l'ordre de convergence on utilise une échelle logarithmique, i.e. on représente $\ln(h)$ sur l'axe des abscisses et $\ln(\text{err})$ sur l'axe des ordonnées. En effet, si $\text{err}=Ch^p$ alors $\ln(\text{err})=\ln(C)+p\ln(h)$. En échelle logarithmique, $p$ représente donc la pente de la ligne droite $\ln(\text{err})$.

Schémas explicites

Schéma de Adam-Bashforth à 1 pas = schéma d'Euler progressif

$$ \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 [4]:
def EE(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[0])]
    for i in range(N):
        uu.append(uu[i]+h*phi(tt[i],uu[i]))
    return uu

Schéma de Adam-Bashforth à 2 pas

$$ \begin{cases} u_0=y_0,\\ u_{1}=y_1,\\ u_{n+1}=u_n+\frac{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 [5]:
def AB2(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[i]) for i in range(2)]
    for i in range(1,N):
        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}=y_1,\\ u_{2}=y_2,\\ u_{n+1}=u_n+\frac{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 [6]:
def AB3(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[i]) for i in range(3)]
    for i in range(2,N):
        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}=y_1,\\ u_{2}=y_2,\\ u_{3}=y_3,\\ u_{n+1}=u_n+\frac{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 [7]:
def AB4(phi,tt,sol_exacte):
    h = tt[1]-tt[0]    
    uu = [sol_exacte(tt[i]) for i in range(4)]
    for i in range(3,N):
        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}=y_1,\\ u_{2}=y_2,\\ u_{3}=y_3,\\ u_{4}=y_4,\\ u_{n+1}=u_n+\frac{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 [8]:
def AB5(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[i]) for i in range(5)]
    for i in range(4,N):
        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}=y_1,\\ u_{n+1}=u_{n-1}+2h\varphi(t_{n},u_{n})& n=1,2,3,4,5,\dots N-1 \end{cases} $$
In [9]:
def N2(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[0])]
    uu.append(sol_exacte(tt[1]))
    for i in range(1,N):
        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}=y_1,\\ u_{2}=y_2,\\ 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 [10]:
def N3(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[i]) for i in range(3)]
    for i in range(2,N):
        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}=y_1,\\ u_{2}=y_2,\\ u_{3}=y_3,\\ 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 [11]:
def N4(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[i]) for i in range(4)]
    for i in range(3,N):
        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 [12]:
def EM(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[0])]
    for i in range(N):
        k1 = h * phi( tt[i], uu[i] )
        uu.append( uu[i]+h*phi(tt[i]+h/2,uu[i]+k1/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 [13]:
def RK4(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[0])]
    for i in range(N):
        k1 = phi( tt[i], uu[i] )
        k2 = phi( tt[i]+h/2 , uu[i]+h*k1/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é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}$$
In [14]:
def RK6_5(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[0])]
    for i in range(N):
        k1 = phi( tt[i]      , uu[i] )
        k2 = phi( tt[i]+h/4  , uu[i]+h*k1/4 )
        k3 = phi( tt[i]+h/4  , uu[i]+h*(k1+k2)/8 )
        k4 = phi( tt[i]+h/2  , uu[i]+h*(-k2+2*k3)/2 )
        k5 = phi( tt[i]+h*3/4, uu[i]+h*(3*k1+9*k4)/16 )
        k6 = phi( tt[i+1]    , uu[i]+h*(-3*k1+2*k2+12*k3-12*k4+8*k5)/7 )
        uu.append( uu[i] + (7*k1+32*k3+12*k4+32*k5+7*k6)*h/90 )
    return uu

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}$$
In [15]:
def RK7_6(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[0])]
    for i in range(N):
        k1 = phi( tt[i]      , uu[i] )
        k2 = phi( tt[i]+h/2  , uu[i]+h*(        k1                                                     )/2 )
        k3 = phi( tt[i]+h*2/3, uu[i]+h*(      2*k1     +4*k2                                           )/9 )
        k4 = phi( tt[i]+h/3  , uu[i]+h*(      7*k1     +8*k2      -3*k3                                )/36 )
        k5 = phi( tt[i]+h*5/6, uu[i]+h*(    -35*k1   -220*k2    +105*k3    +270*k4                     )/144 )
        k6 = phi( tt[i]+h/6  , uu[i]+h*(       -k1   -110*k2     -45*k3    +180*k4     +36*k5          )/360 )
        k7 = phi( tt[i+1]    , uu[i]+h*(-41/260*k1 +22/13*k2 +43/156*k3 -118/39*k4 +32/195*k5 +80/39*k6)     )
        uu.append( uu[i] + (13/200*k1+11/40*k3+11/40*k4+4/25*k5+4/25*k6+13/200*k7)*h )
    return uu

Schémas implicites

Schéma 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 [16]:
def EI(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[0])]
    for i in range(N):
        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}$ un zéro de la fonction $x\mapsto -x+u_n+\frac{h}{2}(\varphi(t_n,u_n)+\varphi(t_{n+1},x))$.

In [17]:
def CN(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[0])]
    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=y_1,\\ 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 [18]:
def AM2(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[i]) for i in range(2)]
    for i in range(1,N):
        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=y_1,\\ u_{2}=y_2,\\ 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 [19]:
def AM3(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[i]) for i in range(3)]
    for i in range(2,N):
        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])[0]
        uu.append(temp)
    return uu

Schéma de AM-4

$$ \begin{cases} u_0=y_0,\\ u_1=y_1,\\ u_{2}=y_2,\\ u_{3}=y_3,\\ 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 [20]:
def AM4(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[i]) for i in range(4)]
    for i in range(3,N):
        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])[0]
        uu.append(temp)
    return uu

Schéma de AM-5

$$ \begin{cases} u_0=y_0,\\ u_1=y_1,\\ u_{2}=y_2,\\ u_{3}=y_3,\\ u_{4}=y_4,\\ 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} $$
In [21]:
def AM5(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[i]) for i in range(5)]
    for i in range(4,N):
        temp = fsolve(lambda x: -x+uu[i]+h*( 475*phi(tt[i+1],x)+1427*phi(tt[i],uu[i])-798*phi(tt[i-1],uu[i-1])+482*phi(tt[i-2],uu[i-2])-173*phi(tt[i-3],uu[i-3])+27*phi(tt[i-4],uu[i-4]))/1440, uu[i])[0]
        uu.append(temp)
    return uu

Schéma MS-2

$$ \begin{cases} u_0=y_0,\\ u_1=y_1,\\ u_{n+1}=u_{n-1}+\frac{h}{3}\Bigl(\varphi(t_{n+1},u_{n+1})+4\varphi(t_n,u_n)+\varphi(t_{n-1},u_{n-1})\Bigr)& n=1,2,\dots N-1 \end{cases} $$
In [22]:
def MS2(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[i]) for i in range(2)]
    for i in range(1,N):
        temp = fsolve(lambda x: -x+uu[i-1]+h*(phi(tt[i+1],x)+4*phi(tt[i],uu[i])+phi(tt[i-1],uu[i-1]) )/3, uu[i])[0]
        uu.append(temp)
    return uu

Schéma BDF2

$$ \begin{cases} u_0=y_0,\\ u_1=y_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 [23]:
def BDF2(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[i]) for i in range(2)]
    for i in range(1,N):
        temp = fsolve(lambda x: -x+4/3*uu[i]-1/3*uu[i-1] + 2/3*h*phi(tt[i+1],x) , uu[i])[0]
        uu.append(temp)
    return uu

Schéma BDF3

$$ \begin{cases} u_0=y_0,\\ u_1=y_1,\\ u_{2}=y_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 [24]:
def BDF3(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[i]) for i in range(3)]
    for i in range(2,N):
        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])[0]
        uu.append(temp)
    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 [25]:
def RK1_M(phi,tt,sol_exacte):
    h  = tt[1]-tt[0]
    uu = [sol_exacte(tt[0])]
    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 [26]:
def heun(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[0])]
    for i in range(N):
        k1 = phi( tt[i], uu[i] )
        k2 = phi( tt[i+1], uu[i] + k1*h  )
        uu.append( uu[i] + (k1+k2)*h/2 )
    return uu

Schéma AM-4 AB-2/3/4/5

In [27]:
def AM4AB2(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[i]) for i in range(4)]
    for i in range(3,N):
        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*( 251*phi(tt[i+1],pred)+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)
    return uu
     
def AM4AB3(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[0])]
    uu.append(sol_exacte(tt[1]))
    uu.append(sol_exacte(tt[2]))
    uu.append(sol_exacte(tt[3]))
    for i in range(3,N):
        k1 = phi( tt[i], uu[i] )
        k2 = phi( tt[i-1], uu[i-1] )
        k3 = phi( tt[i-2], uu[i-2] )
        pred = uu[i] + (23*k1-16*k2+5*k3)*h/12
        uu.append(uu[i]+h*( 251*phi(tt[i+1],pred)+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)
    return uu
     
def AM4AB4(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[i]) for i in range(4)]
    for i in range(3,N):
        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] )
        pred = uu[i] + (55*k1-59*k2+37*k3-9*k4)*h/24
        uu.append(uu[i]+h*( 251*phi(tt[i+1],pred)+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)
    return uu
     
def AM4AB5(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = [sol_exacte(tt[i]) for i in range(5)]   
    for i in range(4,N):
        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] )
        pred = uu[i] + (1901*k1-2774*k2+2616*k3-1274*k4+251*k5)*h/720
        uu.append(uu[i]+h*( 251*phi(tt[i+1],pred)+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)
    return uu

Convergence

On initialise le problème de Cauchy, on définit l'équation différentielle (phi est une fonction python qui contient la fonction mathématique $\varphi(t, y)$ dépendant des variables $t$ et $y$) et enfin on définit la solution exacte:

In [28]:
t0, y0, tfinal = 0, 1, 3

def sol_exacte(t):
    return exp(t)
    #return exp(-t)
    #return sqrt(2.*t+1.)
    #return sqrt(t**2+1.)
    #return 1./sqrt(1.-2.*t)
    
def phi(t,y):
    return y
    #return -y
    #return 1./y 
    #return t/y 
    #return y**3

Pour chaque schéma, on calcule la solution approchée avec différentes valeurs de $h_k=1/N_k$; on sauvegarde les valeurs de $h_k$ dans le vecteur H.

Pour chaque valeur de $h_k$, on calcule le maximum de la valeur absolue de l'erreur et on sauvegarde toutes ces erreurs dans le vecteur err_schema de sort que err_schema[k] contient $e_k=\max_{i=0,...,N_k}|y(t_i)-u_{i}|$.

Nous pouvons utiliser deux méthodes différentes pour stocker ces informations.

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.

In [29]:
H = []

err_ep   = []
err_AB2  = []
err_AB3  = []
err_AB4  = []
err_AB5  = []
err_N2   = []
err_N3   = []
err_N4   = []
err_RK4  = []
err_RK6_5  = []
err_RK7_6  = []

err_er   = []
err_CN   = []
err_AM2  = []
err_AM3  = []
err_AM4  = []
err_AM5  = []
err_BDF2 = []
err_BDF3 = []
err_MS2=[]

err_em   = []
err_heun = []
err_AM4AB2  = []
err_AM4AB3  = []
err_AM4AB4  = []
err_AM4AB5  = []

N = 10
for k in range(6):
    N += 50 #2**(k+3)
    tt = linspace(t0,tfinal,N+1)
    h = tt[1]-tt[0]
    H.append(h)
    yy = array([sol_exacte(t) for t in tt])
    
    # schemas explicites
    uu_ep   = EE(phi,tt,sol_exacte)
    uu_AB2  = AB2(phi,tt,sol_exacte)
    uu_AB3  = AB3(phi,tt,sol_exacte)
    uu_AB4  = AB4(phi,tt,sol_exacte)
    uu_AB5  = AB5(phi,tt,sol_exacte)
    uu_N2   = N2(phi,tt,sol_exacte)
    uu_N3   = N3(phi,tt,sol_exacte)
    uu_N4   = N4(phi,tt,sol_exacte)
    uu_em   = EM(phi,tt,sol_exacte)
    uu_RK4  = RK4(phi,tt,sol_exacte)
    uu_RK6_5  = RK6_5(phi,tt,sol_exacte)
    uu_RK7_6  = RK7_6(phi,tt,sol_exacte)
    # schemas implicites
    uu_er   = EI(phi,tt,sol_exacte)
    uu_CN   = CN(phi,tt,sol_exacte)
    uu_AM2  = AM2(phi,tt,sol_exacte)
    uu_AM3  = AM3(phi,tt,sol_exacte)
    uu_AM4  = AM4(phi,tt,sol_exacte)
    uu_AM5  = AM5(phi,tt,sol_exacte)
    uu_BDF2 = BDF2(phi,tt,sol_exacte)
    uu_BDF3 = BDF3(phi,tt,sol_exacte)
    uu_MS2  = MS2(phi,tt,sol_exacte)
    # schemas predictor-corrector
    uu_heun = heun(phi,tt,sol_exacte)
    uu_AM4AB2 = AM4AB2(phi,tt,sol_exacte)
    uu_AM4AB3 = AM4AB3(phi,tt,sol_exacte)
    uu_AM4AB4 = AM4AB4(phi,tt,sol_exacte)
    uu_AM4AB5 = AM4AB5(phi,tt,sol_exacte)
    
    # erreurs
    err_ep.append( norm(uu_ep-yy,inf))
    err_AB2.append(norm(uu_AB2-yy,inf))
    err_AB3.append(norm(uu_AB3-yy,inf))
    err_AB4.append(norm(uu_AB4-yy,inf))
    err_AB5.append(norm(uu_AB5-yy,inf))
    err_N2.append( norm(uu_N2-yy,inf))
    err_N3.append(norm(uu_N3-yy,inf))
    err_N4.append(norm(uu_N4-yy,inf))
    err_em.append(norm(uu_em-yy,inf))
    err_RK4.append(norm(uu_RK4-yy,inf))
    err_RK6_5.append(norm(uu_RK6_5-yy,inf))
    err_RK7_6.append(norm(uu_RK7_6-yy,inf))
    
    err_er.append(norm(uu_er-yy,inf))
    err_CN.append(norm(uu_CN-yy,inf))
    err_AM2.append(norm(uu_AM2-yy,inf))
    err_AM3.append(norm(uu_AM3-yy,inf))
    err_AM4.append(norm(uu_AM4-yy,inf))
    err_AM5.append(norm(uu_AM5-yy,inf))
    err_BDF2.append(norm(uu_BDF2-yy,inf))
    err_BDF3.append(norm(uu_BDF3-yy,inf))
    err_MS2.append(norm(uu_MS2-yy,inf))
    
    err_heun.append(norm(uu_heun-yy,inf))
    err_AM4AB2.append(norm(uu_AM4AB2-yy,inf))
    err_AM4AB3.append(norm(uu_AM4AB3-yy,inf))
    err_AM4AB4.append(norm(uu_AM4AB4-yy,inf))
    err_AM4AB5.append(norm(uu_AM4AB5-yy,inf))

Pour estimer l'ordre de convergence on estime la pente de la droite qui relie l'erreur au pas $k$ à l'erreur au pas $k+1$ en echelle logarithmique en utilisant la fonction polyfit basée sur la régression linéaire.

In [30]:
# print (f'EE    {polyfit(log(H),log(err_ep),    1)[0]:1.2f}' )
# print (f'AB2   {polyfit(log(H),log(err_AB2),   1)[0]:1.2f}'  )
# print (f'AB3   {polyfit(log(H),log(err_AB3),   1)[0]:1.2f}' )
# print (f'AB4   {polyfit(log(H),log(err_AB4),   1)[0]:1.2f}' )
# print (f'AB5   {polyfit(log(H),log(err_AB5),   1)[0]:1.2f}' )
# print (f'N2    {polyfit(log(H),log(err_N2),    1)[0]:1.2f}' )
# print (f'N3    {polyfit(log(H),log(err_N3),    1)[0]:1.2f}' )
# print (f'N4    {polyfit(log(H),log(err_N4),    1)[0]:1.2f}' )
# print (f'EM    {polyfit(log(H),log(err_em),    1)[0]:1.2f}' )
# print (f'RK4   {polyfit(log(H),log(err_RK4),   1)[0]:1.2f}' )
# print (f'RK6_5 {polyfit(log(H),log(err_RK6_5), 1)[0]:1.2f}' )
# print (f'RK7_6 {polyfit(log(H),log(err_RK7_6), 1)[0]:1.2f}' )
# print('\n')
# print (f'EI   {polyfit(log(H),log(err_er),   1)[0]:1.2f}' )
# print (f'CN   {polyfit(log(H),log(err_CN),   1)[0]:1.2f}' )
# print (f'AM2  {polyfit(log(H),log(err_AM2),  1)[0]:1.2f}' )
# print (f'AM3  {polyfit(log(H),log(err_AM3),  1)[0]:1.2f}' )
# print (f'AM4  {polyfit(log(H),log(err_AM4),  1)[0]:1.2f}' )
# print (f'AM5  {polyfit(log(H),log(err_AM5),  1)[0]:1.2f}' )
# print (f'BDF2 {polyfit(log(H),log(err_BDF2), 1)[0]:1.2f}' ) 
# print (f'BDF3 {polyfit(log(H),log(err_BDF3), 1)[0]:1.2f}' )
# print (f'MS2  {polyfit(log(H),log(err_MS2),  1)[0]:1.2f}' )
# print('\n')
# print (f'Heun    {polyfit(log(H),log(err_heun),   1)[0]:1.2f}' )
# print (f'AM4_AB2 {polyfit(log(H),log(err_AM4AB2), 1)[0]:1.2f}' )
# print (f'AM4_AB3 {polyfit(log(H),log(err_AM4AB3), 1)[0]:1.2f}' )
# print (f'AM4_AB4 {polyfit(log(H),log(err_AM4AB4), 1)[0]:1.2f}' )
# print (f'AM4_AB5 {polyfit(log(H),log(err_AM4AB5), 1)[0]:1.2f}' )

EE    0.97
AB2   1.98
AB3   2.96
AB4   3.94
AB5   4.92
N2    1.99
N3    2.97
N4    3.95
EM    1.98
RK4   3.98
RK6_5 4.95
RK7_6 6.18


EI   1.04
CN   2.00
AM2  2.98
AM3  3.96
AM4  4.95
AM5  5.93
BDF2 1.97
BDF3 2.95
MS2  4.00


Heun    1.98
AM4_AB2 2.95
AM4_AB3 3.94
AM4_AB4 4.93
AM4_AB5 4.76

Si on note

  • $\omega_{[C]}$ l'ordre du corrector (schéma implicite)
  • $\omega_{[P]}$ l'ordre du predictor (schéma explicite)
    alors l'ordre du schéma predictor-corrector est $$\omega_{[PC]}=\min\{\omega_{[C]},\omega_{[P]}+1\}$$

Pour les schémas AM4-ABx, le schéma corrector est d'ordre $p=5$, pour ne pas perdre en ordre convergence il faut choisir un schéma predictor d'ordre $p-1$ (ici AB4).

Pour afficher l'ordre de convergence on utilise une échelle logarithmique : on représente $\ln(h)$ sur l'axe des abscisses et $\ln(\text{err})$ sur l'axe des ordonnées. Le but de cette représentation est clair: si $\text{err}=Ch^p$ alors $\ln(\text{err})=\ln(C)+p\ln(h)$. En échelle logarithmique, $p$ représente donc la pente de la ligne droite $\ln(\text{err})$.

In [38]:
figure(figsize=(24,7))

subplot(1,3,1)
loglog(H,err_ep,    'r-o',label=f'AB-1 = EE  {polyfit(log(H),log(err_ep),    1)[0]:1.2f}')
loglog(H,err_AB2,   'g-+',label=f'AB2   {polyfit(log(H),log(err_AB2),   1)[0]:1.2f}')
loglog(H,err_AB3,   'c-D',label=f'AB3   {polyfit(log(H),log(err_AB3),   1)[0]:1.2f}')
loglog(H,err_AB4,   'y-*',label=f'AB4   {polyfit(log(H),log(err_AB4),   1)[0]:1.2f}')
loglog(H,err_AB5,   'r-.',label=f'AB5   {polyfit(log(H),log(err_AB5),   1)[0]:1.2f}')
loglog(H,err_N2,    'y-*',label=f'N2    {polyfit(log(H),log(err_N2),    1)[0]:1.2f}')
loglog(H,err_N3,    'r-<',label=f'N3    {polyfit(log(H),log(err_N3),    1)[0]:1.2f}')
loglog(H,err_N4,    'b-+',label=f'N4    {polyfit(log(H),log(err_N4),    1)[0]:1.2f}')
loglog(H,err_em,    'c-o',label=f'EM    {polyfit(log(H),log(err_em),    1)[0]:1.2f}')
loglog(H,err_RK4,   'g-o',label=f'RK4_1 {polyfit(log(H),log(err_RK4),   1)[0]:1.2f}')
loglog(H,err_RK6_5, 'b-*',label=f'RK6_5 {polyfit(log(H),log(err_RK6_5), 1)[0]:1.2f}')
loglog(H,err_RK7_6, 'r-D',label=f'RK7_6 {polyfit(log(H),log(err_RK7_6), 1)[0]:1.2f}')
xlabel('$h$')
ylabel('$e$')
title("Schemas explicites")
legend(loc='upper center', bbox_to_anchor=(0.5, -0.05), fancybox=True, shadow=True, ncol=1)
grid(True)

subplot(1,3,2)
loglog(H,err_er,  'r-o',label=f'AM-0 = EI   {polyfit(log(H),log(err_er),   1)[0]:1.2f}')
loglog(H,err_CN,  'b-v',label=f'AM-1 = CN   {polyfit(log(H),log(err_CN),   1)[0]:1.2f}')
loglog(H,err_AM2, 'm->',label=f'AM2  {polyfit(log(H),log(err_AM2),  1)[0]:1.2f}')
loglog(H,err_AM3, 'c-D',label=f'AM3  {polyfit(log(H),log(err_AM3),  1)[0]:1.2f}')
loglog(H,err_AM4, 'g-+',label=f'AM4  {polyfit(log(H),log(err_AM4),  1)[0]:1.2f}')
loglog(H,err_AM5, 'b->',label=f'AM5  {polyfit(log(H),log(err_AM5),  1)[0]:1.2f}')
loglog(H,err_BDF2,'y-*',label=f'BDF2 {polyfit(log(H),log(err_BDF2), 1)[0]:1.2f}')
loglog(H,err_BDF3,'r-<',label=f'BDF3 {polyfit(log(H),log(err_BDF3), 1)[0]:1.2f}')
loglog(H,err_MS2, 'm-o',label=f'MS2  {polyfit(log(H),log(err_MS2),  1)[0]:1.2f}')
xlabel('$h$')
ylabel('$e$')
title("Schemas implicites")
legend(loc='upper center', bbox_to_anchor=(0.5, -0.05), fancybox=True, shadow=True, ncol=1)
grid(True)

subplot(1,3,3)
loglog(H,err_heun,   'y->',label=f'Heun    {polyfit(log(H),log(err_heun),   1)[0]:1.2f}')
loglog(H,err_AM4AB2, 'r-<',label=f'AM4_AB2 {polyfit(log(H),log(err_AM4AB2), 1)[0]:1.2f}')
loglog(H,err_AM4AB3, 'b-v',label=f'AM4_AB3 {polyfit(log(H),log(err_AM4AB3), 1)[0]:1.2f}' )
loglog(H,err_AM4AB4, 'm-*',label=f'AM4_AB4 {polyfit(log(H),log(err_AM4AB4), 1)[0]:1.2f}')
loglog(H,err_AM4AB5, 'g-+',label=f'AM4_AB5 {polyfit(log(H),log(err_AM4AB5), 1)[0]:1.2f}' )
xlabel('$h$')
ylabel('$e$')
title("Schemas predicteur-correcteur")
legend(loc='upper center', bbox_to_anchor=(0.5, -0.05), fancybox=True, shadow=True, ncol=1)
grid(True)

show()

Version compacte

De manière compacte en utilisant une liste des noms des schémas et des dictionnaires:

In [41]:
schemasEXPL = ['EE','AB2','AB3','AB4','AB5','N2','N3','N4', 'EM', 'RK4', 'RK6_5', 'RK7_6']
schemasIMPL = ['EI', 'CN', 'AM2', 'AM3', 'AM4', 'AM5', 'BDF2', 'BDF3', 'MS2', 'RK1_M']
schemasPRCO = ['heun', 'AM4AB2', 'AM4AB3', 'AM4AB4', 'AM4AB5']

schemas = schemasEXPL+schemasIMPL+schemasPRCO

H = []
err = { schemas[s] : [] for s in range(len(schemas)) }

N=10
for k in range(6):
    N += 50
    tt = linspace(t0,tfinal,N+1)
    h = tt[1]-tt[0]
    H.append(h)
    yy = array([sol_exacte(t) for t in tt])
    uu = { s : eval(s)(phi,tt,sol_exacte) for s in schemas }
    for key in uu:
        err[key].append(norm(uu[key]-yy,inf))

ordre = {  key : polyfit(log(H),log(err[key]),1)[0] for key in err  }

    
figure(figsize=(24,7))
markers=['^', 's', 'p', 'h', '8', 'D', '>', '<', '*','+','o', 'x']

subplot(1,3,1)
for idx,s in enumerate(schemasEXPL):
    loglog(H,err[s], label=f'{s}: {ordre[s]:1.2f}',marker=markers[idx])
xlabel('$h$')
ylabel('$e$')
title("Schemas explicites")
legend(loc='upper center', bbox_to_anchor=(0.5, -0.05), fancybox=True, shadow=True, ncol=1)
grid(True)

subplot(1,3,2)
for idx,s in enumerate(schemasIMPL):
    loglog(H,err[s], label=f'{s}: {ordre[s]:1.2f}',marker=markers[idx])
xlabel('$h$')
ylabel('$e$')
title("Schemas implicites")
legend(loc='upper center', bbox_to_anchor=(0.5, -0.05), fancybox=True, shadow=True, ncol=1)
grid(True)

subplot(1,3,3)
for idx,s in enumerate(schemasPRCO):
    loglog(H,err[s], label=f'{s}: {ordre[s]:1.2f}',marker=markers[idx])
xlabel('$h$')
ylabel('$e$')
title("Schemas predicteur-correcteur")
legend(loc='upper center', bbox_to_anchor=(0.5, -0.05), fancybox=True, ncol=1)
grid(True);
In [ ]: