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

Auteur·rice

Gloria Faccanoni

Date de publication

27 mars 2026

Code

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
import sympy as sp
from IPython.display import display, Math, Markdown, Latex

PARTIE I : Étude de 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 (MPL) à $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 ($2p+1$) coefficients donnés.

Méthode explicite : $b_{-1}=0$.
Méthode implicite : $b_{-1}\neq0$.

Théorème de Lax-Richtmyer

Pour $\omega\ge1$, le schéma est convergent d’ordre $\omega$ ssi il est consistant d’ordre $\omega$ et zéro-stable.

Zéro-stabilité (comportement à $h\to0$)

On définit le premier polynôme caractéristique $ \varrho(r) = r^{p+1} - \sum_{j=0}^p a_jr^{p-j}. $

Condition des racines : la méthode est zéro-stable si toutes les racines $r_j$ de $\varrho(r)$ sont dans le disque unité ($|r_j| \le 1$), et si celles sur le cercle ($|r_j| = 1$) sont simples.

Consistance et ordre

On définit

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

Conditions d’ordre :

Ordre $1$ (= consistance) : $\xi(0)=\xi(1)=1$.
Ordre $\omega$ (= consistance d’ordre $\omega$) : $\xi(0)=\xi(1)=\dots=\xi(\omega)=1$.

Première barrière de Dahlquist

Il existe une limite sur l’ordre maximal $\omega$ en fonction du nombre de pas $q$ :

$ \omega_{\max} = \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} $

Intervalle de A-stabilité

Considérons le problème test $y'(t)=-\beta y(t)$, $y(0)=1$ avec $\beta>0$ un nombre réel positif. L’intervalle de stabilité absolue est l’ensemble des valeurs $h(\beta)$ pour lesquelles $u_n\to0$ lorsque $n\to\infty$.

Soit le polynôme caractéristique $ \begin{aligned} p_x(r) = \varrho(r) + x \sigma(r) &= r^{p+1} - \sum_{j=0}^p a_jr^{p-j} + x\sum_{j=-1}^p b_jr^{p-j} \\ &= (1 + x b_{-1})r^{p+1} - \sum_{j=0}^p (a_j - x b_j)r^{p-j}. \end{aligned} $

Remarque : $x = \beta h = -\lambda h = -z$.

Condition des racines strictes : si, pour $x\in]0,x_{\max}[$ donné, toutes les racines $r_j(x)$ de $p_x(r)$ sont dans le disque unité ($|r_j(x)| < 1$), alors la méthode est A-stable pour le pas $h\in]0,\frac{x_{\max}}{\beta}[$.

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
# =============================================================================

def xi(i, q):
    # Définition des symboles pour le cas général
    aa = sp.symbols(f'a_0:{q}')
    bb = sp.symbols(f'b_0:{q}')
    bm1 = sp.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 = 3     # 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}) = {sp.latex(xi(i, q))}$"))

# Affichage sous forme de système 
systeme = [sp.Eq(sp.Symbol(f"\\xi({i})") , sp.Eq(xi(i, q),1)) for i in range(omega + 1)]
display(Markdown(r"\(\begin{cases}"+",\\\\ ".join([sp.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=4$ pas, on affiche $\xi(i)$ pour $i = 0, \dots, 6$

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

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

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

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

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

$\qquad$ $\xi(5) = - a_{1} - 32 a_{2} - 243 a_{3} + 5 \left(b_{1} + 16 b_{2} + 81 b_{3} + b_{-1}\right)$

$\qquad$ $\xi(6) = a_{1} + 64 a_{2} + 729 a_{3} - 6 \left(b_{1} + 32 b_{2} + 243 b_{3} - b_{-1}\right)$

$\begin{cases}\xi(0) = a_{0} + a_{1} + a_{2} + a_{3} = 1,\\ \xi(1) = - a_{1} - 2 a_{2} - 3 a_{3} + b_{0} + b_{1} + b_{2} + b_{3} + b_{-1} = 1,\\ \xi(2) = a_{1} + 4 a_{2} + 9 a_{3} - 2 \left(b_{1} + 2 b_{2} + 3 b_{3} - b_{-1}\right) = 1,\\ \xi(3) = - a_{1} - 8 a_{2} - 27 a_{3} + 3 \left(b_{1} + 4 b_{2} + 9 b_{3} + b_{-1}\right) = 1,\\ \xi(4) = a_{1} + 16 a_{2} + 81 a_{3} - 4 \left(b_{1} + 8 b_{2} + 27 b_{3} - b_{-1}\right) = 1,\\ \xi(5) = - a_{1} - 32 a_{2} - 243 a_{3} + 5 \left(b_{1} + 16 b_{2} + 81 b_{3} + b_{-1}\right) = 1,\\ \xi(6) = a_{1} + 64 a_{2} + 729 a_{3} - 6 \left(b_{1} + 32 b_{2} + 243 b_{3} - b_{-1}\right) = 1\end{cases}$

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

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

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(bbox_to_anchor=(1.05, 1), loc='upper left');

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

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

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(bbox_to_anchor=(1.05, 1), loc='upper left');

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

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

sp.var('r')
rho = r**3+r**2-r-1
sp.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 à $q=2$ pas consistant et zéro-stable peut être d’ordre $\omega=q+2=4$. (Pour cela on sait qu’il sera implicite).

Construire l’unique schéma MPL implicite à $q=2$ pas convergent d’ordre $\omega=q+2=4$.
En étudier sa stabilité absolue.

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

On peut résoudre ce système avec la méthode de triangulation de Gauss :

$ \left(\begin{array}{ccccc|c} 1&0& 1& 0&0&1\\ 0&1&-1& 1&1&1\\ 0&0& 1&-2&2&1\\ 0&0&-1& 3&3&1\\ 0&0& 1&-4&4&1 \end{array}\right) \rightarrow \left(\begin{array}{ccccc|c} 1&0& 1& 0&0&1\\ 0&1&-1& 1&1&1\\ 0&0& 1&-2&2&1\\ 0&0& 0& 1&5&2\\ 0&0& 0&-2&2&0 \end{array}\right) \rightarrow \left(\begin{array}{ccccc|c} 1&0& 1& 0& 0&1\\ 0&1&-1& 1& 1&1\\ 0&0& 1&-2& 2&1\\ 0&0& 0& 1& 5&2\\ 0&0& 0& 0&12&4 \end{array}\right) $

donc $b_{-1}=\frac{1}{3}$, $b_0=2-5b_{-1}=\frac{4}{3}$, $b_1=1-2b_{-1}+2b_1=\frac{1}{3}$, $a_0=1$ et $a_1=0$.

Code

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

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

# Résultat
u_np1 = sp.symbols('u_{n+1}')
u_n = sp.symbols('u_n')
u_nm1 = sp.symbols('u_{n-1}')
phi_np1 = sp.symbols(r'\varphi_{n+1}')
phi_n = sp.symbols(r'\varphi_n')
phi_nm1 = sp.symbols(r'\varphi_{n-1}')
h = sp.symbols('h')
eq = sp.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)

{a_0: 0, a_1: 1, b_0: 4/3, b_1: 1/3, b_{-1}: 1/3}

$\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 les deux racines sont sur le cercle unité et simples.

Étudions maintrenant la A-stabilité de ce schéma. Le polynôme caractéristique est $ \begin{aligned} p_x(r) &= \varrho(r) + x \sigma(r) = (1+x b_{-1}) r^2 - (a_0-x b_0) r - (a_1-x b_1) \\ &= \left(1+\frac{x}{3}\right)r^2 - \frac{4x}{3}r - \left(1-\frac{x}{3}\right) \\ &= \frac{\left(3+x\right)r^2 - 4xr - \left(3-x\right)}{3}. \end{aligned} $

Les racines de ce polynôme sont $ r_{1,2}(x) = \frac{2x\pm\sqrt{9+3x^2}}{3+x} $ et la suite s’écrit $ u_{n+1} = A r_1^n(x) + B r_2^n(x) $ pour des constantes $A$ et $B$ dépendant des conditions initiales.

Si $|r_1(x)|<1$ et $|r_2(x)|<1$ alors $u_n\to0$ lorsque $n\to\infty$. Dans notre cas, $|r_1(x)|>1$ pour tout $x>0$, donc la méthode n’est pas A-stable. Le terme associé à cette racine est désigné dans la littérature sous le nom de « composante parasite ». Ces composantes nuisent à la précision globale (à long terme) de la méthode numérique.

Code

r, x = sp.symbols('r x')
poly = (3 + x)*r**2 - 4*x*r - (3 - x)
roots = sp.solve(poly, r)  # liste de racines symboliques
for i, root in enumerate(roots):
    display(Math(f"r_{i+1}(x) = {sp.latex(root)}"))

# Conversion en fonctions numpy pour évaluation
r1_func = sp.lambdify(x, roots[0], 'numpy')
r2_func = sp.lambdify(x, roots[1], 'numpy')

x_vals = np.linspace(0, 10, 100)
r_1 = r1_func(x_vals)
r_2 = r2_func(x_vals)

# Tracé
plt.figure(figsize=(8,5))

# Zone stable (|r|<1)
plt.fill_between(x_vals, -1, 1, color='gray', alpha=0.3, label=r'Zone stable $|r|<1$')

# Courbes des racines
plt.plot(x_vals, r_1, label=r'|r_1(x)|')
plt.plot(x_vals, r_2, label=r'|r_2(x)|')


plt.ylim(-1.5, 1.5)
plt.xlim(x_vals[0], x_vals[-1])
plt.xlabel('$x$')
plt.ylabel('Racines')
plt.title('Racines du polynôme caractéristique en fonction de $x$')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.grid(True)
plt.show()

$\displaystyle r_1(x) = \frac{2 x - \sqrt{3 x^{2} + 9}}{x + 3}$

$\displaystyle r_2(x) = \frac{2 x + \sqrt{3 x^{2} + 9}}{x + 3}$

Voici une figure montrant la solution numérique du problème de Cauchy $y'(t)=-\beta y(t)$, $y(0)=1$ avec $\beta=1$ sur l’intervalle $[0,10]$ obtenue avec la méthode $\text{MS}_2$. On voit que la solution numérique diverge de la solution exacte $y(t)=e^{-t}$ à mesure que le temps augmente, ce qui est une manifestation de la composante parasite associée à la racine $r_1(x)$. Si vous changez l’intervalle de temps pour $[0,1]$ vous verrez que la solution numérique est très proche de la solution exacte, ce qui illustre le fait que la composante parasite a un effet à long terme.

Code

t_fin = 10

# =============================================================================

h = 0.1

phi = lambda t, y: -y
t0, y0 = 0,1
sol_exacte = lambda t: y0*np.exp(-t)

def MS2(phi,tt,yy_init):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    # initialisation
    q = 2
    uu[:q] = yy_init
    # récurrence
    for n in range(q-1,len(tt)-1):
        k0 = phi( tt[n], uu[n] )
        k1 = phi( tt[n-1], uu[n-1] )
        temp = fsolve(lambda x: -x+uu[n-1]+h*(phi(tt[n+1],x)+4*k0+k1)/3, uu[n])
        uu[n+1] = temp[0]
    return uu

tt = np.arange(t0, t_fin, h)
yy = sol_exacte(tt)
y1 = y0+h*phi(tt[0],y0)
uu = MS2(phi, tt, [y0, y1])
plt.plot(tt, uu, 'o-', label=f'MPL', lw=2)

plt.plot(tt, yy, label='Solution exacte')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.show()

👣 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} $
Pour $\gamma=\frac{1}{2}$, étudier la A-stabilité de ce schéma.

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

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

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


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

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

sol_exacte = sp.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

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

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 = {(np.polyfit(np.log(H),np.log(err_mp), 1)[0])}')
plt.grid(True)

Intervalle de A-stabilité. Pour étudier la A-stabilité de ce schéma, on considère le problème test $y'(t)=-\beta y(t)$, $y(0)=1$ avec $\beta>0$ et on pose $x=h\beta$. Le polynôme caractéristique est $ \begin{aligned} p_x(r) &= \varrho(r) + x \sigma(r) = (1+x b_{-1}) r^2 - (a_0+x b_0) r - (a_1+x b_1) \\ &= \left(1+\frac{x}{4}(\gamma+3)\right)r^2 - (1-\gamma) r - \left(\gamma-\frac{x}{4}(3\gamma+1)\right) \end{aligned} $ Pour $\gamma=\frac{1}{2}$, on trouve $ p_x(r) = \left(1+\frac{7x}{8}\right)r^2 - \frac{1}{2} r - \left(\frac{5x}{8}-\frac{1}{2}\right) $

En traçant la valeur absolue des racines en fonction de $x$, on peut visualiser l’intervalle de stabilité absolue. On trouve que la méthode est A-stable pour tout $x>0$. En effet, les racines sont complexes conjuguées et de module strictement inférieur à 1 pour tout $x>0$. Par conséquent, la méthode est A-stable pour tout $h>0$.

Code

r = sp.Symbol('r')
x = sp.Symbol('x', real=True, positive=True)
gamma = sp.Symbol('gamma')

p = 1
aa = [1-gamma, gamma]
bb = [0, (3*gamma+1)/4]
bm1 = (gamma+3)/4

# Construction du polynôme
rho_expr = r**(p+1) - sum(aa[i]*r**(p-i) for i in range(p+1))
sigma_expr = bm1*r**(p+1) + sum(bb[i]*r**(p-i) for i in range(p+1))
poly = sp.Poly(rho_expr + x * sigma_expr, r)

display(Math(f"p_x(r) = {sp.latex(poly.as_expr())}"))

coeffs = poly.all_coeffs() # Renvoie [A, B, C]
display(Math(rf"A(x) = {sp.latex(coeffs[0])}"))
display(Math(rf"B(x) = {sp.latex(coeffs[1])}"))
display(Math(rf"C(x) = {sp.latex(coeffs[2])}"))

# =============================================================================

# Cas gamma = 1/2
gamma_val = sp.Rational(1, 2)
poly_gamma = sp.Poly(poly.subs(gamma, gamma_val), r)

display(Math(f"$p_x(r)$ pour $\\gamma=1/2$: ${sp.latex(poly_gamma.as_expr())}$"))

display(Math(rf"A(x) = {sp.latex(coeffs[0].subs(gamma, gamma_val))}"))
display(Math(rf"B(x) = {sp.latex(coeffs[1].subs(gamma, gamma_val))}"))
display(Math(rf"C(x) = {sp.latex(coeffs[2].subs(gamma, gamma_val))}"))

roots = sp.solve(poly_gamma, r)

display(Markdown(f"Racines : $r_1(x) = {sp.latex(roots[0])}$, $r_2(x) = {sp.latex(roots[1])}$"))

# On définit une fonction qui calcule le module manuellement pour éviter le sqrt(négatif) dans numpy
def get_abs_expr(root):
    re_part = sp.re(root)
    im_part = sp.im(root)
    return sp.sqrt(re_part**2 + im_part**2)

abs_r1_expr = get_abs_expr(roots[0])
abs_r2_expr = get_abs_expr(roots[1])

r1_func = sp.lambdify(x, abs_r1_expr, modules=['numpy'])
r2_func = sp.lambdify(x, abs_r2_expr, modules=['numpy'])


# Le discriminant s'annule en
x_bifurcation = sp.solve(poly_gamma.discriminant(), x)[0]
display(Markdown(f"Point de bifurcation (racines fusionnent) : $x = {sp.latex(x_bifurcation)}$"))

x_vals = np.linspace(0, 5, 1000)

r_abs_1 = r1_func(x_vals)
r_abs_2 = r2_func(x_vals)


plt.figure(figsize=(10, 6))

plt.axhspan(0, 1, color='gray', alpha=0.1, label=r'Zone stable $|r| \leq 1$')

plt.plot(x_vals, r_abs_1, label=r'$|r_1(x)|$', lw=2)
plt.plot(x_vals, r_abs_2, label=r'$|r_2(x)|$', lw=2, linestyle='--')

# Point où les racines fusionnent (deviennent complexes)
plt.axvline(x_bifurcation, color='orange', linestyle=':', label=f'Bifurcation x = {x_bifurcation}')

plt.title(r'Module des racines ($y\' = -\beta y, \gamma=1/2$)')
plt.xlabel('$x = h\\beta$')
plt.ylabel('$|r(x)|$')
plt.legend()
plt.grid(True, alpha=0.3)
plt.ylim(0, 1.2)
plt.show()

$\displaystyle p_x(r) = \frac{3 \gamma x}{4} - \gamma + r^{2} \left(\frac{\gamma x}{4} + \frac{3 x}{4} + 1\right) + r \left(\gamma - 1\right) + \frac{x}{4}$

$\displaystyle A(x) = \frac{\gamma x}{4} + \frac{3 x}{4} + 1$

$\displaystyle B(x) = \gamma - 1$

$\displaystyle C(x) = \frac{3 \gamma x}{4} - \gamma + \frac{x}{4}$

$\displaystyle p_x(r)$ pour $\gamma=1/2$: $r^{2} \cdot \left(\frac{7 x}{8} + 1\right) - \frac{r}{2} + \frac{5 x}{8} - \frac{1}{2}$

$\displaystyle A(x) = \frac{7 x}{8} + 1$

$\displaystyle B(x) = - \frac{1}{2}$

$\displaystyle C(x) = \frac{5 x}{8} - \frac{1}{2}$

Racines : $r_1(x) = \frac{2 - \sqrt{- 35 x^{2} - 12 x + 36}}{7 x + 8}$, $r_2(x) = \frac{\sqrt{- 35 x^{2} - 12 x + 36} + 2}{7 x + 8}$

Point de bifurcation (racines fusionnent) : $x = \frac{6}{7}$

👣 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} $
Pour les valeurs de $\gamma$ qui donnent la convergence, étudier la A-stabilité de ce schéma.

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

t0     = 0
tfinal = 1
y0     = 2

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


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

sol_exacte = sp.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)

Intervalle de A-stabilité
On étudie la A-stabilité seulement pour les valeurs de $\gamma$ qui donnent la convergence.

Code

r = sp.Symbol('r')
x = sp.Symbol('x', real=True, positive=True)
gamma = sp.Symbol('gamma')

p = 1
aa = [gamma**2/3+gamma/2, gamma**2/3+gamma/2-1]
bb = [-gamma**2/6-gamma/4+2, -gamma**2/6-gamma/4]
bm1 = sp.S(0)

# Construction du polynôme
rho_expr = r**(p+1) - sum(aa[i]*r**(p-i) for i in range(p+1))
sigma_expr = bm1*r**(p+1) + sum(bb[i]*r**(p-i) for i in range(p+1))
poly = sp.Poly(rho_expr + x * sigma_expr, r)

# display(Math(f"p_x(r) = {sp.latex(poly.as_expr())}"))

coeffs = poly.all_coeffs() # Renvoie [A, B, C]
# display(Math(rf"A(x) = {sp.latex(coeffs[0])}"))
# display(Math(rf"B(x) = {sp.latex(coeffs[1])}"))
# display(Math(rf"C(x) = {sp.latex(coeffs[2])}"))

# =============================================================================

for gamma_val in [1, -2]:

    display(Markdown(f"**Cas $\\gamma = {gamma_val}$**"))

    poly_gamma = sp.Poly(poly.subs(gamma, gamma_val), r)

    display(Math(f"$p_x(r) = {sp.latex(poly_gamma.as_expr())}$"))

    display(Math(rf"A(x) = {sp.latex(coeffs[0].subs(gamma, gamma_val))}"))
    display(Math(rf"B(x) = {sp.latex(coeffs[1].subs(gamma, gamma_val))}"))
    display(Math(rf"C(x) = {sp.latex(coeffs[2].subs(gamma, gamma_val))}"))

    roots = sp.solve(poly_gamma, r)

    display(Markdown(f"Racines : $r_1(x) = {sp.latex(roots[0])}$, $r_2(x) = {sp.latex(roots[1])}$"))

    # On définit une fonction qui calcule le module manuellement pour éviter le sqrt(négatif) dans numpy
    def get_abs_expr(root):
        re_part = sp.re(root)
        im_part = sp.im(root)
        return sp.sqrt(re_part**2 + im_part**2)

    abs_r1_expr = get_abs_expr(roots[0])
    abs_r2_expr = get_abs_expr(roots[1])

    r1_func = sp.lambdify(x, abs_r1_expr, modules=['numpy'])
    r2_func = sp.lambdify(x, abs_r2_expr, modules=['numpy'])


    # Le discriminant s'annule en
    x_bifurcation = sp.solve(poly_gamma.discriminant(), x)[0]
    display(Markdown(f"Point de bifurcation (racines fusionnent) : $x = {sp.latex(x_bifurcation)}$"))

    x_vals = np.linspace(0, 5, 1000)

    r_abs_1 = r1_func(x_vals)
    r_abs_2 = r2_func(x_vals)


    plt.figure(figsize=(10, 6))

    plt.axhspan(0, 1, color='gray', alpha=0.1, label=r'Zone stable $|r| \leq 1$')

    plt.plot(x_vals, r_abs_1, label=r'$|r_1(x)|$', lw=2)
    plt.plot(x_vals, r_abs_2, label=r'$|r_2(x)|$', lw=2, linestyle='--')

    plt.title(rf'Module des racines ($y\' = -\beta y, \gamma={gamma_val}$)')
    plt.xlabel('$x = h\\beta$')
    plt.ylabel('$|r(x)|$')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.ylim(0, 1.2)
    plt.show()

Cas $\gamma = 1$

$\displaystyle p_x(r) = r^{2} + r \left(\frac{19 x}{12} - \frac{5}{6}\right) - \frac{5 x}{12} + \frac{1}{6}$

$\displaystyle A(x) = 1$

$\displaystyle B(x) = \frac{19 x}{12} - \frac{5}{6}$

$\displaystyle C(x) = \frac{1}{6} - \frac{5 x}{12}$

Racines : $r_1(x) = - \frac{19 x}{24} - \frac{\sqrt{361 x^{2} - 140 x + 4}}{24} + \frac{5}{12}$, $r_2(x) = - \frac{19 x}{24} + \frac{\sqrt{361 x^{2} - 140 x + 4}}{24} + \frac{5}{12}$

Point de bifurcation (racines fusionnent) : $x = \frac{70}{361} - \frac{24 \sqrt{6}}{361}$

Cas $\gamma = -2$

$\displaystyle p_x(r) = r^{2} + r \left(\frac{11 x}{6} - \frac{1}{3}\right) - \frac{x}{6} + \frac{2}{3}$

$\displaystyle A(x) = 1$

$\displaystyle B(x) = \frac{11 x}{6} - \frac{1}{3}$

$\displaystyle C(x) = \frac{2}{3} - \frac{x}{6}$

Racines : $r_1(x) = - \frac{11 x}{12} - \frac{\sqrt{121 x^{2} - 20 x - 92}}{12} + \frac{1}{6}$, $r_2(x) = - \frac{11 x}{12} + \frac{\sqrt{121 x^{2} - 20 x - 92}}{12} + \frac{1}{6}$

Point de bifurcation (racines fusionnent) : $x = \frac{10}{121} + \frac{12 \sqrt{78}}{121}$

👣 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} $
Pour les valeurs de $\gamma$ qui donnent la convergence, étudier la A-stabilité de ce schéma.

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

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


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

cond1 = a0+a1
cond2 = -(0*a0+1*a1)+(bm1+b0+b1)
display(sp.solve([sp.Eq(cond1,1),sp.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}$

[(3,)]

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

Code

display(Math('a_0='+sp.latex(a0.subs(gamma,3))))
display(Math('a_1='+sp.latex(a1.subs(gamma,3))))
display(Math('b_0='+sp.latex(b0.subs(gamma,3))))
display(Math('b_1='+sp.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  = sp.Symbol('t')
y  = sp.Function('y')
edo = sp.Eq( sp.diff(y(t),t) , phi(t,y(t)) )
display(edo)
solpar = sp.dsolve(edo, ics={y(t0): y0})
display(solpar)

sol_exacte = sp.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)

Intervalle de A-stabilité

On étudie la A-stabilité seulement pour les valeurs de $\gamma$ qui donnent la convergence.

Le schéma étant explicite, il n’est au mieux A-stable que pour certaines valeurs de $x$. En traçant la valeur absolue des racines en fonction de $x$, on trouve que la méthode est A-stable pour $x\in]0,x_{\max}[$ avec $x_{\max}=1$.

Code

r = sp.Symbol('r')
x = sp.Symbol('x', real=True, positive=True)
gamma = sp.Symbol('gamma')

p = 1
aa = [gamma**2-9*gamma/2+11/sp.S(2), gamma**2-7*gamma/2+3/sp.S(2)]
bb = [gamma*(gamma-2)/2, gamma*(gamma-10/sp.S(3))/2]
bm1 = sp.S(0)

# Construction du polynôme
rho_expr = r**(p+1) - sum(aa[i]*r**(p-i) for i in range(p+1))
sigma_expr = bm1*r**(p+1) + sum(bb[i]*r**(p-i) for i in range(p+1))
poly = sp.Poly(rho_expr + x * sigma_expr, r)

# display(Math(f"p_x(r) = {sp.latex(poly.as_expr())}"))

coeffs = poly.all_coeffs() # Renvoie [A, B, C]
# display(Math(rf"A(x) = {sp.latex(coeffs[0])}"))
# display(Math(rf"B(x) = {sp.latex(coeffs[1])}"))
# display(Math(rf"C(x) = {sp.latex(coeffs[2])}"))

# =============================================================================

gamma_val = sp.S(3)

display(Markdown(f"**Cas $\\gamma = {gamma_val}$**"))

poly_gamma = sp.Poly(poly.subs(gamma, gamma_val), r)

display(Math(f"$p_x(r) = {sp.latex(poly_gamma.as_expr())}$"))

display(Math(rf"A(x) = {sp.latex(coeffs[0].subs(gamma, gamma_val))}"))
display(Math(rf"B(x) = {sp.latex(coeffs[1].subs(gamma, gamma_val))}"))
display(Math(rf"C(x) = {sp.latex(coeffs[2].subs(gamma, gamma_val))}"))

roots = sp.solve(poly_gamma, r)

display(Markdown(f"Racines : $r_1(x) = {sp.latex(roots[0])}$, $r_2(x) = {sp.latex(roots[1])}$"))

# On définit une fonction qui calcule le module manuellement pour éviter le sqrt(négatif) dans numpy
def get_abs_expr(root):
    re_part = sp.re(root)
    im_part = sp.im(root)
    return sp.sqrt(re_part**2 + im_part**2)

abs_r1_expr = get_abs_expr(roots[0])
abs_r2_expr = get_abs_expr(roots[1])


r1_func = sp.lambdify(x, abs_r1_expr, modules=['numpy'])
r2_func = sp.lambdify(x, abs_r2_expr, modules=['numpy'])


# On calcule les points où |r_1|=1
from scipy.optimize import fsolve
display(Markdown(f"Point où $|r_1|=1$ : $x = {fsolve(lambda x : r1_func(x) - 1, 1)[0]:.4f}$"))


x_vals = np.linspace(0, 5, 1000)

r_abs_1 = r1_func(x_vals)
r_abs_2 = r2_func(x_vals)

plt.figure(figsize=(10, 6))

plt.axhspan(0, 1, color='gray', alpha=0.1, label=r'Zone stable $|r| \leq 1$')

plt.plot(x_vals, r_abs_1, label=r'$|r_1(x)|$', lw=2)
plt.plot(x_vals, r_abs_2, label=r'$|r_2(x)|$', lw=2, linestyle='--')

plt.title(rf'Module des racines ($y\' = -\beta y, \gamma={gamma_val}$)')
plt.xlabel('$x = h\\beta$')
plt.ylabel('$|r(x)|$')
plt.legend()
plt.grid(True, alpha=0.3)
plt.ylim(0, 1.2)
plt.show()

Cas $\gamma = 3$

$\displaystyle p_x(r) = r^{2} + r \left(\frac{3 x}{2} - 1\right) - \frac{x}{2}$

$\displaystyle A(x) = 1$

$\displaystyle B(x) = \frac{3 x}{2} - 1$

$\displaystyle C(x) = - \frac{x}{2}$

Racines : $r_1(x) = - \frac{3 x}{4} - \frac{\sqrt{9 x^{2} - 4 x + 4}}{4} + \frac{1}{2}$, $r_2(x) = - \frac{3 x}{4} + \frac{\sqrt{9 x^{2} - 4 x + 4}}{4} + \frac{1}{2}$

Point où $|r_1|=1$ : $x = 1.0000$

PARTIE II : Construction de schémas multi-pas linéaires

Exercice : Interpolation, Quadrature et EDO

Soit $f$ une fonction de classe $\mathcal{C}^1([-1,1])$. Écrire le polynôme $p\in\mathbb{R}_2[\tau]$ qui interpole $f$ aux points $-1$, $0$ et $1$.
Construire une méthode de quadrature comme suit: $ \int_{0}^{1}f(\tau)\mathrm{d}\tau \approx \int_{0}^{1}p(\tau)\mathrm{d}\tau. $ NB: on intègre sur $[0,1]$ mais on interpole en $-1$, $0$ et $1$.
À l’aide d’un changement de variable affine entre l’intervalle $[0,1]$ et l’intervalle $[a,b]$, en déduire une formule de quadrature pour l’intégrale $ \int_{a}^{b}f(x) \mathrm{d}x $ lorsque $f$ est une fonction de classe $\mathcal{C}^1([2a-b,b])$.
Remarque: $[2a-b,b]=[a-(b-a),a+(b-a)]$
Considérons le problème de Cauchy: trouver $y \colon [t_0,T]\subset \mathbb{R} \to \mathbb{R}$ tel que $ \begin{cases} y'(t) = \varphi(t,y(t)), &\forall t \in [t_0,T],\\ y(t_0) = y_0, \end{cases} $ dont on suppose l’existence d’une unique solution $y$.
- On subdivise l’intervalle $[t_0;T]$ en $N$ intervalles $[t_{n};t_{n+1}]$ de largeur $h=\dfrac{T-t_0}{N}$ avec $t_n=t_0+nh$ pour $n=0,\dots,N$.
- Utiliser la formule obtenue au point 3 pour approcher l’intégrale $ \int_{t_n}^{t_{n+1}}\varphi(t,y(t))\mathrm{d}t. $
- En déduire un schéma à deux pas implicite pour l’approximation de la solution du problème de Cauchy.

On cherche les coefficients $\alpha$, $\beta$ et $\gamma$ du polynôme $p(\tau)=\alpha+\beta \tau+\gamma \tau^2$ tels que $ \begin{cases} p(-1)=f(-1),\\ p(0) =f(0),\\ p(1) =f(1), \end{cases} \qquad\text{c'est à dire}\qquad \begin{cases} \alpha-\beta+\gamma=f(-1),\\ \alpha=f(0),\\ \alpha+\beta+\gamma=f(1). \end{cases} $ Donc $\alpha=f(0)$, $\beta=\frac{f(1)-f(-1)}{2}$ et $\gamma=\frac{f(1)-2f(0)+f(-1)}{2}$.

Code

f_0 = sp.Symbol('f(0)')
f_1 = sp.Symbol('f(1)')
f_m1 = sp.Symbol('f(-1)')
t = sp.Symbol('t')

p = sp.interpolate([(1,f_1),(0,f_0),(-1,f_m1)], t).factor(t)
display(sp.Eq(sp.Symbol('p'), p))
q = (sp.integrate(p,(t,0,1)).simplify())
display(sp.Eq( sp.Integral(sp.Symbol('p(t)'), (t, 0, 1)) , q ))

$\displaystyle p = \frac{2 f(0) + t^{2} \left(f(-1) - 2 f(0) + f(1)\right) + t \left(- f(-1) + f(1)\right)}{2}$

$\displaystyle \int\limits_{0}^{1} p(t)\, dt = - \frac{f(-1)}{12} + \frac{2 f(0)}{3} + \frac{5 f(1)}{12}$

On en déduit la méthode de quadrature $ \int_{0}^{1}f(\tau)\mathrm{d}\tau\approx\int_{0}^{1}p(\tau)\mathrm{d}\tau=\alpha+\frac{\beta}{2}+\frac{\gamma}{3}=\frac{-f(-1)+8f(0)+5f(1)}{12}. $

Soit $x=m\tau+q$, alors $ \int_{a}^{b} f(x)\mathrm{d}x=m\int_{0}^{1}f(m\tau+q)\mathrm{d}\tau \quad\text{avec}\quad \begin{cases} a=q,\\ b=m+q, \end{cases} \quad\text{i.e.}\quad \begin{cases} m=b-a,\\ q=a, \end{cases} $ d’où le changement de variable $x=(b-a)\tau+a$. On en déduit la formule de quadrature $ \int_{a}^{b}\!\!\!\!f(x)\mathrm{d}x=(b-a)\int_{0}^{1}f\left((b-a)\tau+a\right)\mathrm{d}\tau\approx(b-a)\frac{-f(2a-b)+8f(a)+5f(b)}{12}. $

On pose $a=t_n$ et $b=t_{n+1}$ d’où la formule de quadrature $\int_{t_{n}}^{t_{n+1}}\!\!\!\!f(t)\mathrm{d}t\approx(t_{n+1}-t_{n}) \frac{-f(2t_{n}-t_{n+1})+8f(t_{n})+5f(t_{n+1})}{12}=h \frac{-f(t_{n-1})+8f(t_{n})+5f(t_{n+1})}{12}.$

En utilisant la formule de quadrature pour l’intégration de l’EDO $y'(t)=\varphi(t,y(t))$ entre $t_n$ et $t_{n+1}$ on obtient $ y(t_{n+1})=y(t_n)+\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t\approx h \dfrac{-\varphi(t_{n-1},y(t_{n-1}))+8\varphi(t_{n},y(t_{n}))+5\varphi(t_{n+1},y(t_{n+1}))}{12}. $

Si on note $u_n$ une approximation de la solution $y$ au temps $t_n$, on obtient le schéma à deux pas implicite suivant: $ \begin{cases} u_0=y(t_0)=y_0,\\ u_1\text{ à définir},\\ u_{n+1}=u_n+ h \dfrac{-\varphi(t_{n-1},u_{n-1})+8\varphi(t_{n},u_{n})+5\varphi(t_{n+1},u_{n+1})}{12} & n=1,2,\dots N-1 \end{cases} $ On peut utiliser une pédiction d’Euler explicite pour initialiser $u_1$: $ \begin{cases} u_0=y(t_0)=y_0,\\ u_1=u_0+h\varphi(t_0,u_0),\\ u_{n+1}=u_n+ h \dfrac{-\varphi(t_{n-1},u_{n-1})+8\varphi(t_{n},u_{n})+5\varphi(t_{n+1},u_{n+1})}{12} & n=1,2,\dots N-1 \end{cases} $

Exercice : Interpolation, Quadrature et EDO

Soit $h>0$ et $f\colon [a-h,a+h] \to\mathbb{R}$ une fonction de classe $\mathscr{C}^1([a-h,a+h])$. Écrire le polynôme $p\in\mathbb{R}_1[x]$ qui interpole $f$ aux points $a-h$ et $a$, i.e. l’équation de la droite $p\in\mathbb{R}_1[x]$ qui passe par les deux points $(a-h,f(a-h))$ et $(a,f(a))$.
Construire une méthode de quadrature comme suit : $ \int_{a}^{a+h}f(x)\mathrm{d}x \approx \int_{a}^{a+h}p(x)\mathrm{d}x. $ NB: on intègre sur $[a,a+h]$ mais on interpole en $a-h$ et $a$.
Considérons le problème de Cauchy: trouver $y \colon [t_0,T]\subset \mathbb{R} \to \mathbb{R}$ tel que $ \begin{cases} y'(t) = \varphi(t,y(t)), &\forall t \in [t_0,T],\\ y(t_0) = y_0, \end{cases} $ dont on suppose l’existence d’une unique solution $y$. On subdivise l’intervalle $[t_0;T]$ en $N$ intervalles $[t_{n};t_{n+1}]$ de largeur $h=\dfrac{T-t_0}{N}$ avec $t_n=t_0+nh$ pour $n=0,\dots,N$.
Utiliser la formule obtenue au point 2 pour approcher l’intégrale $ \int_{t_n}^{t_{n+1}}\varphi(t,y(t))\mathrm{d}t. $ En déduire un schéma à deux pas explicite pour l’approximation de la solution du problème de Cauchy.

$p(x)=\dfrac{f(a)-f(a-h)}{a-(a-h)}(x-a)+f(a)=\dfrac{f(a)-f(a-h)}{h}(x-a)+f(a)$.

Code

a = sp.Symbol('a')
h = sp.Symbol('h')
f_a = sp.Symbol('f(a)')
f_amh = sp.Symbol('f(a-h)')

p = sp.interpolate([(a,f_a),(a-h,f_amh)], t).factor(t)
display(sp.Eq(sp.Symbol('p'), p))
q = (sp.integrate(p,(t,a,a+h)).simplify())
display(sp.Eq( sp.Integral(sp.Symbol('p(t)'), (t, a, a+h)) , q ))

$\displaystyle p = \frac{- a f(a) + a f(a-h) + f(a) h + t \left(f(a) - f(a-h)\right)}{h}$

$\displaystyle \int\limits_{a}^{a + h} p(t)\, dt = \frac{h \left(3 f(a) - f(a-h)\right)}{2}$

On en déduit la méthode de quadrature

$\begin{aligned} \int_{a}^{a+h}f(x)\mathrm{d}x &\approx \int_{a}^{a+h}p(x)\mathrm{d}x \\&= \dfrac{f(a)-f(a-h)}{h}\left[ \frac{(x-a)^2}{2} \right]_a^{a+h}+f(a)\left[ x \right]_a^{a+h} \\&= \dfrac{f(a)-f(a-h)}{2h}\left((a+h-a)^2-(a-a)^2 \right)+f(a)(a+h-a) \\&= \dfrac{f(a)-f(a-h)}{2h} h^2 + hf(a) \\&= h \dfrac{3f(a)-f(a-h)}{2}. \end{aligned}$

On pose $a=t_n$ et $a+h=t_{n+1}$ d’où la formule de quadrature $ \int_{t_{n}}^{t_{n+1}}\!\!\!\!f(t)\mathrm{d}t\approx(t_{n+1}-t_{n}) \frac{3f(t_n)-f(2t_{n}-t_{n+1})}{2}=h \frac{3f(t_n)-f(t_{n-1})}{2}. $ En utilisant la formule de quadrature pour l’intégration de l’EDO $y'(t)=\varphi(t,y(t))$ entre $t_n$ et $t_{n+1}$ on obtient $ y(t_{n+1})=y(t_n)+\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t\approx h \dfrac{3\varphi(t_{n},y(t_{n}))-\varphi(t_{n-1},y(t_{n-1}))}{2}. $ Si on note $u_n$ une approximation de la solution $y$ au temps $t_n$, on obtient le schéma à deux pas implicite suivant: $ \begin{cases} u_0=y(t_0)=y_0,\\ u_1\text{ à définir},\\ u_{n+1}=u_n+ h \dfrac{3\varphi(t_{n},u_{n})-\varphi(t_{n-1},u_{n-1})}{2} & n=1,2,\dots N-1 \end{cases} $ On peut utiliser une prédiction d’Euler explicite pour initialiser $u_1$: $ \begin{cases} u_0=y(t_0)=y_0,\\ u_1=u_0+h\varphi(t_0,u_0),\\ u_{n+1}=u_n+ h \dfrac{3\varphi(t_{n},u_{n})-\varphi(t_{n-1},u_{n-1})}{2} & n=1,2,\dots N-1 \end{cases} $

Exercice : Construction d’un schéma

Considérons le problème de Cauchy

$\begin{cases} y'(t) = \varphi(t,y(t)), &\forall t \in I=[t_0,T],\\ y(t_0) = y_0, \end{cases}$

avec $y_0$ une valeur donnée et supposons que l’on ait montré l’existence et l’unicité d’une solution $y$ pour $t\in I$.

On subdivise l’intervalle $I=[t_0;T]$, avec $T<+\infty$, en $N$ intervalles $[t_{n};t_{n+1}]$ de largeur $h=(T-t_0)/N$ avec $t_n=t_0+nh$ pour $n=0,\dots,N$.
Pour chaque nœud $t_n$, on note $y_n=y(t_n)$ et on cherche la valeur inconnue $u_n$ qui approche la valeur exacte $y_n$

Si nous intégrons l’EDO $y'(t)=\varphi(t,y(t))$ entre $t_{n-2}$ et $t_{n+1}$ nous obtenons $ y_{n+1}-y_{n-2}=\int_{t_{n-2}}^{t_{n+1}} \varphi(t,y(t))dt. $ Écrire le schéma obtenu en approchant l’intégrale $\int_{t_{n-2}}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t$ par l’intégrale $\int_{t_{n-2}}^{t_{n+1}} p(t) \mathrm{d}t$ où $p$ est le polynôme interpolant $\varphi$ en $t_{n}$ et $t_{n-1}$ (attention à bien initialiser la suite définie par récurrence pour qu’on puisse effectivement calculer tous les termes). Le schéma obtenu est-il explicite ?

Le polynôme interpolant $\varphi$ en $t_{n}$ et $t_{n-1}$ a équation

$ p(t) = \frac{\varphi(t_n,y_n)-\varphi(t_{n-1},y_{n-1})}{t_n-t_{n-1}}(t-t_n)+ \varphi(t_n,y_n). $

On intègre ce polynôme entre $t_{n-2}$ et $t_{n+1}$:

$\begin{aligned} \int_{t_{n-2}}^{t_{n+1}} p(t) \mathrm{d}t &= \frac{\varphi(t_n,y_n)-\varphi(t_{n-1},y_{n-1})}{h}\left[\frac{(t-t_n)^2}{2}\right]_{t_{n-2}}^{t_{n+1}}+ \varphi(t_n,y_n)\left[t\right]_{t_{n-2}}^{t_{n+1}} \\ &=\frac{\varphi(t_n,y_n)-\varphi(t_{n-1},y_{n-1})}{2h}\left[ (t_{n+1}-t_n)^2-(t_{n-2}-t_n)^2 \right]+ \varphi(t_n,y_n)\left[t_{n+1}-t_{n-2}\right] \\ &=\frac{3}{2}h\left(\varphi(t_n,y_n)+\varphi(t_{n-1},y_{n-1})\right). \end{aligned}$

On obtient le schéma explicite

$\begin{cases} u_0=y_0,\\ u_1=u_0+h\varphi(t_0,u_0), \\ u_{n+1}=u_{n-2}+\frac{3}{2}h\left(\varphi(t_n,u_n)+\varphi(t_{n-1},u_{n-1})\right). \end{cases}$

Code

h = sp.Symbol('h')
t_n = sp.Symbol('t_{n}')
t_nm1 = sp.Symbol('t_{n-1}')
t_np1 = sp.Symbol('t_{n+1}')
t_nm2 = sp.Symbol('t_{n-2}')
f_n = sp.Symbol('f(t_{n})')
f_nm1 = sp.Symbol('f(t_{n-1})')

p = sp.interpolate([(t_n,f_n),(t_nm1,f_nm1)], t).factor(t)
display(sp.Eq(sp.Symbol('p'), p))
q = (sp.integrate(p,(t,t_nm2,t_np1)).subs(t_np1,t_n+h).subs(t_nm1,t_n-h).subs(t_nm2,t_n-2*h)).simplify()
display(sp.Eq( sp.Integral(sp.Symbol('p(t)'), (t, t_nm2,t_np1)) , q ))

$\displaystyle p = \frac{- f(t_{n-1}) t_{n} + f(t_{n}) t_{n-1} + t \left(f(t_{n-1}) - f(t_{n})\right)}{t_{n-1} - t_{n}}$

$\displaystyle \int\limits_{t_{n-2}}^{t_{n+1}} p(t)\, dt = \frac{3 h \left(f(t_{n-1}) + f(t_{n})\right)}{2}$

Exercice : construction d’un schéma

Considérons le problème de Cauchy:

$ \begin{cases} y'(t) = \varphi(t,y(t)), &\forall t \in [t_0,T],\\ y(t_0) = y_0. \end{cases} $

Supposons que l’on ait montré l’existence d’une unique solution $y$.

On subdivise l’intervalle $[t_0,T]$ en $N$ intervalles de longueur $h=(T-t_0)/N=t_{n+1}-t_n$. Pour chaque nœud $t_n=t_0 + nh$ ($1 \le n \le N$) on cherche la valeur inconnue $u_n$ qui approche $y(t_n)$. L’ensemble de $N+1$ valeurs $\{u_0 = y_0, u_1,\dots, u_{N} \}$ représente la solution numérique.

Dans cette exercice on va construire des nouveaux schémas numériques basés sur l’intégration de l’EDO $y'(t)=\varphi(t,y(t))$ entre $t_n$ et $t_{n+2}$: $ y(t_{n+2})=y(t_n)+\int_{t_n}^{t_{n+2}} \varphi(t,y(t))\mathrm{d}t. $

En utilisant la formule de quadrature du point milieu pour approcher le membre de droite écrire un schéma numérique explicite permettant de calculer $u_{n+2}$ à partir de $u_{n+1}$ et $u_{n}$. Notons que ce schéma a besoin de deux valeurs initiales ; on posera alors $u_0=y_0$ et $u_{1}$ sera approché par une prédiction d’Euler progressive.
En utilisant la formule de quadrature de Cavalieri-Simpson pour approcher le membre de droite écrire un schéma numérique implicite permettant de calculer $u_{n+2}$ à partir de $u_{n+1}$ et $u_{n}$. Notons que ce schéma a besoin de deux valeurs initiales ; on posera alors $u_0=y_0$ et $u_{1}$ sera approché par une prédiction d’Euler progressive.
Proposer une modification du schéma au point précédent pour qu’il devienne explicite.

Rappels: une formule de quadrature interpolatoire approche l’intégrale d’une fonction $f$ par l’intègrale dun polynôme qui interpole $f$ en des points donnés: $\int_a^b f(t)\mathrm{d}t\approx\int_a^b p(t)\mathrm{d}t$ + la formule de quadrature du point milieu considère comme points d’interpolation le point $\frac{a+b}{2}$ ainsi $ \int_a^b f(t)\mathrm{d}t\approx\int_a^b f\left(\frac{a+b}{2}\right)\mathrm{d}t = (b-a)f\left(\frac{a+b}{2}\right) $ + la formule de Simpson considère comme points d’interpolation les points $\left\{a,\frac{a+b}{2},b\right\}$ ainsi $ \int_a^b f(t)\mathrm{d}t\approx \frac{b-a}{6}\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right] $

Si on utilise la formule de quadrature du point milieu sur l’intervalle $[t_n;t_{n+2}]$, ie $ \int_{t_n}^{t_{n+2}} \varphi(t,y(t))\mathrm{d}t\approx 2h \varphi\left(t_{n+1},y(t_{n+1})\right) $ on obtient $\begin{cases} u_0=y(t_0)=y_0,\\ u_1\approx y(t_1)\\ u_{n+2}=u_{n}+2h \varphi(t_{n+1},u_{n+1})& n=0,1,\dots N-2 \end{cases}$ où $u_{1}$ est une approximation de $y(t_1)$. Nous pouvons utiliser une prédiction d’Euler progressive pour approcher $u_1$. Nous avons construit ainsi un nouveau schéma explicite à 2 pas: $\begin{cases} u_0=y(t_0)=y_0,\\ u_1=u_0+h\varphi(t_{0},u_{0}),\\ u_{n+2}=u_{n}+2h \varphi(t_{n+1},u_{n+1})& n=0,1,\dots N-2 \end{cases}$

Si on utilise la formule de quadrature de Cavalieri-Simpson sur l’intervalle $[t_n;t_{n+2}]$, ie $ \int_{t_n}^{t_{n+2}} \varphi(t,y(t))\mathrm{d}t\approx \frac{h}{3} \left( \varphi(t_{n},y(t_{n}))+4\varphi(t_{n+1},y(t_{n+1}))+\varphi(t_{n+2},y(t_{n+2})) \right), $ et avec une prédiction d’Euler explicite pour approcher $u_1$, nous obtenons un nouveau schéma implicite à 3 pas: $\begin{cases} u_0=y(t_0)=y_0,\\ u_1=u_0+h\varphi(t_{0},u_{0}),\\ u_{n+2}=u_{n}+\frac{h}{3} \left( \varphi(t_{n},u_{n})+4\varphi(t_{n+1},u_{n+1})+\varphi(t_{n+2},u_{n+2}) \right)& n=0,1,\dots N-2 \end{cases}$

Pour éviter le calcul implicite de $u_{n+2}$, nous pouvons utiliser une technique de type predictor-corrector et remplacer le $u_{n+2}$ dans le terme $\varphi(t_{n+2},u_{n+2})$ par exemple par $\tilde u_{n+2}=u_{n+1}+h \varphi(t_{n+1},u_{n+1})$. Nous avons construit ainsi un nouveau schéma explicite à 3 pas: $\begin{cases} u_0=y(t_0)=y_0,\\ u_1=u_0+h\varphi(t_{0},u_{0}),\\ u_{n+2}=u_{n}+\frac{h}{3} \left( \varphi(t_{n},u_{n})+4\varphi(t_{n+1},u_{n+1})+\varphi(t_{n+2},u_{n+1}+h \varphi(t_{n+1},u_{n+1}) \right)& n=0,1,\dots N-2 \end{cases}$

Exercice : Construction du schéma de Boole-Moulton

La formule d’intégration de Boole (parfois appelée Bode’s rule suite à une erreur historique de typographie) est une méthode de Newton-Cotes fermée. Elle approche l’intégrale d’une fonction $\varphi$ sur l’intervalle $[t_{n-3}, t_{n+1}]$ par l’intégrale du polynôme d’interpolation de degré 4 passant par les points équirépartis : $ (t_{n+1}, \phi_{n+1}), (t_{n}, \phi_{n}), (t_{n-1}, \phi_{n-1}), (t_{n-2}, \phi_{n-2}), (t_{n-3}, \phi_{n-3}). $

Objectif : en utilisant cette formule de quadrature, construire une méthode multi-pas linéaire pour résoudre l’équation différentielle $y'(t) = \varphi(t, y(t))$. Cette méthode est connue sous le nom de méthode de Boole-Moulton (1960).

Suggestion : utilisez la bibliothèque sympy pour calculer le polynôme d’interpolation et intégrer ce dernier afin d’extraire les coefficients exacts du schéma.

Documentation MathWorld : Boole’s Rule

Le principe des méthodes multi-pas linéaires consiste à intégrer l’équation différentielle $y'(t) = \varphi(t, y(t))$ sur un intervalle de largeur $4h$ et à approcher l’intégrale par une formule de quadrature interpolatoire. En intégrant l’EDO entre $t_{n-3}$ et $t_{n+1}$ on obtient: $ y(t_{n+1}) = y(t_{n-3}) + \int_{t_{n-3}}^{t_{n+1}} \varphi(t, y(t)) \mathrm{d}t $

Étape 1 : Interpolation polynomiale. On remplace la fonction inconnue $\varphi(t, y(t))$ sous l’intégrale par son polynôme d’interpolation de Lagrange $P(t)$ s’appuyant sur les points $\{t_{n+1}, t_n, t_{n-1}, t_{n-2}, t_{n-3}\}$. Astuce calculatoire : pour simplifier le calcul des coefficients du polynôme d’interpolation, on effectue le changement de variable affine $t = t_n + kh$ : $ \int_{t_{n-3}}^{t_{n+1}} \varphi(t, y(t)) \mathrm{d}t \approx \int_{t_{n-3}}^{t_{n+1}} P(t) \mathrm{d}t = h\int_{-3}^{1} P(k) \mathrm{d}k $
Étape 3 : Quadrature de Boole. L’intégration du polynôme conduit à la formule de quadrature de Boole : $ \int_{t_{n-3}}^{t_{n+1}} P(t) \mathrm{d}t = \frac{2h}{45} \left( 7\varphi_{n+1} + 32\varphi_{n} + 12\varphi_{n-1} + 32\varphi_{n-2} + 7\varphi_{n-3} \right) $
Étape 4 : Schéma numérique final. En remplaçant l’intégrale dans l’équation initiale, on obtient le schéma de Boole-Moulton : $ u_{n+1} = u_{n-3} + \frac{2h}{45} \left( 7\varphi(t_{n+1}, u_{n+1}) + 32\varphi(t_{n}, u_{n}) + 12\varphi(t_{n-1}, u_{n-1}) + 32\varphi(t_{n-2}, u_{n-2}) + 7\varphi(t_{n-3}, u_{n-3}) \right). $

Code

# 1. Variables symboliques
t, h, t_n = sp.symbols('t h t_n')
phi_np1, phi_n, phi_nm1, phi_nm2, phi_nm3 = sp.symbols('phi_{n+1} phi_n phi_{n-1} phi_{n-2} phi_{n-3}')
u_np1, u_nm3 = sp.symbols('u_{n+1} u_{n-3}')

# 2. Points d'interpolation (t_i, phi_i)
points = [
    (t_n + h,   phi_np1),
    (t_n,       phi_n),
    (t_n - h,   phi_nm1),
    (t_n - 2*h, phi_nm2),
    (t_n - 3*h, phi_nm3)
]

# 3. Calcul du polynôme d'interpolation P(t)
P = sp.interpolate(points, t)

# 4. Calcul de l'intégrale de t_{n-3} à t_{n+1} de P(t)
integrale = sp.integrate(P, (t, t_n - 3*h, t_n + h))
integrale = sp.simplify(integrale)

# 5. Construction et affichage du schéma final
# u_{n+1} = u_{n-3} + Intégrale
schema = sp.Eq(u_np1, u_nm3 + integrale)

print("Schéma de Boole-Moulton obtenu :")
display(schema)

Schéma de Boole-Moulton obtenu :

$\displaystyle u_{n+1} = \frac{2 h \left(32 \phi_{n} + 7 phi_{n+1} + 12 phi_{n-1} + 32 phi_{n-2} + 7 phi_{n-3}\right)}{45} + u_{n-3}$

BONUS : on généralise cette approche pour construire des méthodes multi-pas linéaires en utilisant des polynômes d’interpolation quelconques et en intégrant ces polynômes sur des intervalles quelconques.

Code

def calculer_multipas(indices_points, a_shift, b_shift):

    t, h, t_n = sp.symbols('t h t_n')
    
    # 1. Création des symboles
    t_syms = [sp.Symbol(f't_{{n{f"+{k}" if k>0 else "" if k==0 else k}}}') for k in indices_points]
    phi_syms = [sp.Symbol(fr'\phi_{{n{f"+{k}" if k>0 else "" if k==0 else k}}}') for k in indices_points]
    
    # 2. Interpolation
    subs_t = {ts: t_n + k*h for ts, k in zip(t_syms, indices_points)}
    pts_calcul = [(ts.subs(subs_t), phis) for ts, phis in zip(t_syms, phi_syms)]
    polynome = sp.interpolate(pts_calcul, t)
    
    # 3. Quadrature
    borne_a, borne_b = t_n + a_shift * h, t_n + b_shift * h
    integrale = sp.integrate(polynome, (t, borne_a, borne_b))
    quadrature = sp.simplify(integrale)
    
    # 4. Coefficients du schéma
    expr_over_h = sp.simplify(quadrature / h)
    coeffs = [expr_over_h.coeff(p) for p in phi_syms]
    denoms = [sp.fraction(c)[1] for c in coeffs if c != 0]
    lcm_denom = sp.lcm(denoms) if denoms else 1
    
    termes_somme = [sp.simplify(c * lcm_denom) * p for c, p in zip(coeffs, phi_syms)]
    somme_interne = sp.Add(*termes_somme, evaluate=False)
    
    # 5. Retour des résultats structurés
    return {
        "t": t, "h": h, "t_n": t_n,
        "t_syms": t_syms, "phi_syms": phi_syms,
        "polynome": polynome.simplify(),
        "quadrature": quadrature,
        "lcm_denom": lcm_denom,
        "somme_interne": somme_interne,
        "indices": {"points": indices_points, "a": a_shift, "b": b_shift}
    }


def afficher_etude_schema(nom, data):
    
    format_idx = lambda k: f"+{k}" if k > 0 else f"{k}" if k < 0 else ""

    # Extraction des variables utiles
    idx = data["indices"]
    u_start = sp.Symbol(f'u_{{n{format_idx(idx["a"])}}}')
    u_target = sp.Symbol(f'u_{{n{format_idx(idx["b"])}}}')
    
    # Construction de l'équation du schéma
    schema = sp.Eq(u_target, u_start + sp.Mul(data["h"]/data["lcm_denom"], data["somme_interne"], evaluate=False), evaluate=False)
    
    # Points d'interpolation en LaTeX
    latex_points = ", ".join([rf"({sp.latex(ts)}, {sp.latex(phis)})" for ts, phis in zip(data["t_syms"], data["phi_syms"])])

    # Polynôme d'interpolation en forme canonique
    poly_canonique = sp.expand(data["polynome"])

    # --- Rendu ---
    display(Markdown(f"Étude du schéma : {nom}"))
    
    display(Markdown(r"- Points d'interpolation $(t_i, \phi_i)$"))
    display(Markdown(f"\({latex_points}\)"))
    
    display(Markdown("- Polynôme d'interpolation $P(t)$ (factorisé)"))
    display(poly_canonique)
    
    display(Markdown(f"- Formule de quadrature (Intégration de $t_{{n{format_idx(idx['a'])}}}$ à $t_{{n{format_idx(idx['b'])}}}$)"))
    symb_int = sp.Integral(sp.Symbol('P(t)'), (data["t"], data["t_n"] + idx['a']*data["h"], data["t_n"] + idx['b']*data["h"]))
    display(sp.Eq(symb_int, data["quadrature"]))
    
    display(Markdown("- Schéma EDO final"))
    display(schema)

Code

res_boole = calculer_multipas(indices_points=[1, 0, -1, -2, -3], a_shift=-3, b_shift=1)
afficher_etude_schema("Boole-Moulton", res_boole)

Étude du schéma : Boole-Moulton

Points d’interpolation $(t_i, \phi_i)$

$(t_{n+1}, \phi_{n+1}), (t_{n}, \phi_{n}), (t_{n-1}, \phi_{n-1}), (t_{n-2}, \phi_{n-2}), (t_{n-3}, \phi_{n-3})$

Polynôme d’interpolation $P(t)$ (factorisé)

\(\displaystyle \frac{\phi_{n+1} t}{4 h} - \frac{\phi_{n+1} t_{n}}{4 h} + \frac{11 \phi_{n+1} t^{2}}{24 h^{2}} - \frac{11 \phi_{n+1} t t_{n}}{12 h^{2}} + \frac{11 \phi_{n+1} t_{n}^{2}}{24 h^{2}} + \frac{\phi_{n+1} t^{3}}{4 h^{3}} - \frac{3 \phi_{n+1} t^{2} t_{n}}{4 h^{3}} + \frac{3 \phi_{n+1} t t_{n}^{2}}{4 h^{3}} - \frac{\phi_{n+1} t_{n}^{3}}{4 h^{3}} + \frac{\phi_{n+1} t^{4}}{24 h^{4}} - \frac{\phi_{n+1} t^{3} t_{n}}{6 h^{4}} + \frac{\phi_{n+1} t^{2} t_{n}^{2}}{4 h^{4}} - \frac{\phi_{n+1} t t_{n}^{3}}{6 h^{4}} + \frac{\phi_{n+1} t_{n}^{4}}{24 h^{4}} - \frac{3 \phi_{n-1} t}{2 h} + \frac{3 \phi_{n-1} t_{n}}{2 h} + \frac{\phi_{n-1} t^{2}}{4 h^{2}} - \frac{\phi_{n-1} t t_{n}}{2 h^{2}} + \frac{\phi_{n-1} t_{n}^{2}}{4 h^{2}} + \frac{\phi_{n-1} t^{3}}{h^{3}} - \frac{3 \phi_{n-1} t^{2} t_{n}}{h^{3}} + \frac{3 \phi_{n-1} t t_{n}^{2}}{h^{3}} - \frac{\phi_{n-1} t_{n}^{3}}{h^{3}} + \frac{\phi_{n-1} t^{4}}{4 h^{4}} - \frac{\phi_{n-1} t^{3} t_{n}}{h^{4}} + \frac{3 \phi_{n-1} t^{2} t_{n}^{2}}{2 h^{4}} - \frac{\phi_{n-1} t t_{n}^{3}}{h^{4}} + \frac{\phi_{n-1} t_{n}^{4}}{4 h^{4}} + \frac{\phi_{n-2} t}{2 h} - \frac{\phi_{n-2} t_{n}}{2 h} + \frac{\phi_{n-2} t^{2}}{6 h^{2}} - \frac{\phi_{n-2} t t_{n}}{3 h^{2}} + \frac{\phi_{n-2} t_{n}^{2}}{6 h^{2}} - \frac{\phi_{n-2} t^{3}}{2 h^{3}} + \frac{3 \phi_{n-2} t^{2} t_{n}}{2 h^{3}} - \frac{3 \phi_{n-2} t t_{n}^{2}}{2 h^{3}} + \frac{\phi_{n-2} t_{n}^{3}}{2 h^{3}} - \frac{\phi_{n-2} t^{4}}{6 h^{4}} + \frac{2 \phi_{n-2} t^{3} t_{n}}{3 h^{4}} - \frac{\phi_{n-2} t^{2} t_{n}^{2}}{h^{4}} + \frac{2 \phi_{n-2} t t_{n}^{3}}{3 h^{4}} - \frac{\phi_{n-2} t_{n}^{4}}{6 h^{4}} - \frac{\phi_{n-3} t}{12 h} + \frac{\phi_{n-3} t_{n}}{12 h} - \frac{\phi_{n-3} t^{2}}{24 h^{2}} + \frac{\phi_{n-3} t t_{n}}{12 h^{2}} - \frac{\phi_{n-3} t_{n}^{2}}{24 h^{2}} + \frac{\phi_{n-3} t^{3}}{12 h^{3}} - \frac{\phi_{n-3} t^{2} t_{n}}{4 h^{3}} + \frac{\phi_{n-3} t t_{n}^{2}}{4 h^{3}} - \frac{\phi_{n-3} t_{n}^{3}}{12 h^{3}} + \frac{\phi_{n-3} t^{4}}{24 h^{4}} - \frac{\phi_{n-3} t^{3} t_{n}}{6 h^{4}} + \frac{\phi_{n-3} t^{2} t_{n}^{2}}{4 h^{4}} - \frac{\phi_{n-3} t t_{n}^{3}}{6 h^{4}} + \frac{\phi_{n-3} t_{n}^{4}}{24 h^{4}} + \phi_{n} + \frac{5 \phi_{n} t}{6 h} - \frac{5 \phi_{n} t_{n}}{6 h} - \frac{5 \phi_{n} t^{2}}{6 h^{2}} + \frac{5 \phi_{n} t t_{n}}{3 h^{2}} - \frac{5 \phi_{n} t_{n}^{2}}{6 h^{2}} - \frac{5 \phi_{n} t^{3}}{6 h^{3}} + \frac{5 \phi_{n} t^{2} t_{n}}{2 h^{3}} - \frac{5 \phi_{n} t t_{n}^{2}}{2 h^{3}} + \frac{5 \phi_{n} t_{n}^{3}}{6 h^{3}} - \frac{\phi_{n} t^{4}}{6 h^{4}} + \frac{2 \phi_{n} t^{3} t_{n}}{3 h^{4}} - \frac{\phi_{n} t^{2} t_{n}^{2}}{h^{4}} + \frac{2 \phi_{n} t t_{n}^{3}}{3 h^{4}} - \frac{\phi_{n} t_{n}^{4}}{6 h^{4}}\)

Formule de quadrature (Intégration de $t_{n-3}$ à $t_{n+1}$)

$\displaystyle \int\limits_{- 3 h + t_{n}}^{h + t_{n}} P(t)\, dt = \frac{2 h \left(7 \phi_{n+1} + 12 \phi_{n-1} + 32 \phi_{n-2} + 7 \phi_{n-3} + 32 \phi_{n}\right)}{45}$

Schéma EDO final

$\displaystyle u_{n+1} = \frac{h}{45} \cdot \left(14 \phi_{n+1} + 24 \phi_{n-1} + 64 \phi_{n-2} + 14 \phi_{n-3} + 64 \phi_{n}\right) + u_{n-3}$

PARTIE III : 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\} $

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

p = 1
bm1 = 0
a0 = 1
a1 = 0
a2 = 0
b0 = sp.Rational(23,12)
b1 = sp.Rational(-16,12)
b2 = sp.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

p = 1
bm1 = sp.Rational(9,24)
a0 = 1
a1 = 0
a2 = 0
b0 = sp.Rational(19,24)
b1 = sp.Rational(-5,24)
b2 = sp.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

# 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

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

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

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

# 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



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

Exercice : Construction d’un schéma Predictor-Corrector

Écrire le schéma predictor-corrector basé sur AM-3 AB-3. Attention à bien initialiser la suite récurrente.

AB-1 $u_{n+1}=u_n+h\varphi(t_n,u_n)$
AB-2 $u_{n+1}=u_n+\frac{h}{2}\left(3\varphi(t_n,u_n)-\varphi(t_{n-1},u_{n-1})\right)$
AB-3 $u_{n+1}=u_n+\frac{h}{12}\left(23\varphi(t_n,u_n)-16\varphi(t_{n-1},u_{n-1})+5\varphi(t_{n-2},u_{n-2})\right)$
AM-3 $u_{n+1}=u_n+\frac{h}{24}\left(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})\right)$

donc $\begin{cases} u_0=y_0,\\ u_1=u_0+h\varphi(t_0,u_0)\\ u_{2}=u_1+\frac{h}{2}\left(3\varphi(t_1,u_1)-\varphi(t_0,u_0)\right) \\ \\ \tilde u=u_n+\frac{h}{12}\left(23\varphi(t_n,u_n)-16\varphi(t_{n-1},u_{n-1})+5\varphi(t_{n-2},u_{n-2})\right)\\ u_{n+1}=u_n+\frac{h}{24}\left(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}\right) \end{cases} $ AM-3 est un schéma implicite à 3 pas d’ordre 4,
AB-3 est un schéma explicite à 3 pas d’ordre 3,
le schéma predictor-corrector est donc un schéma à 3 pas d’ordre 4.