Méthode linéaire à pas multiples

Les méthodes linéaires à pas multiples sont des schémas de calcul utilisés en mathématiques pour la résolution numérique des équations différentielles ordinaires.

En général, les méthodes de résolution numérique permettent de calculer de proche en proche les valeurs de la fonction solution de l'équation différentielle à partir d'un point initial dont la valeur est connue. L'écart entre deux points successifs d'un intervalle est appelé le pas. Les méthodes à un seul pas (comme la méthode d'Euler) utilisent uniquement la valeur de la fonction et de sa dérivée au dernier point connu pour calculer la valeur de la fonction au point suivant.

Il existe des méthodes plus avancées telles que les méthodes de Runge-Kutta qui utilisent des points intermédiaires de l'intervalle pour raffiner les calculs, mais ces informations ne sont pas conservées d'une étape à la suivante.

Les méthodes à pas multiples (ou multi-pas) cherchent à gagner en efficacité en conservant certaines informations des étapes précédentes du calcul. En particulier, les méthodes linéaires à pas multiple utilisent une combinaison linéaire des valeurs de la fonction et de sa dérivée aux points précédents[1].

Définitions

Considérons l'équation différentielle ordinaire à condition initiale suivante: $${\displaystyle y'=f(t,y),\quad y(t_{0})=y_{0}.}$$

Une méthode de résolution numérique est une approximation de la valeur de la fonction y(t) à des valeurs discrètes t_i : $${\displaystyle y_{i}\approx y(t_{i})\quad {\text{et}}\quad t_{i}=t_{0}+ih,}$$ où h est le pas de temps (parfois noté Δt) et i est un entier.

Les méthodes à pas multiples utilisent des informations provenant de s étapes précédentes pour calculer la valeur suivante. Pour calculer la valeur à t_n+s on utilisera donc les valeurs à t_n+s, t_n+s–1, t_n+s–2,...t_{n .}

Les méthodes linéaires à pas multiple utilisent une combinaison linéaire des y_i et des f(t_i,y_i) pour calculer la valeur de y à l'étape suivante. On a donc la formule générale suivante :$${\displaystyle {\begin{aligned}&y_{n+s}+a_{s-1}\cdot y_{n+s-1}+a_{s-2}\cdot y_{n+s-2}+\cdots +a_{0}\cdot y_{n}\\&\qquad {}=h\cdot \left(b_{s}\cdot f(t_{n+s},y_{n+s})+b_{s-1}\cdot f(t_{n+s-1},y_{n+s-1})+\cdots +b_{0}\cdot f(t_{n},y_{n})\right)\\&\Leftrightarrow \sum _{j=0}^{s}a_{j}y_{n+j}=h\sum _{j=0}^{s}b_{j}f(t_{n+j},y_{n+j}),\end{aligned}}}$$ avec a_s = 1 .

Le choix des coefficients a_j et b_j, spécifique à la méthode employée, est souvent un compromis entre la réduction de l'erreur par rapport à la vraie solution et la difficulté des calculs.

On peut dès lors distinguer les méthodes explicites et implicites. Si b_s est nul, alors la méthode est explicite. En effet, on peut calculer directement y_n+s, à partir des autres valeurs connues. Inversement si b_s est non nul, la méthode est implicite, puisque la valeur de y_n+s dépend de la valeur de f(t_n+s,y_n+s). On ne pourra donc calculer cette valeur directement ; il faudra résoudre numériquement une équation d'inconnue y_n+s, par exemple par la méthode de Newton.

Parfois, une première méthode explicite à pas multiple est utilisée pour prédire une valeur de y_n+s, qui est ensuite utilisée dans une seconde formule implicite pour en corriger la valeur selon le principe des méthodes prédicteur-correcteur .

Exemples

Soit l'équation différentielle suivante $${\displaystyle y'=f(t,y)=y,\quad y(0)=1.}$$dont la solution exacte est la fonction exponentielle ${\displaystyle y(t)=\mathrm {e} ^{t}}$ .

Méthode d'Euler

Une méthode numérique courante est la méthode d'Euler : $${\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n}).}$$Cette méthode n'utilise qu'un seul pas. Par exemple avec h = 0,5 : $${\displaystyle {\begin{aligned}y_{1}&=y_{0}+hf(t_{0},y_{0})=1+{\tfrac {1}{2}}\cdot 1=1,5,\\y_{2}&=y_{1}+hf(t_{1},y_{1})=1,5+{\tfrac {1}{2}}\cdot 1,5=2,25,\\y_{3}&=y_{2}+hf(t_{2},y_{2})=2,25+{\tfrac {1}{2}}\cdot 2,25=3,375,\\y_{4}&=y_{3}+hf(t_{3},y_{3})=3,375+{\tfrac {1}{2}}\cdot 3,375=5,0625.\end{aligned}}}$$

Méthode d'Adams–Bashforth

Un exemple de schéma à pas double est la méthode d'Adams-Bashforth d'ordre 2.$${\displaystyle y_{n+2}=y_{n+1}+{\tfrac {3}{2}}hf(t_{n+1},y_{n+1})-{\tfrac {1}{2}}hf(t_{n},y_{n}).}$$

Cette méthode prend en compte deux valeurs successives : y_n et y_n+1, pour calculer y_n+2. Pour initialiser le calcul à l'étape 1 (où l'on ne dispose que de la condition initiale y₀) on peut utiliser la méthode d'Euler.

Dans ces conditions, les 4 premières étapes donnent, à 10^-4 près: $${\displaystyle {\begin{aligned}y_{2}&=y_{1}+{\tfrac {3}{2}}hf(t_{1},y_{1})-{\tfrac {1}{2}}hf(t_{0},y_{0})=1,5+{\tfrac {3}{2}}\cdot {\tfrac {1}{2}}\cdot 1,5-{\tfrac {1}{2}}\cdot {\tfrac {1}{2}}\cdot 1=2,375,\\y_{3}&=y_{2}+{\tfrac {3}{2}}hf(t_{2},y_{2})-{\tfrac {1}{2}}hf(t_{1},y_{1})=2,375+{\tfrac {3}{2}}\cdot {\tfrac {1}{2}}\cdot 2,375-{\tfrac {1}{2}}\cdot {\tfrac {1}{2}}\cdot 1,5=3,7812,\\y_{4}&=y_{3}+{\tfrac {3}{2}}hf(t_{3},y_{3})-{\tfrac {1}{2}}hf(t_{2},y_{2})=3,7812+{\tfrac {3}{2}}\cdot {\tfrac {1}{2}}\cdot 3,7812-{\tfrac {1}{2}}\cdot {\tfrac {1}{2}}\cdot 2,375=6,0234.\end{aligned}}}$$

La solution exacte à t₄ = 2 est e^2 = 7,3891... ; la méthode d'Adams-Bashforth est plus précise que celle d'Euler.

Familles de méthodes

Trois familles de schémas linéaires multi-pas sont couramment utilisées : les méthodes d'Adams-Bashforth, les méthodes d'Adams-Moulton et les formules de différences finies rétrogrades.

Méthodes d'Adams-Bashforth

Les schémas d’Adams-Bashforth sont des méthodes explicites. On pose a_s–1= –1, et les autres coefficients a_j sont nuls. Les coefficients b_j sont choisis de telle sorte que la méthode soit d'ordre s, ce qui les détermine de manière unique.

À titre d'exemple, voici les coefficients formules utilisées jusqu'à l'ordre 5[2],[3].

$${\displaystyle {\begin{aligned}y_{n+1}&=y_{n}+hf(t_{n},y_{n})\\y_{n+2}&=y_{n+1}+h\left({\frac {3}{2}}f(t_{n+1},y_{n+1})-{\frac {1}{2}}f(t_{n},y_{n})\right),&\\y_{n+3}&=y_{n+2}+h\left({\frac {23}{12}}f(t_{n+2},y_{n+2})-{\frac {16}{12}}f(t_{n+1},y_{n+1})+{\frac {5}{12}}f(t_{n},y_{n})\right),\\y_{n+4}&=y_{n+3}+h\left({\frac {55}{24}}f(t_{n+3},y_{n+3})-{\frac {59}{24}}f(t_{n+2},y_{n+2})+{\frac {37}{24}}f(t_{n+1},y_{n+1})-{\frac {9}{24}}f(t_{n},y_{n})\right),\\y_{n+5}&=y_{n+4}+h\left({\frac {1901}{720}}f(t_{n+4},y_{n+4})-{\frac {2774}{720}}f(t_{n+3},y_{n+3})+{\frac {2616}{720}}f(t_{n+2},y_{n+2})-{\frac {1274}{720}}f(t_{n+1},y_{n+1})+{\frac {251}{720}}f(t_{n},y_{n})\right).\end{aligned}}}$$

Ces coefficients b_j se calculent en déterminant une interpolation polynomiale P de la fonction f de degré s-1 telle que $${\displaystyle P(t_{n+i})=f(t_{n+i},y_{n+i}),\qquad {\text{pour }}i=0,\ldots ,s-1.}$$Ce polynôme s'obtient par la formule de l'interpolation de Lagrange avec valeurs précédemment calculées de la fonction f$${\displaystyle P(t)=\sum _{j=0}^{s-1}{\frac {(-1)^{s-j-1}f(t_{n+j},y_{n+j})}{j!(s-j-1)!h^{s-1}}}\prod _{i=0 \atop i\neq j}^{s-1}(t-t_{n+i}).}$$Le polynôme P est, par construction, une bonne approximation de la fonction f. Ce qui permet de substituer à f dans l'équation différentielle originale

${\displaystyle y^{\prime }=f(y,t)\approx P(t)}$

Or cette nouvelle équation peut être résolue simplement : il suffit d'intégrer P :$${\displaystyle y_{n+s}=y_{n+s-1}+\int _{t_{n+s-1}}^{t_{n+s}}p(t)\,\mathrm {d} t.}$$On obtient donc ainsi les formules pour les coefficients ${\displaystyle b_{j}}$ de la méthode d'Adams-Bashforth en identifiant après intégration : $${\displaystyle b_{s-j-1}={\frac {(-1)^{j}}{j!(s-j-1)!}}\int _{0}^{1}\prod _{i=0 \atop i\neq j}^{s-1}(u+i)\,\mathrm {d} u,\qquad {\text{pour }}j=0,\ldots ,s-1.}$$La substitution de f par son polynôme interpolateur P engendre une erreur de l'ordre de h^s, ce qui démontre que la méthode d'Adams-Bashforth à s étapes est bien d'ordre s.

Les méthodes d'Adams-Bashforth furent originellement conçues par John Couch Adams pour résoudre une équation différentielle modélisant l'action capillaire due à Francis Bashforth[4],[5].

Méthodes d'Adams-Moulton

La version implicite du même procédé donne les méthodes d'Adams-Moulton. Là encore on pose a_s–1 = -1 et les autres coefficients a_j sont nuls. On utilise également l'interpolation de Lagrange pour approximer f ; en revanche, on s'autorise à ce que le polynôme soit d'ordre s + 1 ; autrement dit, b_s sera non nul.

Ci-dessous, on détaille les valeurs des coefficients pour s = 0, 1, 2, 3, 4 sont :$${\displaystyle {\begin{aligned}y_{n}&=y_{n-1}+hf(t_{n},y_{n}),\\y_{n+1}&=y_{n}+{\frac {1}{2}}h\left(f(t_{n+1},y_{n+1})+f(t_{n},y_{n})\right),\\y_{n+2}&=y_{n+1}+h\left({\frac {5}{12}}f(t_{n+2},y_{n+2})+{\frac {8}{12}}f(t_{n+1},y_{n+1})-{\frac {1}{12}}f(t_{n},y_{n})\right),\\y_{n+3}&=y_{n+2}+h\left({\frac {9}{24}}f(t_{n+3},y_{n+3})+{\frac {19}{24}}f(t_{n+2},y_{n+2})-{\frac {5}{24}}f(t_{n+1},y_{n+1})+{\frac {1}{24}}f(t_{n},y_{n})\right),\\y_{n+4}&=y_{n+3}+h\left({\frac {251}{720}}f(t_{n+4},y_{n+4})+{\frac {646}{720}}f(t_{n+3},y_{n+3})-{\frac {264}{720}}f(t_{n+2},y_{n+2})+{\frac {106}{720}}f(t_{n+1},y_{n+1})-{\frac {19}{720}}f(t_{n},y_{n})\right).\end{aligned}}}$$On remarque que le cas s = 0 correspond à la méthode d'Euler implicite et le cas s = 1 à la méthode des trapèzes.

Le calcul des coefficients en général est similaire à celui vu plus haut pour la méthode Adams-Bashforth ; cependant, le polynôme d'interpolation utilise non seulement les points ${\displaystyle t_{n-1},\dots ,t_{n-s}}$, mais aussi ${\displaystyle t_{n}}$ .

Les coefficients sont donnés par $${\displaystyle b_{s-j}={\frac {(-1)^{j}}{j!(s-j)!}}\int _{0}^{1}\prod _{i=0 \atop i\neq j}^{s}(u+i-1)\,\mathrm {d} u,\qquad {\text{pour }}j=0,\ldots ,s.}$$

Le schéma d’Adams–Moulton sont également dues à John Couch Adams, comme celui d’Adams–Bashforth. Forest Ray Moulton réalisa qu'elles pouvaient être utilisées en tandem avec les méthodes d'Adams-Bashforth comme paire prédicteur-correcteur[6].

Formules de différences finies rétrogrades

Les méthodes de différences finies rétrogrades sont des méthodes implicites dans lesquelles les coefficients b_s–1,..., b₀ sont nuls et les autres coefficients sont choisis de telle sorte que la méthode atteigne l'ordre s. Ces méthodes sont particulièrement utiles pour résoudre des équations différentielles raides.

Propriétés

On s'intéresse ici aux propriétés importantes de ces méthodes multi-pas : leur convergence, l'ordre de leur erreur et leur stabilité.

Cohérence et ordre

On souhaite d'abord vérifier que la méthode est cohérente, c'est à dire si l'équation aux différences finies $${\displaystyle {\begin{aligned}&a_{s}y_{n+s}+a_{s-1}y_{n+s-1}+a_{s-2}y_{n+s-2}+\cdots +a_{0}y_{n}\\&\qquad {}=h\left(b_{s}f(t_{n+s},y_{n+s})+b_{s-1}f(t_{n+s-1},y_{n+s-1})+\cdots +b_{0}f(t_{n},y_{n})\right),\end{aligned}}}$$est une bonne approximation de l'équation différentielle ${\displaystyle y'=f(t,y)}$.

Formellement, une méthode sera cohérente si l'erreur de troncature locale tend vers zéro plus rapidement que la taille du pas h lorsque h tend vers zéro. L'erreur de troncature locale est la différence entre le résultat ${\displaystyle y_{n+s}}$ donné par le calcul, en supposant que toutes les valeurs précédentes ${\displaystyle y_{n+s-1},\ldots ,y_{n}}$ soient exactes, et la solution exacte de l'équation à l'instant ${\displaystyle t_{n+s}}$ . Un calcul par série de Taylor montre qu'une méthode linéaire à pas multiples est cohérente si et seulement si $${\displaystyle \sum _{k=0}^{s-1}a_{k}=-1\quad {\text{et}}\quad \sum _{k=0}^{s}b_{k}=s+\sum _{k=0}^{s-1}ka_{k}.}$$On peut vérifier que c'est bien le cas pour toutes les méthodes vues ci-dessus.

Une fois assurés de la cohérence, on cherche à évaluer l'ordre de l'erreur d'approximation. On dit qu'une méthode est d'ordre p si l'erreur locale est en ${\displaystyle O(h^{p+1})}$ lorsque h tend vers zéro. Pour les méthodes linéaires à pas multiples, il suffit que : $${\displaystyle \sum _{k=0}^{s-1}a_{k}=-1\quad {\text{et}}\quad q\sum _{k=0}^{s}k^{q-1}b_{k}=s^{q}+\sum _{k=0}^{s-1}k^{q}a_{k}{\text{ pour }}q=1,\ldots ,p.}$$Comme vu précédemment, la méthode d'Adams-Bashforth à s étapes est d'ordre s, tandis que celle d'Adams-Moulton est d'ordre s + 1[2].

On résume souvent ces conditions à l'aide des deux polynômes caractéristiques $${\displaystyle \rho (z)=z^{s}+\sum _{k=0}^{s-1}a_{k}z^{k}\quad {\text{et}}\quad \sigma (z)=\sum _{k=0}^{s}b_{k}z^{k}.}$$Les conditions précédentes se pour que la méthode soit d'ordre p est alors$${\displaystyle \rho (e^{h})-h\sigma (e^{h})=O(h^{p+1})\quad {\text{quand }}h\to 0.}$$ En particulier, la méthode est cohérente si elle est au moins d'ordre 1, ce qui est le cas si ${\displaystyle \rho (1)=0}$ et ${\displaystyle \rho '(1)=\sigma (1)}$ .

Stabilité et convergence

La solution numérique d'une méthode à pas multiples dépend non seulement de la condition initiale y₀, mais aussi des s – 1 autres valeurs de départ : y₁,...y_s . La question de la stabilité de la solution par rapport aux variations de ces valeurs initiales se pose alors.

Une méthode linéaire à pas multiple est dite stable en zéro pour une équation différentielle sur un intervalle de temps donné, si une perturbation dans les valeurs de départ d'amplitude ε engendre une erreur de la solution numérique sur cet intervalle de temps ne dépassant pas Kε, où K est une constante indépendante du pas h. Il suffit de vérifier cette condition sur l'équation y' = 0, ce qui explique l'appellation de stabilité en zéro[7].

Une condition nécessaire et suffisante pour qu'une méthode soit stable en zéro est la condition des racines. Il faut et suffit que les racines du polynôme caractéristique ρ soient toutes de module inférieur ou égal à 1, et que les racines de module 1 soient de multiplicité 1, pour que la méthode soit stable en zéro.

Ce théorème d'équivalence de Dahlquist est un théorème analogue au théorème d'équivalence de Lax pour les méthodes aux différences finies.

De plus, pour une méthode stable en zéro d'ordre p, alors l'erreur globale (la différence entre la solution numérique et la solution exacte à en un point donné) est d'ordre O(h^p)[7].

Pour une méthode convergente :

• une méthode dont z = 1 est la seule racine de module 1 sera dite fortement stable
• une méthode ayant plusieurs racines de module 1, sans qu'aucune soit multiple sera dite relativelement stable.

Puisqu'il faut que 1 soit racine du polynome pour que la méthode soit convergente, toute méthode convergente sera nécéssairement soit fortement stable, soit relativement stable.

Pour évaluer les performances des méthodes linéaires multi-étapes sur des équations raides, considérons par exemple l'équation linéaire y' = λ y .

Une méthode linéaire à pas multiples appliquée à cette équation différentielle avec un pas h donne une relation de récurrence linéaire qui fait apparaître un troisième polynôme caractéristique : $${\displaystyle \pi (z;h\lambda )=(1-h\lambda \beta _{s})z^{s}+\sum _{k=0}^{s-1}(\alpha _{k}-h\lambda \beta _{k})z^{k}=\rho (z)-h\lambda \sigma (z).}$$Ce polynôme est appelé le polynôme de stabilité de la méthode multi-étapes. Si toutes ses racines sont de module inférieur à un, alors la solution numérique de la méthode multi-étapes tend vers zéro. On dit que la méthode est absolument stable pour cette valeur de h λ.

La méthode est dite A-stable si elle est absolument stable pour tout h λ de partie réelle négative. La région de stabilité absolue est l'ensemble de tous les h λ pour lesquels la méthode multi-étapes est absolument stable.

Exemple

Soit le schéma d’Adams-Bashforth à 3 étapes :$${\displaystyle y_{n+3}=y_{n+2}+h\left({23 \over 12}f(t_{n+2},y_{n+2})-{4 \over 3}f(t_{n+1},y_{n+1})+{5 \over 12}f(t_{n},y_{n})\right).}$$Le premier polynôme caractéristique ρ est alors $${\displaystyle \rho (z)=z^{3}-z^{2}=z^{2}(z-1)}$$dont les racines sont z = 0 et z = 1.

Les conditions de stabilité sont satisfaites : comme z = 1 est une racine simple, la méthode est fortement stable.

L'autre polynôme caractéristique est $${\displaystyle \sigma (z)={\frac {23}{12}}z^{2}-{\frac {4}{3}}z+{\frac {5}{12}}}$$

Première et deuxième barrières de Dahlquist

Germund Dahlquist établit les deux résultats suivants concernant la convergence et la A-stabilité des schémas linéaires à pas multiples.

Première barrière de Dahlquist

L'ordre de convergence de toute méthode linéaire à pas multiple à q étapes zéro-stable est au plus q +1 si q est impair, et au plus q + 2 si q est pair.

De plus, si la méthode est explicite, alors l'ordre sera au plus q[8].

Deuxième barrière Dahlquist

Si un schéma linéaire à pas multiple est explicite, alors il ne peut être A-stable.

S'il est implicite et A-stable, alors il sera au plus d'ordre 2. La méthode des trapèzes est la méthode d'ordre 2 A-stable ayant la plus petite constante d'erreur[9].

Références

• (en) Cet article est partiellement ou en totalité issu de l’article de Wikipédia en anglais intitulé « Linear multistep method » (voir la liste des auteurs).

1. « Les Méthodes multi-pas :: Cours Tan », sur feelpp.github.io (consulté le 8 octobre 2024)
2. (en) Ernst Hairer, Gerhard Wanner et Syvert P. Nørset, Solving Ordinary Differential Equations I : Nonstiff problems, Springer, 1993 (ISBN 978-3-540-56670-0, DOI 10.1007/978-3-540-78862-1, lire en ligne), III.1
3. John C. Butcher, Numerical methods for ordinary differential equations, Wiley [u.a.], 2003 (ISBN 978-0-471-96758-3), p. 103
4. (en) Francis Bashforth et John Couch Adams, An attempt to test the theories of capillary action by comparing the theoretical and measured forms of drops of fluid, Cambridge University Press, 1883 (lire en ligne)
5. Herman H. Goldstine, A history of numerical analysis from the 16th through the 19th century, Springer, coll. « Studies in the history of mathematics and physical sciences », 1977 (ISBN 978-0-387-90277-7 et 978-3-540-90277-5)
6. Roger Sherman Hoar et F. R. Moulton, « New Methods in Exterior Ballistics. », The American Mathematical Monthly, vol. 34, n^o 6,‎ juin 1927, p. 325 (DOI 10.2307/2299335, lire en ligne, consulté le 8 octobre 2024)
7. Endre Süli et David Francis Mayers, An introduction to numerical analysis, Cambridge University Press, 2003 (ISBN 978-0-521-00794-8 et 978-0-521-81026-5), p.323-340
8. Germund Dahlquist, « Convergence and stability in the numerical integration of ordinary differential equations », MATHEMATICA SCANDINAVICA, vol. 4,‎ 1^er décembre 1956, p. 33 (ISSN 1903-1807 et 0025-5521, DOI 10.7146/math.scand.a-10454, lire en ligne, consulté le 8 octobre 2024)
9. (en) Germund G. Dahlquist, « A special stability problem for linear multistep methods », BIT, vol. 3, n^o 1,‎ mars 1963, p. 27–43 (ISSN 0006-3835 et 1572-9125, DOI 10.1007/BF01963532, lire en ligne, consulté le 8 octobre 2024)