Convergence et A-stabilité des schémas multi-pas linéaires

Auteur·rice

Gloria Faccanoni

Date de publication

29 mars 2026

Code

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

Rappels

Problème de Cauchy de référence $ \begin{cases} y'(t)=\varphi(t,y(t)), & t\in[t_0,T],\\ y(t_0)=y_0, \end{cases} $ avec solution unique $y(t)$ et approximation $u_n \approx y(t_n)$ sur la grille $ t_n = t_0 + n h, \quad n=0,\dots,N, \quad h = \frac{T-t_0}{N}. $

Schémas multi-pas linéaires (MPL)

Notons $\varphi_{n-j} = \varphi(t_{n-j}, u_{n-j})$.
Soit $p\ge0$ et posons $q=p+1$, $q$ étant le nombre de pas du schéma.
Soient les $2q+1$ coefficients : $\{a_j\}_{j=0}^p$, $\{b_j\}_{j=-1}^p$. Une méthode multi-pas linéaire à $q$ pas s’écrit $ \begin{aligned} u_{n+1} &= \sum_{j=0}^p a_ju_{n-j} + h\sum_{j=-1}^p b_j\varphi_{n-j} \\ &= \sum_{j=0}^p \Big(a_ju_{n-j} + b_j\varphi_{n-j} \Big) + hb_{-1}\varphi_{n+1}, \end{aligned} \qquad n=p,\dots,N-1. $ Initialisation : valeurs nécessaires : $u_0,\dots,u_p$ avec $u_0 = y_0$ connu et $u_1,\dots,u_p$ calculés via une autre méthode.

Elle est explicite si $b_{-1}=0$, et implicite sinon.

Généralités : convergence et stabilités

La convergence se décompose en une propriété locale (consistance) et une propriété de propagation des erreurs (zéro-stabilité).

Théorème de Lax–Richtmyer : Consistance + zéro-stabilité $\iff$ convergence

La convergence étudie l’évolution des erreurs lorsque $N\to\infty$ sur un intervalle borné (et donc $h\to 0$). Pour des intervalles très grands (voir infinis) la zéro-stabilité ne suffit plus pour garantir un bon comportement numérique. Il faut alors introduire la notion de A-stabilité.
Stabilité absolue (A-stabilité) : la solution approchée $u_n$ doit tendre vers zéro pour tout pas $h>0$ fixé (ou sous certaines conditions sur $h$) lorsque le schéma est appliqué au toy-model $y'(t)=-\beta y(t)$ et $y(0)=1$ avec $\beta>0$.

Polynômes caractéristiques

Dans ce chapitre, nous allons étudier d’abord la convergence (à travers les notions de consistance et de zéro-stabilité) puis la A-stabilité d’une méthode multi-step.

Dans l’analyse de ces méthodes, trois polynômes jouent un rôle clé. Pour cela, commençons par réécrire le schéma multi-step comme suit :

$ u_{n+1} -\sum_{j=0}^p a_ju_{n-j}= h\left(b_{-1}\varphi_{n+1}+\sum_{j=0}^p b_j\varphi_{n-j}\right). $

Définition.

Soit une méthode MPL. On lui associe les trois polynômes caractéristiques suivants définis dans $\mathbb{C}$ :

premier polynôme caractéristique

$ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j} $

second polynôme caractéristique

$ \sigma(r)=b_{-1}r^{p+1}+\sum_{j=0}^p b_jr^{p-j} $

polynôme de A-stabilité (une famille de polynômes paramétrés par $z$)

$ p_z(r) = \varrho(r) - z \sigma(r) = (1 - z b_{-1})r^{p+1} - \sum_{j=0}^p (a_j + z b_j) r^{p-j} $

Formellement,

pour construire le premier polynôme caractéristique on ne regarde que le LHS et les $u_{n-j}$ deviennent des $r^{p-j}$,
pour construire le second polynôme caractéristique on ne regarde que le RHS divisé par $h$ et les $\varphi_{n-j}$ deviennent des $r^{p-j}$.

On verra que les racines du polynôme $\varrho$ sont liées à la zéro-stabilité tandis que celles du polynôme $p_z$ sont liées à la stabilité absolue (en association avec un choix particulier de $z=h\lambda$).

Exemple : méthodes classiques à $q=1$ pas :

$ \begin{array}{|l|l|c|c|c|c|c|c|} \hline \text{Méthode} & \text{Formule} & p & a_0 & b_{-1} & b_0 & \varrho(r) & \sigma(r) \\ \hline \text{Euler explicite} & u_{n+1}=u_n+h\varphi_n & 0 & 1 & 0 & 1 & r-1 & 1 \\ \text{Euler implicite} & u_{n+1}=u_n+h\varphi_{n+1} & 0 & 1 & 1 & 0 & r-1 & r \\ \text{Crank-Nicolson} & u_{n+1}=u_n+\frac{h}{2}(\varphi_n+\varphi_{n+1}) & 0 & 1 & \frac{1}{2} & \frac{1}{2} & r-1 & \frac{1}{2}(r+1) \\ \hline \end{array} $

Exemple : méthodes à $q>1$ pas :

$ \begin{array}{|l|l|c|l|l|} \hline \text{Méthode} & \text{Formule} & p & \varrho(r) & \sigma(r) \\ \hline \text{AB}_q \text{ ou } \text{AM}_q & u_{n+1} = u_n + h \sum_{j=-1}^{p} b_j \varphi_{n-j} & q-1 & r^{q}(r-1) & \sum_{j=-1}^{p} b_j r^{p-j} \\ \hline \text{N}_q \text{ ou } \text{MS}_q & u_{n+1} = u_{n-1} + h \sum_{j=-1}^{p} b_j \varphi_{n-j} & q-1 & r^{q-1}(r^2-1) & \sum_{j=-1}^{p} b_j r^{p-j} \\ \hline \text{BDF}_2 & u_{n+1} - \frac{4}{3}u_n + \frac{1}{3}u_{n-1} = \frac{2}{3}h \varphi_{n+1} & q-1 & (r-1)(r-\frac{1}{3}) & \frac{2}{3}r^2 \\ \hline \end{array} $

Convergence d’un schéma multi-pas linéaire

Théorème (Dahlquist 1956, Lax–Richtmyer)

Une méthode multi-pas linéaire est convergente si et seulement si elle est consistante, (zero)-stable et les erreurs sur les conditions initiales sont $O(h)$.

La convergence est d’ordre $\omega$ si la méthode est convergente, consistante d’ordre $\omega$, et les erreurs sur les conditions initiales sont $O(h^{\omega})$.

Consistance

Une méthode est dite consistante (ou consistante d’ordre au moins 1) si l’erreur de troncature $\tau(h)=\max_n |\tau_n(h)|$ tend vers zéro quand $h$ tend vers zéro, où $\tau_n(h)$ est l’erreur de troncature locale de la méthode multi-pas définie par

$ \tau_{n+1}(h)=\frac{1}{h}\left(y_{n+1}-\sum_{j=0}^p a_jy_{n-j}-h\sum_{j=0}^p b_j\varphi(t_{n-j},y_{n-j})-hb_{-1}\varphi(t_{n+1},y_{n+1})\right). $

Si, de plus, $\tau(h)=\mathscr{O}(h^\omega)$, on dit que la méthode est consistante d’ordre $\omega$.

Pour évaluer la consistance et son ordre, nous allons introduire une nouvelle quantité $\xi(i)$.

Définition.

Pour $i=0,1,\dots$ introduisons la quantité

$ \xi(i) = \sum_{j=0}^p (-j)^{i}a_j+i\sum_{j=-1}^p (-j)^{i-1}b_j, $

ayant posé $(-j)^0=1$ même pour $j=0$.

Pour bien comprendre cette définition, on peut écrire les premières valeurs de $\xi(i)$ :

$ \begin{aligned} \xi(0) &= \sum_{j=0}^p a_j, \\ \xi(1) &= -\sum_{j=0}^p j a_j + \sum_{j=-1}^p b_j, \\ \xi(2) &= \sum_{j=0}^p j^{2} a_j - 2 \sum_{j=-1}^p j b_j, \\ \xi(3) &= -\sum_{j=0}^p j^{3} a_j + 3 \sum_{j=-1}^p j^{2} b_j, \\ &\quad \vdots \end{aligned} $

Exemple avec un schéma à $q=4$ pas (et donc $p=3$) :

$\xi(0)=\Big(a_0+a_1+a_2+a_3\Big)$
$\xi(1)=-\Big(a_1+2a_2+3a_3\Big)+b_{-1}+\Big(b_1+b_2+b_3\Big)$
$\xi(2)=\Big(a_1+4a_2+9a_3\Big)+2b_{-1}-2\Big(b_1+2b_2+3b_3\Big)$
$\xi(3)=-\Big(a_1+8a_2+27a_3\Big)+3b_{-1}+3\Big(b_1+4b_2+9b_3\Big)$

Il n’y a pas de difficulté à calculer ces quantités pour un schéma donné, cependant il faut faire attention à ne pas se tromper dans les calculs, notamment les signes des termes. Pour vous aider, voici un petit script Python qui calcule ces quantités de façon automatique pour un schéma donné. Il suffit de rentrer les valeurs de $p$ et jusqu’à quel ordre vous voulez calculer les $\xi(i)$.

Code

## ===========================================================================
## Calcul et affichage de xi(i) pour un i donné pour un schema de pas q=p+1
## ===========================================================================
def xi(i,p):
    q = p+1
    aa = sp.symbols(f'a_0:{q}')
    bb = sp.symbols(f'b_0:{q}')
    bm1 = sp.Symbol('b_{-1}')
    p = len(aa)
    sa = sum( [ (-j)**i * aa[j] for j in range(p) ] )
    sb = bm1+sum( [(-j)**(i-1) * bb[j] for j in range(1,p)] )
    ## si j=0 et i=0, sympy calcule (-j)**i =1 et on a bien a_0 dans la somme
    ## mais si j=0 et i=1 sympy ne calcule pas b_0 et il faut l'ajouter à la main
    if i==1: sb += bb[0]
    xi_rhs  = (sa).factor()+(i*sb).factor()                  ## calcul de xi(i)
    display(Markdown(f"$\\xi({i}) = {sp.latex(xi_rhs)}$") )  ## affichage de xi(i)
    ## return xi_rhs                                         

## ==========================================================================
## Calcul et affichage des xi(i) pour i=0,...,omega
## ==========================================================================
def xi_cond(p,omega):
    ## j = 0...p ## coeffs du schema
    ## i = 0, ..., omega ordre du schema
    q = p+1
    display(Markdown(f"**Pour ${p=}$, on a une méthode à $q={p+1}$ pas et on affiche $ξ(i)$ pour $i=0,...,{omega}$**"))
    for i in range(omega+1):
        xi(i,p)

Code

## ==========================================================================    
## Exemple d'utilisation
## Calcul des xi(i) pour i=0,...,omega
## ==========================================================================

xi_cond( p=1, omega=3 )

Pour $p=1$, on a une méthode à $q=2$ pas et on affiche $ξ(i)$ pour $i=0,...,3$

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

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

$\xi(2) = a_{1} - 2 \left(b_{1} - b_{-1}\right)$

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

La proposition suivante donne une condition nécessaire et suffisante pour la consistance d’une méthode multi-pas.

Proposition [ consistance et ordre de consistance ]

La méthode est consistante (ou, plus précisément, consistante d’ordre au moins 1) ssi $\xi(0)=\xi(1)=1$.

La méthode est consistante d’ordre $\omega\ge2$ ssi $\xi(0)=\xi(1)=\dots=\xi(\omega)=1$ et $\xi(\omega+1)\neq 1$.

Remarques :

La condition de consistance d’ordre au moins 1 est une condition nécessaire pour la convergence d’une méthode numérique.
La condition de consistance d’ordre $\omega$ signifie que toutes les valeurs de $\xi(i)$ sont égales à 1 pour $i$ allant de $0$ à $\omega$, et non pas seulement $\xi(\omega)$.

Preuve. On notera $y_{n-j} = y(t_{n-j})$ , $\varphi(t_{n-j}, y_{n-j}) = \varphi(t_{n-j}, y(t_{n-j})) = y'(t_{n-j})$ et $u_{n-j}\approx y(t_{n-j})$ pour $j=0,\dots,p$.

On notera $u^*_{n+1}$ la valeur que prendrait $u_{n+1}$ si on remplaçait $u_{n-j}$ par $y_{n-j}$ et $\varphi(t_{n-j}, u_{n-j})$ par $\varphi(t_{n-j}, y_{n-j})$ dans le schéma numérique, soit $ u^*_{n+1} = \sum_{j=0}^p a_j y_{n-j} + h \sum_{j=-1}^p b_j \varphi(t_{n-j}, y_{n-j}). $

En développant $y_{n-j}$ et $\varphi(t_{n-j},y_{n-j})$ en série de Taylor autour de $t_n$, on trouve pour tout $n\ge p$ $ \begin{aligned} y_{n-j} &= y_n - j h y'(t_n) + O(h^2), \\ \varphi(t_{n-j}, y_{n-j}) &= y'(t_{n-j}) = y'(t_n) + O(h). \end{aligned} $

En remplaçant ces développements dans le schéma numérique $ \begin{aligned} \sum_{j=0}^p a_j y_{n-j} &= \sum_{j=0}^p a_j \left(y_n - j h y'(t_n) + O(h^2)\right) = \left(\sum_{j=0}^p a_j\right) y_n - h \left(\sum_{j=0}^p j a_j\right) y'(t_n) + O(h^2),\\ \sum_{j=-1}^p b_j \varphi(t_{n-j}, y_{n-j}) &= \sum_{j=-1}^p b_j \left(y'(t_n) + O(h)\right) = \left(\sum_{j=-1}^p b_j\right) y'(t_n) + O(h), \end{aligned} $ on trouve $ \begin{aligned} u^*_{n+1} &= \sum_{j=0}^p a_j y_{n-j} + h\sum_{j=-1}^p b_j \varphi(t_{n-j}, y_{n-j}) \\ &= \left(\sum_{j=0}^p a_j\right) y_n - h \left(\sum_{j=0}^p j a_j\right) y'(t_n) + h\left(\sum_{j=-1}^p b_j\right) y'(t_n) + O(h^2) \\ &= \left(\sum_{j=0}^p a_j\right) y_n + h \left(\sum_{j=0}^p (-j) a_j + \sum_{j=-1}^p b_j\right) y'(t_n) + O(h^2) \\ &= \xi(0) y_n + h \xi(1) y'(t_n) + O(h^2). \end{aligned} $

L’erreur de troncature locale est alors donnée par $ \begin{aligned} \tau_{n+1}(h) = \frac{y_{n+1}-u^*_{n+1}}{h} &= \frac{y_{n}+hy'(t_n)}{h} - \frac{\xi(0) y_n + h \xi(1) y'(t_n)}{h} + O(h)\\ &= \left(1-\xi(0)\right) \frac{y_n}{h} +\left(1-\xi(1)\right) y'(t_n) + O(h)\\ \end{aligned} $

On a $\tau_{n+1}(h)=O(h)$ pour tout $n\ge p$ ssi $\xi(0)=\xi(1)=1$. On conclut alors que le schéma est consistant (d’ordre au moins 1) ssi $\xi(0)=\xi(1)=1$.

Ces calculs peuvent être poursuivis pour trouver les conditions de consistance d’ordre $\omega\ge2$. Par exemple, pour l’ordre 2 : $ \begin{aligned} y_{n-j} &= y_n - j h y'(t_n) + \frac{j^2 h^2}{2} y''(t_n) + O(h^3), \\ \varphi(t_{n-j}, y_{n-j}) &= y'(t_{n-j}) = y'(t_n) - j h y''(t_n) + O(h^2). \end{aligned} $

En remplaçant ces développements dans le schéma numérique $ \begin{aligned} \sum_{j=0}^p a_j y_{n-j} &= \sum_{j=0}^p a_j \left(y_n - j h y'(t_n) + \frac{j^2 h^2}{2} y''(t_n) + O(h^3)\right) \\ &= \left(\sum_{j=0}^p a_j\right) y_n - h \left(\sum_{j=0}^p j a_j\right) y'(t_n) + \frac{h^2}{2} \left(\sum_{j=0}^p j^2 a_j\right) y''(t_n) + O(h^3),\\ \sum_{j=-1}^p b_j \varphi(t_{n-j}, y_{n-j}) &= \sum_{j=-1}^p b_j \left(y'(t_{n-j}) = y'(t_n) - j h y''(t_n) + O(h^2)\right) \\ &= \left(\sum_{j=-1}^p b_j\right) y'(t_n) - h \left(\sum_{j=-1}^p j b_j\right) y''(t_n) + O(h^2), \end{aligned} $ on trouve $ \begin{aligned} u^*_{n+1} &= \sum_{j=0}^p a_j y_{n-j} + h\sum_{j=-1}^p b_j \varphi(t_{n-j}, y_{n-j}) \\ &= \left(\sum_{j=0}^p a_j\right) y_n + h \left(\sum_{j=0}^p -j a_j\right) y'(t_n) + \frac{h^2}{2} \left(\sum_{j=0}^p (-j)^2 a_j\right) y''(t_n) \\ &+ h\left(\sum_{j=-1}^p b_j\right) y'(t_n) - h^2\left(\sum_{j=-1}^p (-j) b_j\right) y''(t_n) + O(h^3) \\ &= \xi(0) y_n + \xi(1) h y'(t_n) + \xi(2) \frac{h^2}{2} y''(t_n) + O(h^3). \end{aligned} $

L’erreur de troncature locale est alors donnée par $ \begin{aligned} \tau_{n+1}(h) = \frac{y_{n+1}-u^*_{n+1}}{h} &= \frac{y_{n}+hy'(t_n)+\frac{h^2}{2}y''(t_n)}{h} - \frac{\xi(0) y_n + h \xi(1) y'(t_n) + \frac{h^2}{2} \xi(2) y''(t_n)}{h}+ O(h^2) \\ &= (1-\xi(0)) \frac{y_n}{h} +(1-\xi(1)) y'(t_n) + (1-\xi(2)) \frac{h}{2} y''(t_n) + O(h^2)\\ \end{aligned} $

On conclut que $\tau_{n+1}(h)=O(h^2)$ si et seulement si $\xi(0)=\xi(1)=\xi(2)=1$.

Exercice : montrer que

$ \begin{cases} \xi(0)=1\\ \xi(1)=1 \end{cases} \iff \begin{cases} \varrho(1)=0\\ \varrho'(1)=\sigma(1) \end{cases} $

Remarque : il s’agit de démontrer que le premier système est équivalent au second. Cela ne signifie pas que chaque équation du premier système est équivalente à l’équation correspondante du second système, ce qui n’est pas le cas.

Correction

Rappelons les définitions de $\xi(0)$, $\xi(1)$ et des polynômes $\varrho$ et $\sigma$. Calculons ensuite $\varrho'$ :

$\begin{aligned} \xi(0)&=\sum_{j=0}^p a_j & \xi(1)&=\sum_{j=0}^p -j a_j + \sum_{j=-1}^p b_j \\ \varrho(r)&=r^{p+1}-\sum_{j=0}^p a_jr^{p-j} & \varrho'(r)&=(p+1)r^{p}-\sum_{j=0}^p (p-j)a_jr^{p-j-1} \\ \sigma(r)&=b_{-1}r^{p+1}+\sum_{j=0}^p b_jr^{p-j}=\sum_{j=-1}^p b_jr^{p-j} && \end{aligned}$

Évaluons $\varrho$ et $\varrho'-\sigma$ en $1$:

$\begin{aligned} \varrho(1) &=1-\sum_{j=0}^p a_j = 1-\xi(0) \\ \varrho'(1)-\sigma(1) &=(p+1)-\sum_{j=0}^p (p-j)a_j -\sum_{j=-1}^p b_j= p\left(1-\sum_{j=0}^p a_j\right) +1-\sum_{j=0}^p -ja_j -\sum_{j=-1}^p b_j= p(1-\xi(0))+1-\xi(1) \end{aligned}$

donc

$ \varrho(1)=0 \qquad\iff\qquad 1-\xi(0)=0 $

$ \varrho'(1)-\sigma(1) = 0 \qquad\iff\qquad p(1-\xi(0))+1-\xi(1)=0 $

Bien noter que

$\xi(0)=1\quad\iff\quad\varrho(1)=0$ : la première équation du premier système est effectivement équivalente à la première équation du second système;
$\xi(1)=1\quad\xcancel{\iff}\quad\varrho'(1)=\sigma(1)$ : la deuxième équation du premier système n’est pas équivalente à la deuxième équation du second système ! Ce que l’on peut dire est :
- $\xi(0)=\xi(1)=1\quad\implies\quad\varrho'(1)=\sigma(1)$ : la deuxième équation du second système est une conséquence de tout le premier système ;
- $\begin{cases}\varrho'(1)=\sigma(1)\\ \xi(0)=1\end{cases}\implies\quad \xi(1)=1$ : la deuxième équation du premier système est une conséquence de la deuxième équation du second système et de la première équation du premier système ;
- $\begin{cases}\varrho'(1)=\sigma(1)\\ \xi(1)=1\end{cases}\implies\quad \xi(0)=1$ : la première équation du premier système est une conséquence de la première équation du second système et de la deuxième équation du premier système.

Exercice : étudier la consistance et l’ordre de consistance du schéma multi-pas linéaire $\text{AB}_2$

$ u_{n+1} = u_n + \frac{h}{2} \left( 3\varphi_n - \varphi_{n-1} \right) $

Correction

$p=1$ $\leadsto$ $q=2$ : c’est un schéma de la forme

$ u_{n+1}= \underset{= \, 1}{a_0} u_n + \underset{= \, 0}{a_1} u_{n-1} + h \left( \underset{= \, \frac{3}{2}}{b_0} \varphi_n + \underset{= \, -\frac{1}{2}}{b_1} \varphi_{n-1} \right) + h \underset{= \, 0}{b_{-1}} \varphi_{n+1} $

Notons que $b_{-1}=0$ $\leadsto$ le schéma est explicite.

On vérifie la consistance en calculant les premières deux valeurs de $ (i) $. Pour afficher $\xi(0)$ et $\xi(1)$ dans le cas général, on peut utiliser la fonction précédente :

Code

xi_cond(p=1,omega=1)

Pour $p=1$, on a une méthode à $q=2$ pas et on affiche $ξ(i)$ pour $i=0,...,1$

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

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

Maintenant on remplace les valeurs des coefficients dans les formules de $\xi(i)$ pour obtenir les valeurs de $\xi(0)$ et $\xi(1)$.

$\begin{aligned} \xi(0) &= 1 + 0 = 1 \\ \xi(1) &= -0+\frac{3}{2}-\frac{1}{2}+0 = 1 \end{aligned}$

Le schéma est consistante, c’est à dire qu’il est d’ordre au moins 1.

Pour vérifier si le schéma est consistant d’ordre 2, on doit calculer $\xi(2)$ et vérifier si $\xi(2)=1$. Pour cela, affichons $\xi(i)$ pour $i=2$ dans le cas général :

Code

xi(i=2,p=1) ## xi(2) avec p=1

$\xi(2) = a_{1} - 2 \left(b_{1} - b_{-1}\right)$

On évalue alors $\xi(2)$ :

$ \xi(2) = 0 - 2\left(-\frac{1}{2}-0\right) = 1 $

donc le schéma est consistant d’ordre au moins 2.

Pour vérifier si le schéma est consistant d’ordre au moins 3, on doit calculer $\xi(3)$ et vérifier si $\xi(3)=1$. Pour cela, affichons $\xi(i)$ pour $i=3$ dans le cas général :

Code

xi(i=3,p=1) ## xi(3) avec p=1

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

On évalue alors $\xi(3)$ :

$ \xi(3) = 0 + 3\left(-\frac{1}{2}+0\right) \neq 1 $

donc le schéma n’et pas d’ordre 3. On conclut qu’il est consistant d’ordre exactement 2.

Zéro-stabilité

Proposition [zéro-stabilité et condition des racines]

Soient $r_0, r_1,\dots, r_p$ les $p+1$ racines (dans $\mathbb{C}$) du premier polynôme caractéristique $\varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}$.

La méthode est zéro-stable ssi

$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ r_j\text{ a multiplicité}>1\quad\implies \quad |r_j|<1. \end{cases} $

Rappel : si le schéma est consistent, alors $r=1$ est une racine de $\varrho$.

Explication :

La première condition signifie que, dans le plan complexe, toutes les racines doivent être contenues dans le cercle unitaire, c’est-à-dire que leur module doit être inférieur ou égal à 1. Si une racine est située à l’extérieur du disque unité, la méthode n’est pas zéro-stable.
La seconde condition stipule qu’une racine ayant une multiplicité supérieure à 1 ne peut pas se trouver sur le bord du disque unitaire. Si une racine de multiplicité supérieure à 1 est sur le bord du disque unitaire, alors la méthode est dite semi-stable, mais pas zéro-stable.

Exemple : étudier la zéro-stabilité des méthodes Euler explicite, Euler implicite, Crank-Nicolson, de Adams à $q$ pas, de Nystrom et Milne-Simpson à $q$ pas et des méthodes BDF.

Méthodes à un pas (Euler explicite, Euler implicite, Crank-Nicolson) :

Ces méthodes ont pour polynôme caractéristique :

$ \varrho(r) = r - 1 \quad\leadsto\quad\text{il n'y a qu'une seule racine : } r = 1 $

Par conséquent, toutes ces méthodes sont zéro-stables.
Méthodes de Adams à $q$ pas :

Ces méthodes ont pour polynôme caractéristique :

$ \varrho(r) = r^{p+1} - r^p = r^p(r - 1) \quad\leadsto\quad\text{il y a deux racines : } \begin{cases} r = 1 & \text{a une multiplicité de 1, et $|r|=1$} \\ r = 0 & \text{a une multiplicité de $p$, et $|r|=0<1$} \end{cases} $

Ces méthodes sont donc zéro-stables.
Méthodes de Nystrom et Milne-Simpson à $q$ pas :

Pour ces méthodes, le polynôme caractéristique est :

$ \varrho(r) = r^{p+1} - r^{p-1} = r^{p-1}(r - 1)(r + 1) \quad\leadsto\quad\text{il y a trois racines :} \begin{cases} r = 1 & \text{a une multiplicité de 1, et $|r|=1$} \\ r = 0 & \text{a une multiplicité de $p-1$, et $|r|=0<1$} \\ r = -1 & \text{a une multiplicité de 1, et $|r|=1$} \end{cases} $

Ces méthodes sont également zéro-stables.
Méthodes BDF :

Il est possible de montrer que les méthodes BDF sont zéro-stables si et seulement si $p \leq 5$ (ce qui correspond à $q \leq 6$, comme le montre le calcul effectué avec sympy ci-dessous).

Ordre de convergence maximale d’une méthode MPL

Il est possible de démontrer que la zero stabilité implique une restriction sur l’ordre maximal atteignable pour une méthode multi-pas. On appelle cette contrainte la première barrière de Dahlquist :

⚠ Première barrière de Dahlquist [limitation de l’ordre maximal de convergence]

L’ordre $\omega$ d’une méthode MPL à $q$ pas convergente satisfait les inégalités suivantes :

si $b_{-1}=0$ (la méthode est explicite) alors $\omega\le q$
si $b_{-1}\neq0$ (la méthode est implicite) on doit considérér deux cas
- si $q$ est impair alors $\omega\le q+1$
- si $q$ est pair alors $\omega\le q+2$

Remarque : la première barrière de Dalquist affirme qu’on peut augmenter arbitrairement l’ordre d’un schéma multi-pas en augmentant le nombre de pas. La convergence dépendra ensuite juste du premier polynôme caractéristique (c’est à dire de la zéro-stabilité).

Exercie : vérifier la première barrière de Dahlquist pour les méthodes Euler explicite, Euler implicite, Crank-Nicolson, de Adams-Bashforth et Adams-Moulton avec $q=1,2$ pas et Milne-Simpson avec $q=2$ pas.

Schéma	Type	Ordre maximal prévu par la barrière	Ordre réel $\omega$	Remarque
EE	Explicite	$q=1$	1	-
$\text{AB}_2$	Explicite	$q=2$	2	-
EI	Implicite, $q$ impair	$q+1 = 2$	1	Non optimal : ordre réel $1$, inférieur à l’ordre maximal $2$
CN	Implicite, $q$ impair	$q+1 = 2$	2	-
$\text{AM}_2$	Implicite, $q$ pair	$q+2 = 4$	3	Non optimal : ordre réel $3$, inférieur à l’ordre maximal $4$
$\text{MS}_2$	Implicite, $q$ pair	$q+2 = 4$	4	-

Le script suivant permet de vérifier la première barrière de Dahlquist pour une méthode donnée. Il suffit de rentrer la valeur de $p$ et le script affiche l’ordre maximal de convergence prédit par la barrière de Dahlquist ainsi que les équations à vérifier pour chaque cas.

Code

## j = 0...p ## coeffs du schema
## q = p+1 ## nb de pas
## ω ## ordre de la méthode

p = 3
q = p+1 
ordre_max = q+2 if q%2==0 else q+1

display(Markdown(f"**On considère une méthode à q={q} pas (et donc p={p})**"))
display(Markdown(f"**Elle est d'ordre au plus {ordre_max}. Voici les conditions à satisfaire.**"))

aa = sp.symbols(f'a_0:{q}')
bb = sp.symbols(f'b_0:{q}')
bm1 = sp.Symbol('b_{-1}')

i = 1
sa = sum( [-j*aa[j] for j in range(len(aa))] )
sb = bm1+sum( [bb[j] for j in range(len(bb))] )
display(Markdown(f"- Pour que la méthode soit consistante (c'est-à-dire d'ordre au moins {i}) elle doit vérifier les deux conditions : "))
display(Markdown( r"$\qquad\qquad " + sp.latex(sp.Eq( sum(aa).factor() , 1 ) ) + "$" ))
display(Markdown( r"$\qquad\qquad " + sp.latex(sp.Eq( (sa).factor()+(i*sb).factor() , 1 ) ) + "$" ))

for i in range(2,ordre_max+1):
    sa = sum( [(-j)**i*aa[j] for j in range(len(aa))] )
    sb = bm1+sum( [(-j)**(i-1)*bb[j] for j in range(1,len(bb))] )
    display(Markdown(f"- Pour que la méthode soit d'ordre au moins {i} elle doit vérifier, en plus des précédentes, la condition : "))
    display(Markdown( r"$\qquad\qquad " + sp.latex(sp.Eq( (sa).factor()+(i*sb).factor() , 1 ) ) + "$" ))

On considère une méthode à q=4 pas (et donc p=3)

Elle est d’ordre au plus 6. Voici les conditions à satisfaire.

Pour que la méthode soit consistante (c’est-à-dire d’ordre au moins 1) elle doit vérifier les deux conditions :

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

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

Pour que la méthode soit d’ordre au moins 2 elle doit vérifier, en plus des précédentes, la condition :

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

Pour que la méthode soit d’ordre au moins 3 elle doit vérifier, en plus des précédentes, la condition :

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

Pour que la méthode soit d’ordre au moins 4 elle doit vérifier, en plus des précédentes, la condition :

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

Pour que la méthode soit d’ordre au moins 5 elle doit vérifier, en plus des précédentes, la condition :

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

Pour que la méthode soit d’ordre au moins 6 elle doit vérifier, en plus des précédentes, la condition :

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

Récapitulatif des méthodes multi-pas vues jusqu’ici:

Code

r = sp.symbols('r')

sc = { 
    "EE = AB_1":           {"aa": [1], "bb": [1], "bm1": 0},
    "EI = AM_0 = BDF_1":   {"aa": [1], "bb": [0], "bm1": 1},
    "CN = AM_1":           {"aa": [1], "bb": [sp.Rational(1,2)], "bm1": sp.Rational(1,2)}, 
    "AB_2":                {"aa": [1, 0], "bb": [sp.Rational(3,2), sp.Rational(-1,2)], "bm1": 0}, 
    "AB_3":                {"aa": [1, 0, 0], "bb": [sp.Rational(23,12), sp.Rational(-16,12), sp.Rational(5,12)], "bm1": 0}, 
    "AB_4":                {"aa": [1, 0, 0, 0], "bb": [sp.Rational(55,24), sp.Rational(-59,24), sp.Rational(37,24), sp.Rational(-9,24)], "bm1": 0}, 
    "AB_5":                {"aa": [1, 0, 0, 0, 0], "bb": [sp.Rational(1901,720), sp.Rational(-2774,720), sp.Rational(2616,720), sp.Rational(-1274,720), sp.Rational(251,720)], "bm1": 0}, 
    "AM_2":                {"aa": [1, 0], "bb": [sp.Rational(8,12), sp.Rational(-1,12)], "bm1": sp.Rational(5,12)}, 
    "AM_3":                {"aa": [1, 0, 0], "bb": [sp.Rational(19,24), sp.Rational(-5,24), sp.Rational(1,24)], "bm1": sp.Rational(9,24)}, 
    "AM_4":                {"aa": [1, 0, 0, 0], "bb": [sp.Rational(646,720), sp.Rational(-264,720), sp.Rational(106,720), sp.Rational(-19,720)], "bm1": sp.Rational(251,720)}, 
    "BDF_2":               {"aa": [sp.Rational(4,3), sp.Rational(-1,3)], "bb": [0,0], "bm1": sp.Rational(2,3)}, 
    "BDF_3":               {"aa": [sp.Rational(18,11), sp.Rational(-9,11), sp.Rational(2,11)], "bb": [0,0,0], "bm1": sp.Rational(6,11)}, 
    "BDF_4":               {"aa": [sp.Rational(48,25), sp.Rational(-36,25), sp.Rational(16,25), sp.Rational(-3,25)], "bb": [0,0,0,0], "bm1": sp.Rational(12,25)}, 
    "BDF_5":               {"aa": [sp.Rational(300,137), sp.Rational(-300,137), sp.Rational(200,137), sp.Rational(-75,137), sp.Rational(12,137)], "bb": [0,0,0,0,0], "bm1": sp.Rational(60,137)}, 
    "BDF_6":               {"aa": [sp.Rational(120,49), sp.Rational(-150,49), sp.Rational(400,147), sp.Rational(-75,49), sp.Rational(24,49), sp.Rational(-10,147)], "bb": [0,0,0,0,0,0], "bm1": sp.Rational(20,49)}, 
    "MS_2":                {"aa": [0, 1], "bb": [sp.Rational(4,3), sp.Rational(1,3)], "bm1": sp.Rational(1,3)}
}

for schema, data in sc.items():
    # display(Markdown(f"---"))
    display(Markdown(f"**{schema}**"))
    
    aa = data["aa"]
    bb = data["bb"]
    bm1 = data["bm1"]
    q = len(aa)

    # 1. Type et Pas
    etat = "explicite" if bm1 == 0 else f"implicite"
    display(Markdown(f"- Schéma à $q = {q}$ pas, {etat}."))

    # 2. Calcul de l'Ordre 
    def check_order(p):
        if p == 0: return sum(aa) == 1
        term_a = sum([((-j)**p * aa[j]) for j in range(q)])
        term_b = p * sum([((-j)**(p-1) * bb[j]) for j in range(q)])
        term_bm1 = p * bm1
        return sp.simplify(term_a + term_b + term_bm1) == 1

    ordre = 0
    if check_order(0):
        ordre = 1
        while check_order(ordre + 1):
            ordre += 1
    
    display(Markdown(f"- **Consistance** : {'Oui' if ordre >= 1 else 'Non'}"))
    display(Markdown(f"- **Ordre** : {ordre}"))

    # 3. Polynôme rho et racines
    rho_sp = r**q - sum(aa[j] * r**(q-1-j) for j in range(q))
    display(Markdown(f"- $\\varrho(r) = {sp.latex(sp.expand(rho_sp))}$"))
    
    # Calcul des modules des racines
    racines = sp.solve(rho_sp, r)
    # On définit une fonction qui calcule le module manuellement 
    # pour éviter les problèmes d'évaluation de sp.Abs sur des 
    # expressions symboliques complexes
    def get_abs_expr(root):
        re_part = sp.re(root)
        im_part = sp.im(root)
        return sp.sqrt(re_part**2 + im_part**2)
    # moduli = [sp.Abs(rac).evalf(4) for rac in racines]
    moduli = [get_abs_expr(rac).evalf(4) for rac in racines]
    display(Markdown(f"- Module de ses racines ${{|r_j|}}$ : ${sp.latex(moduli)}$"))

EE = AB_1

Schéma à $q = 1$ pas, explicite.

Consistance : Oui

Ordre : 1

$\varrho(r) = r - 1$

Module de ses racines ${|r_j|}$ : $\left[ 1.0\right]$

EI = AM_0 = BDF_1

Schéma à $q = 1$ pas, implicite.

Consistance : Oui

Ordre : 1

$\varrho(r) = r - 1$

Module de ses racines ${|r_j|}$ : $\left[ 1.0\right]$

CN = AM_1

Schéma à $q = 1$ pas, implicite.

Consistance : Oui

Ordre : 2

$\varrho(r) = r - 1$

Module de ses racines ${|r_j|}$ : $\left[ 1.0\right]$

AB_2

Schéma à $q = 2$ pas, explicite.

Consistance : Oui

Ordre : 2

$\varrho(r) = r^{2} - r$

Module de ses racines ${|r_j|}$ : $\left[ 0, \ 1.0\right]$

AB_3

Schéma à $q = 3$ pas, explicite.

Consistance : Oui

Ordre : 3

$\varrho(r) = r^{3} - r^{2}$

Module de ses racines ${|r_j|}$ : $\left[ 0, \ 1.0\right]$

AB_4

Schéma à $q = 4$ pas, explicite.

Consistance : Oui

Ordre : 4

$\varrho(r) = r^{4} - r^{3}$

Module de ses racines ${|r_j|}$ : $\left[ 0, \ 1.0\right]$

AB_5

Schéma à $q = 5$ pas, explicite.

Consistance : Oui

Ordre : 5

$\varrho(r) = r^{5} - r^{4}$

Module de ses racines ${|r_j|}$ : $\left[ 0, \ 1.0\right]$

AM_2

Schéma à $q = 2$ pas, implicite.

Consistance : Oui

Ordre : 3

$\varrho(r) = r^{2} - r$

Module de ses racines ${|r_j|}$ : $\left[ 0, \ 1.0\right]$

AM_3

Schéma à $q = 3$ pas, implicite.

Consistance : Oui

Ordre : 4

$\varrho(r) = r^{3} - r^{2}$

Module de ses racines ${|r_j|}$ : $\left[ 0, \ 1.0\right]$

AM_4

Schéma à $q = 4$ pas, implicite.

Consistance : Oui

Ordre : 5

$\varrho(r) = r^{4} - r^{3}$

Module de ses racines ${|r_j|}$ : $\left[ 0, \ 1.0\right]$

BDF_2

Schéma à $q = 2$ pas, implicite.

Consistance : Oui

Ordre : 2

$\varrho(r) = r^{2} - \frac{4 r}{3} + \frac{1}{3}$

Module de ses racines ${|r_j|}$ : $\left[ 0.3333, \ 1.0\right]$

BDF_3

Schéma à $q = 3$ pas, implicite.

Consistance : Oui

Ordre : 3

$\varrho(r) = r^{3} - \frac{18 r^{2}}{11} + \frac{9 r}{11} - \frac{2}{11}$

Module de ses racines ${|r_j|}$ : $\left[ 1.0, \ 0.4264, \ 0.4264\right]$

BDF_4

Schéma à $q = 4$ pas, implicite.

Consistance : Oui

Ordre : 4

$\varrho(r) = r^{4} - \frac{48 r^{3}}{25} + \frac{36 r^{2}}{25} - \frac{16 r}{25} + \frac{3}{25}$

Module de ses racines ${|r_j|}$ : $\left[ 1.0, \ 0.5609, \ 0.5609, \ 0.3815\right]$

BDF_5

Schéma à $q = 5$ pas, implicite.

Consistance : Oui

Ordre : 5

$\varrho(r) = r^{5} - \frac{300 r^{4}}{137} + \frac{300 r^{3}}{137} - \frac{200 r^{2}}{137} + \frac{75 r}{137} - \frac{12}{137}$

Module de ses racines ${|r_j|}$ : $\left[ 1.0, \ 0.7087 - 5.0 \cdot 10^{-26} i, \ 0.4176 + 1.0 \cdot 10^{-25} i, \ 0.7087 - 5.0 \cdot 10^{-26} i, \ 0.4176 + 1.0 \cdot 10^{-25} i\right]$

BDF_6

Schéma à $q = 6$ pas, implicite.

Consistance : Oui

Ordre : 6

$\varrho(r) = r^{6} - \frac{120 r^{5}}{49} + \frac{150 r^{4}}{49} - \frac{400 r^{3}}{147} + \frac{75 r^{2}}{49} - \frac{24 r}{49} + \frac{10}{147}$

Module de ses racines ${|r_j|}$ : $\left[ 1.0, \ 0.4061, \ 0.474, \ 0.474, \ 0.8634, \ 0.8634\right]$

MS_2

Schéma à $q = 2$ pas, implicite.

Consistance : Oui

Ordre : 4

$\varrho(r) = r^{2} - 1$

Module de ses racines ${|r_j|}$ : $\left[ 1.0, \ 1.0\right]$

A-stabilité

Rappels

Problème de Dalquist : $y'(t)=\lambda y(t)$, $y(0)=1$ avec $\lambda\in\mathbb{C}$ tel que $\Re(\lambda)<0$. Sa solution exacte est donnée par $y(t)=e^{\lambda t}$, et on a $\lim_{t\to+\infty}|y(t)|=0$.
Schéma A-stable : pour $h$ fixé, un schéma numérique est A-stable si, pour tout $\lambda$ vérifiant $\Re(\lambda) < 0$, la suite numérique $u_n$ satisfait $\lim_{n\to+\infty} |u_n| =0$. Si 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.

Lorsqu’on applique une méthode MPL au problème de Dahlquist, la suite numérique $u_n$ vérifie une relation de récurrence linéaire du type $ \left(1- (h \lambda) b_{-1}\right)u_{n+1} = \sum_{j=0}^p (a_j + (h \lambda) b_j) u_{n-j} $ pour $n\ge p$.

En notant $z$ le produit $h \lambda$, cette relation peut s’écrire sous la forme : $ (1 - z b_{-1}) u_{n+1} = \sum_{j=0}^{p} (a_j + z b_j)\, u_{n-j}. $

Si $q=1$ on a un schéma à un pas et l’équation récurrente s’écrit $ u_{n+1} = q(z) u_n, \qquad q(z) = \frac{a_0 + z b_0}{1 - z b_{-1}}. $ Ainsi, si $|q(z)| < 1$ alors $\lim_{n\to+\infty} |u_n| = 0$.
Si $q>1$ on a un schéma à plusieurs pas. Pour étudier la stabilité de cette équation récurrente on verra ci-dessous qu’il faut étudier les racines du polynôme de stabilité (ou plutôt de la famille de polynômes paramétrée par $z$), défini par $ p_z(r) = (1 - z b_{-1}) r^{p+1} - \sum_{j=0}^p (a_j + z b_j) r^{p-j} = \varrho(r) - z \sigma(r). $

Proposition (condition stricte des racines) :

Une méthode MPL est absolument stable si et seulement si son polynôme de A-stabilité satisfait (éventuellement sous certaines conditions sur $z$) la condition stricte des racines. Autrement dit, les racines du polynôme

$ p_z(r) = \varrho(r) - z \sigma(r) $

doivent être strictement à l’intérieur du disque unité dans $\mathbb{C}$.

CAS RÉEL : lorsqu’on étudie la A-stabilité seulement dans $\mathbb{R}$, on peut se contenter d’étudier les racines du polynôme de stabilité pour $z=-x$ avec $x\in\mathbb{R}^+$. Cela revient à étudier les fonctions $x\mapsto r_j(x)$ racines du polynôme de stabilité. L’intervalle de A-stabilité est alors l’ensemble des valeurs de $x$ pour lesquelles $|r_j(x)|<1$ pour tout $j=0,\dots,p$.

NB : contrairement à la zéro-stabilité (où $|r|=1$ est toléré pour les racines simples), la A-stabilité exige une contraction stricte pour compenser l’amortissement de la solution exacte.

⚠ Deuxième barrière de Dahlquist : limitation sur $h$

Aucune méthode multi-pas explicite n’est inconditionnellement A-stable.
Aucune méthode multi-pas implicite d’ordre supérieur à 2 n’est inconditionnellement A-stable dans $\mathbb{C}$.

En d’autres termes, le 🎅 n’existe pas :

Seules les méthodes multi-pas implicites d’ordre 1 et 2 peuvent être inconditionnellement A-stables dans $\mathbb{C}$.
Méthodes multi-pas explicites ou Méthodes multi-pas implicites d’ordre supérieur à 2 $\leadsto$ toujours conditionnellement A-stables, c’est-à-dire qu’il existe une valeur maximale de $h$ pour laquelle elles restent A-stables dans $\mathbb{C}$.

🔎 Conséquences : augmenter le nombre de pas $q\dots$

améliore l’ordre de consistance
mais impose une restriction sur $h$ pour garantir l’A-stabilité
peut affecter la convergence : la zéro-stabilité devient délicate pour assurer la convergence de la méthode.

Exemples : analyse de l’A-stabilité des schémas d’Adams

Les schémas $\text{AB}_q$ (Adams-Bashforth) sont toujours conditionnellement A-stables car ils sont explicites.
- Par exemple, pour le schéma $\text{AB}_1$ (Euler explicite), on a :
  $ p_z(r) = (r - 1) - z $ La seule racine est $r = 1 + z$. Pour assurer la stabilité, il faut que $\lvert 1 + z \rvert < 1$, ce qui définit un disque ouvert de centre $(-1,0)$ et de rayon $1$ dans le plan complexe. Dans le cas réel, où $\varphi(t,y) = -\beta y$ avec $\beta = -\Re(\lambda)$, la condition se réécrit $0 < \beta h < 1$, ce qui correspond au résultat vu au chapitre précédent.
Les schémas $\text{AM}_q$ (Adams-Moulton) avec $q \le 1$ sont inconditionnellement A-stables, car ils sont implicites et d’ordre $\le2$.
- $\text{AM}_0$ : Euler implicite
  $ p_z(r) = (r - 1) - z r $ La seule racine est $r = \frac{1}{1 + z}$. La condition de stabilité est donc :
  $ \left\lvert \frac{1}{1 + z} \right\rvert < 1 $ Cette condition définit l’extérieur d’un disque de centre $(1,0)$ et de rayon $1$.
  Comme le demi-plan $\Re(z) < 0$ est inclus dans cette région, la méthode est inconditionnellement A-stable.
- $\text{AM}_1$ : Crank-Nicolson
  $ p_z(r) = (r - 1) - z \frac{1}{2} (r + 1) $ Sa seule racine est :
  $ r = \frac{2 + z}{2 - z} $ Pour la stabilité, on impose :
  $ \left\lvert \frac{2 + z}{2 - z} \right\rvert < 1 $ Ce qui revient à exiger $\Re(z) < 0$.
  Comme $\mathbb{C}^-$ est entièrement inclus dans cette région, la méthode est inconditionnellement A-stable.
Les schémas $\text{AM}_q$ (Adams-Moulton) avec $q > 1$ sont conditionnellement A-stables, car ils sont implicites et d’ordre $>2$.

Les images affichées ci-dessous montrent les régions de A-stabilité des quelques méthodes vues en précédence et sont issues du livre
A. Quarteroni, F. Saleri, P. Gervasio Calcul Scientifique: Cours, exercices corrigés et illustrations en Matlab et Octave.

Régions de A-stabilité des méthodes d’Adams-Bashforth et d’Adams-Moulton

Les figures ci-dessous illustrent les régions de stabilité absolue des méthodes :

À gauche : Adams-Bashforth (AB), qui sont des schémas explicites.
À droite : Adams-Moulton (AM), qui sont des schémas implicites.

On observe que :

Ces régions sont bornées et diminuent à mesure que l’ordre augmente.
De plus, toutes les méthodes explicites et les méthodes implicites d’ordre supérieur à 2 ne sont pas inconditionnellement A-stables.
Elles ne conviennent donc pas aux problèmes nécessitant une intégration sur de très longs intervalles de temps, car elles imposent une contrainte sur la taille du pas de temps $h$.

Régions de A-stabilité des méthodes BDF

L’illustration ci-dessous montre les régions de stabilité absolue des méthodes Backward Differentiation Formula (BDF) pour $q \leq 4$ qui sont des schémas implicites.

Contrairement aux méthodes d’Adams :

Les régions de stabilité des BDF sont conditionnellement A-stable sur $\mathbb{C}$ mais incluent toujours l’axe réel négatif, ce qui leur confère une inconditionnelle A-stabilité sur $\mathbb{R}$.
Ces méthodes sont donc particulièrement adaptées aux problèmes nécessitant une intégration sur des intervalles de temps très longs, voire infinis.
Elles sont aussi privilégiées pour la résolution des problèmes raides, car, avec un changement de variable, cela revient à résoudre un problème de Cauchy sur un intervalle de temps grand.

Exercice de synthèse : consistance, zéro-stabilité et A-stabilité des schémas EE, EI, CN

Pour les méthodes d’Euler explicite, Euler implicite et Cranck-Nicolson :

Trouver les coefficients $p\ge0$, $\{a_j\}_{j=0}^p$ et $\{b_j\}_{j=-1}^p$ tels que $ u_{n+1} = \sum_{j=0}^p (a_ju_{n-j}+h b_j\varphi_{n-j})+hb_{-1}\varphi_{n+1}, \qquad n=p,p+1,\dots,N-1$
Calculer les $p+1$ racines $\{r_j\}_{j=0}^p$ du premier polynôme caractéristique $\varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}.$ En déduire que la méthode est zéro-stable, i.e. que $ \begin{cases} |r_j| \le 1 \quad\text{pour tout } j=0,\dots,p \\ r_j \text{ a multiplicité} >1 \quad\implies \quad |r_j|<1.\end{cases} $
Donner l’ordre de consistance maximal théorique puis étudier l’ordre de consistance réel, autrement dit trouver le plus grand $\omega\ge1$ tel que $\xi(0)=\xi(1)=\dots=\xi(\omega)=1$.
Conclure sur la convergence de la méthode et son ordre de convergence.
Étudier la A-stabilité

Il s’agit de trois méthodes à $q=1$ pas, donc $p=0$.

Commençons par l’étude de la convergence en étudiant la consistance (et son ordre) et la zéro-stabilité.

Selon la première barrière de Dahlquist, si la méthode est explicite l’ordre de convergence maximal est 2 (car $q$ impair), si la méthode est explicite l’ordre maximal est 1. Pour pouvoir étudier l’ordre de consistance, commençons par afficher les $\xi(i)$ pour $i=0,1,2$ :

Code

xi_cond(p=0, omega=2)

Pour $p=0$, on a une méthode à $q=1$ pas et on affiche $ξ(i)$ pour $i=0,...,2$

$\xi(0) = a_{0}$

$\xi(1) = b_{0} + b_{-1}$

$\xi(2) = 2 b_{-1}$

Résumons donc les quantités indiquées dans le tableau suivant :

\( \def\arraystretch{2.8} \begin{array}{|c|c|c|c|} \hline & \textbf{EE} & \textbf{EI} & \textbf{CN} \\ \hline \textbf{Méthode : } u_{n+1} = & u_n + h \varphi_n & u_n + h \varphi_{n+1} & u_n + \frac{h}{2} \big(\varphi_{n+1} + \varphi_n\big) \\ \hline \textbf{Coefficients : } p = 0 \text{ et } & a_0 = 1, b_{-1} = 0, b_0 = 1 & a_0 = 1, b_{-1} = 1, b_0 = 0 & a_0 = 1, b_{-1} = \frac{1}{2}, b_0 = \frac{1}{2} \\ \hline \textbf{Polynôme caractéristique : } \varrho & \varrho(r) = r - 1 & \varrho(r) = r - 1 & \varrho(r) = r - 1 \\ \hline \textbf{Racines : } & r_0 = 1 & r_0 = 1 & r_0 = 1 \\ \hline \textbf{Zéro-stabilité : } & \text{oui, } \lvert r_0 \rvert = 1 \text{ et mult. 1} & \text{oui, } \lvert r_0 \rvert = 1 \text{ et mult. 1} & \text{oui, } \lvert r_0 \rvert = 1 \text{ et mult. 1} \\ \hline \textbf{Consistance : } & \text{oui, } \xi(0) = 1, \xi(1) = 1 & \text{oui, } \xi(0) = 1, \xi(1) = 1 & \text{oui, } \xi(0) = 1, \xi(1) = 1 \\ \hline \textbf{Ordre maximal (Dahlquist) : } & \omega \le 1 & \omega \le 2 & \omega \le 2 \\ \hline \textbf{Ordre réel : } \omega & \omega = 1 & \omega = 1 \text{ car } \xi(2) = 2 \neq 1 & \omega = 2 \text{ car } \xi(2) = 1 \\ \hline \end{array} \)

Passons maintenant à l’étude de la A-stabilité. Pour cela, considérons le toy-model : soit $\lambda\in\mathbb{C}$ un nombre complexe tel que $\Re(\lambda)=-\beta<0$ et considérons le problème de Cauchy

$\begin{cases} y'(t)=\lambda y(t), &\text{pour }t>0,\\ y(0)=1. \end{cases}$

Posons $z = \lambda h$ et calculons les racines du polynôme caractéristique $p_z(r) = \varrho(r) - z \sigma(r)$ :

$ \def\arraystretch{2.8} \begin{array}{|c|c|c|c|} \hline & \textbf{EE} & \textbf{EI} & \textbf{CN} \\ \hline \textbf{Schéma appliqué au toy-model} & u_n = (1 + z)^n & u_n = \left(\frac{1}{1 + z}\right)^n & u_n = \left(\frac{2 + z}{2 - z}\right)^n \\ \hline \textbf{Condition pour } \lvert u_n \rvert \xrightarrow[n\to\infty]{}0 & \lvert 1 + z \rvert < 1 & \left\lvert \frac{1}{1 + z} \right\rvert < 1 & \left\lvert \frac{2 + z}{2 - z} \right\rvert < 1 \\ \hline \textbf{A-stabilité ?} & \textbf{Uniquement} \text{ si } -2 < |z| < 0 & \textbf{Oui} \text{ (car } \Re(\lambda) < 0\text{)} & \textbf{Oui} \text{ (car } \Re(\lambda) < 0\text{)} \\ \hline \end{array} $

Conclusion :

Euler explicite : la méthode est consistante d’ordre $\omega = 1$ qui est l’ordre maximal prévu par la barrière de Dahlquist (schéma explicite). Le schéma étant aussi zéro-stable, il est donc convergent d’ordre 1. Concerant la A-stabilité, il est nécessaire que $-2 < |z| < 0$ pour que la méthode soit A-stable. Dans $\mathbb{R}$, cela implique $0 < \beta h < 2$.
Euler implicite : la méthode est consistante d’ordre $\omega = 1$ qui est inférieur à l’ordre maximal prévu par la barrière de Dahlquist (schéma implicite). Le schéma étant aussi zéro-stable, il est donc convergent d’ordre 1. Concernant la A-stabilité, la méthode est inconditionellement A-stable.
Cranck-Nicolson : la méthode est consistante d’ordre $\omega = 2$ qui coïncide avec l’ordre maximal prévu par la barrière de Dahlquist (schéma implicite, nombre de pas impair). Le schéma étant aussi zéro-stable, il est donc convergent d’ordre 2. Concernant la A-stabilité, la méthode est inconditionellement A-stable.

Propriétés des schémas Prédicteur-Correcteur basés sur les méthodes multi-pas

🔴 Proposition: ordre de convergence d’un schéma PC

Si on note

$\omega_{[P]}$ l’ordre du predictor (schéma explicite)
$\omega_{[C]}$ l’ordre du corrector (schéma implicite)

alors l’ordre du schéma predictor-corrector est donné par

$ \omega_{[PC]}=\min\Big\{\omega_{[C]},\omega_{[P]}+1\Big\} $

Ainsi, pour obtenir une méthode d’ordre $\omega$, on peut choisir un prédicteur d’ordre $\omega-1$ et un correcteur d’ordre $\omega$ (par exemple, deux méthodes d’Adam ayant le même nombre de pas).

Zéro-stabilité : on peut démontrer que la zéro-stabilité du prédicteur n’influe pas sur la zéro-stabilité de la méthode prédicteur-correcteur. Il est même possible d’utiliser un prédicteur instable.

A-stabilité : les méthodes predictor-corrector sont soumises à une condition de A-stabilité, qui est typiquement celle du prédicteur. Elles ne sont donc pas adaptées à la résolution de problèmes de Cauchy sur des intervalles non bornés ou de problèmes stiff.

Méthodes cycliques composites : comment surmonter les limites de Dahlquist

On peut surmonter les limites imposées par les barrières de Dahlquist en combinant plusieurs méthodes multi-pas. Par exemple, les deux méthodes suivantes à trois pas

$ \begin{aligned} u_{n+1} &= -\frac{8}{11}u_{n} + \frac{19}{11}u_{n-1} + \frac{h}{30} \big( 30\varphi_{n+1}+57\varphi_n+24\varphi_{n-1}-\varphi_{n-2} \big) \\ u_{n+1} &= \frac{449}{240}u_{n} - \frac{19}{30}u_{n-1} -\frac{361}{240}u_{n-2} + \frac{h}{720}\big(251\varphi_{n+1}+456\varphi_n-1347\varphi_{n-1}-350\varphi_{n-2}\big) \end{aligned} $

sont d’ordre cinq, mais sont instables. Cependant, en les combinant comme suit

$ u_{n+1} = \begin{cases} -\frac{8}{11}u_{n} + \frac{19}{11}u_{n-1} + \frac{h}{30} \big( 30\varphi_{n+1}+57\varphi_n+24\varphi_{n-1}-\varphi_{n-2} \big) &\text{si $n$ est pair}\\ \frac{449}{240}u_{n} - \frac{19}{30}u_{n-1} -\frac{361}{240}u_{n-2} + \frac{h}{720}\big(251\varphi_{n+1}+456\varphi_n-1347\varphi_{n-1}-350\varphi_{n-2}\big) &\text{sinon.} \end{cases} $

on obtient une méthode A-stable à trois pas d’ordre cinq.

Annexe : comment tracer les régions de A-stabilité ?

Soit $z = h\lambda \in \mathbb{C}^-$. La région de stabilité absolue $\mathcal{A}$ est l’ensemble des points $z$ pour lesquels toutes les racines $r_j$ du polynôme caractéristique $p_z(r) = \varrho(r) - z\,\sigma(r)$ vérifient $|r_j| < 1$. L’intervalle de $\mathcal{A}$-stabilité réelle est l’intersection $\mathcal{A} \cap \mathbb{R}^-$.

Il est en général difficile de déterminer directement $\mathcal{A}$ car il faut tester, pour chaque $z$, si le critère $|r_j| < 1$ est satisfait pour toutes les racines. Il est plus aisé de déterminer la frontière $\partial\mathcal{A}$. Sur cette frontière, au moins une racine a un module exactement égal à $1$ : on peut l’écrire $r = e^{i\theta}$, $\theta \in [0, 2\pi]$. En substituant dans $p_z(r) = 0$, on obtient : $ \varrho(e^{i\theta}) - z\,\sigma(e^{i\theta}) = 0 \quad\leadsto\quad z(\theta) = \frac{\varrho(e^{i\theta})}{\sigma(e^{i\theta})}. $ Cette relation définit $\partial\mathcal{A}$ comme une courbe paramétrée dans le plan complexe. Pour en obtenir une approximation numérique, il suffit d’évaluer ce quotient sur une discrétisation de $[0, 2\pi]$ (en vérifiant au préalable que $\sigma(e^{i\theta}) \neq 0$). La courbe ainsi tracée localise les points où une racine est de module $1$, mais ne garantit pas que les autres racines restent dans le disque unité. C’est pourquoi la frontière est complétée par un remplissage coloré (test sur grille) afin d’identifier la véritable zone de stabilité.

Code

# ── Helpers symboliques ──────────────────────────────────────────────────────

def _build_rho_sigma(schema):
    """Construit ρ et σ symboliquement (SymPy)."""
    r = sp.Symbol('r')
    aa, bb, bm1 = schema["aa"], schema["bb"], schema["bm1"]
    p = len(aa) - 1                                      
    rho   = r**(p+1) - sum(aa[j] * r**(p-j) for j in range(p+1))
    sigma = bm1 * r**(p+1) + sum(bb[j] * r**(p-j) for j in range(p+1))
    return rho, sigma


def polyCar(schema):
    """Affiche ρ et σ (méthode explicite : bm1 = 0)."""
    rho, sigma = _build_rho_sigma(schema)
    display(Math(r'\varrho(r) = ' + sp.latex(rho)))
    display(Math(r'\sigma(r)  = ' + sp.latex(sigma)))


def polyCarIMPL(schema):
    """
    Affiche ρ, σ et p_z(r) = ρ(r) − z·σ(r) avec z symbolique.
    Les termes sont regroupés en puissances de r pour plus de lisibilité.
    """
    r, z = sp.symbols('r z')
    rho, sigma = _build_rho_sigma(schema)
    pz = sp.collect(sp.expand(rho - z * sigma), r)   # ← collecte en r
    display(Math(r'\varrho(r) = '                          + sp.latex(rho)))
    display(Math(r'\sigma(r)  = '                          + sp.latex(sigma)))
    display(Math(r'p_z(r) = \varrho(r)-z\,\sigma(r) = '   + sp.latex(pz)))


# ── Fonction de stabilité ────────────────────────────────────────────────────

def stabilityFunction(z, schema):
    """
    Retourne max|r_j| pour le polynôme caractéristique en z.
    Valeur > 1  ⟺  z est hors de la région de stabilité absolue.
    Ne modifie pas schema.
    """
    aa  = np.array(schema["aa"],  dtype=float)
    bb  = np.array(schema["bb"],  dtype=float)
    bm1 = float(schema["bm1"])
    # Coefficients de p_z : [1−z·bm1, −(aa+z·bb)]  (degré décroissant)
    coeffs = np.concatenate([[1 - z * bm1], -(aa + z * bb)])
    return float(np.max(np.abs(np.roots(coeffs))))


# ── Analyse de stabilité absolue ─────────────────────────────────────────────

def Astabilite(schema):
    aa  = np.array(schema["aa"],  dtype=float)
    bb  = np.array(schema["bb"],  dtype=float)
    bm1 = float(schema["bm1"])
    p   = len(aa) - 1

    # ── Frontière paramétrée z(θ) = ρ(e^iθ) / σ(e^iθ) ──────────────────────
    theta   = np.linspace(0, 2 * np.pi, 201)
    c       = np.exp(1j * theta)
    rho_c   = c**(p+1) - np.polyval(aa, c)
    sigma_c = bm1 * c**(p+1) + np.polyval(bb, c)
    mask    = np.abs(sigma_c) > 1e-14
    zz      = np.where(mask, rho_c / sigma_c, np.nan + 1j * np.nan)
    zz_fin  = zz[np.isfinite(zz)]

    # ── Limites de la grille ─────────────────────────────────────────────────
    x_min = min(zz_fin.real.min(), -0.5)
    # ↓ On inclut Re(z) > 0 pour capturer les frontières qui y résident
    #   (méthodes A-stables : ∂A est entièrement dans le demi-plan droit)
    x_max = max(zz_fin.real.max() * 1.1, 0.5)
    y_ext = max(np.abs(zz_fin.imag).max(), 0.5) * 1.1
    Lx, Ly = (x_min, x_max), (-y_ext, y_ext)

    x = np.linspace(Lx[0], Lx[1], 301)
    y = np.linspace(Ly[0], Ly[1], 301)
    X, Y = np.meshgrid(x, y)
    _stab = np.vectorize(lambda xv, yv: stabilityFunction(xv + 1j * yv, schema))
    Z = _stab(X, Y)

    # ── Tracé ────────────────────────────────────────────────────────────────
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4.5))

    # --- Frontière paramétrée ---
    ax1.plot(zz.real, zz.imag, 'r', lw=1.5)
    ax1.axvline(0, color='k', lw=0.8)
    ax1.axhline(0, color='k', lw=0.8)
    ax1.set_aspect('equal')
    ax1.grid(True, alpha=0.4)
    ax1.set_xlabel(r'$\Re(h\lambda)$')
    ax1.set_ylabel(r'$\Im(h\lambda)$')
    ax1.set_title(f"Frontière $\\partial\\mathcal{{A}}$ — {schema['Nom']}")

    # --- Zone de stabilité remplie ---
    ax2.contourf(X, Y, Z, levels=[0, 1], colors=[[1, 0.5, 0.8]])

    # ↓ Correction : on ne trace le contour que si le niveau 1 est dans la plage
    if Z.min() <= 1.0 <= Z.max():
        ax2.contour(X, Y, Z, levels=[1], colors='k', linewidths=1.2)
    else:
        # La frontière est hors de la vue → méthode A-stable sur cette région
        ax2.text(0.97, 0.97, r"$\mathcal{A}$-stable : $\partial\mathcal{A} \notin \mathbb{C}^-$",
                 transform=ax2.transAxes, ha='right', va='top',
                 fontsize=8, color='gray',
                 bbox=dict(boxstyle='round,pad=0.3', fc='white', alpha=0.7))

    ax2.axvline(0, color='k', lw=2)
    ax2.grid(True, alpha=0.4)
    ax2.set_xlabel(r'$\Re(h\lambda)$')
    ax2.set_ylabel(r'$\Im(h\lambda)$')
    ax2.set_title(fr"Zone de stabilité absolue — {schema['Nom']}")

    # --- Intervalle de A-stabilité sur ℝ⁻ ---
    # Recherche sur un domaine large, indépendant de la grille graphique
    X_SEARCH = 7.0
    x_real = np.linspace(-X_SEARCH, 0, 401)
    stab_r = np.array([stabilityFunction(xv + 0j, schema) for xv in x_real])
    stable = x_real[stab_r <= 1.0 + 1e-10]

    if stable.size:
        x_lo = stable.min()
        # Non borné ssi le segment atteint la borne gauche de la recherche
        unbounded = x_lo <= -X_SEARCH + (X_SEARCH / len(x_real))
        label = (r"$\mathcal{A}\cap\mathbb{R}^-=(-\infty,\,0]$"
                 if unbounded
                 else fr"$\mathcal{{A}}\cap\mathbb{{R}}^-\approx[{x_lo:.3f},\,0]$")
        # On ne trace que la partie visible sur le graphique
        visible = stable[stable >= Lx[0]]
        if visible.size:
            ax2.plot(visible, np.zeros_like(visible),
                     color='darkred', lw=9, solid_capstyle='butt',
                     alpha=0.85, label=label)
        ax2.legend(fontsize=8, loc='lower right')

    plt.tight_layout()
    plt.show()

Code

EE = { "Nom":"EE=AB1", "aa":[1], "bb":[1], "bm1":0  }
polyCar(EE)
Astabilite(EE)

$\displaystyle \varrho(r) = r - 1$

$\displaystyle \sigma(r) = 1$

Code

AB2 = { "Nom":"AB2", "aa":[1,0], "bb":[sp.Rational(3,2),-sp.Rational(1,2)], "bm1":0  }
polyCar(AB2)
Astabilite(AB2)

$\displaystyle \varrho(r) = r^{2} - r$

$\displaystyle \sigma(r) = \frac{3 r}{2} - \frac{1}{2}$

Code

AB3 = { "Nom":"AB3", "aa":[1,0,0], "bb":[sp.Rational(23,12),-sp.Rational(16,12), sp.Rational(5,12)], "bm1":0  }
polyCar(AB3)
Astabilite(AB3)

$\displaystyle \varrho(r) = r^{3} - r^{2}$

$\displaystyle \sigma(r) = \frac{23 r^{2}}{12} - \frac{4 r}{3} + \frac{5}{12}$

Code

AB4 = { "Nom":"AB4", "aa":[1,0,0,0], "bb":[sp.Rational(55,24),-sp.Rational(59,24), sp.Rational(37,24), -sp.Rational(9,24)], "bm1":0  }
polyCar(AB4)
Astabilite(AB4)

$\displaystyle \varrho(r) = r^{4} - r^{3}$

$\displaystyle \sigma(r) = \frac{55 r^{3}}{24} - \frac{59 r^{2}}{24} + \frac{37 r}{24} - \frac{3}{8}$

Code

AB5 = { "Nom":"AB5", "aa":[1,0,0,0,0], "bb":[sp.Rational(1901,720), -sp.Rational(1387,360), sp.Rational(109,30), -sp.Rational(637,360), sp.Rational(251,720)], "bm1":0  }
polyCar(AB5)
Astabilite(AB5)

$\displaystyle \varrho(r) = r^{5} - r^{4}$

$\displaystyle \sigma(r) = \frac{1901 r^{4}}{720} - \frac{1387 r^{3}}{360} + \frac{109 r^{2}}{30} - \frac{637 r}{360} + \frac{251}{720}$

Code

AM0 = { "Nom":"AM0=EI", "aa":[1], "bb":[0], "bm1": 1  }
polyCarIMPL(AM0)
Astabilite(AM0)

$\displaystyle \varrho(r) = r - 1$

$\displaystyle \sigma(r) = r$

$\displaystyle p_z(r) = \varrho(r)-z\,\sigma(r) = r \left(1 - z\right) - 1$

Code

AM1 = { "Nom":"AM1=CN",  "aa":[1], "bb":[sp.Rational(1,2)], "bm1": sp.Rational(1,2)  }
polyCarIMPL(AM1)
Astabilite(AM1)

$\displaystyle \varrho(r) = r - 1$

$\displaystyle \sigma(r) = \frac{r}{2} + \frac{1}{2}$

$\displaystyle p_z(r) = \varrho(r)-z\,\sigma(r) = r \left(1 - \frac{z}{2}\right) - \frac{z}{2} - 1$

Code

AM2 = { "Nom":"AM2", "aa":[1,0], "bb":[sp.Rational(8,12),-sp.Rational(1,12)], "bm1":sp.Rational(5,12) }
polyCar(AM2)
Astabilite(AM2)

$\displaystyle \varrho(r) = r^{2} - r$

$\displaystyle \sigma(r) = \frac{5 r^{2}}{12} + \frac{2 r}{3} - \frac{1}{12}$

Code

AM3 = { "Nom":"AM3", "aa":[1,0,0], "bb":[sp.Rational(19,24),-sp.Rational(5,24),sp.Rational(1,24)], "bm1":sp.Rational(9,24)  }
polyCar(AM3)
Astabilite(AM3)

$\displaystyle \varrho(r) = r^{3} - r^{2}$

$\displaystyle \sigma(r) = \frac{3 r^{3}}{8} + \frac{19 r^{2}}{24} - \frac{5 r}{24} + \frac{1}{24}$

Code

AM4 = { "Nom":"AM4", "aa":[1,0,0,0], "bb":[sp.Rational(646,720), sp.Rational(-264,720), sp.Rational(106,720), sp.Rational(-19,720)], "bm1":sp.Rational(251,720)  }
polyCar(AM4)
Astabilite(AM4)

$\displaystyle \varrho(r) = r^{4} - r^{3}$

$\displaystyle \sigma(r) = \frac{251 r^{4}}{720} + \frac{323 r^{3}}{360} - \frac{11 r^{2}}{30} + \frac{53 r}{360} - \frac{19}{720}$

Code

MS2 = { "Nom":"MS2", "aa":[0,1], "bb":[sp.Rational(4,3),sp.Rational(1,3)], "bm1":sp.Rational(1,3) }
polyCar(MS2)
##Astabilite(MS2)
# Jamais A-stable (voir le TD sur les méthodes de multi-step)

$\displaystyle \varrho(r) = r^{2} - 1$

$\displaystyle \sigma(r) = \frac{r^{2}}{3} + \frac{4 r}{3} + \frac{1}{3}$

Code

BDF1 = { "Nom":"BDF1=EI", "aa":[1], "bb":[0], "bm1":1  }
polyCarIMPL(BDF1)
Astabilite(BDF1)

$\displaystyle \varrho(r) = r - 1$

$\displaystyle \sigma(r) = r$

$\displaystyle p_z(r) = \varrho(r)-z\,\sigma(r) = r \left(1 - z\right) - 1$

Code

BDF2 = { "Nom":"BDF2", "aa":[sp.Rational(4,3),-sp.Rational(1,3)], 
                       "bb":[0,0], "bm1":sp.Rational(2,3) }
polyCarIMPL(BDF2)
Astabilite(BDF2)

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

$\displaystyle \sigma(r) = \frac{2 r^{2}}{3}$

$\displaystyle p_z(r) = \varrho(r)-z\,\sigma(r) = r^{2} \cdot \left(1 - \frac{2 z}{3}\right) - \frac{4 r}{3} + \frac{1}{3}$

Code

BDF3 = { "Nom":"BDF3", "aa":[sp.Rational(18,11),-sp.Rational(9,11),sp.Rational(2,11)], 
                       "bb":[0,0,0], "bm1":sp.Rational(6,11)  }
##polyCarIMPL(BDF3)
Astabilite(BDF3)

Code

BDF4 = { "Nom":"BDF4", "aa":[sp.Rational(48,25),-sp.Rational(36,25),sp.Rational(16,25),-sp.Rational(3,25)], 
                       "bb":[0,0,0,0], "bm1":sp.Rational(12,25)  }
## polyCarIMPL(BDF4)
Astabilite(BDF4)

Code

BDF5 = { "Nom":"BDF5", "aa":[sp.Rational(300,137),-sp.Rational(300,137),sp.Rational(200,137),-sp.Rational(75,137),sp.Rational(12,137)], 
                       "bb":[0,0,0,0,0], "bm1":sp.Rational(60,137) }
## polyCar(BDF5)
Astabilite(BDF5)

Code

BDF6 = { "Nom":"BDF6", "aa":[sp.Rational(120,49),-sp.Rational(150,49),sp.Rational(400,147),-sp.Rational(75,49),sp.Rational(24,49),-sp.Rational(10,147)], 
                       "bb":[0,0,0,0,0,0], "bm1":sp.Rational(20,49)  }
## polyCar(BDF6)
Astabilite(BDF6)

Schéma	Type	Ordre maximal prévu par la barrière	Ordre réel \(\omega\)	Remarque
EE	Explicite	\(q=1\)	1	-
\(\text{AB}_2\)	Explicite	\(q=2\)	2	-
EI	Implicite, \(q\) impair	\(q+1 = 2\)	1	Non optimal : ordre réel \(1\), inférieur à l’ordre maximal \(2\)
CN	Implicite, \(q\) impair	\(q+1 = 2\)	2	-
\(\text{AM}_2\)	Implicite, \(q\) pair	\(q+2 = 4\)	3	Non optimal : ordre réel \(3\), inférieur à l’ordre maximal \(4\)
\(\text{MS}_2\)	Implicite, \(q\) pair	\(q+2 = 4\)	4	-