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

Objectifs du chapitre

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

• énoncer le problème de Cauchy et la condition de Lipschitz qui en garantit l'existence et l'unicité, et distinguer en permanence la valeur calculée $y_n$ de la valeur exacte $y(t_n)$;
• construire la méthode d'Euler explicite de trois manières — différence finie, développement de Taylor, quadrature de la forme intégrale — et dire laquelle des trois se généralise;
• écrire et programmer Euler implicite, la méthode du trapèze, Heun et RK4, avec leur tableau de Butcher, et chiffrer ce que coûte un pas de chacune;
• définir l'erreur locale de troncature et l'erreur globale, démontrer que consistance et stabilité entraînent la convergence, et en déduire l'ordre 1 d'Euler;
• lire un ordre de convergence sur un tableau de rapports et sur un graphique log–log, et reconnaître qu'un rapport de $1,90$ est la signature d'un ordre 1, non d'un défaut;
• contrôler le pas par une paire emboîtée, et ramener un système ou une équation d'ordre supérieur — l'oscillateur amorti d'Analyse II, par exemple — à un problème du premier ordre vectoriel.

Le problème de Cauchy

Vous avez rencontré les équations différentielles ordinaires en Analyse II, chapitre 1, du point de vue de la résolution exacte: variables séparables, équations linéaires du premier ordre, équations linéaires du second ordre à coefficients constants. Ce chapitre commence là où celui-là s'arrête, c'est-à-dire presque tout de suite: la très grande majorité des équations différentielles issues de la physique, de la chimie ou de la biologie n'admettent aucune solution en termes de fonctions élémentaires. Pour celles-là — et elles sont la règle — on calcule.

Avant de calculer quoi que ce soit, il faut savoir qu'il y a quelque chose à calculer: une solution, et une seule. C'est l'objet du théorème suivant, dont l'hypothèse porte sur la régularité de $f$ par rapport à sa seconde variable.

Ce théorème est admis, et il faut dire pourquoi. Sa démonstration réécrit (9.1) sous la forme intégrale $y(t)=y_0+{\int }_a^{t}f(s,y(s))\mathrm{d}s$, puis construit la solution comme point fixe de l'opérateur qui à une fonction continue associe le membre de droite; la condition (9.2) fait de cet opérateur une contraction sur un espace de fonctions continues, et le théorème du point fixe de Banach conclut. L'outil est donc le même que celui du de ce cours — le théorème du point fixe — mais appliqué dans un espace de fonctions, ce qui relève d'un cours d'équations différentielles ou d'analyse fonctionnelle. La version locale de ce théorème est énoncée, avec la même justification, au .

Retenez plutôt ce que l'hypothèse fait. La seconde inégalité du théorème 9.1 est le fait décisif pour tout ce chapitre: deux solutions issues de données voisines peuvent s'écarter, mais au plus exponentiellement, avec l'exposant $L$. Le facteur ${\mathrm{e}}^{L(b-a)}$ va réapparaître, identique, dans la majoration de l'erreur numérique — et pour la même raison, car une méthode numérique ne fait rien d'autre que de redémarrer, à chaque pas, depuis une donnée légèrement fausse.

Le problème témoin de ce chapitre

Tout ce chapitre s'appuie sur un problème dont on connaît la solution exacte — c'est la seule manière d'exhiber des erreurs vraies plutôt que des estimations:

dont la solution, que l'on obtient par la méthode des équations linéaires du premier ordre d'Analyse II, est

Vérifions l'hypothèse: $\partial f/\partial y=1$ partout, donc $f$ est lipschitzienne en $y$ avec $L=1$, et le théorème 9.1 s'applique sur $[0,2]$ tout entier. Notons au passage $y^{\prime \prime }=2-\frac{1}{2}{\mathrm{e}}^t$, dont le module maximal sur vaut , atteint en : ces deux constantes, et , sont tout ce dont la théorie de l'erreur aura besoin.

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

Une méthode numérique ne calcule pas une fonction: elle calcule un nombre fini de nombres. On découpe donc $[a,b]$ en $n$ intervalles égaux.

Par convention, $y_0=y_0$ — la valeur calculée au premier nœud est la donnée initiale, exacte. L'erreur initiale $e_0$ est donc nulle, ce qui simplifiera les majorations; nous verrons dans la démonstration du théorème 9.3 que l'hypothèse n'est pas nécessaire, et ce qu'il advient si la donnée elle-même est entachée d'erreur.

La forme (9.5) dit tout de suite quelque chose d'utile: une méthode à un pas est une récurrence, elle n'a besoin d'aucun démarrage particulier, et l'on peut changer $h$ d'un pas à l'autre sans rien réorganiser. Cette dernière propriété est ce qui rendra possible le contrôle automatique du pas, plus loin dans le chapitre.

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

Vous avez déjà croisé cette méthode en Analyse II (chapitre 1), comme illustration géométrique du champ des directions. Nous la reprenons ici parce que la manière dont on l'obtient détermine ce que l'on saura en faire ensuite, et parce que la troisième des trois dérivations qui suivent est celle qui engendre tout le reste du chapitre.

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

Le chapitre 8 a construit l'approximation de la dérivée par différence finie progressive:

En remplaçant $y^{\prime }(t_n)$ par $f(t_n,y(t_n))$ selon l'équation, puis les valeurs exactes par les valeurs calculées, on obtient

C'est la méthode d'Euler explicite, de fonction d'incrément $\Phi (t,y,h)=f(t,y)$. La dérivation est courte et elle rattache immédiatement la méthode au chapitre 8 — le prix à payer est qu'elle n'indique pas comment faire mieux: passer à une différence centrée donnerait un schéma qui a besoin de deux valeurs passées, donc une méthode multipas.

Deuxième dérivation: Taylor

La formule de Taylor avec reste de Lagrange (Analyse I, chapitre 9), appliquée à la solution exacte entre $t_n$ et $t_{n+1}=t_n+h$, donne, si $y$ est deux fois dérivable,

avec ${\xi }_n\in ]t_n,t_{n+1}[$. Euler explicite consiste à jeter le reste. Cette lecture-là est celle qui donne le contrôle de l'erreur: le terme négligé est explicitement $\frac{h^2}{2}{y}^{\prime \prime }({\xi }_n)$, et il se majore. Elle suggère aussi une famille de méthodes — garder le terme en , celui en , etc. — au prix de dérivées de que l'on n'a généralement pas: pour un second membre un peu compliqué, et les ordres supérieurs deviennent inutilisables à la main comme en machine.

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

C'est celle qu'il faut retenir. Intégrons l'équation $y^{\prime }=f(t,y)$ entre $t_n$ et $t_{n+1}$; le théorème fondamental du calcul intégral (Analyse I, chapitre 10) donne l'identité exacte

Aucune approximation n'a encore été faite. Toute méthode à un pas s'obtient maintenant en choisissant une formule de quadrature pour l'intégrale de droite — et le chapitre 8 en a fourni un catalogue.

• Le rectangle à gauche, ${\int }_{t_n}^{t_{n+1}}g\approx hg(t_n)$, redonne exactement Euler explicite (9.6).

Euler explicite en Python

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

Le coût est de une évaluation de $f$ par pas, soit $n$ évaluations pour parcourir $[a,b]$, plus une multiplication et une addition. C'est le minimum absolu: aucune méthode ne peut se contenter de moins d'une évaluation du second membre par pas. Notez que $t$ est recalculé par $t_0+(k+1)h$ et non accumulé par t += h: accumuler $n$ additions flottantes ferait dériver la grille elle-même, d'une quantité bornée par d'après le chapitre 1. Sur pas, la différence est visible.

L'ordre d'Euler, lu dans un tableau

Faisons ce que le chapitre 8 nous a appris à faire: mesurer l'erreur pour une suite de pas et regarder le rapport entre erreurs consécutives.

La colonne des rapports est la seule qui apprenne quelque chose, et elle ne dit pas 2. Elle dit $1,65$, puis $1,90$, puis $1,95$, puis $1,97$: elle tend vers 2 sans jamais l'atteindre.

Euler implicite, et ce que coûte un pas

Reprenons (9.8) et choisissons cette fois le rectangle à droite: ${\int }_{t_n}^{t_{n+1}}g\approx hg(t_{n+1})$. Il vient

C'est la méthode d'Euler implicite (ou rétrograde). L'inconnue $y_{n+1}$ apparaît des deux côtés: (9.9) n'est plus une formule, c'est une équation à résoudre à chaque pas.

Que coûte un pas implicite? Il faut résoudre, en l'inconnue $z$, l'équation

C'est le sujet du chapitre 2, et la réponse y est la méthode de Newton:

initialisée par $z^{(0)}=y_n$, ou mieux par le prédicteur d'Euler explicite $z^{(0)}=y_n+hf(t_n,y_n)$, qui est déjà à de la réponse. Deux ou trois itérations suffisent alors, chacune coûtant une évaluation de et une de . , et il exige en plus une dérivée que l'on n'a pas toujours.

Notons aussi, avec le chapitre 2, que si $hL<1$ l'application $z\mapsto y_n+hf(t_{n+1},z)$ est une contraction de rapport $hL$, donc l'itération de point fixe la plus naïve converge déjà. C'est pourquoi la condition revient partout dans la littérature des méthodes implicites: elle garantit que le pas est bien défini et calculable.

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

Retour à (9.8), avec cette fois la formule du trapèze du chapitre 8:

On obtient la méthode du trapèze, appelée schéma de Crank-Nicolson dans le contexte des équations aux dérivées partielles:

Elle est implicite, comme (9.9), et pour la même raison. Mais elle est d'ordre 2, ce qui compense largement: le chapitre 8 a établi que l'erreur de la formule du trapèze sur un intervalle de longueur $h$ vaut $-\frac{h^3}{12}{g}^{\prime \prime }(\eta )$, d'où une erreur locale en $h^3$ et, par le théorème 9.3 ci-dessous, une erreur globale en .

L'idée de Heun consiste à supprimer le caractère implicite en remplaçant l'inconnue $y_{n+1}$ du membre de droite par une prédiction, fournie par Euler explicite:

Le prédicteur commet une erreur en $O(h^2)$; elle entre dans le correcteur multipliée par $\frac{h}{2}{\partial }_{y}f$, donc contribue en $O(h^3)$ à l'erreur du pas — . C'est là tout le mécanisme prédicteur-correcteur: on a le droit d'être grossier dans la prédiction, parce que l'erreur de prédiction entre dans le résultat affectée d'un facteur supplémentaire.

Runge-Kutta: échantillonner la pente plusieurs fois

Le principe est maintenant visible. Sur le pas $[t_n,t_{n+1}]$, la formule (9.8) demande la moyenne de $f(s,y(s))$; une seule évaluation, au bord gauche, donne l'ordre 1; deux évaluations bien placées donnent l'ordre 2. Il est naturel de continuer.

Les méthodes à deux étages

Avec $s=2$, le tableau se réduit à trois paramètres $c_2$, $a_{21}=c_2$, $b_1$, , et

Un développement de Taylor à deux variables (Analyse II, chapitre 5) montre que la méthode est d'ordre 2 si et seulement si

C'est une famille à un paramètre. Le choix $c_2=1$, $b_1=b_2=\frac{1}{2}$ donne ; le choix , , donne le ; le choix donne la méthode de Ralston, qui minimise la constante d'erreur. Leurs tableaux de Butcher:

RK4, la méthode classique

Avec quatre étages, le choix devenu canonique est celui de Kutta (1901), qui reproduit la formule de Simpson lorsque $f$ ne dépend pas de $y$:

Son tableau de Butcher:

Les poids $\frac{1}{6},\frac{1}{3},\frac{1}{3},\frac{1}{6}$ sont ceux de Simpson, avec le point milieu compté deux fois parce qu'il est évalué deux fois: à partir de , puis à partir de . Cette double évaluation du milieu est ce qui fait passer de l'ordre 3 à l'ordre 4.

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)

Le coût est de quatre évaluations de $f$ par pas, plus une dizaine d'opérations arithmétiques, négligeables devant une évaluation de $f$ dès que celle-ci est un tant soit peu compliquée. Un pas de RK4 coûte donc quatre pas d'Euler.

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

Jusqu'à quatre étages, le nombre d'évaluations et l'ordre coïncident: une évaluation donne l'ordre 1, deux l'ordre 2, trois l'ordre 3, quatre l'ordre 4. Cela s'arrête là. Butcher a montré qu'une méthode de Runge-Kutta explicite à $s$ étages ne peut pas dépasser l'ordre $s$, et surtout que l'ordre $s$ n'est atteignable que pour $s\le 4$: il faut 6 étages pour l'ordre 5, 7 pour l'ordre 6, 9 pour l'ordre 7 et 11 pour l'ordre 8. Le rendement décroît donc brutalement à partir du cinquième étage.

C'est la raison pour laquelle RK4 est resté, plus d'un siècle après Kutta, la méthode par défaut des cours et de beaucoup de programmes: elle est le dernier point où l'on obtient un ordre entier pour le prix d'une évaluation. Au-delà, les méthodes d'ordre élevé existent et sont excellentes — mais elles ne se justifient que pour des tolérances très serrées, où leur ordre finit par l'emporter sur leur coût par pas.

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

Nous avons maintenant un catalogue de méthodes et une observation empirique: leur erreur se comporte comme $h^p$ avec un $p$ qui dépend de la méthode. Il est temps de démontrer pourquoi.

La division par $h$ dans (9.16) est une convention, mais elle est universelle et il faut savoir la lire: l'erreur du pas vaut $h{\tau }_n$, donc une méthode consistante d'ordre $p$ commet une erreur en $h^{p+1}$ à chaque pas et il y a $N=(b-a)/h$ pas. Le compte naïf — « erreurs de taille , donc une erreur totale de » — donne le bon exposant, et c'est déjà rassurant. Mais il est faux comme raisonnement, car il oublie que chaque erreur passée est par les pas suivants, comme l'exemple 9.1 l'a montré. Le théorème de convergence est précisément ce qui répare ce raisonnement.

Le lemme d'accumulation

Tout le travail tient dans une inégalité de récurrence, qui est la version discrète du lemme de Grönwall.

Démonstration. Récurrence sur $n$. Pour $n=0$, le membre de droite de (9.18) vaut ${\epsilon }_0$ et l'inégalité est une égalité. Supposons (9.18) vraie au rang $n$. Alors

ce qui est (9.18) au rang $n+1$. Le cas $A=0$ est immédiat. $\square$

La forme (9.18) est encore un peu opaque; on la rend lisible avec l'inégalité $1+x\le {\mathrm{e}}^x$, valable pour tout réel $x$ (elle exprime la convexité de l'exponentielle, dont la tangente en $0$ est $1+x$). Il vient $(1+A{)}^n\le {\mathrm{e}}^{nA}$, donc

On reconnaît déjà le facteur exponentiel du théorème 9.1: il ne s'agit pas d'une coïncidence, mais de la même propagation, vue une fois sur la solution exacte et une fois sur la solution discrète.

Consistance et stabilité entraînent la convergence

Cette définition mérite d'être commentée avant d'être utilisée. La stabilité ne dit rien sur la qualité de l'approximation: elle dit que la méthode ne fabrique pas d'amplification supplémentaire, que deux exécutions parties de données voisines restent voisines sur un intervalle borné. C'est une condition de continuité de la récurrence par rapport à ses données, exactement comme (9.2) l'est pour l'équation elle-même. Pour Euler explicite, $\Phi =f$ et (9.20) est vérifiée avec $\Lambda =L$; pour Heun, un calcul de deux lignes donne $\Lambda =L+\frac{h}{2}{L}^2\le L+\frac{h_0}{2}{L}^2$; pour RK4, . , ce qui rend cette hypothèse peu contraignante — et c'est justement ce que le chapitre 10 remettra en cause, sous un autre sens du mot «stabilité».

Démonstration. Écrivons côte à côte ce que fait la méthode et ce que fait la solution exacte. La méthode:

La solution exacte, par la définition (9.16) de l'erreur locale de troncature, réarrangée:

Cette seconde ligne est une identité, pas une approximation: ${\tau }_n$ a été défini pour qu'elle le soit. Soustrayons, avec $e_n=y_n-y(t_n)$:

C'est l'équation centrale du chapitre, et elle se lit ainsi: l'erreur nouvelle est l'erreur ancienne, transportée par la méthode, plus l'erreur fraîchement commise. Le crochet est le transport, le terme $h{\tau }_n$ est la fraîche.

Majorons. L'inégalité triangulaire, puis la stabilité (9.20) appliquée au crochet avec $u=y_n$ et $v=y(t_n)$, puis la consistance (9.17):

Nous sommes exactement dans les hypothèses du théorème 9.2, avec ${\epsilon }_n=|e_n|$, $A=h\Lambda >0$ et $B=C{h}^{p+1}$. La conclusion (9.19) donne

Il ne reste qu'à simplifier: $nh=t_n-a$, et le facteur $h$ du dénominateur absorbe une puissance de $h$ du numérateur, laissant $C{h}^p/\Lambda$. On obtient (9.21).

Si $e_0=0$, le premier terme disparaît; et comme $t_n\le b$, l'exponentielle est majorée par ${\mathrm{e}}^{\Lambda (b-a)}$, d'où $|e_n|\le K{h}^p$ avec la constante annoncée, indépendante de et de .

L'ordre 1 d'Euler explicite, en entier

Démonstration. Nous vérifions les deux hypothèses du théorème 9.3, puis nous appliquons (9.21).

Stabilité. Pour Euler explicite, $\Phi (t,y,h)=f(t,y)$, et (9.20) est exactement l'hypothèse de Lipschitz (9.2) sur $f$: la méthode est stable avec $\Lambda =L$, pour tout $h$.

Consistance d'ordre 1. La solution exacte étant ${\mathcal{C}}^2$, la formule de Taylor avec reste de Lagrange (Analyse I, chapitre 9) donne, pour chaque $n$, l'existence d'un ${\xi }_n\in ]t_n,t_{n+1}[$ tel que

Comme $y$ est solution, $y^{\prime }(t_n)=f(t_n,y(t_n))=\Phi (t_n,y(t_n),h)$. En reportant dans la définition (9.16):

Donc $|{\tau }_n(h)|\le \frac{M}{2}h$ pour tout $n$: la méthode est consistante d'ordre $p=1$ avec $C=M/2$.

Accumulation. On pourrait invoquer le théorème 9.3 et conclure. Refaisons-le ici en entier, car c'est le cas où l'on voit tout. L'identité fondamentale s'écrit, avec $\Phi =f$:

d'où, par l'inégalité triangulaire et la condition de Lipschitz (9.2),

C'est le lemme 9.2 avec $A=hL$ et $B=M{h}^2/2$, et ${\epsilon }_0=|e_0|=0$ puisque . Il donne

Enfin $1+hL\le {\mathrm{e}}^{hL}$, donc $(1+hL{)}^n\le {\mathrm{e}}^{nhL}={\mathrm{e}}^{L(t_n-a)}$, ce qui donne (9.22).

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

Tout ce qui précède suppose un pas fixe, décidé à l'avance. Ce n'est presque jamais ce que l'on veut. La solution d'une équation différentielle a rarement la même échelle de temps du début à la fin: une réaction chimique démarre vite puis se calme, un projectile passe par un choc, un système mécanique traverse une phase transitoire avant son régime permanent. Un pas assez fin pour la phase rapide serait du gaspillage partout ailleurs.

Il faut donc estimer l'erreur commise pendant le calcul, sans connaître la solution exacte. L'idée, aujourd'hui universelle, est simple: on calcule le même pas avec deux méthodes d'ordres différents, et l'on prend leur différence comme estimation de l'erreur de la moins précise.

Soient ${\hat{y}}_{n+1}$ le résultat d'une méthode d'ordre $p$ et ${\hat{z}}_{n+1}$ celui d'une méthode d'ordre $p+1$, calculés depuis le même . Alors

donc, en soustrayant,

La différence des deux résultats estime l'erreur du moins bon, à un ordre près. C'est un estimateur a posteriori: il ne coûte aucune connaissance de la solution, seulement un second calcul.

Le second calcul serait cher s'il fallait le refaire de zéro. D'où l'idée décisive, qui est celle des paires emboîtées: construire deux méthodes de Runge-Kutta d'ordres $p$ et $p+1$ qui partagent leurs étages, c'est-à-dire le même tableau $A$ et les mêmes $c_i$, et ne diffèrent que par les poids $\mathbf{b}$. Le tableau de Butcher porte alors deux lignes de poids:

et l'estimation d'erreur est gratuite: elle vaut $h{\sum }_i(b_i-{\hat{b}}_i)k_i$, une simple combinaison des pentes déjà calculées.

La plus simple des paires emboîtées est Heun-Euler, d'ordres 2 et 1: les deux étages $k_1$ et $k_2$ de Heun contiennent déjà le pas d'Euler.

Reste à choisir le nouveau pas. Si l'on demande une tolérance $\text{tol}$ par pas et que l'erreur se comporte comme $\mathrm{est}\approx K{h}^{p+1}$, le pas $h_{\text{nouv}}$ qui atteindrait exactement la tolérance vérifie $\text{tol}=K{h}_{\text{nouv}}^{p+1}$, d'où, en divisant les deux relations:

où $\sigma \approx 0,9$ est un coefficient de sécurité qui rend acceptable le pas proposé plutôt que marginal. L'algorithme complet tient en quelques lignes: on calcule le pas, on compare ${\mathrm{est}}_n$ à $\text{tol}$, on accepte si c'est assez bon et l'on recommence sinon, et dans les deux cas on met $h$ à jour par (9.25) — en bridant le facteur multiplicatif, typiquement entre $0,2$ et $5$, pour éviter les oscillations.

Les paires réellement utilisées sont d'ordre plus élevé. Runge-Kutta-Fehlberg RKF45 utilise 6 étages pour les ordres 4 et 5; Dormand-Prince DOPRI5, également d'ordres 4 et 5 sur 7 étages, est le moteur de la fonction ode45 de Matlab et du solveur par défaut de la plupart des bibliothèques scientifiques. Son septième étage est celui de la propriété first same as last: la dernière pente d'un pas est la première du suivant, ce qui ramène le coût effectif à six évaluations par pas. Bogacki-Shampine 3(2), quatre étages avec la même astuce, est le choix courant pour les tolérances lâches.

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

Tout ce chapitre a été écrit pour une équation scalaire, et il s'applique tel quel à un système — c'est le principal argument en faveur de la forme (9.1).

Aucune formule de ce chapitre ne change. Euler s'écrit ${\mathbf{u}}_{n+1}={\mathbf{u}}_n+h\mathbf{f}(t_n,{\mathbf{u}}_n)$, RK4 calcule quatre vecteurs , et les démonstrations des théorèmes 9.2 à 9.4 restent valables mot pour mot en remplaçant les valeurs absolues par une norme vectorielle (annexe A) et la condition de Lipschitz par . Seul le coût change: une évaluation du second membre coûte maintenant composantes.

La réduction d'une équation d'ordre $m$ à un système du premier ordre est la même qu'en Analyse II. Soit

On pose $u_1=y$, $u_2=y^{\prime }$, …, $u_m=y^{(m-1)}$, et le système devient

Synthèse

• Le problème de Cauchy $y^{\prime }=f(t,y)$, $y(a)=y_0$ possède une solution unique sur $[a,b]$ dès que $f$ est continue et (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 — le même facteur qui bornera l'erreur numérique.

Exercices

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

Références

• Quarteroni, A., Sacco, R. et Saleri, F., Méthodes numériques — algorithmes, analyse et applications, Springer, Milan, chap. 11 (problème de Cauchy, méthodes à un pas, consistance et convergence).
• Hairer, E., Nørsett, S. P. et Wanner, G., Solving Ordinary Differential Equations I — Nonstiff Problems, 2e éd., Springer, Berlin, chap. II (la référence sur les méthodes de Runge-Kutta, les conditions d'ordre et les paires emboîtées).
• Rappaz, J. et Picasso, M., Introduction à l'analyse numérique, Presses polytechniques et universitaires romandes, Lausanne, chap. 8.
• Burden, R. L. et Faires, J. D., Numerical Analysis, Cengage, Boston, chap. 5 (le problème témoin de ce chapitre y est traité en détail).
• Butcher, J. C., Numerical Methods for Ordinary Differential Equations, 3e éd., Wiley, Chichester (les barrières d'ordre et la théorie des arbres).
• Dormand, J. R. et Prince, P. J., «A family of embedded Runge-Kutta formulae», Journal of Computational and Applied Mathematics, vol. 6, n° 1, 1980 (la paire DOPRI5, moteur de la plupart des intégrateurs actuels).