Stabilité absolue, méthodes multipas et problèmes raides

Objectifs du chapitre

À la fin de ce chapitre, vous serez capable de:

• réduire l'étude de la stabilité d'un schéma sur un système linéaire à celle de l'équation-test $y^{\prime }=\lambda y$, en justifiant la réduction par la diagonalisation;
• déterminer le facteur d'amplification $R(z)$ d'une méthode à un pas et en déduire sa région de stabilité absolue dans le plan complexe $z=h\lambda$;
• démontrer que la région d'Euler explicite est le disque $|1+z|\le 1$, que celle d'Euler implicite contient tout le demi-plan gauche, et énoncer ce que signifient l' et la ;

Ce que le chapitre 9 laisse ouvert

Le chapitre 9 s'est achevé sur un théorème rassurant — consistance plus stabilité entraîne convergence — et sur une observation qui ne l'est pas du tout. Reprenons-la. La suspension d'instrument du chapitre 9 ($m=2$ kg, $k=200$ N/m, $c=8$ N·s/m) conduit à l'équation $x^{\prime \prime }+4x^{\prime }+100x=0$, dont la solution décroît vers zéro avec une pseudo-période s. Intégrée par Euler explicite avec — un pas treize fois plus petit que la pseudo-période, ce qui semble prudent —, la solution calculée ne se contente pas d'être imprécise: son amplitude , et atteint mm en s là où la solution exacte est tombée sous le micromètre.

Or Euler explicite est une méthode convergente d'ordre 1, démontrée telle au théorème 9.4. Le théorème n'est pas faux; il est asymptotique. Il affirme que l'erreur est majorée par $Kh$ pour $h\le h_0$, et la constante $h_0$ y est cachée dans les hypothèses. Ce chapitre est l'étude de ce $h_0$: d'où il vient, comment on le calcule, et pourquoi il devient parfois si petit que la méthode, quoique convergente, est inutilisable.

Trois sujets sont donc à instruire, dans cet ordre: l'outil d'analyse (l'équation-test et les régions de stabilité), le phénomène qui rend l'outil indispensable (les problèmes raides), et la famille de méthodes que les bibliothèques leur opposent (les multipas, et les BDF en particulier).

L'équation-test, et pourquoi une seule équation scalaire suffit

Étudier la stabilité d'une méthode sur une équation linéaire scalaire à coefficient constant paraît dérisoire. Ce n'est pourtant pas un appauvrissement arbitraire, et l'argument qui le justifie mérite d'être écrit en entier, car c'est lui qui donne au plan complexe tout son intérêt: le $\lambda$ de (10.1) n'est pas un paramètre libre, c'est une valeur propre.

Considérons le système linéaire à coefficients constants

et supposons $A$ diagonalisable: il existe $V$ inversible et $D=\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_m)$ telles que $A=VD{V}^{-1}$, les étant les valeurs propres de et les colonnes de ses vecteurs propres (Algèbre linéaire; annexe A pour les rappels). Posons le changement d'inconnue

Le système (10.2) devient $m$ équations découplées, chacune de la forme $w_i^{\prime }={\lambda }_i{w}_i$: c'est-à-dire $m$ copies de l'équation-test (10.1). Voilà pour l'équation exacte. Le point décisif est que la méthode numérique survit au même changement de variable, pourvu qu'elle soit linéaire en ses arguments. Prenons Euler explicite appliqué à (10.2):

En multipliant à gauche par $V^{-1}$ et en posant ${\mathbf{w}}_n=V^{-1}{\mathbf{y}}_n$,

puisque $V^{-1}IV=I$ et $V^{-1}AV=D$. La matrice $I+hD$ est diagonale: la $i$-ième composante vérifie , c'est-à-dire . Le même calcul vaut pour toute méthode de Runge–Kutta et pour toute méthode multipas linéaire, parce que toutes s'écrivent, sur un système linéaire, à l'aide de polynômes ou de fractions rationnelles en , et que pour toute telle fonction .

Comme ${\mathbf{y}}_n=V{\mathbf{w}}_n$ et que $V$ est fixée, ${\mathbf{y}}_n$ reste borné si et seulement si chacune des $m$ suites $w_{n,i}$ l'est. D'où la conclusion:

Deux honnêtetés s'imposent. D'abord, la réduction suppose $A$ diagonalisable; si $A$ ne l'est pas, on passe par la forme de Jordan et la conclusion subsiste, au prix de facteurs polynomiaux en $n$. Ensuite, $V$ intervient dans les constantes: si $V$ est très mal conditionnée, «borné» peut vouloir dire «borné par $\kappa (V)$ fois quelque chose», et $\kappa (V)$ peut être énorme (chapitre 4). Les transitoires d'un système non normal ne sont pas capturés par les seules valeurs propres. Cela dit, pour les problèmes que ce chapitre traite, l'analyse par valeurs propres est exacte et suffisante.

Le facteur d'amplification et la région de stabilité absolue

Appliquons une méthode à un pas à l'équation-test. Toutes celles du chapitre 9 ont alors la même structure: $y_{n+1}$ est un multiple de $y_n$, par un nombre qui ne dépend que du produit $h\lambda$.

Le calcul de $R$ est mécanique. Pour Euler explicite, $y_{n+1}=y_n+h\lambda y_n$ donne $R(z)=1+z$. Pour Euler implicite, donne , soit . Pour le trapèze, donne

Pour RK4 le calcul vaut la peine d'être fait, car il livre le polynôme dont on tracera la frontière. Avec $f(t,y)=\lambda y$ et $z=h\lambda$, les quatre pentes du tableau de RK4 (relation 9.15) s'enchaînent:

En reportant dans $y_{n+1}=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)$ et en factorisant :

Ce résultat n'a rien d'un hasard: $R$ est le polynôme de Taylor de degré 4 de ${\mathrm{e}}^z$. C'est la traduction exacte de «RK4 est d'ordre 4», puisque la solution exacte est multipliée par ${\mathrm{e}}^z$ à chaque pas et que l'erreur locale est $R(z)-{\mathrm{e}}^z=-z^5/120+O(z^6)$. Plus généralement, une méthode à un pas est d'ordre si et seulement si .

Le vocabulaire mérite une mise en garde. On dit «la méthode est stable pour ce pas» par abus: la stabilité absolue est une propriété du produit $h\lambda$, et $\lambda$ est une donnée du problème. Une même méthode est stable à $h=1$ sur $y^{\prime }=-y$ et instable à $h=0,01$ sur $y^{\prime }=-10^{4}y$ — c'est le même point dans les deux cas.

Euler explicite: un disque, et rien de plus

Démonstration. La récurrence d'Euler explicite sur $y^{\prime }=\lambda y$ est $y_{n+1}=y_n+h\lambda y_n=(1+h\lambda )y_n$. Une récurrence immédiate donne , d'où

La suite géométrique de raison $r$ tend vers $0$ si $r<1$, est constante en module si $r=1$, et tend vers $+\infty$ si $r>1$ (dès que $y_0\ne 0$). Elle est donc bornée si et seulement si , ce qui est exactement (10.8) avec : l'ensemble des tels que est, par définition de la distance dans , le disque fermé de centre et de rayon .

Reste (10.9). Si $\lambda =-a$ avec $a>0$, alors $z=-ha$ est réel et

la borne inférieure $ha\ge 0$ étant automatique. $\square$

Le résultat appelle trois commentaires, qui sont le cœur du chapitre.

Premier commentaire: la condition n'a rien à voir avec la précision. Elle ne fait intervenir ni la tolérance demandée, ni la régularité de la solution, ni l'ordre de la méthode. C'est une contrainte structurelle, et l'exiger de $h$ n'améliore en rien l'erreur — cela empêche seulement le désastre.

Deuxième commentaire: le disque est minuscule. Il est entièrement contenu dans le demi-plan $\mathrm{Re}z\le 0$, mais il n'en occupe qu'une bulle de rayon 1. Toutes les valeurs propres à partie réelle négative et de grand module en sont exclues dès que $h$ n'est pas assez petit.

Troisième commentaire: sur l'axe réel, il y a deux seuils et non un. Pour $\lambda <0$ réel, $R=1+h\lambda$ est réel, et:

Le régime intermédiaire est stable au sens de la définition et pourtant qualitativement faux: la solution exacte est positive et monotone, la solution calculée change de signe à chaque pas. Nous le verrons à l'œuvre dans un instant.

Euler implicite: pas de condition du tout

Démonstration. La relation (9.9) appliquée à $f(t,y)=\lambda y$ s'écrit $y_{n+1}=y_n+h\lambda y_{n+1}$, soit . Comme , on peut diviser et obtenir , puis, par récurrence, . Le facteur d'amplification est donc et

ce qui est (10.10): le complémentaire du disque ouvert de centre $1$ et de rayon $1$.

Montrons que ce complémentaire contient le demi-plan gauche. Écrivons $z=x+\mathrm{i}y$ avec $x\le 0$. Alors

Comme $x\le 0$, on a $1-x\ge 1$, donc $(1-x{)}^2\ge 1$ et par conséquent

Ainsi $|1-z|\ge 1$, c'est-à-dire $z\in {\mathcal{A}}_{\text{EI}}$, et l'inégalité est stricte dès que $x<0$ ou $y\ne 0$: le module du facteur d'amplification est alors strictement inférieur à 1 et la suite tend vers , comme la solution exacte.

Enfin, si $\lambda$ est réel négatif et $h>0$, alors $z=h\lambda <0$ et

Cette inégalité est vraie pour tout $h>0$ sans exception: aucune borne supérieure sur le pas n'est imposée par la stabilité. $\square$

Le trapèze: exactement le demi-plan gauche

Avec (10.5), la condition $|R(z)|\le 1$ s'écrit $|1+z/2|\le |1-z/2|$. En posant $z=x+\mathrm{i}y$ et en élevant au carré,

La région de stabilité du trapèze est donc exactement le demi-plan $\mathrm{Re}z\le 0$ — ni plus, ni moins. Géométriquement, $|1+z/2|\le |1-z/2|$ dit que $z$ est plus proche de $-2$ que de $+2$, ce qui est la définition du demi-plan délimité par la médiatrice de ces deux points, c'est-à-dire l'axe imaginaire.

RK4: une région bornée, mais quatre fois plus large

Le polynôme (10.6) n'admet pas de description en forme close de sa région, mais celle-ci se calcule sans difficulté en résolvant $|R(z)|=1$. La frontière tracée à la figure 10.2 a été obtenue ainsi, par dichotomie colonne par colonne. Deux nombres méritent d'être retenus:

• l'intersection avec l'axe réel négatif est $z^{\star }=-2,785294$ (à $10^{-6}$ près), racine de $R(x)=+1$. Attention au signe: sur l'axe réel, $R$ ne descend jamais en dessous de $0,270395$ (minimum atteint en ), donc n'a solution réelle et la borne de l'intervalle est l'endroit où jusqu'à ;

L'intervalle réel de RK4 est donc $[-2,785;0]$, contre $[-2;0]$ pour Euler explicite: à peine 39 % plus long, pour quatre évaluations de $f$ par pas au lieu d'une. Rapporté au coût, RK4 est moins économe qu'Euler explicite en matière de stabilité pure. Son avantage est ailleurs, entièrement du côté de la précision — et c'est pourquoi, sur un problème dont le pas est dicté par la stabilité et non par la précision, monter en ordre ne sert à rien.

Notons enfin un détail que la figure montre et qu'on oublie souvent: la région de RK4 déborde légèrement dans le demi-plan droit, jusqu'à $\mathrm{Re}z=0,23742$, atteint en $\mathrm{Im}z=\pm 2,3036$. Il existe donc des problèmes dont la solution exacte croît et que RK4 amortit. C'est une anomalie sans conséquence pratique, mais c'est un bon rappel que la région est ce que le calcul donne, pas ce qu'on s'attend à trouver.

Euler implicite est L-stable, puisque $|R(z)|=1/|1-z|\to 0$. Le trapèze est A-stable mais pas L-stable, car $R(z)\to -1$ quand $|z|\to \infty$: les composantes les plus rapides, au lieu d'être écrasées, sont renvoyées presque intactes avec le signe changé. Nous mesurerons plus loin ce que cela coûte; c'est l'un des pièges les mieux dissimulés du sujet.

Le problème raide témoin, déroulé jusqu'au bout

Passons à l'expérience que ce chapitre doit à tout prix rendre visible. Le problème est celui fixé pour tout le cours:

La solution exacte est positive, monotone décroissante, et elle a pratiquement disparu après $t=0,5$ puisque ${\mathrm{e}}^{-7,5}=5,53\cdot 10^{-4}$. Rien, dans cette fonction, ne suggère la moindre difficulté.

D'après (10.9), Euler explicite est absolument stable si et seulement si

Deux pas encadrent donc ce seuil: $h=0,1$, qui le respecte, et $h=0,2$, qui le viole.

Qu'est-ce qu'un problème raide?

Un problème à une seule valeur propre, comme (10.13), illustre la mécanique mais ne montre pas le phénomène dans toute sa perversité. Sur $y^{\prime }=-15y$, la contrainte $h\le 2/15$ n'est pas absurde: la solution elle-même varie sur l'échelle de temps $1/15$, et il faudrait de toute façon un pas de cet ordre pour la représenter. Le vrai problème raide est celui où les deux échelles divergent.

Un système $2\times 2$ construit pour l'occasion

Construisons un système dont nous connaissons tout. Choisissons d'abord les valeurs propres, ${\lambda }_1=-1$ et ${\lambda }_2=-1001$, et les vecteurs propres ${\mathbf{v}}_1=(1,1{)}^{\mathsf{T}}$ et , puis reconstituons avec

Le produit donne une matrice à coefficients entiers:

La vérification se fait de tête: $A{\mathbf{v}}_1=(-501+500,500-501{)}^{\mathsf{T}}=-{\mathbf{v}}_1$ et . Comme , la solution exacte est

Les deux constantes de temps sont ${\tau }_1=1/1=1$ et ${\tau }_2=1/1001=9,99\cdot 10^{-4}$, dans un rapport . Voici ce que cela donne:

Dès $t=0,02$, la composante rapide a disparu — elle vaut deux milliardièmes — et la solution n'est plus qu'une exponentielle décroissante paisible, de constante de temps 1. Pour la représenter à $1\%$ près, un pas de $h=0,02$ suffirait largement. Or la contrainte d'Euler explicite ne se soucie pas de ce que vaut la composante rapide: elle ne connaît que la valeur propre, et ${\lambda }_2=-1001$ est là pour l'éternité dans la matrice. Donc

et il faut $\lceil 5/0,001998\rceil =2503$ pas pour atteindre $t=5$, dont $2493$ servent à intégrer une composante numériquement nulle. RK4 ne sauve rien: sa limite est $2,785/1001=2,78\cdot 10^{-3}$, soit pas — mais à quatre évaluations par pas, c'est évaluations contre .

Une vérification numérique le montre sans appel. En intégrant (10.16) jusqu'à $t=1$ par Euler explicite:

Entre $h=0,0019$ et $h=0,0021$ — une variation de $10\%$ du pas — le résultat passe de correct à $6\cdot 10^{19}$. Le seuil n'est pas une zone floue, c'est une falaise.

Ne pas confondre raideur et oscillation

La suspension d'instrument du chapitre 9 a pour valeurs propres $\lambda =-2\pm \mathrm{i}\sqrt{96}$, de module exactement 10. Elles ont la même partie réelle, donc $S=1$: le problème n'est pas raide. Et pourtant Euler explicite y est contraint à $h\le 0,04$, comme l'exercice l'a vérifié. La contrainte vient ici de la , non d'une séparation d'échelles: le pas doit résoudre l'oscillation de pseudo-période s, et donne seize pas par période, ce qui est tout juste raisonnable pour la précision aussi. La contrainte de stabilité et la contrainte de précision sont du même ordre, et c'est exactement ce qui fait que le problème n'est pas raide.

La leçon est qu'un petit pas n'est pas en soi le symptôme de la raideur. Le symptôme est l'écart entre le pas que la précision demanderait et celui que la stabilité impose. Sur (10.16), cet écart est d'un facteur 25 environ; sur les problèmes de cinétique chimique ou d'électronique, il atteint couramment $10^6$.

Ce que coûte un pas implicite, et pourquoi il reste moins cher

L'argument contre l'implicite est le coût du pas, et le chapitre 9 l'a chiffré: trois à dix fois un pas explicite. Reprenons-le en détail, puis comparons ce qui doit l'être — non pas le coût d'un pas, mais le coût du calcul entier.

Pour un système ${\mathbf{y}}^{\prime }=\mathbf{f}(t,\mathbf{y})$ de dimension $m$, un pas d'Euler implicite demande de résoudre en $\mathbf{z}\in {\mathbb{R}}^m$ l'équation non linéaire

La méthode de Newton du chapitre 2, dans sa version vectorielle, consiste à itérer

où $J=\partial \mathbf{f}/\partial \mathbf{y}$ est la matrice jacobienne. Chaque itération demande donc: une évaluation de $\mathbf{f}$, une évaluation de $J$, une factorisation $LU$ de $I-hJ$ au coût $\frac{2}{3}{m}^3$ (chapitre 3), et une descente-remontée en . Le point crucial, et il est très utilisé en pratique, est que : on garde la même matrice pendant plusieurs itérations de Newton et même pendant plusieurs pas, tant que et ne changent pas trop. C'est ce qu'on appelle un . Le pas implicite coûte alors une descente-remontée en , plus une évaluation de — et une factorisation amortie sur des dizaines de pas.

Le prédicteur d'Euler explicite ${\mathbf{z}}^{(0)}={\mathbf{y}}_n+h\mathbf{f}(t_n,{\mathbf{y}}_n)$ place le départ à de la solution, si bien que deux itérations suffisent presque toujours. Notons enfin, avec le chapitre 2, que si ( étant la constante de Lipschitz de en ) l'itération de point fixe naïve convergerait déjà — mais sur un problème raide et la condition imposerait , c'est-à-dire exactement la contrainte qu'on cherchait à fuir. C'est un détail que les exposés généraux escamotent et qui décide de tout.

Méthodes multipas: Adams, Adams–Moulton, BDF

Toutes les méthodes du chapitre 9 étaient à un pas: $y_{n+1}$ ne dépendait que de $y_n$. À chaque pas, RK4 jette quatre évaluations de $f$ et n'en garde aucune. L'idée des méthodes multipas est de réutiliser le passé: les valeurs $f_n,f_{n-1},f_{n-2}$ déjà calculées portent de l'information sur la courbure de la solution, et rien n'oblige à la redécouvrir.

Adams–Bashforth: interpoler le passé, extrapoler

On repart de l'identité exacte (9.8):

Au lieu d'approcher l'intégrande par une quadrature sur $[t_n,t_{n+1}]$, on l'interpole (chapitre 6) aux $k$ points passés $t_n,t_{n-1},\dots ,t_{n-k+1}$, où est connue, puis on intègre ce polynôme sur — c'est-à-dire de tous les nœuds. On obtient les méthodes d', explicites:

AB$k$ est d'ordre $k$ et ne coûte qu'une seule évaluation de $f$ par pas, les autres étant conservées d'un pas au précédent. Comparé à RK4 (quatre évaluations pour l'ordre 4), AB4 est quatre fois moins cher par pas. L'erreur locale de AB2 vaut $\frac{5}{12}{h}^3{y}^{\prime \prime \prime }(\xi )$, celle de AB3 .

Adams–Moulton: inclure le point d'arrivée

Si l'on ajoute $t_{n+1}$ aux nœuds d'interpolation, le polynôme est intégré à l'intérieur de son intervalle de nœuds, ce qui gagne un ordre — au prix de l'implicitation, puisque $f_{n+1}$ apparaît. Ce sont les méthodes d'Adams–Moulton:

À nombre de pas égal, AM gagne un ordre sur AB — AM3 est d'ordre 3 avec deux valeurs passées, AB3 en demande trois — et son erreur est bien plus petite: $-\frac{1}{24}{h}^4{y}^{(4)}(\xi )$ pour AM3 contre $\frac{3}{8}{h}^4{y}^{(4)}(\xi )$ pour AB3, soit .

Prédicteur–correcteur

Comme Heun au chapitre 9, on peut résoudre l'équation implicite non pas par Newton mais par une seule itération de point fixe, initialisée par une méthode explicite du même ordre. Le schéma dit PECE (predict, evaluate, correct, evaluate) apparie par exemple AB3 et AM3:

Deux évaluations de $f$ par pas, l'ordre du correcteur, et — cadeau — l'écart $|y_{n+1}-y_{n+1}^{\ast }|$ est un estimateur d'erreur a posteriori gratuit, du même genre que celui des paires emboîtées du chapitre 9. C'est la structure des solveurs d'Adams à pas et ordre variables, dont le plus célèbre est le code LSODE.

Mais il faut dire immédiatement ce que PECE ne fait pas: une seule itération de point fixe ne donne pas la stabilité du correcteur. La région de stabilité de l'appariement est proche de celle du prédicteur explicite, donc bornée et petite. Un prédicteur–correcteur d'Adams n'est pas un schéma pour problème raide, quelle que soit la A-stabilité du correcteur pris isolément. Sur un problème raide, il faut résoudre l'équation implicite pour de bon, par Newton.

BDF: la famille que les solveurs raides utilisent vraiment

Les méthodes d'Adams discrétisent l'intégrale. Les BDF (backward differentiation formulas, formules de différentiation rétrograde) font l'inverse: elles interpolent les $y$ — et non les $f$ — aux points $t_{n+1},t_n,\dots ,t_{n+1-k}$, puis imposent que la du polynôme interpolateur en vaille . Il vient

où $\nabla$ est l'opérateur de différence rétrograde, $\nabla y_{n+1}=y_{n+1}-y_n$. En développant:

BDF$k$ est d'ordre $k$, et n'évalue $f$ qu'au seul point $t_{n+1}$ — ce qui rend le pas de Newton particulièrement simple, avec la même matrice $I-h{\beta }_{k}J$ que pour Euler implicite.

Voici ce qui les distingue de tout le reste. Leurs angles de A($\alpha$)-stabilité, calculés à partir du lieu frontière $\zeta \mapsto \rho (\zeta )/\sigma (\zeta )$:

BDF1 et BDF2 sont A-stables; au-delà, l'angle se referme, lentement d'abord (à l'ordre 4, le secteur couvre encore $\pm 73°$, ce qui suffit pour les valeurs propres réelles négatives et pour tout ce qui n'est pas trop oscillant), puis brutalement. À $k=7$, le polynôme $\rho$ a des racines de module $1,022>1$ et la méthode n'est même plus convergente, quel que soit son ordre. C'est pourquoi tous les codes BDF s'arrêtent à l'ordre 5 ou 6.

Et surtout: BDF$k$ est L-stable pour $k\le 2$ et, plus généralement, $R\to 0$ à l'infini dans son secteur, parce que $\sigma (\zeta )={\beta }_k{\zeta }^k$ est un monôme. Les composantes rapides sont donc écrasées, et non renvoyées comme avec le trapèze. C'est la conjonction de trois propriétés — un secteur large, l'amortissement des modes rapides, et une seule évaluation de par pas — qui fait des BDF le cœur de tous les intégrateurs raides, de et jusqu'au de SciPy.

La condition des racines, la zéro-stabilité, et les barrières de Dahlquist

Il reste une question que les méthodes à un pas ne posaient pas. Une méthode multipas transforme une équation différentielle du premier ordre en une récurrence d'ordre $k$, qui possède $k$ solutions indépendantes alors que l'équation différentielle n'en a qu'une. Les $k-1$ solutions supplémentaires sont des artefacts du schéma: on les appelle racines parasites. Si l'une d'elles croît, elle finit par dominer la solution utile — et cela peut arriver à une méthode parfaitement consistante, d'ordre aussi élevé qu'on veut.

La consistance ne suffit plus

Commençons par rappeler ce que consistance veut dire pour (10.20). En injectant la solution exacte et en développant par Taylor, on trouve que la méthode est consistante (ordre $\ge 1$) si et seulement si

et d'ordre $p$ si et seulement si, pour $q=1,\dots ,p$,

La première condition, $\rho (1)=0$, dit que $\zeta =1$ est toujours racine de $\rho$: c'est la racine principale, celle qui porte la solution qu'on cherche. Les $k-1$ autres sont les parasites, et (10.27) ne dit rien d'elles.

Le nom «zéro-stabilité» vient de ce que la condition porte sur le comportement du schéma quand $h\to 0$, c'est-à-dire sur $\rho$ seul, $\sigma$ disparaissant avec $h$. Elle n'a rien à voir avec la stabilité absolue, qui porte sur un $h$ fixé et fait intervenir les deux polynômes. Les deux mots se ressemblent et ne désignent pas la même chose: la zéro-stabilité est la condition minimale sans laquelle la méthode ne converge pas du tout; la stabilité absolue est une condition supplémentaire sur le pas.

Démonstration. Notons ${\zeta }_1,\dots ,{\zeta }_r$ les racines distinctes de $\rho$, de multiplicités $m_1,\dots ,m_r$ avec (le degré de vaut exactement puisque ).

Étape 1: la forme générale des solutions. L'ensemble $E$ des suites complexes vérifiant ${\sum }_j{\alpha }_j{u}_{n+j}=0$ est un espace vectoriel de dimension $k$, car une telle suite est entièrement déterminée par et que, réciproquement, tout -uplet en engendre une (on résout en , licite car ). Montrons que les suites

appartiennent à $E$. Partons de l'identité, valable pour tout $\zeta$ et tout $n$,

et appliquons-lui $q$ fois l'opérateur $\theta =\zeta \frac{\mathrm{d}}{\mathrm{d}\zeta }$. À gauche, $\theta$ multiplie chaque terme ${\zeta }^{n+j}$ par son exposant, donc ${\theta }^q$ le multiplie par :

À droite, la formule de Leibniz donne une combinaison linéaire de $\rho (\zeta ),{\rho }^{\prime }(\zeta ),\dots ,{\rho }^{(q)}(\zeta )$, à coefficients polynomiaux en $\zeta$ et $n$. Évaluons en $\zeta ={\zeta }_i$, racine de multiplicité : toutes les dérivées jusqu'à l'ordre y sont nulles, donc le membre de droite s'annule et

ce qui est exactement la récurrence appliquée à $u_n=n^q{\zeta }_i^n$ — car $(n+j{)}^q{\zeta }_i^{n+j}$ est bien . Ces suites sont linéairement indépendantes (c'est le déterminant de Vandermonde confluent, non nul), donc elles forment une base de .

Étape 2: (i) implique (ii). Supposons la condition des racines vérifiée et soit $u_n=n^q{\zeta }_i^n$ l'un des éléments de la base. Si $|{\zeta }_i|<1$, alors (l'exponentielle l'emporte sur la puissance), donc la suite est bornée. Si , la condition des racines impose , donc et : bornée. Toute suite de est une combinaison linéaire fixe de ces suites bornées, donc bornée. Pour la constante : l'application linéaire va de dans l'espace des suites bornées; en notant la solution issue du -ième vecteur de la base canonique et (fini par ce qui précède), la linéarité donne , et ne dépend que de .

Étape 3: (ii) implique (i), par contraposition. Supposons la condition des racines violée. Ou bien il existe une racine $\zeta$ avec $|\zeta |>1$, et alors $u_n={\zeta }^n$ appartient à $E$ et $|u_n|=|\zeta {|}^n\to +\infty$. Ou bien il existe une racine de module 1 et de multiplicité , et alors appartient à (étape 1 avec ) et . Dans les deux cas (ii) est en défaut.

Étape 4: la conséquence. Appliquons la méthode au problème trivial $y^{\prime }=0$, $y(0)=0$, dont la solution exacte est $y\equiv 0$ et que toute méthode consistante devrait intégrer sans peine. Ici $f\equiv 0$, donc (10.20) se réduit à la récurrence homogène ${\sum }_j{\alpha }_j{y}_{n+j}=0$. Supposons la condition des racines violée et soit une solution non bornée fournie par l'étape 3. Prenons comme valeurs de départ pour : elles tendent vers les valeurs exactes quand , à la vitesse , ce que toute procédure de démarrage raisonnable garantit. Par linéarité, pour tout . Au dernier nœud , on obtient

avec $|u_N|\ge |\zeta {|}^N$ (cas d'une racine de module $>1$) ou $|u_N|=N=(b-a)/h$ (cas d'une racine double de module 1). Dans le premier cas quand , car l'exponentielle en écrase le facteur ; dans le second, , qui ne tend pas vers . L'erreur ne tend donc pas vers zéro: la méthode n'est pas convergente.

Un contre-exemple qui vaut tous les discours

Considérons la méthode à deux pas

Les conditions d'ordre (10.28) se vérifient à la main avec $\alpha =(-5,4,1)$ et $\beta =(2,4,0)$: elles sont satisfaites pour $q=1,2,3$ et échouent pour $q=4$. , un ordre que la première barrière de Dahlquist autorise pour . Mais

dont la seconde racine vaut $-5$, de module 5. La condition des racines est violée. Voici ce que donne l'intégration de $y^{\prime }=-y$, $y(0)=1$ sur $[0,1]$, avec des valeurs de départ exactes ($y_0=1$, ), alors que la réponse est :

Diviser le pas par deux multiplie l'erreur par $10^{13}$. C'est le comportement exactement inverse de celui d'une méthode convergente, et il s'explique en une ligne: la racine parasite $-5$ engendre la solution $(-5{)}^n$, et au dernier nœud $n=N=1/h$ elle vaut $5^{1/h}$, qui explose quand . L'erreur d'arrondi suffit à l'exciter — même avec des valeurs de départ exactes, la première opération arrondie introduit une composante parasite de l'ordre de , et la ramène au premier plan. Aucune précaution d'implémentation n'y change rien: le défaut est dans la formule.

Les deux barrières

Les deux résultats sont dus à Dahlquist (1956 et 1963); leurs démonstrations reposent sur des arguments d'analyse complexe (fonctions à partie réelle positive, transformation $\zeta \mapsto (1+w)/(1-w)$) qui sortent du cadre de ce cours. Nous les admettons, mais il faut voir ce qu'ils interdisent.

La première barrière est atteinte par la méthode de Milne–Simpson,

d'ordre 4 avec $k=2$ pas, et zéro-stable puisque $\rho (\zeta )={\zeta }^2-1$ a pour racines $\pm 1$, toutes deux simples. C'est le meilleur qu'une méthode à deux pas puisse faire. Et pourtant elle est inutilisable sur un problème dissipatif: on vérifie que sa région de stabilité absolue ne contient aucun point du demi-plan gauche ouvert — elle se réduit au segment $[-\mathrm{i}\sqrt{3},\mathrm{i}\sqrt{3}]$ de l'axe imaginaire. Autrement dit, pour tout réel négatif et tout , l'une des deux racines a un module supérieur à 1. La racine , simple et de module 1, est une racine parasite «à la limite», et la moindre dissipation la fait sortir du disque. Ordre maximal, zéro-stabilité, et stabilité absolue nulle: les trois propriétés sont bien indépendantes.

La seconde barrière est celle qui organise le paysage. Elle dit qu'il n'existe pas de méthode multipas A-stable d'ordre 3, ni d'ordre 4, ni davantage. Ce n'est pas un aveu d'impuissance des constructeurs de méthodes: c'est un théorème. Deux échappatoires seulement:

• Renoncer à un peu d'A-stabilité: c'est la voie des BDF, A($\alpha$)-stables avec $\alpha =86°$, $73°$, $52°$, $18°$ aux ordres 3 à 6. Le secteur suffit dès que les valeurs propres du problème ne sont pas trop proches de l'axe imaginaire, ce qui est le cas de la plupart des problèmes dissipatifs (cinétique chimique, diffusion, circuits).
• Renoncer au multipas: la barrière ne concerne que les méthodes de la forme (10.20). Les méthodes de Runge–Kutta implicites y échappent complètement. Les méthodes de Gauss à $s$ étages sont A-stables et d'ordre $2s$ — donc d'ordre arbitrairement élevé — et celles de Radau IIA à $s$ étages sont L-stables et d'ordre $2s-1$. Leur prix est que le système non linéaire à résoudre à chaque pas a inconnues au lieu de .

C'est ce partage du monde que la figure 10.3 résume: à gauche de la barrière, les deux familles coexistent; à droite, les multipas doivent abandonner quelque chose et les Runge–Kutta implicites non.

Ce que fait une vraie bibliothèque

Le lecteur qui ouvrira un code scientifique rencontrera ce vocabulaire sous une forme légèrement différente, et il vaut la peine de faire la traduction. En Python, l'interface de référence est scipy.integrate.solve_ivp — que nous ne pouvons pas exécuter dans la page, puisque ce cours s'en tient à la bibliothèque standard, mais dont il faut connaître la forme:

from scipy.integrate import solve_ivp

def f(t, y):
    return [-501 * y[0] + 500 * y[1], 500 * y[0] - 501 * y[1]]

sol = solve_ivp(f, (0.0, 5.0), [1.0, 0.0], method="BDF"

Les choix offerts se rangent exactement selon ce chapitre:

Trois remarques pratiques, qui sont autant de corollaires du chapitre. D'abord, RK45 sur un problème raide ne plante pas: il réduit son pas jusqu'à respecter la stabilité, et il rend la bonne réponse — après un temps qui peut être mille fois plus long. Un solveur qui «rame» sans diverger est le symptôme le plus fréquent de la raideur, et c'est exactement ce que la figure 10.1 explique. Ensuite, l'argument jac fournit la jacobienne analytique; sans lui, le solveur l'approche par différences finies, ce qui coûte $m$ évaluations de $\mathbf{f}$ et introduit l'erreur de troncature du chapitre 8 — la donner quand on la connaît est le meilleur rapport travail/gain de tout le sujet. Enfin, rtol et atol contrôlent la précision, pas la stabilité: sur un problème raide intégré par RK45, desserrer rtol ne fait pas gagner de temps, parce que ce n'est pas la précision qui limite le pas.

Synthèse

• Sur un système linéaire à matrice diagonalisable, la diagonalisation découple le schéma en autant de copies de l'équation-test $y^{\prime }=\lambda y$ qu'il y a de valeurs propres. Le pas admissible est donc fixé par la position, dans le plan complexe, des points $h{\lambda }_i$ — et le plus contraignant décide. Pour un problème non linéaire, on applique le même raisonnement à la jacobienne gelée: c'est une heuristique, et elle prédit juste.
• Le facteur d'amplification $R(z)$, $z=h\lambda$, résume une méthode à un pas: polynôme si elle est explicite, fraction rationnelle sinon, et toujours pour une méthode d'ordre . La est le disque pour Euler explicite (théorème 10.1, d'où ), le complémentaire du disque pour Euler implicite (théorème 10.2, sans aucune condition sur ), exactement le demi-plan gauche pour le trapèze, et un domaine borné coupant l'axe réel en pour RK4.

Exercices

Vous pouvez afficher le corrigé directement sous chaque énoncé après avoir cherché la solution.

Références

• Hairer, E. et Wanner, G., Solving Ordinary Differential Equations II — Stiff and Differential-Algebraic Problems, 2e éd., Springer, Berlin, chap. IV et V (la référence sur les régions de stabilité, les BDF et les barrières de Dahlquist).
• Hairer, E., Nørsett, S. P. et Wanner, G., Solving Ordinary Differential Equations I — Nonstiff Problems, 2e éd., Springer, Berlin, chap. III (méthodes multipas, condition des racines, théorème d'équivalence).
• Quarteroni, A., Sacco, R. et Saleri, F., Méthodes numériques — algorithmes, analyse et applications, Springer, Milan, chap. 11 (stabilité absolue, méthodes multipas, problèmes raides).
• Rappaz, J. et Picasso, M., Introduction à l'analyse numérique, Presses polytechniques et universitaires romandes, Lausanne, chap. 6.
• Burden, R. L. et Faires, J. D., Numerical Analysis, Cengage, Boston, chap. 5 (sections sur les méthodes multipas et les équations raides).
• Dahlquist, G., «A special stability problem for linear multistep methods», BIT Numerical Mathematics, vol. 3, 1963 (l'article où la seconde barrière est démontrée).