L'intégration numérique par la méthode de Runge-Kutta d'ordre 4 (RK4)

Les équations de la balistique — intérieure comme extérieure — n'admettent pas de solution analytique en forme close dès qu'on emploie un modèle réaliste : la traînée suit une table expérimentale fonction du nombre de Mach, l'air change de densité avec l'altitude, la poudre brûle selon une loi empirique. On les résout donc numériquement, pas de temps par pas de temps. L'algorithme employé pour cela sur Tireur.org est la méthode de Runge-Kutta d'ordre 4, abrégée RK4, due aux travaux de Carl Runge (1895) et Wilhelm Kutta (1901).

Cette page part du problème lui-même, expose la méthode la plus élémentaire qui le résout — celle d'Euler — puis ce que RK4 y ajoute, et décrit enfin la façon dont l'algorithme est appliqué dans nos outils, et ce qui, chez eux, n'en relève pas. Aucune connaissance préalable en analyse numérique n'est supposée.

Un calculateur balistique répond à une question banale : où est la balle, et à quelle vitesse, une seconde après le départ du coup ? Ce qu'on lui donne pour cela n'est pas une trajectoire, mais un point de départ et une règle :

• l'état initial — position et vitesse à la bouche ;
• une règle de changement — pour un état donné, on sait calculer l'accélération que subit la balle : gravité vers le bas, traînée opposée au vecteur vitesse, éventuellement Coriolis.

Toute la difficulté tient au fait que cette règle ne donne jamais la position directement. Elle ne donne que la tendance : de combien la vitesse est en train de changer, à cet instant précis, dans cet état précis. Et cette tendance change avec l'état — la traînée dépend de la vitesse, or la vitesse est précisément ce que la traînée modifie. Chacune des deux grandeurs est l'entrée de l'autre.

C'est ce que les mathématiques appellent une équation différentielle : une relation qui ne dit pas ce que vaut la quantité cherchée, mais à quelle vitesse elle varie. Elle s'écrit toujours de la même façon :

$$ \frac{dy}{dt} = f(t, y) $$

À gauche, la variation de l'état ; à droite, la règle de changement $f$, qui calcule cette variation à partir de l'instant $t$ et de l'état courant $y$. Cette écriture gouverne tout le reste de la page : $y$ y désignera le paquet « position et vitesse » du projectile, et $f$ le calcul des forces qui s'y appliquent.

Le cas d'école est la chute libre sans air, dont la règle de changement ne dépend de rien : elle vaut $-g$ à tout instant. On sait alors « remonter » de la variation à la quantité elle-même, une fois pour toutes — c'est ce qu'on appelle intégrer l'équation :

$$ \frac{dv}{dt} = -g \quad \Longrightarrow \quad v = -gt \quad \Longrightarrow \quad h = -\frac{1}{2}gt^2 $$

La formule obtenue donne la hauteur à n'importe quel instant, sans calcul intermédiaire. Ce qui met fin à ce confort, c'est que la vraie règle de changement dépend de l'état : $f$ réclame la vitesse pour calculer la traînée, et la traînée modifie la vitesse. Quelques cas de ce genre restent solubles à la main — on s'en servira plus bas comme étalon — mais aucune formule ne survit à un coefficient de traînée lu dans une table de mesures.

Il reste alors une méthode dont personne ne peut nous priver : avancer par petits pas. On connaît l'état maintenant, donc on connaît la tendance maintenant ; on suppose qu'elle reste à peu près valable pendant un très court instant $\Delta t$ ; on en déduit l'état à $t + \Delta t$ ; et on recommence. Mille fois, dix mille fois, jusqu'à la cible. C'est exactement ce que fait un GPS de randonnée qui reconstitue un parcours à partir d'un cap et d'une vitesse relevés chaque seconde.

Cette approche s'appelle l'intégration numérique. Elle ne donne pas une formule, elle donne une liste de points — et c'est bien assez pour dresser une table de chute. Tout le savoir-faire consiste alors à choisir comment franchir chaque petit pas : c'est là que se joue la différence entre la méthode d'Euler et RK4, et c'est l'objet de cette page.

Sur un projectile en vol :

• la gravité exerce une accélération quasi constante vers le bas ;
• la traînée le freine, mais elle dépend de la vitesse, et pas linéairement : le coefficient de traînée est lui-même fonction du nombre de Mach, avec un pic marqué au passage transsonique ;
• la vitesse du son varie avec la température de l'air traversé, si bien qu'une même vitesse ne correspond pas au même nombre de Mach d'un bout à l'autre de la trajectoire.

L'accélération dépend donc de la vitesse, qui dépend elle-même de l'accélération : le système d'équations différentielles est couplé et non linéaire. Des solutions analytiques existent pour des lois de traînée simples — traînée proportionnelle à $v$ ou à $v^2$ — et l'artillerie a longtemps vécu sur des solutions approchées de ce type, dont la méthode de Siacci pour le tir tendu. Aucune ne vaut pour une table de traînée mesurée comme G1 ou G7.

Ces lois de traînée simplifiées gardent malgré tout une utilité : elles servent d'étalon. C'est ce qui permet, plus bas, de chiffrer l'erreur d'une méthode — on l'applique à un problème dont on connaît par ailleurs la réponse rigoureuse, et l'écart se lit directement.

L'intégration numérique, elle, ne demande rien d'autre que de savoir évaluer la règle de changement en un état quelconque : la table G7 peut être aussi biscornue qu'elle veut, on l'interroge point par point. Ce qu'on fait avancer de proche en proche s'appelle l'état du système : position et vitesse en balistique extérieure, auxquelles s'ajoutent pression, température et fraction de poudre brûlée en balistique intérieure.

La façon la plus directe de franchir un pas consiste à prendre la tendance telle qu'elle est au départ et à la tenir pour constante pendant tout le pas. Si $y_n$ est l'état à l'instant $t_n$ et $f(t_n, y_n)$ la règle de changement évaluée sur cet état, alors :

$$ y_{n+1} = y_n + \Delta t \cdot f(t_n, y_n) $$

C'est la méthode d'Euler, décrite par Leonhard Euler au XVIIIᵉ siècle, et elle tient tout entière dans cette ligne. Géométriquement, elle remplace la courbe par sa tangente au point de départ : on quitte la vraie trajectoire au point où l'on est, on file tout droit dans la direction où elle pointait, et on décrète que le point d'arrivée est le nouvel état.

Son défaut se lit dans cette description même. La tendance était juste au départ du pas, et seulement là ; dès qu'on a bougé, elle ne l'est plus. Une balle qui ralentit subit une traînée d'autant plus faible qu'elle va moins vite : en gardant pendant tout le pas la traînée du début, la plus forte, Euler freine trop. L'erreur n'est pas une dispersion aléatoire qui se compenserait, c'est un biais systématique, toujours du même côté, qui s'accumule pas après pas.

Prenons un seul pas de $\Delta t = 0{,}1$ s depuis la bouche. À 800 m/s, la décélération vaut 581 m/s², soit près de 59 g. Euler en conclut :

$$ v = 800 - 0{,}1 \times 581{,}1 = 741{,}9\ \text{m/s} $$

La solution exacte donne 745,8 m/s. En un seul pas, 3,9 m/s d'écart — et toujours par défaut, jamais par excès.

Le remède évident est de raccourcir le pas : la tangente s'écarte moins de la courbe si on la suit moins longtemps. La question est de savoir combien on gagne, et c'est ce que mesure l'ordre de la méthode.

Une méthode est d'ordre $p$ lorsque son erreur décroît comme $\Delta t^p$. Euler est d'ordre 1 : diviser le pas par deux ne divise l'erreur que par deux. RK4 est d'ordre 4 : diviser le pas par deux divise l'erreur par $2^4 = 16$. Sur le banc d'essai ci-dessus, à une seconde de vol (vitesse exacte : 463,40 m/s) :

Chaque ligne divise le pas par deux. La colonne Euler se divise par deux, la colonne RK4 par seize. L'écart entre les deux méthodes ne se referme donc jamais : il s'ouvre à chaque raffinement.

C'est la raison de fond pour laquelle personne n'intègre une trajectoire avec Euler. Non parce que la méthode serait fausse — elle converge bien vers la bonne réponse — mais parce qu'elle paie très cher chaque décimale gagnée.

L'idée de Runge et Kutta est que la faiblesse d'Euler n'est pas d'aller tout droit, mais de choisir sa direction sur la seule information disponible au point de départ. Avant de s'engager, on peut aller sonder la pente ailleurs dans l'intervalle, puis partir dans une direction moyenne mieux informée — sans jamais avoir besoin de connaître la solution.

Pour chaque pas de temps $\Delta t$, RK4 évalue quatre fois la dérivée, réparties sur trois instants de l'intervalle :

1. $k_1$ : la pente au début de l'intervalle, sur l'état courant.
2. $k_2$ : la pente au milieu de l'intervalle, sur un état estimé par un demi-pas dans la direction $k_1$.
3. $k_3$ : à nouveau au milieu de l'intervalle, sur un état estimé cette fois avec $k_2$.
4. $k_4$ : la pente à la fin de l'intervalle, sur un état estimé avec $k_3$.

$k_2$ et $k_3$ sont donc pris au même instant $t + \Delta t/2$, sur deux estimations différentes de l'état — ce sont les évaluations, non les instants, qui sont au nombre de quatre.

L'état suivant $y_{n+1}$ est une moyenne pondérée de ces quatre pentes :

$$ y_{n+1} = y_n + \frac{\Delta t}{6} \left( k_1 + 2k_2 + 2k_3 + k_4 \right) $$

Le poids double accordé à $k_2$ et $k_3$, les évaluations centrales, revient à privilégier le milieu du pas — c'est ce qui fait la précision de la formule. L'erreur commise sur un pas isolé décroît comme $\Delta t^5$ ; c'est son accumulation sur l'ensemble des pas d'une trajectoire qui ramène la méthode à l'ordre 4.

Sur le même pas de 0,1 s que tout à l'heure, les quatre sondages donnent $k_1 = -581{,}1$, $k_2 = -539{,}6$, $k_3 = -542{,}5$ et $k_4 = -505{,}0$ m/s². Leur moyenne pondérée vaut −541,7 m/s² : RK4 avance donc avec une décélération sensiblement plus faible que les 581 m/s² d'Euler, parce qu'il a vu venir l'allègement de la traînée au lieu de le découvrir au pas suivant. Le résultat, 745,8259 m/s, coïncide avec la solution exacte à 5 × 10⁻⁵ m/s près — quand Euler s'en écartait de 3,9 m/s.

L'objection est immédiate : RK4 évalue la dérivée quatre fois par pas, là où Euler ne le fait qu'une. Un pas de RK4 coûte donc quatre pas d'Euler, et c'est le nombre d'évaluations — pas le nombre de pas — qui fait le temps de calcul.

La comparaison honnête se fait donc à budget égal. Reprenons la balle du banc d'essai, cette fois avec la gravité, et regardons où elle est après 1,2 s de vol (environ 690 m) ; la référence est la même intégration menée à pas ridiculement court :

À budget strictement identique — 48 évaluations de la dérivée — Euler place la balle deux mètres trop loin et huit centimètres trop haut ; RK4 tombe sous le millimètre. Et même en lui accordant cinquante fois plus de travail, Euler reste plus grossier que RK4 au pas le plus large. C'est tout l'argument : les quatre évaluations ne sont pas un surcoût, elles sont le meilleur usage possible de quatre évaluations.

Nos solveurs avancent à pas constant, choisi une fois pour toutes : $\Delta t = 0{,}0005$ s pour la trajectoire extérieure, 1 microseconde pour la simulation intérieure, où tout se joue en une à deux millisecondes.

Les intégrateurs à pas adaptatif — Runge-Kutta-Fehlberg, Dormand-Prince, ceux que l'on trouve derrière un ode45 — estiment l'erreur à chaque pas en comparant deux formules d'ordres différents, puis resserrent ou relâchent $\Delta t$ en conséquence. Sur une trajectoire, ils raccourciraient le pas au passage transsonique et l'allongeraient sur la fin de course. À pas fixe, la marge se prend en amont : on retient un pas assez court pour la phase la plus raide et on le garde partout. C'est plus simple, un peu plus coûteux en calcul, et suffisant tant que le pas reste petit devant les temps caractéristiques du problème.

En balistique, l'état $y$ n'est pas un scalaire mais un vecteur d'état. Dans notre calculateur 3-DOF, il compte six composantes :

• trois pour la position du projectile ($x, y, z$) ;
• trois pour sa vitesse ($v_x, v_y, v_z$).

Les dérivées de ces six composantes sont donc, respectivement, la vitesse et l'accélération. Cette dernière rassemble la gravité (projetée selon l'angle de site pour le tir en pente), la traînée — calculée sur la vitesse relative à l'air, ce qui fait entrer le vent dans l'équation — et, si l'option est active, la force de Coriolis. Le même RK4 fait ainsi avancer position et vitesse d'un seul mouvement.

Un modèle 6-DOF suit en plus l'attitude du projectile : son vecteur d'état compte douze variables — trois de position, trois de vitesse, trois d'orientation et trois de vitesse angulaire. Il ne réclame plus un simple coefficient de traînée mais un jeu complet de coefficients aérodynamiques : portance, moment de tangage, amortissement en tangage et en roulis, effet Magnus. Ce jeu complet s'obtient par réduction de tirs en tunnel balistique ou d'essais en soufflerie. Les codes semi-empiriques que les tireurs connaissent de nom n'en fournissent chacun qu'une pièce : McDrag, publié par Robert L. McCoy en 1981, estime le seul coefficient de traînée à incidence nulle, et McGyro (1986) le facteur de stabilité gyroscopique et le pas de rayure nécessaire. Ni l'un ni l'autre ne rend la portance, les amortissements ou l'effet Magnus.

Ils ne sont publiés pour presque aucune balle du commerce : c'est ce qui cantonne le 6-DOF à l'étude de la stabilité et de la dispersion plutôt qu'à la prévision courante de trajectoire, où un modèle point-masse correctement renseigné fait pratiquement aussi bien. L'algorithme d'intégration, lui, ne change pas : RK4 s'applique à un vecteur de douze composantes exactement comme à un vecteur de six.

1. C. Runge — Ueber die numerische Auflösung von Differentialgleichungen, Mathematische Annalen 46 (1895), p. 167-178.
2. W. Kutta — Beitrag zur näherungsweisen Integration totaler Differentialgleichungen, Zeitschrift für Mathematik und Physik 46 (1901), p. 435-453 — la formule classique d'ordre 4.
3. J. C. Butcher — Numerical Methods for Ordinary Differential Equations, Wiley, 3ᵉ éd., 2016.
4. E. Hairer, S. P. Nørsett, G. Wanner — Solving Ordinary Differential Equations I: Nonstiff Problems, Springer, 2ᵉ éd., 1993 (Springer Series in Computational Mathematics, 8).
5. R. L. McCoy — Modern Exterior Ballistics: The Launch and Flight Dynamics of Symmetric Projectiles, Schiffer, 2ᵉ éd., 2012 (modèles point-masse et 6-DOF).
6. R. L. McCoy — MC DRAG — A Computer Program for Estimating the Drag Coefficients of Projectiles, Ballistic Research Laboratory, rapport ARBRL-TR-02293, février 1981.

• Balistique extérieure — les forces que le solveur intègre à chaque pas de temps
• Balistique intérieure — le modèle 0D dont RK4 intègre les équations
• Rappels de probabilités et statistiques — l'autre page de fondamentaux mathématiques du wiki
• Calculateur balistique — le solveur 3-DOF en ligne
• Consolidated Competition Ballistics (PDF) — la formulation mathématique complète du solveur

Page vérifiée. Le code des trois outils cités a été relu et le solveur 3-DOF ré-exécuté le 2026-08-12 : fiche.