TD2 - Exercices sur les schémas multi-pas linéaires et les schémas prédicteur-correcteur

Auteur·rice

Gloria Faccanoni

Date de publication

28 janvier 2026

PARTIE I : Schémas 👣 multi-pas 👣 linéaires

Aide-mémoire

Définitions

Soit $p\ge0$ un entier et notons $\varphi_{n-j}\stackrel{\text{déf}}{=}\varphi(t_{n-j},u_{n-j})$.

Une méthode multi-pas linéaire à $q=p+1\ge1$ étapes (=pas) s’écrit

$ u_{n+1} = \sum_{j=0}^p \Big( a_ju_{n-j} + h b_j\varphi_{n-j} \Big) + hb_{-1}\varphi(t_{n+1},u_{n+1}), \qquad n = p,p+1,\dots,N-1 $

où les $\{a_k\}_{j=0}^p$ et $\{b_k\}_{j=-1}^p$ sont des coefficients donnés.

En écrivant explicitement la série, on a

$ u_{n+1} = a_0u_{n}+a_1u_{n-1}+\cdots+a_pu_{n-p} + h b_0\varphi_{n}+h b_1\varphi_{n-1}+\cdots+h b_{p}\varphi_{n-p} + hb_{-1}\varphi_{n+1}. $

Théorème de Lax-Richtmyer

Soit $\omega\ge1$. Le schéma est convergent (d’ordre $\omega$) ssi il est consistant (d’ordre $\omega$) et zéro-stable.

Zéro-stabilité

Pour étudier la zéro-stabilité, on introduit le premier polynôme caractéristique $ \varrho(r) = r^{p+1} - \sum_{j=0}^p a_jr^{p-j}. $

La méthode est zéro-stable ssi $ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ r_j\text{ a multiplicité}>1\quad\implies \quad |r_j|<1, \end{cases} $

autrement dit ssi les racines de $\varrho$ sont dans le cercle unité complexe et, si une racine est de multiplicité supérieure à 1, alors elle est strictement à l’intérieur du cercle unité.

Consistance (et ordre de consistance)

Pour étudier la consistance (et son ordre de consistance), introduisons la quantité

$ \xi(i) = \sum_{j=0}^p (-j)^{i}a_j+i\sum_{j=-1}^p (-j)^{i-1}b_j, \qquad i=0,1,\dots $

ayant posé $(-j)^0=1$ même pour $j=0$.

La méthode est consistante (d’ordre au moins 1) ssi $\xi(0)=\xi(1)=1$.
La méthode est consistante d’ordre $\omega$ ssi $\xi(0)=\xi(1)=\dots=\xi(\omega)=1$.

Première barrière de Dahlquist

L’ordre $\omega$ d’une méthode multi-pas à $q$ étapes convergente satisfait les inégalités suivantes :

$ \omega \le \begin{cases} q & \text{si $b_{-1}=0$ (la méthode est explicite),} \\ q+1 & \text{si $b_{-1}\neq0$ (la méthode est implicite) et }q\text{ est impair,} \\ q+2 & \text{si $b_{-1}\neq0$ (la méthode est implicite) et }q\text{ est pair.} \end{cases} $

A-stabilité

Soit $\beta>0$ un nombre réel positif et considérons le problème de Cauchy

$\begin{cases} y'(t)=-\beta y(t), &\text{pour }t>0,\\ y(0)=1 \end{cases}$

Sa solution est $y(t)=e^{-\beta t}$ et l’on a $y(t)\xrightarrow[t\to+\infty]{}0.$

Si, sous d’éventuelles conditions sur $h$, on a $|u_n|\xrightarrow[n\to+\infty]{}0$ alors on dit que le schéma est A-stable.

Le calcul des quantités $\xi(i)$ n’est pas difficile, mais il est fastidieux. On peut le faire à la main pour des méthodes avec un petit nombre de pas, mais pour des méthodes ayant plusieurs pas il est préférable de le faire avec un module de calcul formel. Voici un exemple avec sympy :

Code

# =============================================================================
# FONCTION qui calcule xi(i) pour un schéma à q pas
# =============================================================================

from sympy import symbols, Symbol, latex, factor, solve, Eq
from IPython.display import display, Math, Markdown, Latex

def xi(i, q):
    # Définition des symboles pour le cas général
    aa = symbols(f'a_0:{q}')
    bb = symbols(f'b_0:{q}')
    bm1 = Symbol('b_{-1}')
    p = len(aa)
    sa = sum((-j)**i * aa[j] for j in range(p))
    sb = bm1 + sum((-j)**(i-1) * bb[j] for j in range(1, p))
    # Cas particulier : si j=0 et i=0, il calcule (-j)**i =1 et on a bien a_0 dans la somme
    # mais si j=0 et i=1 il ne calcule pas b_0 et il faut l'ajouter à la main
    if i == 1:
        sb += bb[0]
    return (sa).factor() + (i * sb).factor()




# =============================================================================
# On doit seulement définir p
# =============================================================================

p = 2     # p+1 coeffs du schema d'indices 0 à p

# =============================================================================
# EXEMPLE D'APPLICATION 1 : calcul de xi(i) pour i = 0, ..., omega
# =============================================================================

# q et omega sont déduits de p, pour omega on considère le cas d'un schéma implicite
q = p + 1
omega = q + 2 if q%2==0 else q + 1 # i = 0, ..., omega

# Affichage des xi(i) séparément
display(Markdown(f"- Pour une méthode à $q={p + 1}$ pas, on affiche $\\xi(i)$ pour $i = 0, \\dots, {omega}$"))
for i in range(omega + 1):
    display(Markdown("$\\qquad$ "+f"$\\xi({i}) = {latex(xi(i, q))}$"))


# Affichage sous forme de système 
systeme = [Eq(Symbol(f"\\xi({i})") , Eq(xi(i, q),1)) for i in range(omega + 1)]
display(Markdown(r"\(\begin{cases}"+",\\\\ ".join([latex(eq) for eq in systeme])+r"\end{cases}\)"))
    

# =================================================================================================================
# EXEMPLE D'APPLICATION 2 : on peut même résoudre le système {xi(i)=1}_{i=0}^{\omega} 
# pour trouver les coefficients qui vérifient la condition de consistance d'ordre le plus élevé
# =================================================================================================================

# display(Markdown(f"- On résout le système $\\xi(i) = 1$ pour $i = 0, \\dots, {omega}$"))
# systeme = [xi(i, q)-1 for i in range(omega + 1)]
# sol = solve(systeme)
# display(Markdown("$\\qquad$ "+"Les coefficients qui vérifient la condition de consistance d'ordre le plus élevé sont :"))
# for k, v in sol.items():
#     display(Markdown("$\\qquad$ "+f"${k} = {v}$"))

Pour une méthode à $q=3$ pas, on affiche $\xi(i)$ pour $i = 0, \dots, 4$

$\qquad$ $\xi(0) = a_{0} + a_{1} + a_{2}$

$\qquad$ $\xi(1) = - a_{1} - 2 a_{2} + b_{0} + b_{1} + b_{2} + b_{-1}$

$\qquad$ $\xi(2) = a_{1} + 4 a_{2} - 2 \left(b_{1} + 2 b_{2} - b_{-1}\right)$

$\qquad$ $\xi(3) = - a_{1} - 8 a_{2} + 3 \left(b_{1} + 4 b_{2} + b_{-1}\right)$

$\qquad$ $\xi(4) = a_{1} + 16 a_{2} - 4 \left(b_{1} + 8 b_{2} - b_{-1}\right)$

$\begin{cases}\xi(0) = a_{0} + a_{1} + a_{2} = 1,\\ \xi(1) = - a_{1} - 2 a_{2} + b_{0} + b_{1} + b_{2} + b_{-1} = 1,\\ \xi(2) = a_{1} + 4 a_{2} - 2 \left(b_{1} + 2 b_{2} - b_{-1}\right) = 1,\\ \xi(3) = - a_{1} - 8 a_{2} + 3 \left(b_{1} + 4 b_{2} + b_{-1}\right) = 1,\\ \xi(4) = a_{1} + 16 a_{2} - 4 \left(b_{1} + 8 b_{2} - b_{-1}\right) = 1\end{cases}$

Exercices

👣 Exercice : consistance d’un schéma à deux pas explicite

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. La méthode multipas suivante est-elle consistante? $ u_{n+1}=u_{n-1}+ \dfrac{h}{4}(3\varphi_{n}-\varphi_{n-1}) $

Correction

C’est une méthode à q=2 pas explicite :

$p=1$
$b_{-1}=0$
$a_0=0$ et $a_1=1$
$b_0=\frac{3}{4}$ et $b_1=-\frac{1}{4}$

La méthode n’est pas consistante car

$ \begin{cases} \xi(0)= a_{0} + a_{1} =1, \\ \xi(1)= - a_{1} + b_{0} + b_{1} + b_{-1} =-1+\frac{3}{4}-\frac{1}{4}\neq1 \end{cases} $

Voyons sur un exemple ce qui se passe lorsqu’on utilise un schéma non consistant. On compare la solution exacte et la solution approchée du problème de Cauchy $y'(t)=2(y(t)-\cos(t))-\sin(t)$, $y(0)=1$ sur l’intervalle $[0,2\pi]$ dont la solution exacte est $y(t) = cos(t)$. Comme visible sur la figure suivante, la solution numérique ne correspond pas à la solution exacte et on peut vérifier qu’un comportement similaire se produit également si on raffine la discrétisation. Il s’agit d’une situation typique de manque de consistance : l’application d’une méthode non consistante produit l’allure d’une autre fonction plutôt que celle reproduisant la solution exacte.

Code

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt

t0, y0 = 0,1 

# Solution exacte
t = sp.symbols('t')
y = sp.Function('y')
edo = sp.Eq(sp.diff(y(t),t) , 2*(y(t)-sp.cos(t))-sp.sin(t) )
# display(edo)
exacte_sp = sp.dsolve( edo , y(t),  ics={y(t0):y0})
# display(exacte_sp)


# Solution approchée
def schema(tt,phi,sol_exacte):
    uu = np.zeros_like(tt)
    q = 2
    uu[:q] = sol_exacte(tt[:q])
    h = tt[1]-tt[0]
    for n in range(q-1,len(tt)-1):
        k0 = phi(tt[n],uu[n])
        k1 = phi(tt[n-1],uu[n-1])
        uu[n+1] = uu[n-1] + h/4*(3*k0-k1)
    return uu


# Comparaison
sol_exacte = sp.lambdify(t, exacte_sp.rhs, 'numpy')
tfin = 2*np.pi
phi =  lambda t,y: 2*(y-np.cos(t))-np.sin(t)
N = 100
tt = np.linspace(t0,tfin,1+N)
yy = sol_exacte(tt)
uu = schema(tt,phi,sol_exacte)
plt.plot(tt,yy,label='exacte')
plt.plot(tt,uu,label='schema')
plt.legend();

👣 Exercice : convergence d’un schéma à 2 pas implicite

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. La méthode multipas suivante est-elle convergente? $ u_{n+1}=-4u_{n}+5u_{n-1}+ h(4\varphi_{n+1}+2\varphi_{n}) $

Correction

C’est une méthode à 2 pas implicite :

$p=1$
$b_{-1}=4$
$a_0=-4$ et $a_1=5$
$b_0=2$ et $b_1=0$

La méthode est consistante car

$ \begin{cases} \displaystyle \xi(0)=a_{0} + a_{1}=-4+5=1, \\ \displaystyle \xi(1)=- a_{1} + b_{0} + b_{1} + b_{-1}=-5+4+2+0=1. \end{cases} $

Le premier polynôme caractéristique est

$ \begin{aligned} \varrho(r) &=r^{p+1}-\sum_{j=0}^p a_jr^{n-j}=r^2-\big(a_0 r^1+a_1r^0\big) \\ &=r^2-(-4r+5)=r^2+4r-5=(r-1)(r+5) \end{aligned} $

dont les racines sont $r_0=1$ et $r_1=-5$. La méthode n’est pas zéro-stable car $|r_1|>1$, donc la méthode ne converge pas.

Voyons sur un exemple ce qui se passe lorsqu’on utilise un schéma consistant mais non zéro-stable. On compare la solution exacte et la solution approchée du problème de Cauchy $y'(t)=2(y(t)-\cos(t))-\sin(t)$, $y(0)=1$ sur l’intervalle $[0,2\pi]$ dont la solution exacte est $y(t) = cos(t)$.

La cosistance est une propriété de précision locale, tandis que la zéro-stabilité et la convergence sont des propriétés de précision globale. La zéro-stabilité assure que la solution numérique ne diverge pas lorsque $h$ tend vers $0$.

Code

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

t0, y0 = 0,1 

# Solution exacte
t = sp.symbols('t')
y = sp.Function('y')
edo = sp.Eq(sp.diff(y(t),t) , 2*(y(t)-sp.cos(t))-sp.sin(t) )
# display(edo)
exacte_sp = sp.dsolve( edo , y(t),  ics={y(t0):y0})
# display(exacte_sp)

# Solution approchée
def schema(tt,phi,sol_exacte):
    uu = np.zeros_like(tt)
    q = 2
    uu[:q] = sol_exacte(tt[:q])
    h = tt[1]-tt[0]
    for n in range(q-1,len(tt)-1):
        eq = lambda x : -x -4*uu[n]+5*uu[n-1] + h*(2*phi(tt[n],uu[n])+phi(tt[n+1],x))
        uu[n+1] =  fsolve(eq,uu[n])[0]
    return uu


# Comparaison
sol_exacte = sp.lambdify(t, exacte_sp.rhs, 'numpy')
tfin = 2*np.pi
phi =  lambda t,y: 2*(y-np.cos(t))-np.sin(t)
N = 20
tt = np.linspace(t0,tfin,1+N)
yy = sol_exacte(tt)
uu = schema(tt,phi,sol_exacte)
plt.plot(tt,yy,label='exacte')
plt.plot(tt,uu,label='schema')
plt.ylim(-3,3)
plt.legend();

👣 Exercice : convergence d’un schéma à 3 pas explicite

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. La méthode multipas suivante est-elle convergente? $ u_{n+1}=-u_{n}+u_{n-1}+u_{n-2}+ 4h\varphi_{n} $

Correction

C’est une méthode à 3 pas explicite :

$p=2$
$b_{-1}=0$
$a_0=-1$, $a_1=1$ et $a_2=1$
$b_0=4$, $b_1=0$ et $b_2=0$

La méthode est consistante car

$ \begin{cases} \xi(0) = a_{0} + a_{1} + a_{2} = -1+1+1=1, \\ \xi(1) = - a_{1} - 2 a_{2} + b_{0} + b_{1} + b_{2} + b_{-1} = -1-2+4+0+0+0=-3+4=1. \end{cases} $

Le premier polynôme caractéristique est

$ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{n-j}=r^3-\big(-r^2+r+1\big)=r^3+r^2-r-1=(r-1)(r+1)^2 $

dont les racines sont $r_0=1$ de multiplicité $2$ et $r_1=-1$. La méthode n’est pas zéro-stable, donc la méthode ne converge pas.

Code

import sympy as sym
sym.var('r')
rho = r**3+r**2-r-1
sym.factor(rho)

$\displaystyle \left(r - 1\right) \left(r + 1\right)^{2}$

👣 Exercice : schéma à deux pas d’ordre 4

La première barrière de Dahlquist affirme qu’un schéma implicite à $q=2$ pas consistante et zéro-stable peut être d’ordre $\omega=q+2=4$.

Construire un schéma implicite à $q=2$ étapes consistante et zéro-stable d’ordre $\omega=q+2=4$.

Correction

Un schéma implicite à $q=2$ étapes s’écrit

$ u_{n+1}=a_0u_n+a_1u_{n-1}+h\big(b_{-1}\varphi_{n+1}+b_0\varphi_n+b_1\varphi_{n-1}\big) $

La méthode est d’ordre 4 ssi

$ \begin{cases} \xi(0) = a_{0} + a_{1} = 1,\\ \xi(1) = - a_{1} + b_{0} + b_{1} + b_{-1} = 1,\\ \xi(2) = a_{1} - 2 \left(b_{1} - b_{-1}\right) = 1,\\ \xi(3) = - a_{1} + 3 \left(b_{1} + b_{-1}\right) = 1,\\ \xi(4) = a_{1} - 4 \left(b_{1} - b_{-1}\right) = 1 \end{cases} $

Sous forme matricielle on trouve

$ \begin{pmatrix} 1&0& 1& 0&0\\ 0&1&-1& 1&1\\ 0&0& 1&-2&2\\ 0&0&-1& 3&3\\ 0&0& 1&-4&4 \end{pmatrix} \begin{pmatrix} a_0\\b_0\\a_1\\b_1\\b_{-1} \end{pmatrix} {}= \begin{pmatrix} 1\\1\\1\\1\\1 \end{pmatrix} $

Code

import sympy as sym
sym.init_printing()

a_0 = sym.symbols('a_0')
a_1 = sym.symbols('a_1')
b_m1 = sym.symbols('b_{-1}')
b_0 = sym.symbols('b_0')
b_1 = sym.symbols('b_1')

# p=1 => q=2 et omega=4
eq1 = sym.Eq(a_0 + a_1, 1) # eq1 = a_0+a_1-1
eq2 = sym.Eq(-a_1 + b_m1 + b_0 + b_1, 1) # eq2 = -a_1+b_{-1}+b_0+b_1-1
eq3 = sym.Eq(a_1 + 2 * b_m1 - 2 * b_1, 1) # eq3 = a_1+2*b_{-1}-2*b_1-1
eq4 = sym.Eq(-a_1 + 3 * b_m1 + 3 * b_1, 1) # eq4 = -a_1+3*b_{-1}+3*b_1-1
eq5 = sym.Eq(a_1 + 4 * b_m1 - 4 * b_1, 1) # eq5 = a_1+4*b_{-1}-4*b_1-1
sol = sym.solve([eq1, eq2, eq3, eq4, eq5]) # résolution du système d'équations
display(sol)

# Résultat
u_np1 = sym.symbols('u_{n+1}')
u_n = sym.symbols('u_n')
u_nm1 = sym.symbols('u_{n-1}')
phi_np1 = sym.symbols(r'\varphi_{n+1}')
phi_n = sym.symbols(r'\varphi_n')
phi_nm1 = sym.symbols(r'\varphi_{n-1}')
h = sym.symbols('h')
eq = sym.Eq(u_np1, (a_0 * u_n + a_1 * u_nm1 + h * (b_m1 * phi_np1 + b_0 * phi_n + b_1 * phi_nm1)).subs(sol))
display(eq)

$\displaystyle \left\{ a_{0} : 0, \ a_{1} : 1, \ b_{0} : \frac{4}{3}, \ b_{1} : \frac{1}{3}, \ b_{-1} : \frac{1}{3}\right\}$

$\displaystyle u_{n+1} = h \left(\frac{4 \varphi_{n}}{3} + \frac{\varphi_{n+1}}{3} + \frac{\varphi_{n-1}}{3}\right) + u_{n-1}$

On reconnait la méthode $\text{MS}_2$.

Le premier polynôme caractéristique est

$ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}=r^2-a_0r-a_1=(r-1)(r+1) $

dont les racines sont $r_0=1$ et $r_1=-1$. La méthode est zéro-stable car

$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \varrho'(r_j)\neq0 \text{ si } |r_j|=1, \end{cases} $

👣 Exercice : convergence d’un schéma à 2 pas implicite

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. Soit la méthode multipas $ u_{n+1}=(1-\gamma)u_{n}+\gamma u_{n-1}+ \frac{h}{4}((\gamma+3)\varphi_{n+1}+(3\gamma+1)\varphi_{n-1}). $

Pour quelles valeurs de $\gamma$ est-elle consistante?
Pour quelles valeurs de $\gamma$ est-elle zéro-stable?
Pour quelles valeurs de $\gamma$ est-elle convergente?
Pour quelles valeurs de $\gamma$ est-elle convergente d’ordre $2$? et d’ordre $3$?
Pour $\gamma=\frac{1}{2}$ vérifier empiriquement l’ordre de convergence sur le problème de Cauchy $ \begin{cases} y'(t) = -4y(t)+t^2, &\forall t \in I=[0,1],\\ y(0) = 1 \end{cases} $

Correction

C’est une méthode à 2 pas implicite :

$p=1$
$b_{-1}=\frac{\gamma+3}{4}$
$a_0=1-\gamma$ et $a_1=\gamma$
$b_0=0$ et $b_1=\frac{3\gamma+1}{4}$

Consistance La méthode est consistante pour tout $\gamma$ car

$ \begin{cases} \xi(0) = a_{0} + a_{1}=1-\gamma+\gamma=1 \text{ pour tout }\gamma, \\ \xi(1) = - a_{1} + b_{0} + b_{1} + b_{-1}=-\gamma+\big(0+\frac{3\gamma+1}{4}+\frac{\gamma+3}{4}\big)=1 \text{ pour tout }\gamma \end{cases} $

Zéro-stabilité Le premier polynôme caractéristique est

$ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}=r^2-(1-\gamma)r-\gamma=(r-1)(r+\gamma) $

dont les racines sont $r_0=1$ et $r_1=-\gamma$. La méthode est zéro-stable ssi

$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \text{si } r_j\text{ a multiplicité}>1 \quad\implies \quad |r_j|<1 \end{cases} $

donc ssi $-1\le -\gamma\le1$ et $-\gamma\neq1$, i.e. ssi $-1< \gamma\le 1$.

Convergence La méthode est convergente ssi elle est consistante et zéro-stable donc ssi $\boxed{-1< \gamma\le 1}$.

Ordre de convergence La méthode est d’ordre au moins 2 pour toute valeur de $\gamma$ pour laquelle elle est convergente, car

$ \xi(2) = a_{1} - 2 \left(b_{1} - b_{-1}\right) = \gamma-2\left(\frac{3\gamma+1}{4}-\frac{\gamma+3}{4}\right)=1 \text{ pour tout }\gamma. $

Pour que la méthode soit d’ordre 3 il faut que $\xi(3)=1$. On a

$ \xi(3) = - a_{1} + 3 \left(b_{1} + b_{-1}\right) = -\gamma+3\left(\frac{3\gamma+1}{4}+\frac{\gamma+3}{4}\right)= 2\gamma+3 = 1 \quad\iff\quad \gamma=-1. $

Puisque cette valeur n’appartient pas au domaine de zéro-stabilité, on conclut que la méthode est convergente d’ordre $2$ pour tout $-1< \gamma\le 1$ et non convergente sinon.

Vérification empirique On définit la solution exacte (calculée en utilisant le module sympy) pour estimer les erreurs, puis on définit le schéma et enfin on calcule les erreurs pour différentes valeurs de $h$.

Code

%reset -f
%matplotlib inline

# variables globales
t0     = 0
tfinal = 1
y0     = 1

phi = lambda t,y : -4*y+t**2

Code

import sympy as sym
sym.init_printing()

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

solpar = sym.dsolve(edo, ics={y(t0):y0})
display(solpar)

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

$\displaystyle \frac{d}{d t} y{\left(t \right)} = t^{2} - 4 y{\left(t \right)}$

$\displaystyle y{\left(t \right)} = \frac{t^{2}}{4} - \frac{t}{8} + \frac{1}{32} + \frac{31 e^{- 4 t}}{32}$

On calcule la solution approchée pour différentes valeurs de $N$:

Code

import numpy as np
from scipy.optimize import fsolve

def multipas(phi, tt, yy, gamma):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    # initialisation
    q = 2
    uu[:q] = yy[:q]
    # recurrence
    for n in range(q-1,len(tt) - 1):
        temp = fsolve(lambda x: -x + (1-gamma)*uu[n]+gamma*uu[n-1] + h/4 * ( (gamma+3)*phi(tt[n+1],x)+(3*gamma+1)*phi(tt[n-1],uu[n-1]) ), uu[n])
        uu[n+1] = temp[0]
    return uu

Code

import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 12})

H = []
err_mp = []
gamma = 1/2

for N in range(10,100,20):
    tt = np.linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = sol_exacte(tt)
    uu_mp = multipas(phi, tt, yy, gamma)
    H.append(h)
    err_mp.append(  np.linalg.norm(uu_mp-yy, np.inf) )
         
         

plt.figure(figsize=(8,5))
plt.plot(np.log(H), np.log(err_mp), 'r-o')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.title(f'Ordre de convergence = pente de la droite = {(np.polyfit(np.log(H),np.log(err_mp), 1)[0])}')
plt.grid(True)

👣 Exercice : convergence d’un schéma à 2 pas explicite

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. Soit la méthode multipas $ u_{n+1}= \left(\frac{1}{3}\gamma^2+\frac{1}{2}\gamma\right)u_{n}+ \left(\frac{2}{3}\gamma^2+\frac{1}{2}\gamma-1\right)u_{n-1}+ h\left[\left(-\frac{1}{6}\gamma^2-\frac{1}{4}\gamma+2\right)\varphi_{n}+ \left(-\frac{1}{6}\gamma^2-\frac{1}{4}\gamma\right)\varphi_{n-1}\right]. $

Pour quelles valeurs de $\gamma$ est-elle consistante ?
Est-elle convergente pour ces valeurs ?
Quel ordre de convergente a-t-elle ?
Pour les valeurs de $\gamma$ qui donnent la convergence, vérifier empiriquement l’ordre de convergence sur le problème de Cauchy $ \begin{cases} y'(t) = -10\big(y(t)-1\big)^2, &\forall t \in I=[0,1],\\ y(0) = 2. \end{cases} $

Correction

C’est une méthode à $q=2$ pas explicite :

$p=1$
$b_{-1}=0$
$a_0=\frac{1}{3}\gamma^2+\frac{1}{2}\gamma$ et $a_1=\frac{2}{3}\gamma^2+\frac{1}{2}\gamma-1$
$b_0=-\frac{1}{6}\gamma^2-\frac{1}{4}\gamma+2$ et $b_1=-\frac{1}{6}\gamma^2-\frac{1}{4}\gamma$

Consistance
La méthode est consistante ssi

$ \begin{cases} \xi(0) = a_{0} + a_{1} =1, \\ \xi(1) = - a_{1} + b_{0} + b_{1} + b_{-1}=1 \end{cases} \quad\iff\quad \begin{cases} \gamma^2+\gamma-1=1, \\ -\gamma^2-\gamma+3=1 \end{cases} \quad\iff\quad \gamma^2+\gamma-2=0 $

donc ssi $\gamma=-2$ ou $\gamma=1$. On obtient alors les valeurs des coefficients suivantes :

$ \begin{array}{|c|c|c|c|c|c|} \hline \gamma & a_0 & a_1 & b_0 & b_1 & b_{-1}\\ \hline -2 & \frac{1}{3} & \frac{2}{3} & \frac{11}{6} & -\frac{1}{6} & 0\\ 1 & \frac{5}{6} & \frac{1}{6} & \frac{19}{12} & -\frac{5}{12} & 0\\ \hline \end{array} $

Zéro-stabilité
On étudie la zéro-stabilité seulement pour les valeurs de $\gamma$ qui donnent la consistance. Le premier polynôme caractéristique est

$ \varrho(r)= \begin{cases} r^2-\frac{1}{3}r-\frac{2}{3}&\text{si }\gamma=-2, \\ r^2-\frac{5}{6}r-\frac{1}{6}&\text{si }\gamma=1, \end{cases} $

dont les racines sont

$ \begin{cases} r_0=1,\ r_1=-\frac{2}{3}&\text{si }\gamma=-2, \\ r_0=1,\ r_1=-\frac{1}{6}&\text{si }\gamma=1. \end{cases} $

La méthode est zéro-stable ssi
$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \text{si } r_j\text{ a multiplicité}>1 \quad\implies \quad |r_j|<1 \end{cases} $

Cette condition est vérifiée pour les deux cas.

Convergence
La méthode est convergente ssi elle est consistante et zéro-stable donc ssi $\gamma=-2$ ou $\gamma=1$.

Ordre
La méthode est d’ordre au moins 2 pour les deux choix de $\gamma$ car

$ \xi(2) = a_{1} - 2 \left(b_{1} - b_{-1}\right) = \begin{cases} \frac{2}{3}+2\frac{1}{6}=1&\text{si }\gamma=-2, \\ \frac{1}{6}+2\frac{5}{12}=1&\text{si }\gamma=1, \end{cases} $

La méthode n’est d’ordre 3 pour aucun des deux choix de $\gamma$ car
$ \xi(3) = - a_{1} + 3 \left(b_{1} + b_{-1}\right) = \begin{cases} {}-\frac{2}{3}+3\frac{-1}{6}\neq1&\text{si }\gamma=-2, \\ {}-\frac{1}{6}+3\frac{-5}{12}\neq1&\text{si }\gamma=1, \end{cases} $

Solution exacte
On définit la solution exacte (calculée en utilisant le module sympy) pour estimer les erreurs:

Code

%reset -f
%matplotlib inline


import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 16})

from scipy.optimize import fsolve

import sympy as sym
sym.init_printing()

from IPython.display import display, Math

Code

t0     = 0
tfinal = 1
y0     = 2

phi = lambda t,y : -10*(y-1)**2

Code

t  = sym.Symbol('t')
y  = sym.Function('y')
edo= sym.Eq( sym.diff(y(t),t) , phi(t,y(t)) )
display(edo)
solpar = sym.dsolve(edo, ics={y(t0): y0})
display(solpar)

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

$\displaystyle \frac{d}{d t} y{\left(t \right)} = - 10 \left(y{\left(t \right)} - 1\right)^{2}$

$\displaystyle y{\left(t \right)} = \frac{- t - \frac{1}{5}}{- t - \frac{1}{10}}$

On calcule la solution approchée pour différentes valeurs de $N$:

Code

def multipas(gamma, phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[:2] = sol_exacte(tt[:2])
    for n in range(1,len(tt) - 1):
        uu[n+1] = ( (gamma**2/3+gamma/2)*uu[n] + (2*gamma**2/3+gamma/2-1)*uu[n-1] 
                  +h*(-gamma**2/6-gamma/4+2)*phi(tt[n],uu[n])+h*(-gamma**2/6-gamma/4)*phi(tt[n-1],uu[n-1]) )                 
    return uu


H = []
err_mp1 = []
err_mp2 = []

for N in range(20,120,20):
    tt = np.linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = sol_exacte(tt)
    uu_mp1 = multipas( 1, phi, tt, sol_exacte)
    uu_mp2 = multipas(-2, phi, tt, sol_exacte)
    H.append(h)
    err_mp1.append(max([abs(uu_mp1[n] - yy[n]) for n in range(len(yy))]))
    err_mp2.append(max([abs(uu_mp2[n] - yy[n]) for n in range(len(yy))]))

ordre_gamma_1 = np.polyfit(np.log(H),np.log(err_mp1), 1)[0]
ordre_gamma_m2 = np.polyfit(np.log(H),np.log(err_mp2), 1)[0]

plt.figure(figsize=(8,5))
plt.plot(np.log(H), np.log(err_mp1), 'r-o', label=f'Multipas b=1 : ordre={ordre_gamma_1:.2f}')
plt.plot(np.log(H), np.log(err_mp2), 'b-d', label=f'Multipas b=-2: ordre={ordre_gamma_m2:.2f}')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.legend(bbox_to_anchor=(1.04, 1), loc='upper left')
plt.grid(True)

👣 Exercice : convergence d’un schéma à 2 pas explicite

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. Soit la méthode multipas

$ u_{n+1} = \left(\gamma^2-\frac{9}{2}\gamma+\frac{11}{2}\right)u_{n} + \left(\gamma^2-\frac{7}{2}\gamma+\frac{3}{2}\right)u_{n-1} + h\left(\frac{\gamma}{2}(\gamma-2)\varphi_{n} + \frac{\gamma}{2}\left(\gamma-\frac{10}{3}\right)\varphi_{n-1}\right). $

Pour quelles valeurs de $\gamma$ est-elle consistante ?
Est-elle convergente pour ces valeurs ?
Quel ordre de convergente a-t-elle ?
Pour les valeurs de $\gamma$ qui donnent la convergence, vérifier empiriquement l’ordre de convergence sur le problème de Cauchy $ \begin{cases} y'(t) = y(t)(1-y(t)), &\forall t \in I=[0,1],\\ y(0) = 2. \end{cases} $

Correction

C’est une méthode à $q=2$ pas explicite :

$p=1$
$b_{-1}=0$
$a_0=\gamma^2-\frac{9}{2}\gamma+\frac{11}{2}$ et $a_1=\gamma^2-\frac{7}{2}\gamma+\frac{3}{2}$
$b_0=\frac{\gamma}{2}(\gamma-2)$ et $b_1=\frac{\gamma}{2}\left(\gamma-\frac{10}{3}\right)$

La méthode est consistante ssi $\xi(0)=\xi(1)=1$:

$ \begin{cases} \xi(0) = a_{0} + a_{1} =1, \\ \xi(1) = - a_{1} + b_{0} + b_{1} + b_{-1}=1 \end{cases} \quad\iff\quad \begin{cases} 2\gamma^2-8\gamma+7=1, \\ \frac{5}{6}\gamma-\frac{9}{6}=1 \end{cases} \quad\iff\quad \begin{cases} \gamma^2-4\gamma+3=0, \\ 5\gamma=15 \end{cases} \quad\iff\quad \begin{cases} (\gamma-1)(\gamma-3)=0, \\ \gamma=3. \end{cases} $

On conclut que la méthode est consistante ssi $\gamma=3$.

Vérifions ces calculs ci-dessous :

Code

%reset -f
%matplotlib inline


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

from IPython.display import display, Math

Code

gamma = sym.Symbol(r'\gamma')
p = 1
bm1 = 0
a0 = gamma**2-sym.Rational(9,2)*gamma+sym.Rational(11,2)
a1 = gamma**2-sym.Rational(7,2)*gamma+sym.Rational(3,2)
b0 = sym.Rational(1,2)*gamma*(gamma-2)
b1 = sym.Rational(1,2)*gamma*(gamma-sym.Rational(10,3))


display(Math(r"a_0 = "+sym.latex(a0)))
display(Math(r"a_1 = "+sym.latex(a1)))
display(Math(r"b_{-1} = "+sym.latex(bm1)))
display(Math(r"b_0 = "+sym.latex(b0)))
display(Math(r"b_1 = "+sym.latex(b1)))

cond1 = a0+a1
cond2 = -(0*a0+1*a1)+(bm1+b0+b1)
display(sym.solve([sym.Eq(cond1,1),sym.Eq(cond2,1)],gamma))

$\displaystyle a_0 = \gamma^{2} - \frac{9 \gamma}{2} + \frac{11}{2}$

$\displaystyle a_1 = \gamma^{2} - \frac{7 \gamma}{2} + \frac{3}{2}$

$\displaystyle b_{-1} = 0$

$\displaystyle b_0 = \frac{\gamma \left(\gamma - 2\right)}{2}$

$\displaystyle b_1 = \frac{\gamma \left(\gamma - \frac{10}{3}\right)}{2}$

$\displaystyle \left[ \left( 3,\right)\right]$

Dans la suite on considerera $\gamma=3$.

Code

display(Math('a_0='+sym.latex(a0.subs(gamma,3))))
display(Math('a_1='+sym.latex(a1.subs(gamma,3))))
display(Math('b_0='+sym.latex(b0.subs(gamma,3))))
display(Math('b_1='+sym.latex(b1.subs(gamma,3))))

$\displaystyle a_0=1$

$\displaystyle a_1=0$

$\displaystyle b_0=\frac{3}{2}$

$\displaystyle b_1=- \frac{1}{2}$

Le premier polynôme caractéristique est

$ \varrho(r) = r^2-a_0r-a_1 \stackrel{\gamma=3}{=} r^2-r $

dont les racines sont

$ r_0=1,\ r_1=0 $

La méthode est zéro-stable ssi

$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \text{si } r_j\text{ a multiplicité}>1 \quad\implies \quad |r_j|<1 \end{cases} $

ce qui est bien vérifié.

La méthode est convergente ssi elle est consistante et zéro-stable donc ssi $\gamma=3$.

La méthode est d’ordre au moins 2 car $\xi(2)=1$ pour $\gamma=3$ :

$ \xi(2) = a_{1} - 2 \left(b_{1} - b_{-1}\right) \stackrel{\gamma=3}{=} 1 $

La méthode n’est pas d’ordre 3 car $\xi(3)\neq1$ :

$ \xi(3) = - a_{1} + 3 \left(b_{1} + b_{-1}\right) \stackrel{\gamma=3}{\neq} 1 $

On définit la solution exacte (calculée en utilisant le module sympy) pour estimer les erreurs:

Code

t0     = 0
tfinal = 1
y0     = 2

phi = lambda t,y : y*(1-y)

Code

t  = sym.Symbol('t')
y  = sym.Function('y')
edo = sym.Eq( sym.diff(y(t),t) , phi(t,y(t)) )
display(edo)
solpar = sym.dsolve(edo, ics={y(t0): y0})
display(solpar)

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

$\displaystyle \frac{d}{d t} y{\left(t \right)} = \left(1 - y{\left(t \right)}\right) y{\left(t \right)}$

$\displaystyle y{\left(t \right)} = \frac{1}{1 - \frac{e^{- t}}{2}}$

On calcule la solution approchée pour différentes valeurs de $N$:

Code

def multipas(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[:2] = sol_exacte(tt[:2])
    for n in range(1,len(tt) - 1):
        uu[n+1] = uu[n] + h*(3/2*phi(tt[n],uu[n])-1/2*phi(tt[n-1],uu[n-1]))
    return uu


H = []
err_mp = []

for N in range(10,100,20):
    tt = np.linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = sol_exacte(tt)
    uu_mp = multipas(phi, tt, sol_exacte)
    H.append(h)
    err_mp.append(max([abs(uu_mp[i] - yy[i]) for i in range(len(yy))]))

plt.figure(figsize=(8,5))
plt.plot(np.log(H), np.log(err_mp), 'r-o', label=r'Multipas $\gamma=3$')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.legend(bbox_to_anchor=(1.04, 1), loc='upper left')
plt.title(f'Ordre = {np.polyfit(np.log(H),np.log(err_mp), 1)[0]:g}')
plt.grid(True)

PARTIE II : Schémas Prédicteur-Correcteur basés sur deux schémas multi-pas linéaires

Aide-mémoire

Principe de la méthode predictor-corrector (PC)

On choisi deux méthodes multi-pas linéaires :

une méthode explicite à $\tilde q = \tilde p+1$ pas (de coefficients $\{\tilde a_j\}$ et $\{\tilde b_j\}$) : $ u_{n+1} = \displaystyle\sum_{j=0}^{\tilde p} \Big(\tilde a_ju_{n-j} + h \tilde b_j\varphi_{n-j}\Big) \quad\text{ pour $n=p,p+1,\dots,N-1$} $
une méthode implicite à $q = p+1$ pas (de coefficients $\{a_j\}$ et $\{b_j\}$) : $ u_{n+1} = \displaystyle\sum_{j=0}^p \Big(a_ju_{n-j} + h b_j\varphi_{n-j}\Big) + hb_{-1}\varphi_{n+1} \quad\text{ pour $n=p,p+1,\dots,N-1$.} $

On construit un schéma prédicteur-correcteur (explicite) de la manière suivante :

$\begin{cases} {\color{blue}{\tilde u_{n+1}}} = \displaystyle\sum_{j=0}^{\tilde p} \Big(\tilde a_ju_{n-j} + h \tilde b_j\varphi_{n-j}\Big) \\ u_{n+1}=\displaystyle\sum_{j=0}^p \Big(a_ju_{n-j} + h b_j\varphi_{n-j} \Big) + hb_{-1}\varphi(t_{n+1},{\color{blue}{\tilde u_{n+1}}}) \end{cases}$

Ordre : si on note $\omega_{[P]}$ l’ordre du predictor (schéma explicite) et $\omega_{[C]}$ l’ordre du corrector (schéma implicite), alors l’ordre du schéma predictor-corrector est donné par

$ \omega_{[PC]}=\min\Big\{\omega_{[C]},\omega_{[P]}+1\Big\} $

Exercices

Exercice : convergence d’un schéma predictor-corrector basé sur deux schémas 👣

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. Soit les méthodes multipas $ \begin{aligned} u_{n+1}&=u_{n}+\frac{h}{12}\left(23\varphi_{n}-16\varphi_{n-1}+5\varphi_{n-2}\right)&\text{[P]} \\ u_{n+1}&=u_{n}+\frac{h}{24}\left(9\varphi_{n+1}+19\varphi_{n}-5\varphi_{n-1}+\varphi_{n-2}\right)&\text{[C]} \end{aligned} $

Étudier consistance, ordre et zéro-stabilité de la méthode P
Étudier consistance, ordre et zéro-stabilité de la méthode C
Soit la méthode PC suivante: $ \begin{cases} \tilde u=u_{n}+\frac{h}{12}\left(23\varphi_{n}-16\varphi_{n-1}+5\varphi_{n-2}\right) \\ u_{n+1}=u_{n}+\frac{h}{24}\left(9\varphi(t_{n+1},\tilde u)+19\varphi_{n}-5\varphi_{n-1}+\varphi_{n-2}\right) \end{cases} $ Calculer empiriquement l’ordre de convergence de la méthode PC sur le problème de Cauchy $ \begin{cases} y'(t) = \frac{1}{1+t^2}-2y^2(t), &\forall t \in I=[0,1],\\ y(0) = 0, \end{cases} $ dont la solution exacte est $y(t)=\frac{t}{1+t^2}$.

Correction schéma P

P est une méthode à $q=3$ pas explicite : - $p=2$ - $b_{-1}=0$ - $a_0=1$ et $a_1=a_2=0$ - $b_0=\frac{23}{12}$, $b_1=\frac{-16}{12}$ et $b_2=\frac{5}{12}$ - La méthode est consistante d’ordre $\omega=3$ car $ \begin{cases} \displaystyle\sum_{j=0}^p a_j=1, \\ \displaystyle\sum_{j=0}^p -ja_j+\sum_{j=-1}^p b_j=-a_1-2a_2+\big(b_{-1}+b_0+b_1+b_2\big)=1 \\ \displaystyle\sum_{j=0}^p (-j)^{2}a_j+2\sum_{j=-1}^p (-j)^{1}b_j=a_1+4a_2-2\big(-b_{-1}+b_1+2b_2\big)=1 \\ \displaystyle\sum_{j=0}^p (-j)^{3}a_j+3\sum_{j=-1}^p (-j)^{2}b_j=-a_1-8a_2+3\big(b_{-1}+b_1+4b_2\big)=1 \\ \displaystyle\sum_{j=0}^p (-j)^{4}a_j+4\sum_{j=-1}^p (-j)^{3}b_j=a_1+16a_2-4\big(-b_{-1}+b_1+8b_2\big)\neq1 \end{cases} $ Vérifions ces calculs ci-dessous:

Code

import sympy as sym
sym.init_printing()
from IPython.display import display, Math
 
p = 1
bm1 = 0
a0 = 1
a1 = 0
a2 = 0
b0 = sym.Rational(23,12)
b1 = sym.Rational(-16,12)
b2 = sym.Rational(5,12)

aa = [a0, a1, a2]
bb = [b0, b1, b2]

cond = []
cond.append(sum(aa)) 
cond.append(sum([-j*aa[j]+bb[j] for j in range(len(aa))])+bm1)
if cond==[1,1]:
    display(Math(r"\omega\ge1"))
    for n in range(2,9):
        cond.append(sum( [ (-j)**n*aa[j]+n*(-j)**(n-1)*bb[j] for j in range(len(aa)) ])+n*bm1)
        if cond[n]==1:
            display(Math(r"\omega\ge"+str(n)))
        else:
            display(Math(r"\omega\le"+str(n)))
            break

$\displaystyle \omega\ge1$

$\displaystyle \omega\ge2$

$\displaystyle \omega\ge3$

$\displaystyle \omega\le4$

Le premier polynôme caractéristique est $ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}=r^3-a_0r^2-a_1r-a_2 {}= r^3-r^2=r^2(r-1) $ dont les racines sont $ r_0=1,\ r_1=r_2=0 $ La méthode est zéro-stable ssi $ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \varrho'(r_j)\neq0 \text{ si } |r_j|=1, \end{cases} $ ce qui est bien vérifié.

La méthode est convergente car consistante et zéro-stable.

Correction schéma C

C est une méthode à $q=3$ pas implicite :

$p=2$
$b_{-1}=\frac{9}{24}$
$a_0=1$ et $a_1=a_2=0$
$b_0=\frac{19}{24}$, $b_1=\frac{-5}{24}$ et $b_2=\frac{1}{24}$
La méthode est consistante d’ordre $\omega=4$ car $ \begin{cases} \displaystyle\sum_{j=0}^p a_j=1, \\ \displaystyle\sum_{j=0}^p -ja_j+\sum_{j=-1}^p b_j=-a_1-2a_2+\big(b_{-1}+b_0+b_1+b_2\big)=1 \\ \displaystyle\sum_{j=0}^p (-j)^{2}a_j+2\sum_{j=-1}^p (-j)^{1}b_j=a_1+4a_2-2\big(-b_{-1}+b_1+2b_2\big)=1 \\ \displaystyle\sum_{j=0}^p (-j)^{3}a_j+3\sum_{j=-1}^p (-j)^{2}b_j=-a_1-8a_2+3\big(b_{-1}+b_1+4b_2\big)=1 \\ \displaystyle\sum_{j=0}^p (-j)^{4}a_j+4\sum_{j=-1}^p (-j)^{3}b_j=a_1+16a_2-4\big(-b_{-1}+b_1+8b_2\big)=1 \\ \displaystyle\sum_{j=0}^p (-j)^{5}a_j+5\sum_{j=-1}^p (-j)^{4}b_j=-a_1+32a_2+5\big(b_{-1}+b_1+16b_2\big)\neq1 \end{cases} $ Vérifions ces calculs ci-dessous:

Code

import sympy as sym
sym.init_printing()
from IPython.display import display, Math
 
p = 1
bm1 = sym.Rational(9,24)
a0 = 1
a1 = 0
a2 = 0
b0 = sym.Rational(19,24)
b1 = sym.Rational(-5,24)
b2 = sym.Rational(1,24)

aa = [a0,a1,a2]
bb = [b0,b1,b2]

cond = []
cond.append(sum(aa)) 
cond.append(sum([-j*aa[j]+bb[j] for j in range(len(aa))])+bm1)
if cond==[1,1]:
    display(Math(r"\omega\ge1"))
    for n in range(2,9):
        cond.append(sum( [ (-j)**n*aa[j]+n*(-j)**(n-1)*bb[j] for j in range(len(aa)) ])+n*bm1)
        if cond[n]==1:
            display(Math(r"\omega\ge"+str(n)))
        else:
            display(Math(r"\omega\le"+str(n)))
            break

$\displaystyle \omega\ge1$

$\displaystyle \omega\ge2$

$\displaystyle \omega\ge3$

$\displaystyle \omega\ge4$

$\displaystyle \omega\le5$

Le premier polynôme caractéristique est $ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}=r^3-a_0r^2-a_1r-a_2 {}= r^3-r^2=r^2(r-1) $ dont les racines sont $ r_0=1,\ r_1=r_2=0 $ La méthode est zéro-stable ssi $ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \varrho'(r_j)\neq0 \text{ si } |r_j|=1, \end{cases} $ ce qui est bien vérifié.

La méthode est convergente car consistante et zéro-stable.

Correction schéma PC

Il est d’ordre $4$ car le prédictor est d’ordre $3$ et le corrector d’ordre $4$.

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

Code

%reset -f
%matplotlib inline

# variables globales
t0     = 0
tfinal = 1
y0     = 0

phi = lambda t,y : 1/(1+t**2)-2*y**2

sol_exacte = lambda t : t/(1+t**2)

On calcule la solution approchée pour différentes valeurs de $N$ pour les trois schémas:

Code

import numpy as np
from scipy.optimize import fsolve

def multipas_P(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[:3] = sol_exacte(tt[:3])
    for n in range(2,len(tt) - 1):
        k_0 = phi(tt[n],uu[n])
        k_1 = phi(tt[n-1],uu[n-1])
        k_2 = phi(tt[n-2],uu[n-2])
        uu[n+1] = uu[n]+h/12* ( 23*k_0-16*k_1+5*k_2 )
    return uu

def multipas_C(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[:3] = sol_exacte(tt[:3])
    for n in range(2,len(tt) - 1):
        k_0 = phi(tt[n],uu[n])
        k_1 = phi(tt[n-1],uu[n-1])
        k_2 = phi(tt[n-2],uu[n-2])
        temp = fsolve ( lambda x : -x+uu[n]+h/24* ( 9*phi(tt[n+1],x) + 19*k_0-5*k_1+k_2 ), uu[n])
        uu[n+1] = temp[0]
    return uu

def multipas_PC(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[:3] = sol_exacte(tt[:3])
    for n in range(2,len(tt) - 1):
        k_0 = phi(tt[n],uu[n])
        k_1 = phi(tt[n-1],uu[n-1])
        k_2 = phi(tt[n-2],uu[n-2])
        u_tilde = uu[n]+h/12 * (23*k_0-16*k_1+5*k_2)
        uu[n+1] = uu[n]+h/24* ( 9*phi(tt[n+1],u_tilde) + 19*k_0-5*k_1+k_2 )
    return uu

Code

import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 16})


H = []
err_p = []
err_c = []
err_pc = []
N = 10
for k in range(7):
    N += 20
    tt = np.linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = sol_exacte(tt)
    uu_p = multipas_P(phi, tt, sol_exacte)
    uu_c = multipas_C(phi, tt, sol_exacte)
    uu_pc = multipas_PC(phi, tt, sol_exacte)
    H.append(h)
    err_p.append(np.linalg.norm(uu_p-yy, np.inf))
    err_c.append(np.linalg.norm(uu_c-yy, np.inf))
    err_pc.append(np.linalg.norm(uu_pc-yy, np.inf))
    
print (f'Multipas P ordre  = {np.polyfit(np.log(H),np.log(err_p) , 1)[0]:1.2f}' )
print (f'Multipas C ordre  = {np.polyfit(np.log(H),np.log(err_c) , 1)[0]:1.2f}' )
print (f'Multipas PC ordre = {np.polyfit(np.log(H),np.log(err_pc), 1)[0]:1.2f}')

plt.figure(figsize=(8,5))
plt.plot(np.log(H), np.log(err_p), 'r-o', label='Multipas P')
plt.plot(np.log(H), np.log(err_c), 'b-o', label='Multipas C')
plt.plot(np.log(H), np.log(err_pc), 'c-o', label='Multipas PC')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.legend(bbox_to_anchor=(1.04, 1), loc='upper left')
plt.grid(True);

Multipas P ordre  = 2.99
Multipas C ordre  = 3.86
Multipas PC ordre = 4.00

Exercice : convergence d’un schéma predictor-corrector basé sur deux schémas 👣

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. Soit les méthodes multipas $ \begin{aligned} u_{n+1}&=\frac{4}{3}u_{n}-\frac{1}{3}u_{n-1}+h\left(b_0\varphi_{n}+b_1\varphi_{n-1}\right) &\text{[P]} \\ u_{n+1}&=\frac{4}{3}u_{n}-\frac{1}{3}u_{n-1}+h\frac{2}{3}\varphi_{n+1} &\text{[C]} \end{aligned} $

Étudier consistance, ordre et zéro-stabilité des deux méthodes.
Écrire et étudier consistance, ordre et zéro-stabilité de la méthode prédictor-corrector associée.
Calculer empiriquement l’ordre de convergence des trois méthodes sur le problème de Cauchy $ \begin{cases} y'(t) = \frac{1}{1+t^2}-2y^2(t), &\forall t \in I=[0,1],\\ y(0) = 0, \end{cases} $ dont la solution exacte est $y(t)=\frac{t}{1+t^2}$.

Correction

Les deux schémas multipas linéaires à $q=2$ pas sont de la forme $ u_{n+1} = a_0 u_n + a_1 u_{n-1} + h \left( b_0 \varphi_n + b_1 \varphi_{n-1} \right) + h b_{-1} \varphi_{n+1} $ avec $ \begin{array}{c|cccccc} & a_0 & a_1 & b_0 & b_1 & b_{-1} \\ \hline \text{[P]} & \frac{4}{3} & -\frac{1}{3} & b_0 & b_1 & 0 \\ \text{[C]} & \frac{4}{3} & -\frac{1}{3} & 0 & 0 & \frac{2}{3} \end{array} $

Les deux schémas ont le même premier polynôme caractéristique : $ \varrho(r)=r^2-a_0r-a_1=(r-1)\left(r-\frac{1}{3}\right). $ Ses racines sont $r_0=1$ et $r_1=\frac{1}{3}$. Comme elles appartiennent au disque unité et sont simples, les deux méthodes sont zéro-stables.

D’après la première barrière de Dalquist, un schéma linéaire à $q=2$ pas consistant et zéro-stable est au mieux

d’ordre $\omega = q+2 = 4$ s’il est implicite,
d’ordre $\omega = q = 2$ s’il est explicite.

Étudions l’ordre de convergence en calculant $\xi(i)$ pour $i=0,\dots,4$ (les deux méthodes ayant le même nombre de pas, on a la même expression pour $\xi(i)$) : $ \begin{cases} \xi(0) = a_{0} + a_{1} \\ \xi(1) = - a_{1} + b_{0} + b_{1} + b_{-1} \\ \xi(2) = a_{1} - 2 \left(b_{1} - b_{-1}\right) \\ \xi(3) = - a_{1} + 3 \left(b_{1} + b_{-1}\right) \\ \xi(4) = a_{1} - 4 \left(b_{1} - b_{-1}\right) \end{cases} $

Pour le schéma [C], on a $\xi(0)=\xi(1)=\xi(2)=1$ et $\xi(3)\neq1$, donc la méthode est d’ordre 2 : $ \begin{cases} \xi(0) = a_{0} + a_{1} = \frac{4}{3} -\frac{1}{3}=1\\ \xi(1) = - a_{1} + b_{0} + b_{1} + b_{-1} = \frac{1}{3}+\frac{2}{3}=1\\ \xi(2) = a_{1} - 2 \left(b_{1} - b_{-1}\right) = -\frac{1}{3}+2\frac{2}{3}=1\\ \xi(3) = - a_{1} + 3 \left(b_{1} + b_{-1}\right) = \frac{1}{3}+3\frac{2}{3}=\frac{7}{3}\neq1 \end{cases} $

Pour le schéma [P], nous allons chercher $b_0$ et $b_1$ pour obtenir un schéma d’ordre 2. Pour cela on doit résoudre le système linéaire : $ \begin{cases} \xi(0) = a_{0} + a_{1} = \frac{4}{3} -\frac{1}{3} = 1 \\ \xi(1) = - a_{1} + b_{0} + b_{1} + b_{-1} = \frac{1}{3}+ b_{0} + b_{1}=1\\ \xi(2) = a_{1} - 2 \left(b_{1} - b_{-1}\right) = -\frac{1}{3}-2b_1=1 \end{cases} $ On trouve $b_0=\frac{4}{3}$ et $b_1=-\frac{2}{3}$.

Le schéma predictor-corrector s’écrit $ \begin{cases} u_0 = y_0,\\ u_1 \approx y_1 ,\\ \tilde u = \frac{4}{3}u_{n}-\frac{1}{3}u_{n-1}+h\left(\frac{4}{3}\varphi_{n}-\frac{2}{3}\varphi_{n-1}\right) &\text{[P]} \\ u_{n+1}=\frac{4}{3}u_{n}-\frac{1}{3}u_{n-1}+h\frac{2}{3}\varphi(t_{n+1},\tilde u) \end{cases} $ Comme le prédictor est d’ordre 2 et le corrector d’ordre 2 aussi, le schéma prédictor-corrector est d’ordre 2 (car $\min(\omega_{[C]},\omega_{[P]}+1)=2$).

Code

%reset -f

import numpy as np
from scipy.optimize import fsolve


# variables globales
t0     = 0
tfinal = 1
y0     = 0

phi = lambda t,y : 1/(1+t**2)-2*y**2
sol_exacte = lambda t : t/(1+t**2)


# ==================================================
# Multipas
# ==================================================

q = 2
aa = [4/3,-1/3] # P et C
bb = [4/3,-2/3] # P
bm1 = 2/3 # C

def multipas_P(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[:q] = sol_exacte(tt[:q])
    for n in range(q-1,len(tt) - 1):
        k_0 = phi(tt[n],uu[n])
        k_1 = phi(tt[n-1],uu[n-1])
        uu[n+1] = aa[0]*uu[n] + aa[1]*uu[n-1] + h*(bb[0]*k_0+bb[1]*k_1)
    return uu

def multipas_C(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[:q] = sol_exacte(tt[:q])
    for n in range(q-1,len(tt) - 1):
        k_0 = phi(tt[n],uu[n])
        k_1 = phi(tt[n-1],uu[n-1])
        temp = fsolve ( lambda x : -x + aa[0]*uu[n] + aa[1]*uu[n-1] + h*bm1*phi(tt[n+1],x) , uu[n])
        uu[n+1] = temp[0]
    return uu

def multipas_PC(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[:q] = sol_exacte(tt[:q])
    for n in range(q-1,len(tt) - 1):
        k_0 = phi(tt[n],uu[n])
        k_1 = phi(tt[n-1],uu[n-1])
        u_tilde = aa[0]*uu[n] + aa[1]*uu[n-1] + h*(bb[0]*k_0+bb[1]*k_1)
        uu[n+1] = aa[0]*uu[n] + aa[1]*uu[n-1] + h*bm1*phi(tt[n+1],u_tilde)
    return uu



import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 16})


H = []
err_p = []
err_c = []
err_pc = []
N = 10
for k in range(7):
    N += 20
    tt = np.linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = sol_exacte(tt)
    uu_p = multipas_P(phi, tt, sol_exacte)
    uu_c = multipas_C(phi, tt, sol_exacte)
    uu_pc = multipas_PC(phi, tt, sol_exacte)
    H.append(h)
    err_p.append(np.linalg.norm(uu_p-yy, np.inf))
    err_c.append(np.linalg.norm(uu_c-yy, np.inf))
    err_pc.append(np.linalg.norm(uu_pc-yy, np.inf))
    

plt.figure(figsize=(8,5))
plt.plot(np.log(H), np.log(err_p), 'r-o', label=f'Multipas P ordre  = {np.polyfit(np.log(H),np.log(err_p) , 1)[0]:1.2f}')
plt.plot(np.log(H), np.log(err_c), 'b-o', label=f'Multipas C ordre  = {np.polyfit(np.log(H),np.log(err_c) , 1)[0]:1.2f}' )
plt.plot(np.log(H), np.log(err_pc), 'c-o', label=f'Multipas PC ordre = {np.polyfit(np.log(H),np.log(err_pc), 1)[0]:1.2f}')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.legend(bbox_to_anchor=(1.04, 1), loc='upper left')
plt.grid(True);