Code
import matplotlib.pyplot as plt
import numpy as np
plt.rcParams.update({'font.size': 12})
from scipy.optimize import fsolve
import sympy as sym
from IPython.display import MarkdownPropriétés d’un schéma sur un intervalle borné
Quelles sont les propriétés d’un schéma numérique qui nous permettent de garantir que les solutions approchées que nous obtenons sont de « bonne qualité » ?
Nous avons déjà vu que les problèmes de Cauchy doivent être mathématiquement bien posés avant d’essayer d’en approcher la solution numériquement. Nous avons aussi vu que certains problèmes de Cauchy sont plus ou moins « difficiles » à résoudre numériquement : cela fait intervenir les notions de
- problème numériquement bien posé, i.e. pas trop sensible aux perturbations
- raideur, liée à la présence de plusieurs échelles temporelles dans le problème.
Dans ce chapitre, nous allons introduire les propriétés d’un schéma numérique qui nous permettent de garantir que les solutions approchées que nous obtenons sont de « bonne qualité ». Nous allons notamment introduire les notions de convergence, consistance et stabilité d’un schéma numérique, et montrer comment elles sont liées entre elles.
Convergence et ordre de convergence
Définition: schéma convergent et ordre de convergence
Une méthode numérique est convergente si
\( |e_n| := |y_n-u_n| \xrightarrow[h\to0]{}0 \qquad\forall n=0,\dots,N \)
Si \(|e_n| = \mathcal{O}(h^{\omega})\) pour \(\omega > 0\), on dit que la convergence de la méthode est d’ordre \(\omega\).
Rappel : \(f(x) = \mathcal{O}(g(x))\) lorsque \(x\to \infty\) signifie que il existe \(C>0\) tel que \(|f(x)|\le C\times |g(x)|\) asymptotiquement.
Pour étudier la convergence d’une méthode numérique, on exprime l’erreur \(e_{n+1}\) au pas \(t_{n+1}\) sous la forme :
\(\begin{aligned} \boxed{e_{n+1}} & := y_{n+1} - u_{n+1} \\ &= (y_{n+1} - u_{n+1}^*) + (u_{n+1}^* - u_{n+1}), \\ &=\boxed{h\tau_{n+1}(h)} + \boxed{u_{n+1}^* - u_{n+1}}, \qquad \boxed{\tau_{n+1}(h) := \frac{y_{n+1} - u_{n+1}^*}{h}} \end{aligned}\)
\(u_{n+1}^*\) désignant la valeur obtenue en appliquant le schéma numérique à la solution exacte au temps \(t_{n+1}\).
Le terme \(h\tau_{n+1}(h) = y_{n+1} - u_{n+1}^*\) représente l’erreur engendrée par une seule itération du schéma (erreur locale de consistance).
Le terme \(u_{n+1}^* - u_{n+1}\) représente la propagation de \(t_{n}\) à \(t_{n+1}\) de l’erreur accumulée au temps précédent \(t_{n}\).
Si les deux termes \((y_{n+1} - u_{n+1}^*)\) et \((u_{n+1}^* - u_{n+1})\) tendent vers zéro lorsque \(h \to 0\), alors la méthode est convergente.
Exemple : convergence de la méthode d’Euler explicite
Considérons le problème de Cauchy \(y'(t)=\varphi(t,y(y))\) pour \(t\in[t_0;T]\) et \(y(t_0)=y_0\). Supposons que la fonction \(\varphi\) soit lipschitzienne de constante \(L\) par rapport à sa deuxième variable. Nous allons montrer que la méthode d’Euler explicite
\( \begin{cases} u_0=y_0,\\ u_{n+1}=u_n+h\varphi(t_n,u_n),&n=0,\dots,N-1 \end{cases} \)
est convergente d’ordre 1.
Pour cela, on étudie séparément l’erreur de consistance et l’accumulation de ces erreurs. On en déduit ensuite la convergence.
Terme \(h\tau_{n+1}(h) := y_{n+1} - u_{n+1}^*\) (erreur de consistance).
Il représente l’erreur engendrée par une seule itération de la méthode d’Euler explicite : \(\begin{aligned}
y_{n+1}-u_{n+1}^*
&=y_{n+1}-\Big(y_{n} + h\varphi(t_{n}, y_{n})\Big)
\\
&=y(t_{n+1})-y(t_{n}) - h y'(t_n).
\end{aligned}\)
En supposant que la dérivée seconde de \(y\) existe et est continue, on sait qu’il existe un point \(\eta_n\) de l’intervalle \(]t_n;t_{n+1}[\) tel que \( y(t_{n+1})=y(t_n)+hy'(t_n)+\frac{h^2}{2}y''(\eta_n) \) (développement de Taylor de \(y\) au voisinage de \(t_n\) évalué en \(t_{n+1}\) avec reste de Lagrange). On a alors \( y_{n+1}-u_{n+1}^* =y(t_{n+1})-y(t_{n}) - h y'(t_n)=\frac{h^2}{2}y''(\eta_n) \le \frac{M}{2}h^2 \xrightarrow[h\to0]{}0 \) où \(M\equiv\max_{t\in [t_0,T]}|y''(t)|\).
L’erreur de troncature locale de la méthode d’Euler explicite est donc de la forme \( \tau_{n+1}(h)=\frac{y_{n+1} -u_{n+1}^*}{h}=\mathcal{O}(h). \)
Terme \(u_{n+1}^* - u_{n+1}\) (accumulation des erreurs).
Il représente la propagation de \(t_{n}\) à \(t_{n+1}\) de l’erreur accumulée au temps précédent \(t_{n}\). On a
\(\begin{aligned} u_{n+1}^* - u_{n+1} &= \Big(y_n + h \varphi(t_n, y_n)\Big)-\Big(u_n + h \varphi(t_n, u_n)\Big) \\ &=e_n+h\Big( \varphi(t_n,y_n)-\varphi(t_n,u_n) \Big). \end{aligned}\)
Comme \(\varphi\) est lipschitzienne par rapport à sa deuxième variable, on a
\(|\varphi(t_n,y_n)-\varphi(t_n,u_n)|\le L|y_n-u_n|=L|e_n|.\)
En factorisant on trouve alors
\( |u_{n+1}^{*} - u_{n+1}|\le (1 + hL)|e_{n}|. \)
Convergence.
Comme \(e_0 = 0\), les relations précédentes donnent
\(\begin{aligned} |e_n| & \le |y_n - u_n^*| + |u_n^* - u_n| \\ &\le h|\tau_n(h)| + (1+hL)|e_{n-1}| \\ &\le h|\tau_n(h)| + (1+hL)\left( h|\tau_{n-1}(h)| + (1+hL)|e_{n-2}| \right)| \\ &\le \left( 1+(1+hL)+\dots+(1+hL)^{n-1} \right)h\tau(h) \qquad\text{avec } \tau(h) = \max_{n=0,\dots,N} |\tau_n(h)| \\ &= \left(\sum_{i=0}^{n-1}(1+hL)^i\right)h\tau(h) \\ &= \frac{(1+hL)^n-1}{hL}h\tau(h) \\ &\le \frac{(e^{hL})^n-1}{hL}h\tau(h) \qquad\text{car }(1+x)\le e^x \\ &= \frac{(e^{hL})^{(t_n-t_0)/h}-1}{L}\tau(h)\qquad\text{car } t_n-t_0=nh \\ &= \frac{e^{L(t_n-t_0)}-1}{L}\tau(h) \\ &\le \frac{e^{L(t_n-t_0)}-1}{L}\frac{M}{2}h=\mathcal{O}(h) \end{aligned}\)
On peut conclure que la méthode d’Euler explicite est convergente d’ordre 1.
Remarque : l’estimation de convergence est obtenue en supposant seulement \(\varphi\) lipschitzienne. On peut établir une meilleure estimation si \(\partial_y\varphi\) existe et est non négative pour tout \(t \in [t_0;T]\) et tout \(y\in\mathbb{R}\). En effet dans ce cas \(\begin{aligned} u_n^* - u_n &= \left(y_{n-1} + h\varphi(t_{n-1}, y_{n-1})\right) - \left( u_{n-1} + h\varphi(t_{n-1},u_{n-1})\right) \\ &= e_{n-1}+h\left( \varphi(t_{n-1},y_{n-1})-\varphi(t_{n-1},u_{n-1}) \right) \\ &= e_{n-1}+h\left( e_{n-1}\partial_y\varphi(t_{n-1},\eta_n) \right) \\ &= \left( 1+h\partial_y\varphi(t_{n-1},\eta_n) \right)e_{n-1} \end{aligned}\) où \(\eta_n\) appartient à l’intervalle dont les extrémités sont \(y_{n-1}\) et \(u_{n-1}\). Ainsi, si \( 0<h<\frac{2}{\max\limits_{t\in[t_0,T]}\partial_y\varphi(t,y(t))} \) alors \( |u_n^* - u_n |\le |e_{n-1}|. \) On en déduit \(|e_n| \le |y_n - u_n^*| + |e_{n-1}| \le nh\tau(h) + |e_0|\) et donc \( |e_n|\le M \frac{h}{2}(t_n-t_0). \) La restriction sur le pas de discrétisation \(h\) est une condition de stabilité, comme on le verra dans la suite.
Exercice
Calculer l’ordre de consistance du schéma d’Euler implicite. Montrer ensuite que la méthode est convergente.
On considère le même problème de Cauchy que précédemment : \(y'(t)=\varphi(t,y(y))\) pour \(t\in[t_0;T]\) et \(y(t_0)=y_0\) où \(\varphi\) est supposée être lipschitzienne de constante \(L\) par rapport à sa deuxième variable : \( |\varphi(t,y_1)-\varphi(t,y_2)|\le L|y_1-y_2| \qquad\forall t\in[t_0;T], \forall y_1,y_2\in\mathbb{R}. \) On suppose que la dérivée seconde de \(y\) existe et est continue.
Le schéma d’Euler implicite est défini par : \( \begin{cases} u_0=y_0,\\ u_{n+1}=u_n+h\varphi(t_{n+1},u_{n+1}),&n=0,\dots,N-1 \end{cases} \)
Soit \(u_{n+1}^* = y_n + h\varphi(t_{n+1}, y_{n+1})\). Calculons l’erreur de troncature locale \(\tau_{n+1}(h)\) : \(\begin{aligned} \tau_{n+1}(h) &\equiv \frac{y_{n+1} - u_{n+1}^*}{h} \\ &= \frac{y(t_{n+1}) - \left(y(t_n) + h\varphi(t_{n+1}, y(t_{n+1}))\right)}{h} \\ &= \frac{y(t_{n+1}) - y(t_n)}{h} - \varphi(t_{n+1}, y(t_{n+1})) \end{aligned}\) Par développement de Taylor de \(y\) au voisinage de \(t_{n+1}\) évalué en \(t_{n}\), il existe un point \(\eta_n\) de l’intervalle \(]t_n;t_{n+1}[\) tel que \(\begin{aligned} y(t_n) & = y(t_{n+1}) - hy'(t_{n+1}) + \frac{h^2}{2}y''(\eta_n) \\ & = y(t_{n+1}) - h\varphi(t_{n+1}, y(t_{n+1})) + \frac{h^2}{2}y''(\eta_n) \end{aligned}\) Comme \(y'(t) = \varphi(t, y(t))\), on obtient alors \( y(t_{n+1}) = y(t_n) + h\varphi(t_{n+1}, y(t_{n+1})) - \frac{h^2}{2}y''(\eta_n). \) On en déduit que l’erreur de troncature locale de la méthode d’Euler implicite est de la forme \(\begin{aligned} \tau_{n+1}(h) & = \frac{y(t_{n}) + h\varphi(t_{n+1}, y(t_{n+1})) - \frac{h^2}{2}y''(\eta_n) - y(t_n)}{h} - \varphi(t_{n+1}, y(t_{n+1})) \\ & = - \frac{h}{2}y''(\eta_n) \end{aligned}\) On a alors, comme pour le schéma d’Euler explicite, \( |\tau_{n+1}(h)| \le \frac{M}{2}h \xrightarrow[h\to0]{}0 \) où \(M\equiv\max_{t\in [t_0,T]}|y''(t)|\) et le schéma d’Euler implicite est consistant d’ordre 1.
La solution exacte vérifie donc \( y(t_{n+1}) = y(t_n) + h\varphi(t_{n+1}, y(t_{n+1})) +h\tau_{n+1}(h). \) En soustrayant l’expression du schéma d’Euler implicite à celle de la solution exacte, on trouve \(\begin{aligned} e_{n+1} & = y(t_{n+1}) - u_{n+1} \\ & = \Big(y(t_n) + h\varphi(t_{n+1}, y(t_{n+1})) +h\tau_{n+1}(h)\Big) - \Big(u_n + h\varphi(t_{n+1}, u_{n+1})\Big) \\ & = e_n + h\Big( \varphi(t_{n+1}, y(t_{n+1})) - \varphi(t_{n+1}, u_{n+1}) \Big) + h\tau_{n+1}(h) \end{aligned} \) Comme \(\varphi\) est lipschitzienne par rapport à sa deuxième variable, \( |\varphi(t_{n+1}, y(t_{n+1})) - \varphi(t_{n+1}, u_{n+1})| \le L|y(t_{n+1}) - u_{n+1}| = L|e_{n+1}|. \) En factorisant on trouve alors \( (1-hL)|e_{n+1}| \le |e_n| + h|\tau_{n+1}(h)|. \) Pour \(h\) suffisamment petit, on a \(1-hL > 0\) et donc \( |e_{n+1}| \le \frac{1}{1-hL} \left(|e_n| + h|\tau_{n+1}(h)|\right). \) Comme \(e_0 = 0\), les relations précédentes donnent \(\begin{aligned} |e_n| & \le |y(t_n) - u_n^*| + |u_n^* - u_n| \\ &\le h|\tau_n(h)| + \frac{1}{1-hL} \left(|e_{n-1}| + h|\tau_{n}(h)|\right) \\ &\le \left(1 + \frac{1}{1-hL} + \dots + \left(\frac{1}{1-hL}\right)^{n-1}\right)h\tau(h) \qquad\text{avec } \tau(h) = \max_{n=0,\dots,N} |\tau_n(h)| µ \\ &= \left(\sum_{i=0}^{n-1}\left(\frac{1}{1-hL}\right)^i\right)h\tau(h) \\ &= \frac{\left(\frac{1}{1-hL}\right)^n-1}{\frac{1}{1-hL}-1}h\tau(h) \\ &= \frac{(1-hL)^{-n}-1}{\frac{hL}{1-hL}}h\tau(h)\\ &= \frac{(1-hL)^{-n}-1}{L}h\tau(h) \qquad\text{car }t_n-t_0=nh \\ &= \frac{e^{L(t_n-t_0)}-1}{L}h\tau(h) \\ &= \frac{e^{L(t_n-t_0)}-1}{L}\frac{M}{2}h=\mathcal{O}(h) \end{aligned}\)
Théorème d’équivalence de Lax-Ritchmyer : consistance + zéro-stabilité = convergence
Un schéma consistant est convergent ssi il est zéro-stable.
Consistance et ordre de consistance
Définition : erreurs de troncature
L’erreur de troncature locale est donnée par :
\( \tau_{n+1}(h) \equiv \frac{y_{n+1} - u_{n+1}^*}{h}. \)
Elle représente (à un facteur \(1/h\) près) l’erreur introduite en insérant la solution exacte dans le schéma numérique.
L’erreur de troncature globale (ou simplement erreur de troncature) est définie comme :
\( \tau(h) = \max_{n=0,\dots,N} |\tau_n(h)|. \)
Elle correspond à la borne supérieure de l’erreur de troncature locale sur l’ensemble des pas de temps.
Définition : Consistance et ordre de consistance
- Une méthode numérique est dite consistante si \(\lim_{h\to0} \tau(h) = 0\).
- Elle est consistante d’ordre \(\omega\) si l’erreur de troncature vérifie \(\tau(h) = \mathcal{O}(h^\omega)\), pour un certain \(\omega \geq 1\).
Remarque : La propriété de consistance est une condition nécessaire pour la convergence. En effet, sans consistance, la méthode introduirait à chaque itération une erreur qui ne tendrait pas vers zéro lorsque \(h\) diminue. L’accumulation de ces erreurs empêcherait alors l’erreur globale de s’annuler lorsque \(h \to 0\).
Remarque : l’ordre de convergence coïncide en général avec l’ordre de son erreur de troncature.
Zéro-stabilité
Prenons le problème de Cauchy, supposé mathématiquement et numériquement bien posé :
\( \begin{cases} y'(t) = \varphi(t,y(t)), &\forall t \in I=]t_0,T[,\\ y(t_0) = y_0. \end{cases} \)
Que se passe-t-il lorsque l’on fixe le temps final \(T\) et que l’on fait tendre le pas \(h\) vers \(0\) ?
La zéro-stabilité garantit que, sur un intervalle borné, de petites perturbations des données entraînent des perturbations contrôlées de la solution numérique lorsque \(h\to0\). Elle garantit que la méthode numérique est peu sensible aux petites perturbations des données. Si la méthode numérique n’est pas zéro-stable, les erreurs d’arrondi sur \(y_0\) et le calcul de \(\varphi(t_n,y_n)\) rendent la solution calculée inutilisable.
Plus précisément, une méthode numérique visant à approcher le problème de Cauchy sur un intervalle fini \(I=[t_0,T]\) est dite zéro-stable si les conditions suivantes sont satisfaites :
\( \exists h_0>0,\ \exists C>0,\ \exists \varepsilon_0>0\ \text{tels que}\ \forall h\in]0,h_0],\ \forall \varepsilon\in]0,\varepsilon_0],\ \text{si } |\varrho_n| \le \varepsilon,\ 0 \le n \le N, \)
alors
\( |z_n - u_n| \le C\varepsilon,\ 0 \le n \le N \)
où :
- \(C\) est une constante indépendante de \(h\),
- \(u_n\) est la solution obtenue par le schéma numérique sans perturbation,
- \(z_n\) est la solution obtenue par le même schéma appliqué au problème perturbé,
- \(\varrho_n\) représente la perturbation au \(n\)-ième pas,
- \(\varepsilon_0\) est la perturbation maximale admissible.
Il est important de noter que \(\varepsilon_0\) doit être suffisamment petit pour que le problème perturbé conserve une solution unique sur l’intervalle d’intégration considéré.
Dans le cas d’une méthode consistante à un pas, on peut démontrer que la zéro-stabilité découle du fait que \(\varphi\) est lipschitzienne par rapport à sa seconde variable. Dans ce contexte, la constante \(C\) dépend de \(\exp((T-t_0)L)\), où \(L\) est la constante de Lipschitz. Toutefois, cette propriété n’est pas nécessairement vérifiée pour d’autres familles de méthodes numériques.
Exemple : zéro-stabilité de la méthode d’Euler explicite
Considérons \(u_n\) comme l’approximation obtenue avec la méthode d’Euler explicite pour le problème initial :
\( \begin{cases} u_0 = y_0, \\ u_{n+1} = u_n + h\varphi(t_n, u_n), & n = 0, \dots, N - 1 \end{cases} \)
et \(z_n\) comme l’approximation obtenue avec la méthode d’Euler explicite pour le problème perturbé :
\( \begin{cases} z_0 = y_0 + \varrho_0, \\ z_{n+1} = z_n + h\left(\varphi(t_n, z_n) + \varrho_n\right), & n = 0, \dots, N - 1 \end{cases} \)
Nous avons alors :
\(\begin{aligned} |u_n - z_n| &\le |u_{n-1} - z_{n-1}| + h|\varphi(t_{n-1}, u_{n-1}) - \varphi(t_{n-1}, z_{n-1})| + h|\varrho_n| \\ &\le (1 + hL)|u_{n-1} - z_{n-1}| + h\varepsilon \\ &\le (1 + hL)^2|u_{n-2} - z_{n-2}| + (1 + (1 + hL))h\varepsilon \\ &\vdots \\ &\le (1 + hL)^n|\varrho_0| + \left(1 + (1 + hL) + \dots + (1 + hL)^{n-1}\right)h\varepsilon \\ &= (1 + hL)^n|\varrho_0| + \left(\sum_{i=0}^{n-1}(1 + hL)^i\right)h\varepsilon \\ &= (1 + hL)^n|\varrho_0| + \frac{(1 + hL)^n - 1}{L}\varepsilon \\ &\le (e^{hL})^n|\varrho_0| + \frac{(e^{hL})^n - 1}{hL}h\varepsilon \qquad\text{car } (1 + x) \le e^x \\ &= e^{nhL}|\varrho_0| + \frac{e^{nhL} - 1}{L}\varepsilon \\ &\le \left(e^{nhL} + \frac{e^{nhL} - 1}{L}\right)\varepsilon \\ &= \frac{(1 + L)e^{L(T - t_0)} - 1}{L}\varepsilon \qquad\text{car } t_n - t_0 = nh \end{aligned}\)
Ainsi, nous avons montré que la méthode d’Euler explicite est zéro-stable.