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é
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
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
# =============================================================================# On doit seulement définir pp =1# p+1 coeffs du schema d'indices 0 à p# =============================================================================# =============================================================================# FONCTION qui calcule xi(i) pour un schéma à q pas# =============================================================================from sympy import symbols, Symbol, latex, factor, solvefrom IPython.display import display, Math, Markdown, Latexdef xi(i, q):# Définition des symboles pour le cas général aa = symbols(f'a_0:{</span>q<span class="sc">}') bb = symbols(f'b_0:{</span>q<span class="sc">}') bm1 = Symbol('b_{-1}') p =len(aa) sa =sum((-j)**i * aa[j] for j inrange(p)) sb = bm1 +sum((-j)**(i-1) * bb[j] for j inrange(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 mainif i ==1: sb += bb[0]return (sa).factor() + (i * sb).factor()# =============================================================================# 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 impliciteq = p +1omega = q +2if q%2==0else q +1# i = 0, ..., omega# Affichage des xi(i)display(Markdown(f"- Pour une méthode à $q={</span>p <span class="op">+</span> <span class="dv">1</span><span class="sc">}$ pas, on affiche $\\xi(i)$ pour $i = 0, \\dots, {</span>omega<span class="sc">}$"))for i inrange(omega +1): display(Markdown("$\\qquad$ "+f"$\\xi({</span>i<span class="sc">}) = {</span>latex(xi(i, q))<span class="sc">}$"))# =================================================================================================================# 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, {</span>omega<span class="sc">}$"))systeme = [xi(i, q)-1for i inrange(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"${</span>k<span class="sc">} = {</span>v<span class="sc">}$"))
Pour une méthode à \(q=2\) pas, on affiche \(\xi(i)\) pour \(i = 0, \dots, 4\)
On résout le système \(\xi(i) = 1\) pour \(i = 0, \dots, 4\)
\(\qquad\) Les coefficients qui vérifient la condition de consistance d’ordre le plus élevé sont :
\(\qquad\)\(a_0 = 0\)
\(\qquad\)\(a_1 = 1\)
\(\qquad\)\(b_0 = 4/3\)
\(\qquad\)\(b_1 = 1/3\)
\(\qquad\)\(b_{-1} = 1/3\)
1.2 Exercices
1.2.1 👣 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})
\)
1.2.1.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}
\displaystyle\sum_{j=0}^p a_j=1,
\\
\displaystyle-\sum_{j=0}^p ja_j+\sum_{j=-1}^p b_j=-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.
1.2.2 👣 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})
\)
1.2.2.1 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)=\sum_{j=0}^p a_j=-4+5=1,
\\
\displaystyle \xi(1)=-\sum_{j=0}^p ja_j+\sum_{j=-1}^p b_j=-\big(0\times a_0+1\times a_1\big)+\big(b_{-1}+b_0+b_1\big)=-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\).
1.2.3 👣 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}
\)
1.2.3.1 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}
\displaystyle\sum_{j=0}^p a_j=-1+1+1=1,
\\
\displaystyle-\sum_{j=0}^p ja_j+\sum_{j=-1}^p b_j=-(0\times(-1)+1\times1+2\times1)+(0+4+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^2-1)=(r-1)^2(r+1)
\) dont les racines sont \(r_0=1\) de multiplicité \(2\) et \(r_1=-1\). La méthode n’est pas zéro-stable car \(\varrho'(r_0)\neq0\), donc la méthode ne converge pas.
Code
import sympy as symsym.var('r')rho = r**3+r**2-r-1sym.factor(rho)
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}
\)
1.2.5 👣 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 la 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}
\)
1.2.5.1 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}
\displaystyle\sum_{j=0}^p a_j=1-\gamma+\gamma=1 \text{ pour tout }\gamma,
\\
\displaystyle-\sum_{j=0}^p ja_j+\sum_{j=-1}^p b_j=-\big(0\times(1-\gamma)+1\times \gamma\big)+\big(\frac{\gamma+3}{4}+0+\frac{3\gamma+1}{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
\\
\varrho'(r_j)\neq0 \text{ si } |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 \(-1< \gamma\le 1\).
Ordre de convergence La méthode est d’ordre au moins 2 pour tout \(\gamma\) car \(
\sum_{j=0}^p (-j)^{2}a_j+2\sum_{j=-1}^p (-j)^{1}b_j
{}=
\big(0^2a_0+1^2a_1\big)-2\big(-1b_{-1}+0b_0+1b_1\big)
{}=
\gamma-2\big(-\frac{\gamma+3}{4}+\frac{3\gamma+1}{4}\big)
{}=1
\)
Pour que la méthode soit d’ordre 3 il faut que \(
1
{}=
\sum_{j=0}^p (-j)^{3}a_j+3\sum_{j=-1}^p (-j)^{2}b_j
{}=
{}-\big(0^3a_0+1^3a_1\big)+3\big((-1)^2b_{-1}+0^2b_0+1^2b_1\big)
{}=
{}-\gamma+3\big(\frac{\gamma+3}{4}+\frac{3\gamma+1}{4}\big)
{}=
2\gamma+3
\)
donc ssi \(\gamma=-1\). Cependant, cette valeur n’appartient pas au domaine de zéro-stabilité.
On conclut que la méthode est convergente d’ordre \(2\) pour \(-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\).
import matplotlib.pyplot as pltplt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})H = []err_mp = []N =10gamma =1/2for k inrange(7): N +=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 = {</span>(np.polyfit(np.log(H),np.log(err_mp), <span class="dv">1</span>)[<span class="dv">0</span>])<span class="sc">}')plt.grid(True)
1.2.6 👣 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].
\) 1. Pour quelles valeurs de \(\gamma\) est-elle consistant e? 1. Est-elle convergente pour ces valeurs? 1. Quel ordre de convergente a-t-elle? 1. 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}
\)
1.2.6.1 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}
\displaystyle\sum_{j=0}^p a_j=\gamma^2+\gamma-1=1,
\\
\displaystyle-\sum_{j=0}^p ja_j+\sum_{j=-1}^p b_j=-\big(0a_0+1a_1\big)+\big(b_{-1}+b_0+b_1\big)=\gamma^2+\gamma-1=1
\end{cases}
\) donc ssi \(\gamma=-2\) ou \(\gamma=1\).
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-a_0r-a_1=
\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
\\
\varrho'(r_j)\neq0 \text{ si } |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 \(
\sum_{j=0}^p (-j)^{2}a_j+2\sum_{j=-1}^p (-j)^{1}b_j
{}=
\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 pas d’ordre 3 pour aucun des deux choix de \(\gamma\) car \(
\sum_{j=0}^p (-j)^{3}a_j+3\sum_{j=-1}^p (-j)^{2}b_j
{}=
\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 inlineimport numpy as npimport matplotlib.pyplot as pltplt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})from scipy.optimize import fsolveimport sympy as symsym.init_printing()from IPython.display import display, Math
1.2.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(\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}
\)
1.2.7.1 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)=\displaystyle\sum_{j=0}^p a_j=2\gamma^2-8\gamma+7=1,
\\
\xi(1)=\displaystyle-\sum_{j=0}^p ja_j+\sum_{j=-1}^p b_j=-\big(0a_0+1a_1\big)+\big(b_{-1}+b_0+b_1\big)=\frac{5}{6}\gamma-\frac{19}{6}=1
\end{cases}
\) donc ssi \(\gamma=3\). Vérifions ces calculs ci-dessous:
Code
%reset -f%matplotlib inlineimport numpy as npimport matplotlib.pyplot as pltplt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})from scipy.optimize import fsolveimport sympy as symsym.init_printing()from IPython.display import display, Math
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
\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{multiplicité de $r_j$ supérieur à 1} \implies |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\): \(
\xi(2)=
\sum_{j=0}^p (-j)^{2}a_j+2\sum_{j=-1}^p (-j)^{1}b_j
{}=
\big(0^2a_0+1^2a_1\big)+2\big(-(-1)b_{-1}+0b_0+1b_1\big)
\stackrel{\gamma=3}{=}
1
\)
La méthode n’est pas d’ordre 3 car \(\xi(3)\neq1\): \(
\xi(3)=
\sum_{j=0}^p (-j)^{3}a_j+3\sum_{j=-1}^p (-j)^{2}b_j
{}=
{}-a_1+3\big(b_{-1}+b_1\big)
\stackrel{\gamma=3}{\neq}1
\)
On définit la solution exacte (calculée en utilisant le module sympy) pour estimer les erreurs:
2 PARTIE II : Schémas Prédicteur-Correcteur basés sur deux schémas multi-pas linéaires
2.1 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 :
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
2.2.1 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}\).
2.2.1.1 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:
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.
2.2.1.2 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:
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.
2.2.1.3 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:
Multipas P ordre = 2.99
Multipas C ordre = 3.86
Multipas PC ordre = 4.00
2.2.2 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}\).
2.2.2.1 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\)).
