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 décrit le principe de l'algorithme, puis la façon dont il est appliqué dans nos outils — et ce qui, chez eux, n'en relève pas.

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.

L'intégration numérique contourne l'obstacle : on découpe le temps en petites tranches (le pas de temps $\Delta t$) et on recalcule l'état de proche en proche — position et vitesse en balistique extérieure, auxquelles s'ajoutent pression, température et fraction de poudre brûlée en balistique intérieure.

L'algorithme ne se contente pas de prolonger la courbe en ligne droite à partir de la pente initiale — ce serait la méthode d'Euler, d'ordre 1, immédiate à écrire mais qui accumule l'erreur trop vite pour un tir à longue distance.

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 globale décroît comme $\Delta t^4$ : diviser le pas de temps par deux divise l'erreur d'intégration par $2^4 = 16$. L'erreur commise sur un pas isolé est d'un ordre plus élevée encore ($\Delta t^5$) ; c'est son accumulation sur l'ensemble des pas d'une trajectoire qui ramène à l'ordre 4.

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. Ces coefficients se mesurent au tunnel balistique ou en soufflerie, ou s'estiment par des codes comme McDrag et McGyro, dus à Robert L. McCoy.

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 — Über die numerische Auflösung von Differentialgleichungen, Mathematische Annalen 46, 1895.
2. W. Kutta — Beitrag zur näherungsweisen Integration totaler Differentialgleichungen, 1901 (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.
5. R. L. McCoy — Modern Exterior Ballistics, Schiffer, 2ᵉ éd., 2012 (modèles point-masse et 6-DOF, codes McDrag et McGyro).

• 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