Vous recherchez Runge-Kutta 8ème ordre en C / C ++

Je voudrais utiliser la méthode Runge-Kutta de 8e ordre (89) dans une application de mécanique / astrodynamique céleste, écrite en C ++, à l'aide d'une machine Windows. Par conséquent, je me demande si quelqu'un connaît une bonne bibliothèque / implémentation documentée et gratuite à utiliser? C'est correct s'il est écrit en C, tant qu'il n'y a pas de problèmes de compilation à prévoir.

Jusqu'à présent, j'ai trouvé cette bibliothèque (mymathlib) . Le code semble correct, mais je n'ai trouvé aucune information sur les licences.

Pouvez-vous m'aider en révélant certaines des alternatives que vous pourriez connaître et qui conviendraient à mon problème?

EDIT:
Je vois qu'il n'y a pas vraiment autant de codes sources C / C ++ disponibles que je m'y attendais. Par conséquent, une version Matlab / Octave serait également acceptable (doit toujours être gratuite).

La bibliothèque scientifique GNU (GSL) (C) et Boost Odeint (C ++) proposent des méthodes Runge-Kutta de 8e ordre.

Les deux sont open source, et sous linux et mac, ils devraient être directement disponibles depuis le gestionnaire de paquets. Sous Windows, il vous sera probablement plus facile d'utiliser Boost plutôt que GSL.

GSL est publié sous licence GPL, et Boost Odeint sous licence Boost.

Edit: Ok, Boost Odeint n'a PAS la méthode Runge-Kutta 89, seulement le 78, mais il fournit une recette pour faire des steppers Runge-Kutta arbitraires.

Cependant, les méthodes de 8e ordre sont assez élevées et sont probablement excessives pour votre problème.

Prince-Dormand fait référence à un type spécifique de Runge-Kutta et n'est pas directement lié à la commande, mais le plus courant est 45. Matlabs ode45, qui est leur algorithme ODE recommandé, est une implémentation Prince-Dormand 45. Il s'agit du même algorithme que celui implémenté dans Boost Odeint Runge_Kutta_Dopri5 .

Si vous faites de la mécanique céleste sur de longues échelles de temps, l'utilisation d'un intégrateur Runge-Kutta classique ne préservera pas l'énergie. Dans ce cas, l'utilisation d'un intégrateur symplectique serait probablement préférable. Boost.odeint implémente également un schéma Runge-Kutta symplectique de 4ème ordre qui fonctionnerait mieux pendant de longs intervalles de temps. GSL ne met en œuvre aucune méthode symplectique, pour autant que je sache.

résumant certains points:

1. S'il s'agit d'une intégration à long terme d'un modèle non dissapatif, un intégrateur symplectique est ce que vous recherchez.
2. Sinon, puisqu'il s'agit d'une équation de mouvement, les méthodes Runge-Kutta Nystrom seront plus efficaces qu'une transformation en un système de premier ordre. Il existe des méthodes RKN d'ordre élevé en raison de DP. Il y a quelques implémentations, comme ici dans Julia, elles sont documentées et en voici une MATLAB .
3. Des méthodes Runge-Kutta d'ordre élevé ne sont nécessaires que si vous voulez une solution de haute précision. S'il s'agit de tolérances inférieures, un RK de 5e ordre sera probablement plus rapide (pour la même erreur). La meilleure chose à faire si vous devez résoudre ce problème souvent est de tester un tas de méthodes différentes. Dans cet ensemble de repères sur les problèmes à 3 corps, nous voyons que (pour la même erreur) les méthodes RK d'ordre élevé ne sont qu'une amélioration marginale de la vitesse, mais comme erreur -> 0, vous pouvez voir que l'amélioration va déjà à> 5x contre Dormand -Prince 45 ( DP5) lorsque vous regardez à 4 chiffres de précision (les tolérances sont cependant beaucoup plus faibles pour cela. Les tolérances ne sont qu'un stade approximatif dans tous les problèmes). Au fur et à mesure que vous réduisez les tolérances, l'amélioration d'une méthode RK d'ordre élevé augmente, mais vous devrez peut-être commencer à utiliser des nombres plus précis.
4. L'algorithme d'ordre 7/8 de Dormand-Prince a un tableau du 8e ordre différent de la méthode DP853 de Hairer's dop853et DifferentialEquations.jl DP8(qui sont les mêmes). La dernière méthode 853 ne peut pas être implémentée dans la version tableau standard d'une méthode Runge-Kutta car son estimateur d'erreur n'est pas standard. Mais cette méthode est beaucoup plus efficace et je ne recommanderais même pas d'utiliser les anciennes méthodes Fehlberg 7/8 ou DP 7/8.
5. Pour les méthodes RK d'ordre élevé, les méthodes Verner «efficaces» sont l'étalon-or. Cela apparaît dans les repères que j'ai liés. Vous pouvez les coder dans Boost vous-même, ou utiliser l'un des 2 packages qui les implémentent si vous les souhaitez plus facilement (Mathematica ou DifferentialEquations.jl).

Je voudrais ajouter que même si ce que Geoff Oxberry suggère pour une intégration à long terme (en utilisant des intégrateurs symplectiques) est vrai, dans certains cas, cela ne fonctionnera pas. Plus précisément, si vous avez des forces dissipatives, votre système ne conserve plus l'énergie, et vous ne pouvez donc pas recourir à des intégrateurs symplectiques dans ce cas. La personne qui posait la question parlait d'orbites terrestres basses, et ces orbites présentent une grande quantité de traînée atmosphérique, c'est-à-dire une force dissipative qui empêche d'utiliser de tels intégrateurs symplectiques.

Dans ce cas spécifique (et pour les cas où vous ne pouvez pas utiliser / n'avez pas accès à / ne souhaitez pas utiliser d'intégrateurs symplectiques), je recommanderais l'utilisation de l'intégrateur Bulirsch-Stoer si vous avez besoin de précision et d'efficacité sur de longues périodes. Il fonctionne bien par expérience et est également recommandé par les recettes numériques (Press et al., 2007).