RK4 versus integrateur leapfrog

Bonjour
Voilà mon problème : je souhaite résoudre les équations de mouvement pour des particules soumises à un champ de potentiel donné et indépendant du temps. Celui-ci a une formule analytique mais je ne peux pas utiliser de méthode implicite car son expression contient des trucs assez moches (intégrale elliptique entre autre, et j'aimerais bien que mon code marche même quand je change profondément la formule du champs).

J'utilise actuellement (dans GEANT4, j'ai aussi besoin des interactions particules/matière) RK4 avec pas adaptatif, hélas les simulations sont longues. Je me demandais alors si je ne peux pas gagner du temps en utilisant un intégrateur plus efficace (c'est-à-dire même erreur sur la conservation de l’énergie pour un calcul plus rapide). De ce que j'ai compris les intégrateurs type leapfrog sont assez efficaces dans certain cas, mais je ne suis pas sûr de bien saisir dans quel cas ils sont plus efficaces que RK4 ?
J'ai notamment fait un très rapidement un benchmark (sur scilab) pour des particules dans un potentiel harmonique 2D. Les méthodes RK4 et leapfrog d'ordre 4 (les deux sans pas adaptatif) semblent plus ou moins équivalentes.

Mon hypothèse, naïve, est pour l'instant que les méthodes leapfrog sont plus à même de maintenir les solutions proches de l’équilibre mais que pour les trajectoires chaotiques elles ne sont pas spécialement plus efficaces que les méthodes RK.

Dans ma situation le potentiel est indépendant du temps (les particules qui bougent ne modifient pas le potentiel, elles sont toutes indépendantes et donc peuvent être résolues de manière individuelle, elle n'interagissent pas entre elles, d'aucune manière). J'essaye de couvrir tout l’espace des phases, et donc l’écrasante majorité des trajectoires sont chaotiques.
(et n'importe quelle remarque sur d'autres méthodes potentiellement efficaces est bienvenue).
(je précise aussi que mon bagage mathématique me rend relativement obscur les publications théorique sur les intégrateurs numériques, manque de connaissance, de vocabulaire etc.)

Réponses

  • Quand on résout numériquement un système Hamiltonien, on a intérêt à utiliser des intégrateurs qui préservent exactement certains aspects géométriques. Les intégrateurs de type leapfrog (appelé aussi Störmer–Verlet), et plus généralement les méthodes de Runge-Kutta Symplectiques, préservent exactement un Hamiltonien proche de ton Hamiltonien de départ sur des temps exponentiellement longs. Du coup elles sont beaucoup plus fiables que le basique RK4 à pas variable, et notamment pour les trajectoires chaotiques (alors que tu pensais l'inverse).

    En anglais mais très accessible il y cet ancien article Hairer et al. http://www.math.kit.edu/ianm3/lehre/geonumint2009s/media/gni_by_stoermer-verlet.pdf Plus actuel, cet article de 2014 d'Hairer http://www.unige.ch/~hairer/preprints/hairer-pisa.pdf Il explique qu'on peut distinguer trois cas, et que ce que j'ai dit plus haut est bien compris théoriquement dans les deux premiers cas, et observé mais pas encore démontré dans le dernier cas (celui où le produit $h\omega$ n'est pas petit, $h$ étant le pas de la méthode symplectique, et $\omega$ la plus haute fréquence du système).

    Et sinon je te conseille vraiment de ne pas rester sous scilab, il te faut un truc rapide en C++ ou Julia http://docs.juliadiffeq.org/latest/solvers/dynamical_solve.html et pour un background https://hal.archives-ouvertes.fr/cel-01830248/document
  • Merci pour la réponse !
    Et oui c'est bien parce que ce sont des intégrateurs symplectiques que je m’intéresse à eux.
    J'utilise scilab juste pour un benchmark rapide (en fait je l'ai aussi fait en C++ mais avec les mêmes résultats). Le code utile à la fin est en c++ (GEANT4 est un framework C++ Monte Carlo pour les simulations en physique des particules).

    Mais mon interrogation initiale provenait du fait que sur mon benchmark RK4 était le meilleurs (d'un point de vu conservation de l’énergie vs temps de calcul). Mais en fait il y avait des erreurs... (plus d'un ans que ce truc traîne dans un coin à cause d'une erreur de code...).

    Par contre j'ai l'impression qu'il y a une taille de pas (extrêmement petite) à partir de laquelle la différence est moins marquée (ou alors sur des temps de simulation extrêmement longs).
    Globalement mon probleme est résolu donc. Merci !

    Question bonus.
    Actuellement j'utilise l'algo PEFRL (ici : https://arxiv.org/pdf/cond-mat/0110585.pdf ) qui est plus efficace que le Yoshida ordre 4 https://en.wikipedia.org/wiki/Leapfrog_integration#cite_note-Yoshida1990-6 . Y a-t-il des endroits où on peut trouver des genres d'inventaires plus ou moins exhaustifs de méthodes avec leur "efficacité" ? Y a-t-il un intérêt à utiliser des ordres 6, 8 ou plus ?
  • Je ne suis pas expert dans le domaine, mais bon à la base une méthode d'ordre plus élevé permet de faire moins de pas pour atteindre une limite $T$ donnée.

    Cet article récent https://ergodic.org.uk/~bl/Data/Preprint/mist.pdf fait un survol des gros logiciels utilisés en dynamique moléculaire et propose un cadre pour développer des méthodes numériques en C++ qui s'interfacent bien. Il y a des listes d'intégrateurs, des graphes de comparaison, etc. Cela pourrait bien t'être utile.

    J'avais utilisé à une époque le fameux DOP853 (méthode d'ordre 8 à pas variable avec continuous output, c'est à dire que ça te fournit une interpolation de bonne qualité pour tout temps entre tes pas) et que j'avais trouvé très fiable. Il est en C dans cette archive tar de E.Hairer http://www.unige.ch/~hairer/prog/nonstiff/cprog.tar
Connectez-vous ou Inscrivez-vous pour répondre.