Équations différentielles: méthodes à un pas

def euler(f, t0, y0, h, n):
    """n pas d'Euler explicite; renvoie la liste des couples (t_k, y_k)."""
    points = [(t0, y0)]
    t, y = t0, y0
    for k in range(n):
        y = y + h * f(t, y)
        t = t0 + (k + 1) * h
        points.append((t, y))
    return points

Exemple 9.1 · Quatre pas sur le problème témoin

Prenons (9.3) avec $h = 0, 5$ , donc $n = 4$ pas sur $[0, 2]$ . Le second membre vaut $f (t, y) = y - t^{2} + 1$ , et la récurrence (9.6) s'écrit $y_{k + 1} = y_{k} + 0, 5 (y_{k} - t_{k}^{2} + 1)$ .

$k$	$t_{k}$	$f (t_{k}, y_{k})$	$y_{k}$ (calculé)

Deux observations, toutes deux importantes. D'abord, l'erreur croît de pas en pas et ne fait que croître: elle passe de $0, 176$ à $0, 868$ , c'est-à-dire qu'elle est multipliée par près de 5 en quatre pas. Ensuite, elle ne croît pas linéairement: les accroissements successifs valent $0, 176$ , $0, 215$ , $0, 243$ et $0, 234$ . L'erreur du pas $k$ n'est donc pas simplement la somme des erreurs commises à chaque pas — chaque erreur passée est elle-même transportée et amplifiée par les pas suivants. C'est exactement ce que le théorème 9.3 quantifiera, et c'est là que le facteur $e^{L (b - a)}$ du théorème 9.1 va revenir.

$h$	$n$	$y_{n}$ en $t = 2$	erreur	rapport
0,5	4	4,437500	0,86797	—
0,25	8	4,779652	0,52582	1,65
0,1	20	5,063500	0,24197	—
0,05	40	5,178006	0,12747	1,90
0,025	80	5,239977	0,06550	1,95
0,0125	160	5,272264	0,03321	1,97

Exemple 9.2 · Euler implicite sur le problème témoin, et une surprise

Le second membre $f (t, y) = y - t^{2} + 1$ est affine en $y$ , donc l'équation (9.9) se résout exactement, sans Newton:

y_{n + 1} = y_{n} + h (y_{n + 1} - t_{n + 1}^{2} + 1) ⟺ y_{n + 1} = \frac{y _{n} + h ( 1 - t _{n + 1}^{2} )}{1 - h},

formule licite tant que $h \neq = 1$ . Les résultats, mis en regard de ceux d'Euler explicite:

$h$	Euler implicite	erreur	Euler explicite	erreur
0,5	8,500000	3,19453	4,437500	0,86797
0,25	6,252820	0,94735	4,779652	0,52582
0,1	5,609895	0,30442	5,063500	0,24197
0,05	5,448386	0,14291	5,178006	0,12747

Sur ce problème, la méthode implicite est nettement moins précise que l'explicite, tout en coûtant plus cher. Il faut le dire, et l'expliquer plutôt que de le cacher. Sur le modèle linéaire $y^{'} = λ y$ , un pas explicite multiplie par $1 + hλ$ et un pas implicite par $1/ (1 - hλ)$ , alors que la solution exacte est multipliée par $e^{hλ}$ . Ici $λ = 1$ et : on compare et à . L'explicite sous-estime de , l'implicite surestime de — deux fois et demie plus. Les deux méthodes sont bien d'ordre 1 — mais il faut descendre assez bas pour le voir: l'erreur de l'implicite passe de à puis et pour ; ; et , soit des rapports de ; et qui tendent vers 2 . Entre et le rapport vaut encore : à ce pas-là, la méthode implicite est loin de son régime asymptotique, et citer son ordre n'apprendrait rien sur son erreur.

Alors pourquoi la méthode implicite existe-t-elle? Parce que la comparaison ci-dessus suppose que l'on choisit $h$ pour la précision. Il existe des problèmes où le pas est imposé par une tout autre contrainte, et où l'explicite ne produit alors plus rien d'utilisable du tout. C'est le sujet du chapitre 10, et l'exemple qui clôt ce chapitre en donnera un avant-goût mesuré.

Exemple 9.3 · Trois méthodes d'ordre 2 sur le problème témoin

Le point milieu est obtenu de la même façon, en approchant l'intégrale de (9.8) par $h g (t_{n} + h /2)$ et en prédisant la valeur au milieu par un demi-pas d'Euler:

y_{n + 1} = y_{n} + h f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} f (t_{n}, y_{n})) . (9.12)

Voici les erreurs en $t = 2$ , chacune avec le rapport que le code a imprimé:

$h$	point milieu	rapport	Heun	rapport
0,5	$9, 057 \cdot 1 0^{- 2}$	—	$3, 892 \cdot 1 0^{- 1}$	—
0,25	$2, 357 \cdot 1 0^{- 2}$

$h$	trapèze (implicite)	rapport	RK4	rapport
0,5	$1, 635 \cdot 1 0^{- 1}$	—	$3, 867 \cdot 1 0^{- 3}$	—
0,25	$3, 905 \cdot 1 0^{- 2}$

Les trois méthodes d'ordre 2 ont bien un rapport qui tend vers $4 = 2^{2}$ , et RK4 un rapport qui tend vers $16 = 2^{4}$ — par en dessous pour les unes, par au-dessus pour le trapèze, ce qui rappelle que le signe du terme correctif dépend du signe de $C_{2} / C_{1}$ et n'a pas de règle générale. Les , elles, diffèrent d'un facteur 4 entre le point milieu et Heun, alors que les deux ont exactement le même coût de deux évaluations par pas: à ordre égal, le choix de la méthode se joue sur la constante, et c'est pourquoi les tables de méthodes de Runge-Kutta sont si longues.

def rk4(f, t0, y0, h, n):
    """n pas de Runge-Kutta d'ordre 4; renvoie la valeur finale."""
    t, y = t0, y0
    for k in range(n):
        k1 = f(t, y)
        k2 = f(t + h / 2.0, y + h / 2.0 * k1)
        k3 = f(t + h / 2.0, y + h / 2.0 * k2)
        k4 = f(t + h, y + h * k3)
        y = y + h / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
        t = t0 + (k + 1) * h
    return y

Exemple 9.5 · Le même budget, trois méthodes

Comparons à coût égal, c'est-à-dire à nombre égal d'évaluations du second membre, qui est la seule mesure de coût honnête ici (le chapitre 8 comptait des opérations; ici, une évaluation de $f$ peut coûter mille opérations dans un vrai modèle).

méthode	$h$	évaluations de $f$	erreur en $t = 2$
Euler	0,025	80	$6, 550 \cdot 1 0^{- 2}$
Heun	0,0625	64	$7, 495 \cdot 1 0^{- 3}$
RK4	0,5	16	$3, 867 \cdot 1 0^{- 3}$
RK4	0,25	32	$2, 622 \cdot 1 0^{- 4}$

RK4 au pas $h = 0, 5$ , avec ses quatre pas et ses 16 évaluations, est plus précis qu'Euler au pas $h = 0, 025$ , qui en demande 80. En cherchant le pas qui donnerait à Euler l'erreur de RK4 à $h = 0, 5$ , on trouve $h \approx 0, 00144$ , soit environ 1 400 évaluations contre 16: un facteur 87. Et l'écart grandit avec l'exigence, puisque le pas requis varie comme $tol$ pour Euler et comme $tol^{1/4}$ pour RK4.

C'est l'argument décisif en faveur des méthodes d'ordre élevé, et il ne devient visible qu'en comptant les évaluations plutôt que les pas.

Exemple 9.6 · La borne théorique, confrontée aux erreurs mesurées

Sur le problème témoin, $L = 1$ , $M = ∣ 2 - \frac{1}{2} e^{2} ∣ = 1, 694528$ et $b - a = 2$ . La borne (9.22) en $t = 2$ s'écrit donc

∣ e_{N} ∣ \leq \frac{1 , 694 528}{2} (e^{2} - 1) h = 5, 4132 h .

$h$	borne (9.22)	erreur mesurée	borne / mesure
0,5	2,7066	0,86797	3,12
0,25	1,3533	0,52582	2,57
0,1	0,5413	0,24197	2,24
0,05	0,2707	0,12747	2,12
0,025	0,1353	0,06550	2,07

La borne est valide dans tous les cas — elle ne serait pas un théorème sinon — et elle est pessimiste d'un facteur 2 à 3, qui se stabilise autour de 2. C'est le comportement typique: la majoration a remplacé $y^{''} (ξ_{n})$ par son maximum sur tout l'intervalle alors que $y^{''}$ change de signe en $t = ln 4 = 1, 386$ , et elle a supposé que toutes les erreurs se cumulaient au pire. Retenez les deux lectures: la borne donne l' de façon fiable et gratuite, et elle donne la de façon très grossière. Pour la taille, on mesure.

Exemple 9.7 · L'estimation contre l'erreur vraie

Premier pas de (9.3), depuis $(0; 0, 5)$ , pour plusieurs pas $h$ . On compare l'estimation $est = \frac{h}{2} ∣ k_{2} - k_{1} ∣$ à l'erreur vraie d'Euler, $∣ y_{0} + h k_{1} - y (h)∣$ , que l'on connaît ici puisqu'on a la solution exacte.

$h$	estimation	erreur vraie d'Euler	estimation / vraie
0,5	0,125000	0,175639	0,712
0,25	0,039063	0,045487	0,859
0,1	0,007000	0,007415	0,944
0,05	0,001813	0,001864	0,972
0,025	0,000461	0,000467	0,986

Le rapport tend vers 1, comme (9.24) le promet — mais seulement asymptotiquement: à $h = 0, 5$ , l'estimation sous-estime l'erreur vraie de $29 %$ . Un estimateur d'erreur n'est fiable que dans le régime où il est asymptotique, c'est-à-dire précisément là où l'erreur est déjà petite. C'est pourquoi tout code sérieux applique un coefficient de sécurité et refuse, comme ci-dessous, de faire croître le pas de plus d'un facteur fixé.

Exemple 9.9 · L'oscillateur amorti d'Analyse II, intégré numériquement

Reprenons exactement la suspension d'instrument de l'exemple 1.13 d'Analyse II: masse $m = 2$ kg, raideur $k = 200$ N/m, amortissement $c = 8$ N·s/m, écartée de $x_{0} = 10$ mm et lâchée sans vitesse. L'équation est

\overset{x}{¨} + 2 ζ ω_{0} \overset{x}{˙} + ω_{0}^{2} x = 0, ω_{0} = k / m = 10 rad/s, ζ = \frac{c}{2 k m} = 0, 2,

en régime pseudo-périodique, de pseudo-pulsation $ω_{d} = ω_{0} 1 - ζ^{2} = 9, 797959$ rad/s. La réduction (9.27) avec donne

u^{'} = f (t, u) = (v - 2 ζ ω_{0} v - ω_{0}^{2} x), u (0) = (0, 010 0) .

from math import sqrt

M, K, C = 2.0, 200.0, 8.0
W0 = sqrt(K / M)
ZETA = C / (2.0 * sqrt(K * M))

def oscillateur(t, U):
    return [U[1], -2.0 * ZETA * W0 * U[1] - W0

Analyse II donne la solution exacte $x (t) = x_{0} e^{- ζ ω_{0} t} (cos ω_{d} t + \frac{ζ ω _{0}}{ω _{d}} sin ω_{d} t)$ , d'où mm. Intégrons de à s:

méthode	$h$	$x (1)$ calculé, en mm	erreur, en m
Euler	0,005	−1,686382	$3, 25 \cdot 1 0^{- 4}$
RK4	0,01	−1,360931	$1, 04 \cdot 1 0^{- 8}$
RK4	0,005	−1,360921	$6, 58 \cdot 1 0^{- 10}$

Les rapports des trois erreurs de RK4 valent $15, 81$ puis $15, 92$ : l'ordre 4 se lit sur un système exactement comme sur une équation scalaire. Euler, lui, se trompe de $24 %$ sur le déplacement avec un pas pourtant deux fois plus fin que celui de la première ligne de RK4.

Et une remarque qui appartient au chapitre suivant. Si l'on porte le pas d'Euler à $h = 0, 05$ — dix fois plus grand, ce qui paraît raisonnable devant la pseudo-période $T_{d} = 0, 641$ s — la solution calculée ne se contente pas d'être imprécise: son amplitude croît au lieu de décroître, et atteint $95$ mm à $t = 5$ s alors que la solution exacte est tombée sous $1 0^{- 3}$ mm. Le seuil est net: le calcul reste borné pour et diverge au-delà. Ce n'est plus une question de précision mais de , et cette distinction est le sujet entier du chapitre 10.

Le problème de Cauchy $y^{'} = f (t, y)$ , $y (a) = y_{0}$ possède une solution unique sur $[a, b]$ dès que $f$ est continue et lipschitzienne en $y$ (théorème 9.1, admis: sa démonstration est un point fixe dans un espace de fonctions). Deux solutions issues de données distantes de $δ$ restent distantes d'au plus $e^{L (b - a)} δ$ — le même facteur qui bornera l'erreur numérique.
Sur la grille $t_{k} = a + k h$ , une méthode à un pas produit $y_{k + 1} = y_{k} + h Φ (t_{k}, y_{k}, h)$ . Toutes ces méthodes s'obtiennent en appliquant une formule de quadrature du chapitre 8 à l'identité exacte — rectangle à gauche pour Euler explicite, à droite pour Euler implicite, trapèze pour Crank-Nicolson, Simpson pour RK4.
Une méthode implicite demande la résolution d'une équation non linéaire par pas, donc une itération de Newton (chapitre 2) et trois à dix fois le coût d'un pas explicite. Sur le problème témoin, Euler implicite est même moins précis qu'Euler explicite; son intérêt est ailleurs, et c'est le chapitre 10 qui le dira.
Une méthode de Runge-Kutta à $s$ étages échantillonne la pente $s$ fois par pas, avec les poids du tableau de Butcher. L'ordre égale le nombre d'étages jusqu'à $s = 4$ et pas au-delà: il faut 6 étages pour l'ordre 5 et 11 pour l'ordre 8. À budget égal d'évaluations, RK4 écrase Euler — 16 évaluations de RK4 valent environ 1 400 évaluations d'Euler sur le problème témoin.
Consistance d'ordre $p$ plus stabilité entraîne convergence d'ordre $p$ (théorème 9.3), avec la borne $∣ e_{n} ∣ \leq e^{Λ (t_{n} - a)} ∣ e_{0} ∣ + \frac{C h ^{p}}{Λ} (e^{Λ (t_{n} - a)} - 1)$ . L'équation qui porte la démonstration est : l'erreur nouvelle est l'ancienne transportée, plus la fraîche. Pour Euler, et la borne est , vérifiée et pessimiste d'un facteur 2 environ.
Un ordre est asymptotique. Le rapport d'erreurs d'Euler vaut $1, 65$ , $1, 90$ , $1, 95$ , $1, 97$ et tend vers 2 sans l'atteindre, l'écart à 2 étant lui-même proportionnel à $h$ . Un rapport de $1, 90$ confirme l'ordre 1; l'arrondir à 2 serait effacer la seule information que la mesure apporte au théorème.
Le pas est contrôlé par une paire emboîtée: deux méthodes d'ordres $p$ et $p + 1$ partageant leurs étages, dont la différence estime l'erreur a posteriori et pilote $h_{nouv} = σ h (tol / est)^{1/ (p + 1)}$ . Un système ou une équation d'ordre $m$ se traite sans changer une formule, en posant .

Exercice 9.2 · Retrouver l'ordre à partir d'un tableau

Une méthode que l'on ne vous nomme pas, appliquée au problème témoin $y^{'} = y - t^{2} + 1$ , $y (0) = 0, 5$ , donne les erreurs suivantes en $t = 1$ . Ce sont de vraies mesures, produites en exécutant la méthode et en comparant à $y (1) = 2, 640859086$ .

$h$	0,2	0,1	0,05	0,025
erreur	$7, 692 \cdot 1 0^{- 3}$	$1, 981 \cdot 1 0^{- 3}$	$5, 017 \cdot 1 0^{- 4}$

Calculer la colonne des rapports et en déduire l'ordre.
Estimer l'ordre par la formule $p \approx lo g_{2} (E (h) / E (h /2))$ pour chaque couple, et commenter la tendance.
Estimer la constante $K$ telle que $E (h) \approx K h^{p}$ , et prédire l'erreur pour .

Solution

1. Les rapports valent

\frac{7 , 692 \cdot 1 0 ^{- 3}}{1 , 981 \cdot 1 0 ^{- 3}} = 3, 884, \frac{1 , 981 \cdot 1 0 ^{- 3}}{5 , 017 \cdot 1 0 ^{- 4}} = 3, 948, \frac{5 , 017 \cdot 1 0 ^{- 4}}{1 , 262 \cdot 1 0 ^{- 4}} = 3, 976.

0	0,0	1,500000	0,500000	0,500000	0
1	0,5	2,000000	1,250000	1,425639	0,175639
2	1,0	2,250000	2,250000	2,640859	0,390859
3	1,5	2,125000	3,375000	4,009155	0,634155
4	2,0	—	4,437500	5,305472	0,867972

$k$	$t_{k}$	$y_{k}$	$y (t_{k})$	$∣ e_{k} ∣$
0	0,00	1,000000	1,000000	0
1	0,25	0,750000	0,807602	0,057602
2	0,50	0,625000	0,713061	0,088061
3	0,75	0,593750	0,694733	0,100983
4	1,00	0,632813	0,735759	0,102946

Objectifs du chapitre

Le problème de Cauchy

Le problème témoin de ce chapitre

Discrétiser: la grille, et deux suites qu'il ne faut jamais confondre

La méthode d'Euler explicite, dérivée trois fois

Première dérivation: le quotient de différences

Deuxième dérivation: Taylor

Troisième dérivation: la quadrature de la forme intégrale

Euler explicite en Python

L'ordre d'Euler, lu dans un tableau

Euler implicite, et ce que coûte un pas

Le trapèze, Crank-Nicolson, et Heun comme prédicteur-correcteur

Runge-Kutta: échantillonner la pente plusieurs fois

Les méthodes à deux étages

RK4, la méthode classique

Le nombre d'étages et l'ordre cessent de coïncider

Erreur locale, erreur globale, et le théorème qui relie les deux

Le lemme d'accumulation

Consistance et stabilité entraînent la convergence

L'ordre 1 d'Euler explicite, en entier

Contrôler le pas: les paires emboîtées

Systèmes, et équations d'ordre supérieur

Synthèse

Exercices

Références

Euler explicite

Le point milieu

Heun, et l'estimation d'erreur gratuite

RK4, et le rapport des erreurs