Propriétés d’un schéma : convergence, consistance, stabilité(s)

La 🏺boîte de Pandore🏺 des approximations

Selon Hésiode, lorsque Prométhée déroba le feu du ciel, Zeus, le roi des dieux, se vengea en offrant Pandore à Épiméthée, le frère de Prométhée. Pandore, poussée par la curiosité, ouvrit une jarre qui lui avait été confiée, libérant ainsi la maladie, la mort et d’innombrables maux dans le monde. En refermant précipitamment le récipient, elle découvrit qu’une seule chose y demeurait : l’espoir. C’est de cette légende qu’est née l’expression « ouvrir la boîte de Pandore », symbolisant le fait de déclencher une série de problèmes imprévus et incontrôlables.

Pandora

De manière similaire, le domaine des approximations numériques peut être perçu comme une boîte de Pandore. En s’aventurant dans ce champ, on se confronte à une multitude de défis, tels que la stabilité, les problèmes raides et les systèmes mal conditionnés. Chacun de ces aspects représente un défi à surmonter pour obtenir des résultats fiables et précis.

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

from scipy.optimize import fsolve

import sympy as sym
sym.init_printing()

from IPython.display import Markdown

Convergence et ordre de convergence

Définition: schéma convergent et ordre de convergence

Une méthode numérique est convergente si

\( |y_n-u_n|\le C(h)\xrightarrow[h\to0]{}0 \qquad\forall n=0,\dots,N \)

Si \(C(h) = \mathcal{O}(h^p)\) pour \(p > 0\), on dit que la convergence de la méthode est d’ordre \(p\).

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} e_{n+1} &\equiv y_{n+1} - u_{n+1} \\ &= (y_{n+1} - u_{n+1}^*) + (u_{n+1}^* - u_{n+1}), \\ &=h\tau_{n+1}(h) + (u_{n+1}^* - u_{n+1}), \qquad \tau_{n+1}(h) \equiv \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 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)|\).

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) \\ &= \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) \\ &= \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{align*} 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{align*}\) 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.

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 \(p\) si l’erreur de troncature vérifie \(\tau(h) = \mathcal{O}(h^p)\), pour un certain \(p \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.

Deux notions de stabilité : zéro-stabilité et A-stabilité

De manière générale, un schéma numérique est considéré comme stable s’il permet de maîtriser la solution lorsqu’on perturbe les données. Il existe plusieurs notions de stabilité, nous allons en détailler deux : la zéro-stabilité et l’A-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} \)

Deux questions naturelles se posent alors :

• Que se passe-t-il lorsque l’on fixe le temps final \(T\) et que l’on fait tendre le pas \(h\) vers \(0\) ?
  
• Que se passe-t-il lorsque l’on fixe le pas \(h>0\) et que l’on fait tendre \(T\) vers l’infini ?

Dans ces deux cas, le nombre de nœuds tend vers l’infini, mais les enjeux diffèrent : dans le premier cas, on s’intéresse à l’erreur en chaque point, tandis que dans le second, il s’agit du comportement asymptotique de la solution approchée.

Ces deux problématiques font intervenir deux notions de stabilité :

• Zéro-stabilité
  Elle concerne la résolution du problème de Cauchy sur des intervalles bornés (le nombre \(N\) de sous-intervalles ne tend vers l’infini que lorsque \(h\) tend vers zéro).
  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.

• A-stabilité
  Dans certaines situations, le problème de Cauchy doit être résolu sur des intervalles de temps très longs, voire infinis. Ici, même pour un \(h\) fixé, le nombre \(N\) de sous-intervalles tend vers l’infini. On recherche alors des méthodes capables d’approcher la solution sur des intervalles de temps arbitrairement grands.
  Idéalement, on souhaiterait utiliser des pas de temps \(h\) suffisamment grands (une méthode est dite inconditionnellement A-stable si l’on peut choisir \(h>0\) sans restriction ; sinon, une condition de stabilité sur \(h\) déterminera le temps de calcul).

Étude de la zéro-stabilité (dans \(\mathbb{R}\))

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

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{align} |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{align}\)

Ainsi, nous avons montré que la méthode d’Euler explicite est zéro-stable.

Étude de la A-stabilité (dans \(\mathbb{R}\)) pour les schémas à un pas

Dans la section précédente, nous avons étudié la résolution du problème de Cauchy sur des intervalles bornés. Dans ce contexte, le nombre \(N\) de sous-intervalles ne tend vers l’infini que lorsque \(h\) tend vers zéro.
Cependant, il existe de nombreuses situations où le problème de Cauchy doit être résolu sur des intervalles de temps très longs, voire infinis. Dans ces cas-là, même pour un \(h\) fixé, \(N\) tend vers l’infini, rendant la notion de zéro-stabilité inapplicable puisque \(T - t_0 \to +\infty\).
On cherche donc des méthodes capables d’approcher la solution sur des intervalles de temps arbitrairement grands, y compris avec des pas de temps \(h\) relativement grands.

Définition : A-stabilité dans \(\mathbb{R}\)

Soit \(\beta > 0\) un réel positif et considérons le problème de Cauchy (modèle simplifié = toy-model) :

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

où \(y_0 \neq 0\) est une valeur initiale donnée. La solution exacte est \(y(t) = y_0 e^{-\beta t}\), d’où :

\(\lim_{t \to +\infty} y(t) = 0.\)

Considérons un pas de temps fixé \(h > 0\) et définissons \(t_n = nh\) pour \(n \in \mathbb{N}\). Notons \(u_n \approx y(t_n)\) une approximation de la solution \(y\) au temps \(t_n\).

On dit que le schéma numérique est A-stable si, pour tout \(\beta>0\), la suite satisfait

\(\lim_{n\to+\infty} |u_n| =0.\)

Lorsque cette propriété est vérifiée pour tout pas de temps \(h > 0\), le schéma est dit inconditionnellement A-stable (ou absolument stable). Si elle n’est vérifiée que sous certaines restrictions sur \(h\), on dit que le schéma est A-stable sous condition.

Remarques :

• Des conclusions similaires restent valables lorsque \(\beta\) est un nombre complexe ou une fonction positive de \(t\) (ce cas sera traité dans le chapitre suivant).
  
• À première vue, l’exigence d’une stabilité absolue pour une méthode numérique appliquée à un autre problème que le toy-model peut sembler superflue. Cependant, il est possible de démontrer que si une méthode est absolument stable pour le toy-model, alors, lorsqu’elle est appliquée à un problème plus général, l’erreur de perturbation (c’est-à-dire la différence absolue entre la solution perturbée et la solution non perturbée) reste uniformément bornée par rapport à \(h\).
  En d’autres termes, une méthode absolument stable selon la définition appliquée au toy-model garantit un contrôle des perturbations même dans des contextes plus généraux.

Exemples et exercices sur la A-stabilité pour les schémas à un pas et à deux pas

Nous allons maintenant calculer les conditions de A-stabilité (et de positivité) de différentes méthodes numériques d’approximation à un pas puis à deux pas. On verra dans le chapitre suivant comment généraliser ces résultats à des suites de plus de pas.

\( \def\arraystretch{1.8} \begin{array}{|c|c|c|c|c|c|} \hline \text{Acronyme} & \text{Schéma} & \text{Type} & \text{Ordre} & \text{A-stabilité} & \text{Positivité} \\ \hline \text{EE} & u_{n+1} = u_n + h\varphi_n & \text{Explicite} & 1 & 0 < h < \frac{2}{\beta} & 0 < h < \frac{1}{\beta} \\ \text{EI} & u_{n+1} = u_n + h\varphi_{n+1} & \text{Implicite} & 1 & \forall h > 0 & \forall h > 0 \\ \text{EM} & u_{n+1} = u_n + h\varphi\left(t_n + \frac{h}{2}, u_n + \frac{h}{2}\varphi_n\right) & \text{Explicite} & 2 & 0 < h < \frac{2}{\beta} & 0 < h < \frac{2}{\beta} \\ \text{CN} & u_{n+1} = u_n + \frac{h}{2}\left(\varphi_n + \varphi_{n+1}\right) & \text{Implicite} & 2 & \forall h > 0 & 0 < h < \frac{2}{\beta} \\ \text{H} & u_{n+1} = u_n + \frac{h}{2}\left(\varphi_n + \varphi(t_{n+1}, u_n + h\varphi_n)\right) & \text{Explicite} & 2 & 0 < h < \frac{2}{\beta} & 0 < h < \frac{2}{\beta} \\ \text{AB2} & u_{n+1} = u_n + \frac{h}{2}\left(3\varphi_n - \varphi_{n-1}\right) & \text{Explicite} & 2 & 0 < h < \frac{1}{\beta} & 0 < h < \frac{1}{\beta} \\ \text{AM2} & u_{n+1} = u_n + \frac{h}{12}\left(8\varphi_n - \varphi_{n-1} + 5\varphi_{n+1}\right) & \text{Implicite} & 2 & 0 < h < \frac{6}{\beta} & 0 < h < \frac{6}{\beta} \\ \text{BDF2} & u_{n+1} = \frac{4}{3}u_n - \frac{1}{3}u_{n-1} + \frac{2}{3}h\varphi_{n+1} & \text{Implicite} & 2 & \forall h > 0 & \forall h > 0 \\ \hline \end{array} \)

Schémas à un pas

AvertissementRemarques sur les suites à 1 pas

Considérons une relation de récurrence de la forme

\( u_{n+1} = q u_n. \)

Par récurrence, on a \( u_n = q^n u_0. \)

La suite tend vers 0 si et seulement si \(|q|<1\).

Exemple : étude de l’A-stabilité du schéma d’Euler explicite

Le schéma d’Euler progressif appliqué au toy-model (\(\varphi(t,y)=-\beta y\)) s’écrit

\( u_{n+1}=(1-\beta h)u_n, \qquad n=0,1,2,\dots \)

et par suite

\( u_{n}=(1-\beta h)^nu_0, \qquad n=0,1,2,\dots \)

Par conséquente, \(\lim\limits_{n\to+\infty}u_n=0\) si et seulement si

\( \left|1-\beta h\right|<1, \)

ce qui a pour effet de limiter \(h\) à

\( \boxed{h<\frac{2}{\beta}.} \)

Cette condition de A-stabilité limite le pas \(h\) d’avance en \(t\) lorsqu’on utilise le schéma d’Euler progressif.

Remarques :

• Cette convergence est monotone ssi \(1-\beta h>0\), i.e. ssi \(h<\frac{1}{\beta}\) (bien plus restrictif que la condition de A-stabilité).
• Notons que si \(1-\beta h>1\) alors \(u_n\) tend vers \(+\infty\) lorsque \(t\) tend vers l’infini et si \(1-\beta h<-1\) alors \(u_n\) tend vers l’infini en alternant de signe lorsque \(t\) tend vers l’infini. Nous dirons dans ces cas que le schéma d’Euler progressif est instable.
• Dans le cas du système, le schéma d’Euler progressif est A-stable ssi \(h<\frac{2}{\lambda_\text{max}}\).

Exercice : étudier l’A-stabilité du schéma d’Euler implicite

\( \begin{cases} u_0 = y_0,\\ u_{n+1} = u_{n}+h \varphi_{n+1}, & n=0,1,\dots,N-1. \end{cases} \)

Le schéma d’Euler rétrograde devient dans le cadre de notre exemple \( (1+\beta h)u_{n+1}=u_n, \qquad n=0,1,2,\dots \) et par suite \( u_{n}=\frac{1}{(1+\beta h)^{n}}u_0, \qquad n=0,1,2,\dots \) Dans ce cas nous voyons que pour tout \(h>0\) nous avons \(\lim_{n\to\infty}u_n=0\), le schéma d’Euler rétrograde est donc toujours stable, sans limitations sur \(h\).

Exercice : étudier l’A-stabilité du schéma d’Euler modifié

\( \begin{cases} u_0 = y_0,\\ u_{n+1} = u_{n}+h \varphi\left(t_n+\frac{h}{2},u_n+\frac{h}{2}\varphi_n\right), & n=0,1,\dots,N-1. \end{cases} \)

Le schéma d’Euler modifié pour notre exemple devient \( u_{n+1}=u_n+h\left( -\beta\left( u_n+\frac{h}{2}(-\beta u_n) \right) \right) =\left(1-\beta h +\frac{(\beta h )^2}{2} \right) u_{n} \) Par induction on obtient \( u_{n}=\left(1-\beta h +\frac{(\beta h )^2}{2} \right)^nu_0. \) Par conséquent, \(\lim\limits_{n\to+\infty}u_n=0\) si et seulement si \( \left| 1 - \beta h +\frac{1}{2}(\beta h)^2 \right|<1. \) Notons \(x\) le produit \(\beta h\) et \(q\) le polynôme \(q(x)=\frac{1}{2}x^2-x+1\) dont le graphe est représenté en figure:

xx = np.linspace(0,3,101)
yy = 1-xx+0.5*xx**2
plt.plot(xx,yy,[0,2,2],[1,1,0],'m:')
plt.axis([0,3,0,3])
plt.xlabel(r'$h\beta$')
plt.ylabel(r'$q(h\beta)$');

Nous avons \(|q(x)|<1\) si et seulement si \(0<x<2\). La relation \(\lim\limits_{n\to+\infty}u_n=0\) est donc satisfaite si et seulement si \( \boxed{ h <\frac{2}{\beta}. } \) Cette condition de stabilité limite le pas \(h\) d’avance en \(t\) lorsqu’on utilise le schéma d’Euler modifié.

Notons qu’on a la même condition de A-stabilité que le schéma d’Euler explicite. Cependant, \(q(x)>0\) pour tout \(x\), donc cette convergence est monotone (tandis qu’avec le schéma d’Euler explicite la convergence n’est monotone que si \(h<\frac{1}{\beta}\)).

Exercice : étudier l’A-stabilité du schéma de Cranck-Nicolson

\( \begin{cases} u_0 = y_0,\\ u_{n+1} = u_{n}+\frac{h}{2}\left(\varphi_n+\varphi_{n+1}\right), & n=0,1,\dots,N-1. \end{cases} \)

Le schéma de Crank-Nicolson appliqué à notre exemple s’écrit \( \left(1+\frac{\beta h}{2}\right)u_{n+1} = \left(1-\frac{\beta h}{2}\right) u_{n} \) et par suite \( u_{n}=\left( \frac{2-\beta h}{2+\beta h} \right)^{n}u_0, \qquad n=0,1,2,\dots \) Par conséquent, \(\lim\limits_{n\to+\infty}u_n=0\) si et seulement si \( \left|\frac{2-\beta h}{2+\beta h}\right|<1. \) Notons \(x\) le produit \(\beta h>0\) et \(q\) la fonction \(q(x)=\frac{2-x}{2+x}=1-2\frac{x}{2+x}\). Nous avons \(0<\frac{x}{2+x}<1\) pour tout \(x\in\mathbb{R}_+^*\), donc \(|q(x)|<1\) pour tout \(x\in\mathbb{R}_+^*\). La relation \(\lim\limits_{n\to+\infty}u_n=0\) est donc satisfaite pour tout \(h>0\): le schéma de Crank-Nicolson est donc toujours stable, sans limitations sur \(h\).
Cependant, si on cherche \(0<q<1\) on a (pour \(x>0\)) \( 1-2\frac{x}{2+x}>0 \iff x<2. \)

xx = np.linspace(0,10,101)
yy = 1-2*xx/(2+xx)
plt.plot(xx,yy)
plt.plot([0,10],[1,1],'r--',[0,10],[-1,-1],'r--')
plt.plot([0,10],[0,0],'m:',[2,2],[-1.5,1.5],'m:')
plt.axis([0,10,-1.5,1.5])
plt.xlabel(r'$h\beta$')
plt.ylabel(r'$q(h\beta)$');

Exercice : étudier l’A-stabilité du schéma de Heun

Le schéma de Heun pour notre exemple devient \( u_{n+1} =u_n+\frac{h}{2}\left( -\beta u_n-\beta\left( u_n+h(-\beta u_n) \right) \right) =\left(1-\beta h +\frac{(\beta h )^2}{2} \right) u_{n} \) ce qui coïncide avec la méthode d’Euler modifiée. Par conséquent, \(\lim\limits_{n\to+\infty}u_n=0\) si et seulement si \( \boxed{ h <\frac{2}{\beta} } \) et cette convergence est monotone.
Cette condition de stabilité limite le pas \(h\) d’avance en \(t\) lorsqu’on utilise le schéma de Heun.

Exercice : étudier l’A-stabilité du schéma de Simpson implicite

Le schéma de Simpson implicite appliqué à notre exemple s’écrit \( \left(1+\beta\frac{h}{6}\right)u_{n+1} = \left(1-\frac{5}{6}\beta h+\frac{1}{3}(\beta h)^2\right) u_{n} \) et par suite \( u_{n}=\left( \frac{ 1-\frac{5}{6}\beta h+\frac{1}{3}(\beta h)^2 }{ 1+\beta\frac{h}{6} } \right)^{n}u_0, \qquad n=0,1,2,\dots \) Par conséquent, \(\lim\limits_{n\to+\infty}u_n=0\) si et seulement si \( \left|\frac{ 6-5\beta h+2(\beta h)^2 }{ 6+\beta h }\right|<1. \) Notons \(x\) le produit \(\beta h>0\) et \(q\) la fonction \(q(x)=\frac{6-5x+2x^2}{6+x}\) dont le graphe est représenté en figure.

xx = np.linspace(0,4,101)
yy = (6-5*xx+2*xx**2)/(6+xx)
plt.plot(xx,yy,[0,3,3],[1,1,0],'m:')
plt.axis([0,4,0,3])
plt.xlabel(r'$h\beta$')
plt.ylabel(r'$q(h\beta)$');

Nous avons \(q(x)>0\) pour tout \(x\in\mathbb{R}_+^*\) et \(q(x)<1\) si et seulement si \(x<3\), donc \(|q(x)|<1\) pour tout \(x\in]0;3[\). Par conséquent, \(\lim\limits_{n\to+\infty}u_n=0\) si et seulement si \( h <\frac{3}{\beta}. \) Cette condition de stabilité limite le pas \(h\) d’avance en \(t\).

Exercice : étudier l’A-stabilité du schéma de Simpson explicite

Le schéma de Simpson explicite pour notre exemple devient \( u_{n+1} = u_n+\frac{h}{6}\left( -\beta u_n -4 \beta\left( u_n+\frac{h}{2}(-\beta u_n) \right) -\beta\left( u_n+h(-\beta u_n) \right) \right) =\left(1-\beta h +\frac{3}{2}(\beta h )^2 \right) u_{n} \) Par induction on obtient \( u_{n}=\left(1-\beta h +\frac{3}{2}(\beta h )^2 \right)^nu_0. \) Par conséquent, \(\lim\limits_{n\to+\infty}u_n=0\) si et seulement si \( \left|1-\beta h +\frac{3}{2}(\beta h )^2\right|<1. \) Notons \(x\) le produit \(\beta h\) et \(q\) le polynôme \(q(x)=\frac{3}{2}x^2-x+1\) dont le graphe est représenté en figure.

xx = np.linspace(0,3,101)
yy = 1-xx+1.5*xx**2
plt.plot(xx,yy,[0,2/3,2/3],[1,1,0],'m:')
plt.axis([0,3,0,3])
plt.xlabel(r'$h\beta$')
plt.ylabel(r'$q(h\beta)$');

Nous avons \(|q(x)|<1\) si et seulement si \(0<x<\frac{2}{3}\). La relation \(\lim\limits_{n\to+\infty}u_n=0\) est donc satisfaite si et seulement si \( h <\frac{2}{3\beta}. \) Cette condition de stabilité limite le pas \(h\) d’avance en \(t\).

Schémas à deux pas

AvertissementRemarques sur les suites à 2 pas

Considérons une relation de récurrence de la forme

\( u_{n+1} = \gamma u_n + \eta u_{n-1}. \)

Elle se réécrit sous la forme matricielle \( \begin{pmatrix} u_{n+1}\\ u_n \end{pmatrix}= \begin{pmatrix}\gamma & \eta\\ 1 & 0 \end{pmatrix}\begin{pmatrix} u_n\\ u_{n-1} \end{pmatrix} \) ainsi \( \begin{pmatrix} u_{n+1}\\ u_n \end{pmatrix}= \begin{pmatrix}\gamma & \eta\\ 1 & 0 \end{pmatrix}^n \begin{pmatrix} u_1\\ u_0 \end{pmatrix} \)

La suite \(u_n\) tend vers 0 si et seulement si les valeurs propres de la matrice de transition sont de module strictement inférieur à 1.

Les valeurs propres de cette matrice sont les racines du polynôme caractéristique

\( r^2-\gamma r-\eta=0 \)

et notons \(\Delta = \gamma^2 + 4\eta\) son discriminant.

3 cas se présentent alors :

1. Le polynôme a deux racines réelles distinctes \(r_{1,2} = \frac{\gamma\pm\sqrt{\Delta}}{2}\). La matrice de transition est diagonalisable ainsi \( \begin{pmatrix} \gamma & \eta\\ 1 & 0 \end{pmatrix}^n = \mathbb{P} \begin{pmatrix} r_1^n & 0\\ 0 & r_2^n \end{pmatrix} \mathbb{P}^{-1}. \) La suite \(u_n\) s’écrit \(u_n = A r_1^n + B r_2^n\), où les constantes \(A\) et \(B\) sont déterminées par les conditions initiales \(u_0\) et \(u_1\). La suite tend vers 0 si et seulement si \(|r_1|<1\) et \(|r_2|<1\).

2. Le polynôme a une racine double \(r = \frac{\gamma}{2}\), ce qui correspond au cas où le discriminant \(\Delta\) est nul. Alors \( \begin{pmatrix} \gamma & \eta\\ 1 & 0 \end{pmatrix}^n = \mathbb{P} \begin{pmatrix} r^n & nr^{n-1}\\ 0 & r^n \end{pmatrix} \mathbb{P}^{-1}. \) La suite s’écrit alors \(u_n = (A n + B) r^n\), où les constantes \(A\) et \(B\) sont déterminées par les conditions initiales \(u_0\) et \(u_1\). La suite tend vers 0 si et seulement si \(|r|<1\).

3. Le polynôme a deux racines complexes conjugués \(r_{1,2} = \frac{\gamma}{2}\pm i \frac{\sqrt{-\Delta}}{2} = \varrho e^{i\theta}\), ce qui correspond au cas où le discriminant \(\Delta\) est négatif. Alors \( \begin{pmatrix} \gamma & \eta\\ 1 & 0 \end{pmatrix}^n = \mathbb{P} \begin{pmatrix} \varrho^n e^{i\theta n} & 0\\ 0 & \varrho^n e^{-i\theta n} \end{pmatrix} \mathbb{P}^{-1} = \mathbb{Q} \varrho^n \begin{pmatrix} \cos(n\theta) & -\sin(n\theta)\\ \sin(n\theta) & \cos(n\theta) \end{pmatrix} \mathbb{Q}^{-1} \) La suite s’écrit alors \(u_n = \varrho^n\big(C\cos(n\theta)+D\sin(n\theta)\big)\) avec \(\varrho = \sqrt{-\eta}\), et \(\cos(\theta) = \frac{\gamma}{2\varrho}\), où les constantes \(C\) et \(D\) sont déterminées par les conditions initiales \(u_0\) et \(u_1\). La suite tend vers 0 si et seulement si \(\varrho<1\).

Exemple : étude de l’A-stabilité du schéma d’Adam-Bashforth à 2 pas.

Le schéma d’Adam-Bashforth à 2 pas s’écrit

\( \begin{cases} u_0 = y_0,\\ u_1 \approx y_1,\\ u_{n+1} = u_{n}+\frac{h}{2}\left(3\varphi_n-\varphi_{n-1}\right), & n=1,2,\dots,N-1. \end{cases} \)

Appliqué au problème modèle \(y'(t) = -\beta y(t)\) s’écrit

\( u_{n+1} = u_n + \frac{h}{2}\left( -3\beta u_n + \beta u_{n-1} \right) = \left(1 - \frac{3\beta h}{2}\right) u_n + \frac{\beta h}{2} u_{n-1}. \)

C’est bien une suite à 2 pas de la forme

\( u_{n+1} = \gamma u_n + \eta u_{n-1} \)

avec \(\gamma = 1 - \frac{3\beta h}{2}\) et \(\eta = \frac{\beta h}{2}\).

Dans notre cas, \(\Delta = \frac{9}{4}(\beta h)^2-(\beta h)+1\) qui est toujours strictement positif ainsi on est dans le premier cas :

\( u_n = A r_1^n + B r_2^n \)

avec

\( r_{1,2} = \frac{1 - \frac{3\beta h}{2} \pm \sqrt{\left(1 - \frac{3\beta h}{2}\right)^2 - 2\beta h}}{2} \)

Pour que la méthode soit A-stable, il faut que les racines satisfassent \(|r_i| < 1\) pour \(i=1,2\). Si on note \(x=\beta h\) et on trace le graphe des fonctions \(r_1(x)\) et \(r_2(x)\) on trouve que la condition \(|r_i|<1\) est satisfaite si et seulement si \(x<1\), c’est-à-dire \(h < \frac{1}{\beta}\).

gamma = lambda x : 1 - 3/2*x
eta   = lambda x : 1/2*x
Delta = lambda x : gamma(x)**2 + 4*eta(x)
r1 = lambda x: (gamma(x) - np.sqrt(Delta(x)))/2 
r2 = lambda x: (gamma(x) + np.sqrt(Delta(x)))/2

xx = np.linspace(0,2,101)
yy1 = r1(xx)
yy2 = r2(xx)
plt.plot(xx,yy1, 'g', label=r'$r_1(h\beta)$')
plt.plot(xx,yy2, 'r', label=r'$r_2(h\beta)$')
plt.hlines(1,0,2, label=r'$1$', linestyles='dashed')
plt.hlines(-1,0,2, label=r'$-1$', linestyles='dashed')
plt.legend()
plt.xlabel(r'$h\beta$');

Exercice : étudier l’A-stabilité du schéma d’Adam-Multon à 2 pas.

Le schéma d’Adam-Bashforth à 2 pas s’écrit

\( \begin{cases} u_0 = y_0,\\ u_1 \approx y_1,\\ u_{n+1} = u_{n}+\frac{h}{12}\left(8\varphi_n-\varphi_{n-1}+5\varphi_{n+1}\right), & n=1,2,\dots,N-1. \end{cases} \)

Appliqué au problème modèle \(y'(t) = -\beta y(t)\) s’écrit

\( u_{n+1} = u_n + \frac{h}{12}\left( 8(-\beta u_n) - (-\beta u_{n-1}) + 5(-\beta u_{n+1}) \right). \)

C’est bien une suite à 2 pas de la forme

\( u_{n+1} = \gamma u_n + \eta u_{n-1} \)

avec \(\gamma = \frac{12-8\beta h}{12+5\beta h}\) et \(\eta = \frac{\beta h}{12+5\beta h}\).

Dans notre cas, \(\Delta = \left(\frac{12-8\beta h}{12+5\beta h}\right)^2 + 4 \frac{\beta h}{12+5\beta h}\) qui est toujours strictement positif ainsi on est dans le premier cas :

\( u_n = A r_1^n + B r_2^n \)

avec

\( r_{1,2} = \frac{\frac{12-8\beta h}{12+5\beta h} \pm \sqrt{\left(\frac{12-8\beta h}{12+5\beta h}\right)^2 + 4 \frac{\beta h}{12+5\beta h}}}{2} \)

Pour que la méthode soit A-stable, il faut que les racines satisfassent \(|r_i| < 1\) pour \(i=1,2\). Si on note \(x=\beta h\) et on trace le graphe des fonctions \(r_1(x)\) et \(r_2(x)\) on trouve que la condition \(|r_i|<1\) est satisfaite si et seulement si \(x<6\), c’est-à-dire \(h < \frac{6}{\beta}\).

gamma = lambda x : (12-8*x)/(12+5*x)
eta   = lambda x : x/(12+5*x)
Delta = lambda x : gamma(x)**2 + 4*eta(x)
r1 = lambda x: (gamma(x) - np.sqrt(Delta(x)))/2 
r2 = lambda x: (gamma(x) + np.sqrt(Delta(x)))/2

xx = np.linspace(0,7,101)
yy1 = r1(xx)
yy2 = r2(xx)
plt.plot(xx,yy1, 'g', label=r'$r_1(h\beta)$')
plt.plot(xx,yy2, 'r', label=r'$r_2(h\beta)$')
plt.hlines(1,0,7, label=r'$1$', linestyles='dashed')
plt.hlines(-1,0,7, label=r'$-1$', linestyles='dashed')
plt.legend()
plt.xlabel(r'$h\beta$');

Exercice : étudier l’A-stabilité du schéma BDF à 2 pas.

Le schéma BDF à 2 pas s’écrit

\( \begin{cases} u_0 = y_0,\\ u_1 \approx y_1,\\ u_{n+1} = \frac{4}{3} u_{n} -\frac{1}{3}u_{n-1} + \frac{2}{3}h \varphi_{n+1}, & n=1,2,\dots,N-1. \end{cases} \)

Appliqué au problème modèle \(y'(t) = -\beta y(t)\) s’écrit

\( u_{n+1} = \frac{4}{3} u_{n} -\frac{1}{3}u_{n-1} + \frac{2}{3}h (-\beta u_{n+1}). \)

C’est bien une suite à 2 pas de la forme

\( u_{n+1} = \gamma u_n + \eta u_{n-1} \)

avec \(\gamma = \frac{4}{3+2\beta h}\) et \(\eta = -\frac{1}{3+2\beta h}\).

Dans notre cas, \(\Delta = \gamma^2 + 4\eta = \frac{4-2\beta h}{(3+2\beta h)^2}\) qui change de signe selon les valeurs de \(\beta h\).

gamma = lambda x :  4/(3+2*x)
eta   = lambda x : -1/(3+2*x)
Delta = lambda x : gamma(x)**2 + 4*eta(x)

xx = np.linspace(0,7,101)
dd = Delta(xx)
plt.plot(xx,dd, 'b', label=r'$\Delta(h\beta)$')
plt.hlines(0,0,7, label=r'$0$', linestyles='dashed')
plt.vlines(2/13,-0.5,0.5, label=r'$h\beta=1$', linestyles='dashed')
plt.legend()

• Si \(\beta h < \frac{1}{2}\), alors \(\Delta > 0\) et on est dans le premier cas avec deux racines réelles distinctes et \(|r_1|<1\) et \(|r_2|<1\).
• Si \(\beta h = \frac{1}{2}\), alors \(\Delta = 0\) et on est dans le deuxième cas avec une racine double et \(|r|<1\).
• Si \(\beta h > \frac{1}{2}\), alors \(\Delta < 0\) et on est dans le troisième cas avec deux racines complexes conjugées et \(\rho<1\).

# CAS DELTA > 0 : x<1/2

r1 = lambda x: (gamma(x) - np.sqrt(Delta(x)))/2 
r2 = lambda x: (gamma(x) + np.sqrt(Delta(x)))/2

xx = np.linspace(0,1/2,101)
yy1 = r1(xx)
yy2 = r2(xx)
plt.plot(xx,yy1, 'g', label=r'$r_1(h\beta)$')
plt.plot(xx,yy2, 'r', label=r'$r_2(h\beta)$')
plt.hlines(1,0,1/2, label=r'$1$', linestyles='dashed')
plt.hlines(-1,0,1/2, label=r'$-1$', linestyles='dashed')
plt.legend()
plt.xlabel(r'$h\beta$');

# CAS DELTA > 0 : x=1/2

r = gamma(1/2)/2 

print(f"Pour hβ=1/2, on a r1=r2={r}")

Pour hβ=1/2, on a r1=r2=0.5

# CAS DELTA < 0 : x>1/2
rho = lambda x: np.sqrt(-Delta(x))/2

xx = np.linspace(1/2+1.e-5,7,101)
yy1 = rho(xx)
plt.plot(xx,yy1, 'm', label=r'$\rho(h\beta)$')
plt.hlines(1,1/2,7, label=r'$1$', linestyles='dashed')
plt.hlines(-1,1/2,7, label=r'$-1$', linestyles='dashed')
plt.legend()
plt.xlabel(r'$h\beta$');    

Vérification empirique des conditions de A-stabilité

Exercice : tester numériquement l’A-stabilité des schémas classiques à un pas

On considère le toy-model avec \(\beta=1\) et on teste numériquement l’A-stabilité des schémas d’Euler explicite, d’Euler implicite, d’Euler modifié, de Crank-Nicolson et de Heun.

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

Calculer la solution exacte. Calculer ensuite la solution approchée par les schémas pour différentes valeurs de \(h\). Vérifier que les conditions de stabilité sont respectées.

Solution exacte

L’unique solution du problème de Cauchy est la fonction \(y(t)=e^{-t}\) définie pour tout \(t\) \(\in\) \(\mathbb{R}\).

phi = lambda t,y : -y
sol_exacte = lambda t : np.exp(-t)

# INITIALISATION
t0, y0, tfinal = 0 , 1 , 12

# PREPARATIONS DES CAS TEST pour différentes valeurs de h
H = [ 2**(k-2) for k in range(5) ]
T = [ np.arange(t0,tfinal,h) for h in H] # liste des grilles

yy = sol_exacte(T[0]) # on affichera la sol exacte sur la grille la plus fine

Solution approchée par la méthode d’Euler explicite

La méthode d’Euler explicite pour cette EDO s’écrit

\( u_{n+1}=(1-h)u_n. \)

En procédant par récurrence sur \(n\), on obtient

\( u_{n+1}=(1-h)^{n+1}. \)

De la formule \(u_{n+1}=(1-h)^{n+1}\) on déduit que

• si \(0<h<1\) alors la solution numérique est stable et convergente,
• si \(h=1\) alors la solution numérique est stationnaire \(u_n=0\) pour tout \(n\in\mathbb{N}^*\),
• si \(1<h<2\) alors la solution numérique oscille mais est encore convergente,
• si \(h=2\) alors la solution numérique oscille, plus précisément on a \(u_{2n}=1\) et \(u_{2n+1}=-1\) pour tout \(n\in\mathbb{N}^*\),
• si \(h>2\) alors la solution numérique oscille et diverge.

En effet, la méthode est A-stable si et seulement si \(|1 - h| < 1\).

Remarque : la suite obtenue est une suite géométrique de raison \(q=1-h\). On sait qu’une telle suite

• diverge si \(|q|>1\) ou \(q=-1\),
• est stationnaire si \(q=1\),
• converge vers \(0\) si \(|q|<1\).
• converge de façon monotone vers \(0\) si \(0<q<1\).

Voyons graphiquement ce que cela donne avec différentes valeurs de \(h>0\) :

• si \(h=4\) alors \(t_n=4n\) et \(u_{n}=\left(-4\right)^{n}\) tandis que \(y(t_{n})=e^{-4n}\),
• si \(h=2\) alors \(t_n=2n\) et \(u_{n}=\left(-1\right)^{n}\) tandis que \(y(t_{n})=e^{-2n}\),
• si \(h=1\) alors \(t_n=n\) et \(u_{n}=0\) tandis que \(y(t_{n})=e^{-n}\),
• si \(h=\frac{1}{2}\) alors \(t_n=\frac{n}{2}\) et \(u_{n}=\left(\dfrac{1}{2}\right)^{n}\) tandis que \(y(t_{n})=e^{-n/2}\),
• si \(h=\frac{1}{4}\) alors \(t_n=\frac{n}{4}\) et \(u_{n}=\left(\dfrac{3}{4}\right)^{n}\) tandis que \(y(t_{n})=e^{-n/4}\).

Ci-dessous sont tracées sur l’intervalle \([0;10]\), les courbes représentatives de la solution exacte et de la solution calculée par la méthode d’Euler explicite. En faisant varier le pas \(h\) nous pouvons constater que si \(h>1\) l’erreur commise entre la solution exacte et la solution calculée est amplifiée d’un pas à l’autre.

Remarque : la première itérée a la même pente quel que soit le pas \(h\) (se rappeler la construction géométrique de la méthode d’Euler).

# SCHEMA
def EE(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)  
    uu[0] = y0
    for i in range(len(tt)-1):
        uu[i+1] = uu[i]+h*phi(tt[i],uu[i])
    return uu


# CALCUL DE LA SOLUTION SUR CHAQUE GRILLE
U = [ EE(phi,tt,y0) for tt in T ] 
     
# AFFICHAGE
plt.figure(1, figsize=(15, 10))
plt.axis([t0, tfinal, -1, 1])

plt.plot(T[0],yy,'-',label='$y(t)=e^{-t}$')
for k in range(5):
    plt.plot(T[k],U[k],'-o',label=f'h = {H[k]}')
plt.legend()
plt.grid();

Solution approchée par la méthode d’Euler implicite

La méthode d’Euler implicite pour cette EDO s’écrit

\( u_{n+1}=\frac{1}{1+h}u_n. \)

En procédant par récurrence sur \(n\), on obtient

\( u_{n+1}=\frac{1}{(1+h)^{n+1}}. \)

De la formule \(u_{n+1}=(1+h)^{-(n+1)}\) on déduit que la solution numérique est stable et convergente pour tout \(h\). En effet, la méthode est inconditionnellement A-stable.

Remarque : la suite obtenue est une suite géométrique de raison \(q=\frac{1}{1+h}\in]0;1[\).

# SCHEMA
def EI(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for i in range(len(tt)-1):
        uu[i+1] = fsolve(lambda x:-x+uu[i]+h*phi(tt[i+1],x),uu[i])[0]   
    return uu


# CALCUL DE LA SOLUTION SUR CHAQUE GRILLE
U = [ EI(phi,tt,y0) for tt in T ] 
     
# AFFICHAGE
plt.figure(1, figsize=(15, 10))
plt.axis([t0, tfinal, -1, 1])

plt.plot(T[0],yy,'-',label='$y(t)=e^{-t}$')
for k in range(5):
    plt.plot(T[k],U[k],'-o',label=f'h = {H[k]}')
plt.legend()
plt.grid();

Solution approchée par la méthode d’Euler modifié

La méthode d’Euler modifié pour cette EDO s’écrit

\( u_{n+1}=(1-h+\frac{1}{2}h^2)u_n. \)

En procédant par récurrence sur \(n\), on obtient

\( u_{n+1}=(1-h+\frac{1}{2}h^2)^{n+1}. \)

La suite obtenue est une suite géométrique de raison \(q=1-h+\frac{1}{2}h^2\):

• si \(0<h<2\) alors la solution numérique est stable et convergente,
• si \(h=2\) alors la solution numérique est stationnaire \(u_n=u_0\) pour tout \(n\in\mathbb{N}^*\),
• si \(h>2\) alors la solution numérique diverge vers \(+\infty\).

En effet, la méthode est A-stable si et seulement si \(|1-h+\frac{1}{2}h^2| < 1\). De plus, comme \(1-h+\frac{1}{2}h^2>0\) pour tout \(h>0\), la suite est positive.

# SCHEMA
def EM(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for i in range(len(tt)-1):
        utilde = uu[i]+h/2*phi(tt[i],uu[i])
        uu[i+1] = uu[i]+h*phi(tt[i]+h/2,utilde)
    return uu


# CALCUL DE LA SOLUTION SUR CHAQUE GRILLE
U = [ EM(phi,tt,y0) for tt in T ] 
     
# AFFICHAGE
plt.figure(1, figsize=(15, 10))
plt.axis([t0, tfinal, -1, 2])

plt.plot(T[0],yy,'-',label='$y(t)=e^{-t}$')
for k in range(5):
    plt.plot(T[k],U[k],'-o',label=f'h = {H[k]}')
plt.legend()
plt.grid();

Solution approchée par la méthode de Crank-Nicolson

La méthode de Crank-Nicolson pour cette EDO s’écrit

\( u_{n+1}=\frac{2-h}{2+h}u_n. \)

En procédant par récurrence sur \(n\), on obtient

\( u_{n+1}=\left(\frac{2-h}{2+h}\right)^{n+1}. \)

La suite obtenue est une suite géométrique de raison \(q=\frac{2-h}{2+h}\in]-1;1[\). On déduit que la solution numérique est stable et convergente pour tout \(h\). En effet, la méthode est inconditionnellement A-stable.

De plus, la suite est positive ssi \(0<q<1\) donc ssi \(h<2\).

# SCHEMA
def CN(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for i in range(len(tt)-1):
        uu[i+1] = fsolve(lambda x:-x+uu[i]+h/2*phi(tt[i],uu[i])+h/2*phi(tt[i+1],x),uu[i])[0]
    return uu   


# CALCUL DE LA SOLUTION SUR CHAQUE GRILLE
U = [ CN(phi,tt,y0) for tt in T ] 
     
# AFFICHAGE
plt.figure(1, figsize=(15, 10))
plt.axis([t0, tfinal, -1, 1])

plt.plot(T[0],yy,'-',label='$y(t)=e^{-t}$')
for k in range(5):
    plt.plot(T[k],U[k],'-o',label=f'h = {H[k]}')
plt.legend()
plt.grid();

Solution approchée par la méthode de Heun

La méthode de Heun pour cette EDO s’écrit

\( u_{n+1} = \left( 1-h+\frac{1}{2}h^2 \right) u_n. \)

On retrouve le schéma d’Euler modifié.

Annexe : A-stabilité d’un système

Considérons un système de \(n\) EDO de la forme

\( \mathbf{y}'(t) = -\mathbb{A}\mathbf{y}(t) \)

où \(\mathbb{A}\) est une matrice constante ayant \(n\) valeurs propres positives \(\lambda_i\), \(i = 1, 2, \dots, n\).
La solution générale s’écrit

\( \mathbf{y}(t) = \sum_i c_i \mathbf{v}_i e^{-\lambda_i t} \)

où \(\mathbf{v}_i\) est le \(i\)ᵉ vecteur propre associé à la valeur propre \(\lambda_i\). On a donc :

\( \lim_{t \to +\infty} y_i(t) = 0 \quad \forall i = 1, 2, \dots, n. \)

Systèmes raides (stiff systems) : un système d’équations différentielles est dit raide lorsque les valeurs propres \(\lambda_i\) de la matrice \(\mathbb{A}\) présentent des magnitudes très différentes. En d’autres termes, si certaines valeurs propres sont beaucoup plus grandes que d’autres, alors le système évolue selon des échelles de temps très variées. Cette différence de magnitudes crée un défi pour les méthodes numériques : pour garantir la A-stabilité, il est parfois nécessaire de choisir un pas de temps \(h\) suffisamment petit, proportionnel au plus grand \(\lambda_i\). Cela peut rendre les calculs très coûteux lorsque les valeurs propres sont très disparates. C’est pour cela qu’on préfère des méthodes inconditionnellement A-stables pour approcher la solution d’un système raide, car elles permettent d’utiliser des pas de temps plus grands tout en maintenant la A-stabilité, indépendamment du plus grand \(\lambda_i\).