TP 4 - Étude de la convergence

Auteur·rice

Gloria FACCANONI

Date de publication

8 avril 2024

Code

%reset -f
%matplotlib inline
%autosave 300

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

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

from scipy.optimize import fsolve

import sympy as sym
sym.init_printing()

Autosaving every 300 seconds
Python version 3.10.12 (main, Nov 20 2023, 15:14:05) [GCC 11.4.0]

1 Exercice

Considérons le problème de Cauchy :

trouver la fonction \(y \colon I\subset \mathbb{R} \to \mathbb{R}\) telle que

\( \begin{cases} y'(t) = y(t), &\forall t \in I=[0,1], \\ y(0) = 1. \end{cases} \)

On sait que la solution de ce problème de Cauchy est tout simplement \(y(t)=e^{t}\).

Compléter l’implémentation et estimer l’ordre de convergence des schémas indiqués.

Résultats attendus:

Les méthodes multipas explicites \(\text{AB}_q\) pas et \(\text{N}_q\) à \(q\) convergent à l’ordre \(q\).
Les méthodes multipas implicites \(\text{AM}_q\) et \(\text{MS}_q\) à \(q\) pas convergent à l’ordre \(q+1\) si \(q\) est impair et \(q+2\) si \(q\) est pair.
Les méthodes multipas implicites \(\text{BDF}_q\) à \(q\) pas convergent à l’ordre \(q\).
Les méthodes predictor-corrector où le schéma predictor (explicite) a ordre \(\omega_E\) et le corrector (implicite) a ordre \(\omega_I\) convergent à l’ordre \(\min(\omega_E+1,\omega_I)\).
Les méthodes Runge-Kutta à \(s\) étages ont ordre au mieux \(s\) si elles sont explicites et \(2s\) si elles sont implicites (bornes supérieures).

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

2 Schémas explicites

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

Code

def EE(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    # initialisation
    uu[0] = sol_exacte(tt[0])
    # recurrence
    for n in range(len(tt)-1):
        uu[n+1] = uu[n]+h*phi(tt[n],uu[n])
    return uu

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

Code

def AB2(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    # initialisation
    q = 2
    uu[:q] = sol_exacte(tt[:q])
    # recurrence
    for n in range(q-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

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

Code

def AB3(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    # initialisation
    q = 3
    uu[:q] = sol_exacte(tt[:q])
    # recurrence
    for n in range(q-1,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

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

Code

def AB4(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    # initialisation
    q = 4
    uu[:q] = sol_exacte(tt[:q])
    # recurrence
    for n in range(q-1,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

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

Code

def AB5(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    # initialisation
    q = 5
    uu[:q] = sol_exacte(tt[:q])
    # recurrence
    for n in range(q-1,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

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

Code

def N2(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    # initialisation
    q = 2
    uu[:q] = sol_exacte(tt[:q])
    # recurrence
    for n in range(q-1,len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        uu[n+1] = uu[n-1] + 2*h*k1
    return uu

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

Code

def N3(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    # initialisation
    q = 3
    uu[:q] = sol_exacte(tt[:q])
    # recurrence
    for n in range(q-1,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

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

Code

def N4(phi,tt,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    # initialisation
    q = 4
    uu[:q] = sol_exacte(tt[:q])
    # recurrence
    for n in range(q-1,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

2.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,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    # initialisation
    uu[0] = sol_exacte(tt[0])
    # recurrence
    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

2.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,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    # initialisation
    uu[0] = sol_exacte(tt[0])
    # 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 , 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

2.11 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,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    # initialisation
    uu[0] = sol_exacte(tt[0])
    # 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

2.12 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,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    # initialisation
    uu[0] = sol_exacte(tt[0])
    # 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

3 Schémas implicites

3.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,sol_exacte):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    # initialisation
    uu[0] = sol_exacte(tt[0])
    # recurrence
    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

3.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} \)