Propriétés d’un schéma sur un intervalle non borné

Auteur·rice

Gloria Faccanoni

Date de publication

29 mars 2026

Revenons sur la question fondamentale de la qualité d’un schéma numérique, déjà abordée dans le chapitre précédent : quelles propriétés doit satisfaire un schéma numérique pour que les solutions approchées soient de « bonne qualité » ?

Sur un intervalle de temps borné, nous avons établi que la consistance et la zéro-stabilité sont des conditions nécessaires et suffisantes pour assurer la convergence de la solution approchée vers la solution exacte du problème de Cauchy. Autrement dit, lorsque le pas de temps \(h\) tend vers zéro, la solution numérique converge vers la solution exacte.

Toutefois, ces deux propriétés s’avèrent insuffisantes dès lors que l’on souhaite résoudre le problème de Cauchy sur des intervalles de temps très longs, voire infinis. En effet, que se passe-t-il lorsque le pas \(h > 0\) est fixé et que l’horizon \(T\) tend vers l’infini ?

Dans les deux situations (\(h \to 0\) avec \(T\) fixé, ou \(T \to +\infty\) avec \(h\) fixé), le nombre de nœuds tend vers l’infini. Cependant, les enjeux sont fondamentalement différents :

dans le premier cas, on cherche à contrôler l’erreur locale en chaque point ;
dans le second, c’est le comportement asymptotique de la solution approchée qui est en jeu.

Lorsque \(T \to +\infty\) avec \(h\) fixé, le nombre \(N\) de sous-intervalles tend vers l’infini indépendamment de \(h\), ce qui rend la notion de zéro-stabilité inapplicable. On cherche alors des méthodes capables d’approcher fidèlement la solution sur des intervalles de temps arbitrairement grands, idéalement sans imposer de restriction trop sévère sur le pas \(h\).

Code

import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
from scipy.optimize import fsolve
from IPython.display import Markdown, Math
plt.rcParams.update({'font.size': 12})

Définition de A-stabilité

Pour analyser la stabilité d’un schéma numérique, on étudie son comportement sur un problème test. Dans les études théoriques, on considère généralement le cas où \(\lambda\) est un nombre complexe avec une partie réelle négative (\(\Re(\lambda) < 0\)). Cependant, dans les exercices, nous allons nous concentrer sur le cas où \(\lambda\) est un nombre réel négatif et, pour éviter toute ambiguïté, on posera \(\beta = -\lambda > 0\).

Problème test de Dahlquist

\( \begin{cases} y'(t) = \lambda y(t), \quad t > 0 \\ y(0) = 1 \end{cases} \) avec \(\Re(\lambda) < 0\).

La solution exacte est \(y(t) = e^{\lambda t}\) ainsi \(|y(t)| \to 0\) lorsque \(t \to +\infty\).

Remarque : dans le cas réel, on peut réécrire le problème de Dahlquist sous la forme \(y' = -\beta y\) avec \(\beta > 0\). La solution exacte est alors \(y(t) = e^{-\beta t}\) et \(y(t) \to 0\) lorsque \(t \to +\infty\).

L’enjeu de la A-stabilité est de déterminer si le schéma numérique reproduit ce comportement de décroissance pour un pas de temps \(h\) donné.

Exemple : comportement du schéma d’Euler explicite sur le problème de Dahlquist

Le schéma d’Euler progressif \(u_{n+1} = u_n + h\varphi(t_n, u_n)\) 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 \)

soit encore

\( u_{n+1} = q(h\beta) u_n \qquad\text{avec } q(h\beta) = 1 - \beta h. \)

et par suite

\( u_{n}=[q(\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|q(\beta h)\right|<1. \)

Si on note \(x=\beta h>0\) et on trace la fonction \(x\mapsto q(x)\), on cherche les valeurs de \(x\) pour lesquelles le graphe de la fonction est compris entre les droites \(y=-1\) et \(y=1\).

Le schéma d’Euler progressif a le même comportement que la solution exacte si et seulement si

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

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

Code

x = np.linspace(0, 3, 101)
q = 1 - x
plt.figure(figsize=(8, 4))
plt.grid()
plt.plot(x, q, label='$q(x)$')
plt.ylim(-2, 2) # On zoome sur la zone intéressante [-1, 1]
plt.xlim(x.min(), x.max())
plt.axhline(y=1, color='black', linestyle='-', alpha=0.5)
plt.axhline(y=-1, color='black', linestyle='-', alpha=0.5)
plt.fill_between(x, -1, 1, color='gray', alpha=0.1, label='Zone de A-Stabilité')
plt.axhline(0, color='magenta', linestyle=':', label='$y=0$')
plt.scatter(1, 0, color='black', label='$q(1) = 0$')
plt.xlabel('$x$')
plt.ylabel('$q(x)$')
plt.title('EE : $q(x) = 1 - x$')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.show()

Remarques :

Si \(\left|1-\beta h\right|\ge 1\), alors \(u_n\) ne tend pas vers 0 lorsque \(t\) tend vers l’infini. Nous dirons dans ces cas que le schéma d’Euler progressif est instable. Plus précisément :
- si \(1-\beta h>1\) alors \(u_n\) tend vers \(+\infty\) lorsque \(t\) tend vers l’infini ;
- si \(1-\beta h<-1\) alors \(u_n\) tend vers l’infini en alternant de signe lorsque \(t\) tend vers l’infini ;
- si \(1-\beta h=1\) alors \(u_n\) est constant ;
- si \(1-\beta h=-1\) alors \(u_n\) alterne de signe.
Pour \(h<\frac{2}{\beta}\), la suite \(u_n\) tend vers 0 lorsque \(t\) tend vers l’infini. 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é).

Code

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

BETA = 2
phi = lambda t,y: -BETA*y
y0 = 1
fig, axx = plt.subplots(nrows=1, ncols=5, figsize=(20, 3))
for ax, x in zip(axx, (2.5,2,1.5,1,0.5)):
    tt = np.arange(0, 10, x/BETA)
    uu = EE(tt, y0, phi)
    ax.grid()
    ax.plot(tt, uu, '.-', label='EE')
    ax.plot(tt, np.exp(-BETA*tt), label='Solution exacte')
    ax.set_xlabel('$t$')
    ax.set_ylabel('$u(t)$')
    ax.set_title(f'EE avec $h\\beta={(tt[1]-tt[0])*BETA:g}$ ')
    ax.legend()
plt.tight_layout()

Définition : A-stabilité (stabilité absolue) d’un schéma numérique

Soit \(h>0\) un pas de temps fixé et \(t_n = nh\). On note \(u_n \approx y(t_n)\) la solution numérique et \(y_n\) la solution exacte du problème de Dahlquist à l’instant \(t_n\).

On introduit la variable de stabilité \(z = h\lambda \in \mathbb{C}^\).

Le schéma est dit

inconditionnellement A-stable lorsque \(u_n\to0\) pour tout \(z\in\mathbb{C}^-\), c’est-à-dire pour tout \(\lambda\) avec \(\Re(\lambda) < 0\) et pour tout pas de temps \(h > 0\) (aucune restriction sur \(h\)) ;
conditionnellement A-stable lorsque \(u_n\to0\) seulement si \(0 < h < h_{\max}(\lambda)\).

NB : la stabilité numérique dépend uniquement du produit \(z\), et non de \(h\) ou \(\lambda\) séparément.

Remarque : dans le cas réel, la variable de stabilité s’écrit \(x = h\beta = -h\lambda=-z\) et le schéma est A-stable si \(u_n \to 0\) pour tout \(x \in \mathbb{R}^+\), c’est-à-dire pour tout \(\beta > 0\) et pour tout pas de temps \(h > 0\) (aucune restriction sur \(h\)). Conditionnellement A-stable si \(u_n \to 0\) seulement si \(0 < h < h_{\max}(\beta)\). Dans l’exemple du schéma d’Euler explicite, on trouve que \(h_{\max}(\beta) = \frac{2}{\beta}\).

Définition : Région et Intervalle de stabilité absolue

La région de stabilité absolue \(\mathcal{A}\) est l’ensemble des points du plan complexe où le schéma est stable : \(\mathcal{A} = \{ z \in \mathbb{C}^- \mid |u_n| \to 0 \text{ quand } n \to +\infty \}\)

Si \(\mathcal{A}=\mathbb{C}^-\), le schéma est inconditionnellement A-stable.
Sinon, il est A-stable sous condition, ce qui impose une contrainte sur le pas : \(0 < h < h_{\max}(\lambda)\).

L’intervalle de stabilité est la trace de cette région sur l’axe réel : \(\mathcal{I} = \mathcal{A} \cap \mathbb{R}^- = ]-h_{\max}\beta, 0].\)

L’étude de la A-stabilité revient à analyser une suite de récurrence linéaire dont les coefficients dépendent de \(x\) (ou \(z\)).

Méthodes à un pas : la relation est de type \(u_{n+1} = q(x)u_n\). La condition de stabilité est directe : il suffit que le facteur d’amplification vérifie \(|q(x)| < 1\).
Méthodes multi-pas (MPL) : pour un schéma à \(q\) étapes, la récurrence lie \(u_{n+1}\) à plusieurs états précédents (\(u_n, u_{n-1}, \dots, u_{n-q+1}\)). On obtient une récurrence linéaire à coefficients constants dépendant de \(x\) (ou \(z\)) dont la solution générale s’écrit : \(u_n = \sum_{j=0}^{p} C_j r_j^n(x)\) avec \(C_j\) des constantes qui dépendent des conditions initiales et les racines \(r_j(x)\) qui dépendent de \(x\). L’analyse devient plus délicate : la limite de \(u_n\) dépend de la position des racines du polynôme de stabilité par rapport au cercle unité. Si toutes les racines sont strictement à l’intérieur du disque unité (\(|r_j(x)| < 1\)), alors la méthode est absolument stable pour ce paramètre \(x\).

Remarques : à 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.

Nous allons 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 \textbf{Méthodes à 1 pas} & & & & & \\ \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} \\ \hline \textbf{Méthodes à 2 pas} & & & & & \\ \hline \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} \)

Lien entre A-stabilité et Problèmes Stiff

La A-stabilité, bien qu’elle soit une propriété définie par le comportement asymptotique (\(t \to \infty\)), est l’outil fondamental pour approcher la solution des problèmes raides sur des intervalles de temps borné. Nous pouvons voir cette équivalence comme un changement d’échelle temporelle.

La raideur comme dilatation du temps : considérons le problème de Cauchy “stiff” sur l’intervalle borné \(I = [0, 1]\) avec \(\gamma \gg 1\) :

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

Pour analyser le comportement d’un schéma numérique à l’échelle de la variation réelle de la solution, effectuons le changement de variable (dilatation du temps) suivant : \(\tau = \gamma t\) et posons \(\tilde y(\tau) = y(t=\frac{\tau}{\gamma})\).

Par dérivation composée, on obtient \(\tilde y'(\tau) = y'(t) \cdot t'(\tau) = - \gamma y(\frac{\tau}{\gamma}) \cdot \frac{1}{\gamma} = -\tilde y(\tau)\). Le problème se reformule alors en termes de \(\tilde y\) et \(\tau\) :

\( \begin{cases} \tilde y'(\tau) = -\tilde y(\tau), & \tau \in [0, \gamma], \\ \tilde y(0) = 1. \end{cases} \)

Du problème local au comportement asymptotique : ce changement de variable révèle une équivalence qui lie la raideur à la stabilité.

Le problème initial est un problème stiff sur un petit intervalle \([0, 1]\). La “pente” initiale est énorme (\(-\gamma\)) ;
le problème équivalent est un problème non-stiff (pente \(\beta=-1\)) mais sur un intervalle “presque infini” \([0, \gamma]\). Par exemple, si \(\gamma = 10^6\), intégrer jusqu’à \(t=1\) revient à intégrer jusqu’à \(\tau = 10^6\).

Pour intégrer numériquement sur un intervalle aussi grand sans que les erreurs de discrétisation ne s’accumulent ou n’explosent, le schéma doit simuler correctement le comportement de la solution pour \(\tau \to \infty\). C’est précisément la définition de la A-stabilité.

En résumé : la raideur d’une EDO contraint le schéma numérique à “résoudre” des dynamiques internes extrêmement rapides, même si l’intervalle de temps global est court. Sans A-stabilité inconditionnelle, le pas de temps \(h\) est dicté par la stabilité (il doit être inversement proportionnel à \(\gamma\)), ce qui rend le calcul prohibitivement long. Avec la A-stabilité, on s’affranchit de cette contrainte. Le pas de temps peut être choisi librement en fonction de la précision souhaitée, permettant ainsi de traverser efficacement les zones de forte raideur.

Étude de la A-stabilité pour un schéma à un pas

Un schéma à un pas est un schéma numérique dans lequel la valeur de \(u_{n+1}\) dépend uniquement de \(u_n\) (et éventuellement de \(t_n\)). Ces méthodes, appliquées au toy-model, conduisent à des relations de la forme \(u_{n+1} = q(\beta h) u_n\), où \(q\) est une fonction de la variable \(x = \beta h\). La condition de A-stabilité se traduit alors par l’exigence que \(|q(x)| < 1\) pour tout \(x > 0\).

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

Exercice : étudier l’A-stabilité du schéma d’Euler implicite \(u_{n+1} = u_{n}+h \varphi_{n+1}\).

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, \) soit encore \( u_n = [q(\beta h)]^n u_0 \qquad\text{avec } q(\beta h) = \frac{1}{1+\beta h}. \)

Dans ce cas \(q(\beta h)\in]0,1[\) pour tout \(h>0\), ainsi \(\lim_{n\to\infty}u_n=0\) pour tout \(h>0\) : le schéma d’Euler rétrograde est donc toujours stable, sans limitations sur \(h\).

Notons \(x=\beta h>0\) et traçons la fonction \(x\mapsto q(x)\). On trouve que \(q(x)\) est strictement décroissante, \(q(0)=1\) et tend vers 0 lorsque \(x\) tend vers l’infini.

Code

x = np.linspace(0, 3, 101)
q = 1/(1 + x)
plt.figure(figsize=(8, 5))
plt.grid()
plt.plot(x, q, label='$q(x)$')
plt.ylim(-2, 2) # On zoome sur la zone intéressante [-1, 1]
plt.xlim(x.min(), x.max())
plt.axhline(y=1, color='black', linestyle='-', alpha=0.5)
plt.axhline(y=-1, color='black', linestyle='-', alpha=0.5)
plt.fill_between(x, -1, 1, color='gray', alpha=0.1, label='Zone de A-Stabilité')
plt.axhline(0, color='magenta', linestyle=':', label='$y=0$')
# plt.scatter(1, 0, color='black', label='$q(1) = 0$')
plt.xlabel(r'$x=\beta h$')
plt.ylabel('$q(x)$')
plt.title('EI : $q(x) = \\frac{1}{1 + x}$')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.show()

Exercice : étudier l’A-stabilité du schéma de Cranck-Nicolson \(u_{n+1} = u_{n}+\left(\varphi_n+\varphi_{n+1}\right)\times h/2\).

Code

beta = sp.symbols(r'\beta', real=True, positive=True)
h = sp.symbols('h', real=True, positive=True)
tn, un, unp1 = sp.symbols('t_n u_n u_{n+1}', real=True)
phi = lambda t,y: -beta*y
k1 = phi(tn, un)
k2 = phi(tn + h, unp1)
sol = sp.solve(  sp.Eq(unp1, un + h/2*(k1+k2)), unp1)[0]
sp.Eq(unp1, sol)

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

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, \) soit encore \( u_n = [q(\beta h)]^n u_0 \qquad\text{avec } q(\beta h) = \frac{2-\beta h}{2+\beta h}. \) Par conséquent, \(\lim\limits_{n\to+\infty}u_n=0\) si et seulement si \( \left|q(\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} \xrightarrow[x\to+\infty]{}-1^+\). 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(x)<1\) on a (pour \(x>0\)) \( 1-2\frac{x}{2+x}>0 \iff x<2. \)

Code

x = np.linspace(0, 10, 101)
q = 1-2*x/(2+x)
plt.figure(figsize=(8, 5))
plt.grid()
plt.plot(x, q, label='$q(x)$')
plt.ylim(-2, 2) # On zoome sur la zone intéressante [-1, 1]
plt.xlim(x.min(), x.max())
plt.axhline(y=1, color='black', linestyle='-', alpha=0.5)
plt.axhline(y=-1, color='black', linestyle='-', alpha=0.5)
plt.fill_between(x, -1, 1, color='gray', alpha=0.1, label='Zone de A-Stabilité')
plt.axhline(0, color='magenta', linestyle=':', label='$y=0$')
plt.scatter(2, 0, color='black', label='$q(2) = 0$')
plt.xlabel(r'$x=\beta h$')
plt.ylabel('$q(x)$')
plt.title('CN : $q(x) = 1 - 2x/(2+x)$')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.show()

Exercice : étudier l’A-stabilité du schéma d’Euler modifié \(u_{n+1} = u_{n}+h \varphi\left(t_n+\frac{h}{2},u_n+\frac{h}{2}\varphi_n\right)\).

Code

beta = sp.symbols(r'\beta', real=True, positive=True)
h = sp.symbols('h', real=True, positive=True)
tn, un, unp1 = sp.symbols('t_n u_n u_{n+1}', real=True)
phi = lambda t,y: -beta*y
k1 = phi(tn, un)
k2 = phi(tn + h/2, un+h/2*k1)
sol = sp.solve(  sp.Eq(unp1, un + h*k2), unp1)[0]
sp.Eq(unp1, sol)

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

Pour \(\varphi(t,y)=-\beta y\), on a \( \varphi(t_n,u_n)=-\beta u_n \qquad\text{et}\qquad \varphi\left(t_n+\frac{h}{2},u_n+\frac{h}{2}\varphi_n\right)=-\beta\left( u_n+\frac{h}{2}(-\beta u_n) \right). \) ainsi 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 \) soit encore \( u_n = [q(\beta h)]^n u_0 \qquad\text{avec } q(\beta h) = 1 - \beta h + \frac{1}{2}(\beta h)^2. \)

Par conséquent, \(\lim\limits_{n\to+\infty}u_n=0\) si et seulement si \( \left| q(\beta h) \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 ci-dessous. :

Code

x = np.linspace(0, 3, 101)
q = 1-x+0.5*x**2
plt.figure(figsize=(8, 5))
plt.grid()
plt.plot(x, q, label='$q(x)$')
plt.ylim(-2, 2) # On zoome sur la zone intéressante [-1, 1]
plt.xlim(x.min(), x.max())
plt.axhline(y=1, color='black', linestyle='-', alpha=0.5)
plt.axhline(y=-1, color='black', linestyle='-', alpha=0.5)
plt.fill_between(x, -1, 1, color='gray', alpha=0.1, label='Zone de A-Stabilité')
plt.axhline(0, color='magenta', linestyle=':', label='$y=0$')
# plt.scatter(2, 0, color='black', label='$q(2) = 0$')
plt.xlabel(r'$x=\beta h$')
plt.ylabel('$q(x)$')
plt.title('EM : $q(x) = 1 - x + x^2/2$')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.show()

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 Heun \(u_{n+1} = u_n + \left(\varphi(t_n, u_n) + \varphi(t_{n+1}, u_{n}+h\varphi(t_n, u_n))\right)\times h/2\).

Code

beta = sp.symbols(r'\beta', real=True, positive=True)
h = sp.symbols('h', real=True, positive=True)
tn, un, unp1 = sp.symbols('t_n u_n u_{n+1}', real=True)
phi = lambda t,y: -beta*y
k1 = phi(tn, un)
k2 = phi(tn + h/2, un+h*k1)
sol = sp.solve(  sp.Eq(unp1, un + h/2*(k1+k2)), unp1)[0]
sp.Eq(unp1, sol)

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

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 \(u_{n+1} = u_n + \left(\varphi(t_n, u_n) + 4\varphi(t_{n}+\frac{h}{2}, u_{n}+\frac{h}{2}\varphi_n) + \varphi(t_{n+1}, u_{n+1})\right)\times h/2\).

Code

beta = sp.symbols(r'\beta', real=True, positive=True)
h = sp.symbols('h', real=True, positive=True)
tn, un, unp1 = sp.symbols('t_n u_n u_{n+1}', real=True)
phi = lambda t,y: -beta*y
k1 = phi(tn, un)
k2 = phi(tn + h/2, un + h/2*k1)
k3 = phi(tn + h, unp1)
sol = sp.solve(  sp.Eq(unp1, un + h/6*(k1 + 4*k2 + k3)), unp1)[0]
sp.Eq(unp1, sol)

\(\displaystyle u_{n+1} = \frac{u_{n} \left(2 \beta^{2} h^{2} - 5 \beta h + 6\right)}{\beta h + 6}\)

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{ 6-5\beta h+2(\beta h)^2 }{ 6+\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{ 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.

Code

x = np.linspace(0, 4, 101)
q = (6-5*x+2*x**2)/(6+x)
plt.figure(figsize=(8, 5))
plt.grid()
plt.plot(x, q, label='$q(x)$')
plt.ylim(-2, 2) # On zoome sur la zone intéressante [-1, 1]
plt.xlim(x.min(), x.max())
plt.axhline(y=1, color='black', linestyle='-', alpha=0.5)
plt.axhline(y=-1, color='black', linestyle='-', alpha=0.5)
plt.fill_between(x, -1, 1, color='gray', alpha=0.1, label='Zone de A-Stabilité')
plt.axhline(0, color='magenta', linestyle=':', label='$y=0$')
# plt.scatter(2, 0, color='black', label='$q(2) = 0$')
plt.xlabel(r'$x=\beta h$')
plt.ylabel('$q(x)$')
plt.title('SI : $q(x) = \\frac{6-5x+2x^2}{6+x}$')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.show()

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 \(u_{n+1} = u_n + \left(\varphi(t_n, u_n) + 4\varphi(t_{n}+\frac{h}{2}, u_{n}+\frac{h}{2}\varphi_n) + \varphi(t_{n+1}, u_{n}+h\varphi(t_n, u_n))\right)\times h/2\).

Code

beta = sp.symbols(r'\beta', real=True, positive=True)
h = sp.symbols('h', real=True, positive=True)
tn, un, unp1 = sp.symbols('t_n u_n u_{n+1}', real=True)
phi = lambda t,y: -beta*y
k1 = phi(tn, un)
k2 = phi(tn + h/2, un + h/2*k1)
k3 = phi(tn + h, un+h*k1)
sol = sp.solve(  sp.Eq(unp1, un + h/6*(k1 + 4*k2 + k3)), unp1)[0]
sp.Eq(unp1, sol)

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

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.

Code

x = np.linspace(0, 1, 101)
q = 1-x+3*x**2/2
plt.figure(figsize=(8, 5))
plt.grid()
plt.plot(x, q, label='$q(x)$')
plt.ylim(-2, 2) # On zoome sur la zone intéressante [-1, 1]
plt.xlim(x.min(), x.max())
plt.axhline(y=1, color='black', linestyle='-', alpha=0.5)
plt.axhline(y=-1, color='black', linestyle='-', alpha=0.5)
plt.fill_between(x, -1, 1, color='gray', alpha=0.1, label='Zone de A-Stabilité')
plt.axhline(0, color='magenta', linestyle=':', label='$y=0$')
# plt.scatter(2, 0, color='black', label='$q(2) = 0$')
plt.xlabel(r'$x=\beta h$')
plt.ylabel('$q(x)$')
plt.title('SI : $q(x) = 1 - x + \\frac{3x^2}{2}$')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.show()

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

Cartographie des zones de stabilité pour les schéma EE, EM, Heun (explicites) et EI, CN, SI (implicites)

Étude de la A-stabilité pour un schéma à deux pas

Les schémas à deux pas sont des méthodes numériques dans lesquelles la valeur de \(u_{n+1}\) dépend à la fois de \(u_n\) et de \(u_{n-1}\). Ces méthodes, appliquées au toy-model, conduisent à des relations de la forme \(u_{n+1} = a(\beta h) u_n + b(\beta h) u_{n-1},\) où \(a\) et \(b\) sont des fonctions de la variable \(x = \beta h\). La condition de A-stabilité se traduit alors par l’exigence que les racines de l’équation caractéristique associée à la relation de récurrence soient de module strictement inférieur à 1.

Remarques 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 :

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

Dans tous les cas, la suite \(u_n\) tend vers 0 si et seulement si les racines du polynôme caractéristique sont de module strictement inférieur à 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}\).

Code

beta = sp.symbols(r'\beta', real=True, positive=True)
h = sp.symbols('h', real=True, positive=True)
tn, unm1, un, unp1 = sp.symbols('t_n u_{n-1} u_n u_{n+1}', real=True)
phi = lambda t,y: -beta*y
k1 = phi(tn, un)
k2 = phi(tn -h, unm1)
sol = sp.solve(  sp.Eq(unp1, un + h/2*(3*k1-k2) ), unp1)[0]
sp.Eq(unp1, sol)

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

Code

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)
plt.figure(figsize=(8, 5))
plt.grid()
yy1 = r1(xx)
yy2 = r2(xx)
plt.plot(xx,yy1, label=r'$r_1(x)$')
plt.plot(xx,yy2,  label=r'$r_2(x)$')
plt.ylim(-2, 2) # On zoome sur la zone intéressante [-1, 1]
plt.xlim(xx.min(), xx.max())
plt.axhline(y=1, color='black', linestyle='-', alpha=0.5)
plt.axhline(y=-1, color='black', linestyle='-', alpha=0.5)
plt.fill_between(xx, -1, 1, color='gray', alpha=0.1, label='Zone de A-Stabilité')
plt.xlabel(r'$x=\beta h$')
plt.title('AB2')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.show()

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

Code

beta = sp.symbols(r'\beta', real=True, positive=True)
h = sp.symbols('h', real=True, positive=True)
tn, unm1, un, unp1 = sp.symbols('t_n u_{n-1} u_n u_{n+1}', real=True)
phi = lambda t,y: -beta*y
k1 = phi(tn, un)
k2 = phi(tn-h, unm1)
km1 = phi(tn+h, unp1)
sol = sp.solve(  sp.Eq(unp1, un + h/12*(8*k1-k2+5*km1) ), unp1)[0]
sp.Eq(unp1, sol)

\(\displaystyle u_{n+1} = \frac{- 8 \beta h u_{n} + \beta h u_{n-1} + 12 u_{n}}{5 \beta h + 12}\)

Le schéma d’Adam-Moulton à 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} \)

Code

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)
plt.figure(figsize=(8, 5))
plt.grid()
yy1 = r1(xx)
yy2 = r2(xx)
plt.plot(xx,yy1, label=r'$r_1(x)$')
plt.plot(xx,yy2,  label=r'$r_2(x)$')
plt.ylim(-2, 2) # On zoome sur la zone intéressante [-1, 1]
plt.xlim(xx.min(), xx.max())
plt.axhline(y=1, color='black', linestyle='-', alpha=0.5)
plt.axhline(y=-1, color='black', linestyle='-', alpha=0.5)
plt.fill_between(xx, -1, 1, color='gray', alpha=0.1, label='Zone de A-Stabilité')
plt.xlabel(r'$x=\beta h$')
plt.title('AM2')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.show()

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

Code

beta = sp.symbols(r'\beta', real=True, positive=True)
h = sp.symbols('h', real=True, positive=True)
tn, unm1, un, unp1 = sp.symbols('t_n u_{n-1} u_n u_{n+1}', real=True)
phi = lambda t,y: -beta*y
km1 = phi(tn+h, unp1)
sol = sp.solve(  sp.Eq(unp1, sp.Rational(4,3)*un-sp.Rational(1,3)*unm1+sp.Rational(2,3)*h*km1 ), unp1)[0]
sp.Eq(unp1, sol)

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

Les racines du polynôme caractéristique sont données par \(r_{1,2} = \frac{\gamma \pm \sqrt{\gamma^2 + 4\eta}}{2} = \frac{\frac{4}{3+2\beta h} \pm \sqrt{\left(\frac{4}{3+2\beta h}\right)^2 - \frac{4}{3+2\beta h}}}{2}\)

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

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

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 de module strictement inférieur à 1 pour tout \(x>0\). Par conséquent, la méthode est A-stable pour tout \(h>0\).

Code

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

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

poly = sp.Poly(r**2-gamma(x)*r-eta(x), r)
display(Math(f"p_x(r) = {sp.latex(poly.as_expr())}"))

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

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

roots = sp.solve(poly, 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.discriminant(), x)[0]
display(Markdown(f"Point de bifurcation (racines complexes) : $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$)")
plt.xlabel('$x = h\\beta$')
plt.ylabel('$|r(x)|$')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.grid(True, alpha=0.3)
plt.ylim(0, 1.2)
plt.show()

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

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

Point de bifurcation (racines fusionnent) : \(x = \frac{1}{2}\)

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.

Code

def affichage(schema):
    # CALCUL DE LA SOLUTION SUR CHAQUE GRILLE
    U = [ schema(phi,tt,y0) for tt in T ] 
        
    # AFFICHAGE
    plt.figure(1, figsize=(8, 5))
    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(bbox_to_anchor=(1.05, 1), loc='upper left')
    plt.grid();

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

Code

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

Code

# SCHEMA
def EE(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_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


affichage(EE)

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

Code

# SCHEMA
def EI(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_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



affichage(EI)

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.

Code

# SCHEMA
def EM(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_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


affichage(EM)

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

Code

# SCHEMA
def CN(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_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   


affichage(CN)

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