Résoudre les EDO numériquement : pas, précision et stabilité

Une trajectoire lisse à l’écran peut résulter d’une succession de pas numériques finis. Le pas influe sur l’erreur à chaque étape, mais aussi sur le risque de transformer une décroissance en croissance ou de faire gagner de l’énergie à un oscillateur. Les notes sur le double pendule et l’attracteur de Rössler présentent leurs modèles et leur intégration par RK4 ; le laboratoire orbital utilise velocity Verlet. Pour juger ces méthodes, on peut commencer par une équation de décroissance et un oscillateur harmonique, dont les solutions exactes sont connues.

Ce qu’un pas approxime

Un problème à valeur initiale donne un taux de variation et un état de départ :

$$y'(t)=f(t,y(t)),\qquad y(t_0)=y_0.$$

L’état $y$ peut être un scalaire ou un vecteur regroupant, par exemple, position et vitesse. Avec un pas $h>0$, on calcule des approximations $y_n\approx y(t_n)$ aux instants $t_n=t_0+nh$. Un pas exact vérifie

$$y(t_{n+1})=y(t_n)+\int_{t_n}^{t_{n+1}}f(t,y(t))\,dt.$$

La difficulté tient à la présence de $y(t)$ dans l’intégrale : cet état est encore inconnu. Une méthode numérique approxime l’incrément à partir d’états connus et de plusieurs évaluations de la pente. Interpoler les résultats pour l’affichage rend la courbe plus lisse sans améliorer les états calculés. L’erreur étudiée ici provient surtout de la discrétisation du temps ; l’arrondi est une autre question, traitée dans Nombres à virgule flottante et calcul stable.

Euler explicite : erreur locale et erreur globale

La section 6.2 de Fundamentals of Numerical Computation, de Driscoll et Braun, donne la méthode d’Euler explicite : conserver la pente du point de départ pendant tout le pas.

$$y_{n+1}=y_n+h f(t_n,y_n).$$

Prenons une équation dont on peut vérifier directement la solution :

$$y'=-2y,\qquad y(0)=1,\qquad y(t)=e^{-2t}.$$

Pour $h=0.1$, la récurrence devient $y_{n+1}=0.8y_n$, d’où $y_n=0.8^n$.

Pas $n$	Instant $t_n$	Approximation d’Euler $y_n$	Valeur exacte $e^{-2t_n}$
0	0.0	1.000000000	1.000000000
1	0.1	0.800000000	0.818730753
2	0.2	0.640000000	0.670320046
10	1.0	0.107374182	0.135335283

L’erreur locale sur un pas mesure l’erreur commise en partant de l’état exact. Pour une solution suffisamment régulière, le développement de Taylor donne

$$y(t_n+h)-\bigl[y(t_n)+hf(t_n,y(t_n))\bigr] =\frac{h^2}{2}y''(t_n)+O(h^3).$$

L’erreur sur un pas d’Euler est donc $O(h^2)$. À l’état initial, $y''(0)=4$, et le terme dominant vaut $2h^2$. Les erreurs absolues effectives sont $0.018730753$ pour $h=0.1$ et $0.004837418$ pour $h=0.05$ : leur réduction se rapproche d’un facteur quatre. Le manuel appelle erreur de troncature locale cette erreur divisée par $h$ ; la quantité ainsi définie est $O(h)$. Il faut vérifier cette convention avant de comparer les ordres.

L’erreur globale compare la valeur issue des itérations, $y_n$, à $y(t_n)$ et comprend la propagation des erreurs antérieures. Le théorème de convergence du manuel suppose des bornes uniformes pour l’erreur locale par unité de pas et pour la dérivée de la fonction d’incrément par rapport à l’état. Sur un intervalle fini fixé, avec un second membre et une solution suffisamment réguliers et sous ces conditions, l’erreur globale d’Euler est $O(h)$. Atteindre le même instant final demande environ $1/h$ pas ; les erreurs précédentes modifient aussi les pentes suivantes. L’erreur locale du dernier pas ne représente donc pas l’ensemble du calcul.

À $T=1$, les erreurs globales absolues pour $h=0.1$ et $0.05$ valent environ $0.027961101$ et $0.013758629$, soit un rapport de $2.032$. Pour une méthode du premier ordre, diviser le pas par deux tend à diviser l’erreur globale par deux.

RK4 classique : davantage de pentes pour gagner en précision

La méthode classique de Runge–Kutta d’ordre quatre, au §6.4, évalue quatre pentes par pas. Ici, les $k_i$ désignent des pentes, sans facteur $h$ :

$$\begin{aligned} k_1&=f(t_n,y_n),\\ k_2&=f(t_n+h/2,y_n+hk_1/2),\\ k_3&=f(t_n+h/2,y_n+hk_2/2),\\ k_4&=f(t_n+h,y_n+hk_3),\\ y_{n+1}&=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4). \end{aligned}$$

Sur un intervalle fini fixé, si le second membre et la solution sont suffisamment réguliers et si les conditions de propagation de l’erreur précédentes sont remplies, l’erreur sur un pas est $O(h^5)$ et l’erreur globale $O(h^4)$. Dans le régime asymptotique de convergence, tant que l’arrondi ne domine pas, diviser le pas par deux devrait ramener l’erreur globale à environ $1/16$ de sa valeur précédente.

Pour $y'=-2y$, substituer les quatre pentes donne le multiplicateur d’un pas :

$$R(z)=1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24},\qquad z=-2h.$$

Il coïncide avec les cinq premiers termes du développement de $e^z$. On peut donc vérifier indépendamment l’approximation à $T=1$ par $R(-2h)^{1/h}$. Le programme ci-dessous donne des erreurs absolues finales de $4.265\times10^{-6}$, $2.452\times10^{-7}$ et $1.470\times10^{-8}$ pour $h=0.1,0.05,0.025$. Les rapports successifs, $17.396$ et $16.682$, se rapprochent de 16. Une méthode d’ordre élevé reste soumise à ses conditions de stabilité.

Stabilité absolue : quand la décroissance devient croissance

La section 11.3 définit la stabilité absolue à l’aide de $y'=\lambda y$ : à pas $h$ fixé, la solution numérique reste-t-elle bornée lorsque le nombre de pas tend vers l’infini ? En posant $z=h\lambda$, Euler explicite donne

$$y_{n+1}=(1+z)y_n,\qquad y_n=(1+z)^n y_0.$$

Sa région de stabilité est $|1+z|\le1$, le disque fermé de rayon 1 centré en $-1$ dans le plan complexe. Pour obtenir une décroissance vers zéro, il faut l’inégalité stricte $|1+z|<1$. Pour $\lambda<0$ réel, elle devient

$$0<h<\frac{2}{|\lambda|}.$$

Pour $\lambda=-2$, la frontière est $h=1$. Le multiplicateur y vaut $-1$ : une solution non nulle change de signe à chaque pas sans décroître. Au-delà, son amplitude augmente. Avec $h=1.1$ et $y_0=1$, on obtient

$$1,\;-1.2,\;1.44,\;-1.728,\;2.0736.$$

La solution exacte reste positive et décroît. Même dans la région de stabilité, Euler alterne les signes pour $0.5<h<1$ : la stabilité ne garantit ni la précision quantitative ni la positivité. Pour $\lambda$ complexe, il faut vérifier $|1+h\lambda|$ et non remplacer $\lambda$ par son module dans la restriction réelle.

La région de RK4 sur cette équation test est $|R(z)|\le1$ et possède elle aussi une frontière. Stabilité et ordre de l’erreur globale décrivent des limites différentes : converger à un instant final fixé quand $h\to0$ ne permet pas de prolonger indéfiniment une simulation à pas fixé.

Raideur et Euler implicite

La section 11.4 explique la raideur par la coexistence de plusieurs échelles de temps : des modes à décroissance rapide imposent un petit pas aux méthodes explicites, alors que l’évolution lente qui nous intéresse demande une longue observation. Par exemple,

$$x'=-x,\qquad w'=-1000w,\qquad x(0)=w(0)=1.$$

Les temps caractéristiques de décroissance sont 1 et $0.001$. La composante lente continue d’évoluer lorsque la composante rapide est déjà minuscule ; Euler explicite exige pourtant encore $h<0.002$ pour amortir ce mode rapide. Les erreurs d’arrondi et de discrétisation peuvent le réexciter. La raideur dépend de l’équation, de la méthode et de la précision recherchée ; l’allure plus ou moins abrupte d’une courbe ne suffit pas à l’identifier.

Euler implicite évalue la pente à l’état suivant, encore inconnu :

$$y_{n+1}=y_n+h f(t_{n+1},y_{n+1}).$$

Pour l’équation test, on obtient en réarrangeant

$$y_{n+1}=\frac{y_n}{1-h\lambda}.$$

Le multiplicateur donne la région de stabilité $|1-z|\ge1$, qui contient tout le demi-plan gauche : la méthode est A-stable. Pour $\lambda<0$ réel, tout $h>0$ amortit le mode. Avec $h=0.1$, le multiplicateur de la composante rapide vaut par exemple $1/101$. On peut ainsi franchir un transitoire rapide ; pour le résoudre en détail, il faut toujours réduire le pas. L’erreur globale d’Euler implicite reste du premier ordre. La stabilité n’augmente pas l’ordre de précision.

La section 6.7 du manuel, consacrée à l’implémentation des méthodes implicites, explique pourquoi l’état suivant, encore inconnu, demande généralement une recherche de racine. Pour une équation non linéaire, chaque pas d’Euler implicite exige de résoudre $G(u)=0$, avec

$$G(u)=u-y_n-hf(t_{n+1},u).$$

Une correction de Newton résout

$$\bigl[I-hJ_f(t_{n+1},u^{(k)})\bigr]\Delta u=-G(u^{(k)}), \qquad u^{(k+1)}=u^{(k)}+\Delta u.$$

$J_f$ est la matrice jacobienne de $f$ par rapport à l’état, et $I$ la matrice identité. Calculer ou approximer la jacobienne, résoudre des systèmes linéaires et répéter les itérations a un coût. Il faut vérifier la convergence et maintenir l’erreur de résolution sous le budget d’erreur de discrétisation temporelle ; une résolution implicite qui échoue ne constitue pas un pas acceptable.

Énergie de l’oscillateur à long terme : Euler et méthodes symplectiques

Prenons un oscillateur harmonique de masse et de pulsation unitaires :

$$q'=p,\qquad p'=-q,\qquad (q(0),p(0))=(1,0).$$

Sa solution exacte est $q(t)=\cos t$, $p(t)=-\sin t$, et son énergie $E=(q^2+p^2)/2=1/2$ est constante. Euler explicite met à jour les deux variables à partir de leurs anciennes valeurs :

$$q_{n+1}=q_n+hp_n,\qquad p_{n+1}=p_n-hq_n.$$

En additionnant les deux carrés, les termes croisés s’annulent et l’on obtient

$$E_{n+1}=(1+h^2)E_n.$$

Tout pas positif fixé fait augmenter l’énergie continuellement. Pour $h=0.1$ et $T=100$, il y a 1000 pas, et $E=\tfrac12(1.01)^{1000}\approx10479.577819$. Euler implicite produit l’effet inverse : il divise l’énergie par $1+h^2$ à chaque pas, pour atteindre environ $0.000024$ au même instant final. La trajectoire reste bornée, mais subit un amortissement artificiel.

La deuxième leçon du cours Geometric Numerical Integration de Hairer donne Euler symplectique. Pour ce système hamiltonien séparable, on met d’abord à jour la quantité de mouvement, puis la position avec la nouvelle quantité de mouvement :

$$p_{n+1}=p_n-hq_n,\qquad q_{n+1}=q_n+hp_{n+1}.$$

« Symplectique » désigne la préservation de la structure géométrique du flot hamiltonien. Dans cet exemple à deux dimensions, la matrice de mise à jour a un déterminant égal à 1 et conserve l’aire dans l’espace des phases. La méthode reste du premier ordre et ne conserve pas l’énergie d’origine à chaque pas. Pour cet oscillateur, la substitution dans les formules montre qu’elle conserve exactement une énergie modifiée :

$$H_h=\frac12(q^2+p^2-hqp).$$

Pour $0<h<2$ fixé, cette forme quadratique est définie positive. Avec $|qp|\le E$ et la valeur initiale $H_h=1/2$, on obtient

$$\frac{1/2}{1+h/2}\le E_n\le\frac{1/2}{1-h/2}.$$

Pour $h=0.1$, l’énergie oscille donc entre des bornes d’environ $0.476190$ et $0.526316$. Ces bornes sont propres à cet oscillateur linéaire et ne se transposent pas directement à tous les systèmes hamiltoniens. L’invariant exige aussi un pas constant ; faire varier $h$ invalide cet argument. Une énergie bornée peut s’accompagner d’une erreur de phase qui s’accumule.

Le théorème 2 de cette même leçon donne la méthode velocity Verlet, d’ordre deux : un demi-pas sur la quantité de mouvement, un pas entier sur la position, puis un autre demi-pas sur la quantité de mouvement. Pour $q''=a(q)$,

$$\begin{aligned} p_{n+1/2}&=p_n+\frac h2a(q_n),\\ q_{n+1}&=q_n+hp_{n+1/2},\\ p_{n+1}&=p_{n+1/2}+\frac h2a(q_{n+1}). \end{aligned}$$

Elle est également symplectique pour les systèmes hamiltoniens séparables ; le laboratoire orbital l’utilise pour la gravitation centrale. L’exemple ci-dessous prend $a(q)=-q$ pour comparer directement l’énergie à long terme et l’erreur sur l’état final.

Une comparaison exécutable

Enregistrez ce code dans ode_compare.py, puis lancez python3 ode_compare.py. Il utilise uniquement la bibliothèque standard de Python. Les fonctions Euler et RK4 acceptent des séquences d’état de toute longueur. f(t, y) doit renvoyer le même nombre de composantes de dérivée ; sinon, le pas lève une exception ValueError. Les deux fonctions implicites sont les mises à jour résolues explicitement pour ces exemples linéaires ; les fonctions Euler symplectique et Verlet sont propres à cet oscillateur.

from math import cos, exp, hypot, sin


def slope(f, t, y):
    values = tuple(f(t, y))
    if len(values) != len(y):
        raise ValueError("f(t, y) must have the same length as y")
    return values


def euler(f, t, y, h):
    return tuple(a + h * b for a, b in zip(y, slope(f, t, y)))


def rk4(f, t, y, h):
    def shift(k, scale):
        return tuple(a + scale * b for a, b in zip(y, k))
    k1 = slope(f, t, y)
    k2 = slope(f, t + h / 2, shift(k1, h / 2))
    k3 = slope(f, t + h / 2, shift(k2, h / 2))
    k4 = slope(f, t + h, shift(k3, h))
    return tuple(a + h * (b + 2 * c + 2 * d + e) / 6
                 for a, b, c, d, e in zip(y, k1, k2, k3, k4))


def decay(t, y):
    return (-2 * y[0],)


def decay_implicit(f, t, y, h):
    return (y[0] / (1 + 2 * h),)


def oscillator(t, y):
    q, p = y
    return (p, -q)


def oscillator_implicit(f, t, y, h):
    q, p = y
    return ((q + h * p) / (1 + h * h),
            (p - h * q) / (1 + h * h))


def symplectic_euler(f, t, y, h):
    q, p = y
    p = p - h * q
    return (q + h * p, p)


def verlet(f, t, y, h):
    q, p = y
    half_p = p - h * q / 2
    q = q + h * half_p
    return (q, half_p - h * q / 2)


print("local Euler: h, one-step error")
for h in (0.1, 0.05):
    print(f"{h:.3f} {abs(1 - 2 * h - exp(-2 * h)):.9f}")

print("decay at T=1: method, h, y, absolute error")
for name, step in (("Euler", euler), ("implicit", decay_implicit), ("RK4", rk4)):
    errors = []
    for n in (10, 20, 40):
        h, y = 1 / n, (1.0,)
        for j in range(n):
            y = step(decay, j * h, y, h)
        error = abs(y[0] - exp(-2))
        errors.append(error)
        print(f"{name:8s} {h:.3f} {y[0]:.9f} {error:.3e}")
    print("error ratios: " + " ".join(f"{a / b:.3f}"
          for a, b in zip(errors, errors[1:])))

print("Euler with h=1.1: y0 through y4")
y, values = (1.0,), [1.0]
for j in range(4):
    y = euler(decay, j * 1.1, y, 1.1)
    values.append(y[0])
print(" ".join(f"{value:.4f}" for value in values))

print("oscillator at T=100, h=0.1: method, E_min, E_max, E_final, state_error")
for name, step in (("Euler", euler), ("implicit", oscillator_implicit),
                   ("RK4", rk4), ("symplectic", symplectic_euler), ("Verlet", verlet)):
    y, energies = (1.0, 0.0), [0.5]
    for j in range(1000):
        y = step(oscillator, j * 0.1, y, 0.1)
        energies.append((y[0] ** 2 + y[1] ** 2) / 2)
    error = hypot(y[0] - cos(100), y[1] + sin(100))
    print(f"{name:10s} {min(energies):.6f} {max(energies):.6f} "
          f"{energies[-1]:.6f} {error:.3e}")

Sortie :

local Euler: h, one-step error
0.100 0.018730753
0.050 0.004837418
decay at T=1: method, h, y, absolute error
Euler    0.100 0.107374182 2.796e-02
Euler    0.050 0.121576655 1.376e-02
Euler    0.025 0.128512157 6.823e-03
error ratios: 2.032 2.016
implicit 0.100 0.161505583 2.617e-02
implicit 0.050 0.148643628 1.331e-02
implicit 0.025 0.142045682 6.710e-03
error ratios: 1.966 1.983
RK4      0.100 0.135339548 4.265e-06
RK4      0.050 0.135335528 2.452e-07
RK4      0.025 0.135335298 1.470e-08
error ratios: 17.396 16.682
Euler with h=1.1: y0 through y4
1.0000 -1.2000 1.4400 -1.7280 2.0736
oscillator at T=100, h=0.1: method, E_min, E_max, E_final, state_error
Euler      0.500000 10479.577819 10479.577819 1.438e+02
implicit   0.000024 0.500000 0.000024 9.935e-01
RK4        0.499993 0.500000 0.499993 8.332e-05
symplectic 0.476190 0.526316 0.521321 5.665e-02
Verlet     0.498750 0.500000 0.499724 4.222e-02

Les minima et maxima d’énergie comprennent l’état initial et tous les points de la grille. state_error est la distance euclidienne entre l’état final et $(\cos100,-\sin100)$. Sur cet intervalle fini, RK4 présente une faible erreur ; les méthodes symplectiques ont une variation d’énergie bornée, mais leur état diffère encore de la solution exacte. Pour choisir une méthode, examinez séparément la précision finale, la phase et les invariants à long terme, plutôt que de classer les méthodes par leur seul ordre.

Pas adaptatif et vérification de la solution

La section 6.5 présente l’estimation d’erreur par formules emboîtées : deux approximations d’ordres différents sont calculées sur le même pas, et leur différence estime l’erreur locale de la formule d’ordre inférieur. Si l’estimation dépasse la cible, le pas est rejeté puis retenté avec un $h$ plus petit ; il n’est accepté que si l’estimation est assez faible. Si l’estimation $\eta$ varie approximativement comme $h^{r+1}$, on peut ajuster le pas par

$$h_{\mathrm{new}}=s h\left(\frac{\mathrm{tol}}{\eta}\right)^{1/(r+1)},\qquad 0<s<1,$$

tout en limitant les facteurs d’augmentation et de réduction. Ici, $r$ est l’ordre de la formule dont on estime l’erreur et $s$ un facteur de sécurité. L’exposant doit correspondre à l’ordre de l’estimateur : il s’agit d’un modèle de contrôle local, pas d’une règle universelle pour tous les solveurs.

La documentation SciPy de solve_ivp décrit RK45 comme une méthode emboîtée 5(4) : la formule d’ordre quatre estime l’erreur, et celle d’ordre cinq fait avancer la solution. Elle diffère du RK4 classique ci-dessus. L’échelle de l’erreur locale est atol + rtol * abs(y) ; choisissez un atol adapté près de zéro, avec des valeurs distinctes pour les composantes d’échelles différentes. RK45 convient aux problèmes non raides, Radau ou BDF aux problèmes raides. Ces paramètres de contrôle local ne garantissent pas l’erreur à l’instant final.

Pour vérifier une solution, gardez la même cible de comparaison :

1. Diviser le pas par deux au même instant final. Pour une méthode d’ordre $r$, calculez $y_h,y_{h/2},y_{h/4}$ dans le régime asymptotique. Le rapport $\|y_h-y_{h/2}\|/\|y_{h/2}-y_{h/4}\|$ devrait tendre vers $2^r$ ; on estime l’erreur de $y_{h/2}$ par $\|y_{h/2}-y_h\|/(2^r-1)$. Pour la grille la plus fine, l’erreur de $y_{h/4}$ s’estime par $\|y_{h/4}-y_{h/2}\|/(2^r-1)$. Il faut la stabilité, une régularité suffisante, un terme d’erreur dominant non nul et un arrondi qui ne domine pas encore. Pour des trajectoires entières, comparez les mêmes instants.
2. Vérifier les invariants du modèle. Pour cet oscillateur, examinez l’énergie avec la phase ou l’erreur sur l’état ; pour la gravitation centrale, l’énergie et le moment cinétique. Des invariants proches ne prouvent pas la précision de la trajectoire, et un modèle dissipatif ne doit pas être contraint à conserver l’énergie.
3. Utiliser une référence indépendante. Préférez une solution exacte lorsqu’elle existe. Sinon, choisissez un solveur adapté à l’équation, resserrez ses tolérances, puis resserrez-les encore pour vérifier que la référence se stabilise. Pour comparer solve_ivp à un programme à pas fixé, t_eval choisit les instants de sauvegarde des résultats, pas les pas internes. Conservez les mêmes équations, paramètres, conditions initiales et instants de comparaison.

Le chaos amplifie les petites différences d’état. Pour le double pendule et l’attracteur de Rössler, vérifiez d’abord la convergence en pas sur un court intervalle fini, puis comparez les invariants pertinents ou les caractéristiques statistiques à plus long terme. Une trajectoire qui finit par diverger ne permet pas, à elle seule, de distinguer la sensibilité aux conditions initiales de l’erreur de l’intégrateur.

Mes lectures · Les notes ouvertes sont mémorisées dans ce navigateur.