Erreur de précision ou incompréhension physique

Bonjour,

Je suis actuellement en train de tester la programmation fonctionnelle réactive et pour cela, je réalise un mini moteur physique. L'idée de base étant de simuler un problème à N corps.

Actuellement, si je met deux objet immobile à une distance raisonnable l'un de l'autre et que je lance la simulation, les deux objets se rapprochent petit à petit, de plus en plus vite (jusque là, tout me semble normal), mais lorsqu'ils se rencontrent (il n'y a pas de collision), ils sont tous les deux fortement accélérés. Au final, à partir de deux objets immobiles et à une certaine distance, j'obtiens deux objets avec un grande vitesse et une distance assez vite beaucoup plus grande que celle d'origine.

J'en déduis donc qu'on a une (forte) augmentation de l'énergie interne du système, ce qui n'est, il me semble, pas physiquement très réaliste.

À supposer que mon interprétation est bonne, on peux supposer que le problème viens de la méthode de résolution de l'équation différentielle du système. A priori, c'est la méthode d'Euler que j'utilise : le calcul des intégrales se fait en ajoutant la valeur à intégrer multiplié par la différence de temps. Je sais que la méthode d'Euler peut propager assez facilement des erreurs, est-ce que c'est bien ça qui se passe? Enfin, est-ce qu'utiliser une méthode telle que Runge-Kutta permet de résoudre ce problème ou est-ce que ça ne ferais que le réduire?

Salut Berdes !

Est ce que tu pourrais être plus précis, car j'ai du mal à comprendre ce que tu fais. Ton programme résout quelles équations? Tu utilises la méthode d'Euler pour le calcul de quelle intégrale ? Si j'ai bien compris, tu résous numériquement une équa diff, donc c'est la dérivée que tu approximes non ?

Avant de passer à Runge-Kutta (c'est pas très compliqué à implémenter ), tu peux déjà faire plusieurs test en changeant ton pas de temps et voir si ça te donne les mêmes résultats où pas.

Enfin, je ne comprends pas la phrase "Une distance assez vite beaucoup plus grande que celle d'origine", il me semblait que tes objets se rapprochent.

Je résous le problème à N corps. Donc je calcul l'accélération des objets à partir des positions, à partir de là, la vitesse et la position s'en déduisent en intégrant.

J'ai testé avec un pas de temps de 0.01s et j'ai toujours le problème. J'ai testé avec deux objets de 1000T chacun séparés de 1m et ils se rejoignent en une centaine de seconde. Au moment où ils se rejoignent, ils recommencent à s'éloigner (normal), mais beaucoup plus vite qu'ils se rapprochaient et se retrouvent à la position d'origine de l'autre en même pas une seconde.

En gros, ça donne :

1
2
3
4
5
6

  0s   A                  B
 10s     A              B
...
 95s            AB
 97s            BA
 98s      B            A

Et ensuite, les deux objet s'éloignent très vite.

Ça ne m'étonne pas plus que ça. L'énergie de ce genre de système est souvent en 1/r, donc les objets vont acquérir beaucoup d'énergie en se rapprochant. Ce faisant, ils risquent, une fois trop proche, d'avoir assez d'énergie pour s'éloigner (ce qui donnerait avec une résolution avec précision infinie un ressort (dissipatif ?)). Avec les erreurs de précision, le ressort devient aberrant.

Si c'est effectivement ça, on peux soit augmenter la précision (avec un pas plus petit), soit ajouter un piège (les objets deviennent collés ou détruits lorsque la distance qui les sépare est trop faible).

Oui, c'est vrai que si les objets sont parfaitement superposés, on a une accélération infinie, donc le problème viens peut-être de là.

En revanche, ça m'étonne plus que le phénomène se reproduise (à moindre mesure) si les objets ne rapproche pas trop. J'ai enregistré un exemple : https://youtu.be/BQniVRGFPII . On voit bien dessus que les objets partent avec une légère vitesse et une certaine distance et qu'ils finissent avec une vitesse supérieur et une distance plus grande aussi, ce qui veux dire une plus grande énergie interne.

Quand j'ai vu ça, j'ai tout de suite pensé à l'assistance gravitationnelle, ou effet fronde. Les astrophysicien utilisent la gravité des planètes pour augmenter la vitesse des machines qu'ils envoient dans l'espace.

Pour l'énergie, je ne sais pas trop. As-tu la possibilité de la tracé ? C'est juste l'énergie cinétique ($mv^2/2$) pour chacun des corps et l'énergie gravitationnelle ($-Gm_1m_2/r$).

Édit : l'énergie gravitationnelle écrite ainsi ne doit être compté qu'une fois (la formulation précédente n'était pas bonne).

Oui, y'a clairement un problème avec l'effet « ressort ». Dans ta vidéo, très légèrement avant la sixième seconde, y'a un instant d'hésitation à partir duquel les objets devraient recommencer à se rapprocher l'un de l'autre, puis se ré-éloigner moins loin, etc. Là, la trajectoire continue.

Vu la configuration de départ, je suis quasiment sûr que les deux objets devraient avoir une trajectoire elliptique.

Je vais essayer de voir si ça change quelque chose de faire la simulation avec RK4 plutôt qu'Euler.

Au niveau de l'affichage de l'énergie du système, l'énergie cinétique me pose pas problème. En revanche l'énergie potentielle lié à la gravité me pose problème. Je n'ai pas de soucis pour la calculer dans le cas d'un champs gravitationnel constant, mais c'est très loin d'être le cas dans mon exemple. J'ai trouvé cette formule sur Wikipedia : , mais je trouve ça étrange d'avoir une énergie négative (après, ça permet de diminuer l'énergie cinétique quand l'énergie potentielle augmente, donc ça s'explique aussi).

mais je trouve ça étrange d'avoir une énergie négative

Un système tend à minimiser son énergie. Avec un terme inversement proportionnel à la distance, l'énergie gravitationnelle devient $+/-\infty$ pour une distance nulle. La distance nulle doit être un minimum d’énergie gravitationnelle, donc l'énergie gravitationnelle doit être négative. On ne le dira jamais assez : vérifier ce qui se passe aux limites, c'est le bien.

Sinon, la formule est exacte. Pour chaque corps, il faut sommer l'énergie gravitationnelle due à chacun des autres corps. Ce qui, dans le cas à 2 corps, donne juste -Gm₁m₂/r (pour tout le système, donc à ne compter qu'une fois ; si on veut le compter une fois par corps, il faut mettre le 1/2).

Désolé pour les temps de réponse, je passe vraiment pas beaucoup de temps sur ce projet, d'autant plus que j'ai mal choisi ça lib graphique (OpenGL, c'est bien trop bas niveau), du coup je perds aussi du temps à changer ce genre de choses.

Du coup, j'ai quand même eu des résultats sur l'énergie du système. Pour vous faire une idée des valeurs utilisés, les deux objets font chacun $10^8 kg$, ils démarrent à $1 m$ l'un de l'autre et avec une vitesse de $0.02 m.s^{-1}$.

L'énergie au départ est d'environ $-6273 J$, varie très peu au début mais quand les deux objets commencent à vraiment se rapprocher, elle diminue assez vite jusqu'à atteindre $-6800 J$ et elle augmente de manière fulgurante lorsque les objets recommencent à s'éloigner pour monter jusqu'à $10500 J$ et elle recommence à diminuer ensuite de plus en plus lentement en tendant vers un peu moins de $9780 J$.

Il y a donc bien un net soucis puisque cette énergie devrait être constante. Le seul point positif, c'est qu'elle varie de manière cohérente avec les observations : on a bien plus d'énergie après que les objets se soient croisés qu'avant.

La prochaine étape est de simuler avec RK4 et de voir s'il y a une amélioration.

Il y a donc bien un net soucis puisque cette énergie devrait être constante.

Pour ça, ça dépend de ce que tu as calculé, si c'est juste l'énergie gravitationnelle, elle n'a aucune raison d'être constante. C'est la somme énergie gravitationnelle + énergie cinétique qui l'est.

Par ailleurs je plussoie fortement Ryuuken.

@adri1 : je te rassure, c'est bien la somme $E_c + E_p$ que j'ai calculé.

@Ryuuken : je n'avais jamais entendu parlé des intégrateurs symplectiques. Je vais regarder de côté. Par contre, je ne sais pas si la modification sera si facile à faire que ça : j'utilise le Haskell et (en partie) la programmation fonctionnelle réactive qui est assez différent de ce qui se fait habituellement dans un langage impératif.

J'ai un peu avancé dans mes tests et j'ai réussi à corriger en grosse partie le problème.

La première chose que j'ai faite a été de simplifier/nettoyer mon code histoire d'avoir une meilleur maîtrise de ce que je fais.

La deuxième a été d'utiliser la méthode d'Euler semi-implicite. Le truc, c'est que cette méthode correspond exactement à ce que je pensais que mon code faisait à la base. Et en le modifiant encore un peu, j'en suis arrivé à la conclusion qu'il ne fait pas exactement ce que je pense. En fait, chaque intégrale, il y a un délai. J'ai donc le calcul de l'accélération qui se fait avec la position courante (pas de délai), mais la vitesse qui se fait sur l'accélération précédente et la position qui se fait aussi sur la vitesse précédente.

Du coup, au lieu que le changement de position s'impacte sur la nouvelle position dès le pas suivant, ça met deux pas actuellement. Je vais regarder comment faire pour éviter d'avoir ce problème, mais ça va nécessiter que je comprenne d'abord pourquoi l'intégration fait un délai.

Au final, j'ai supprimé le délai que j'avais mis de base (je pensais que l'intégrale n'en avais pas, donc je l'ai mis pour éviter une boucle infini) et j'ai aussi modifié le code pour que chaque frame fasse non plus 1 itération du système, mais 2000, ce qui m'a permis de diminuer le pas d'itération d'un facteur 1000. J'ai aussi passé les calculs de float à double, histoire de diminuer au maximum les erreurs liés aux approximations qui pourrait s'accumuler. Actuellement, à chaque croisement, l'énergie du système varie de l'ordre de 0.2%, ce qui est quand même une très nette amélioration des résultats.

Edit : c'est bon, j'ai repris les fonctions de la librairie que j'utilise pour faire une intégrale sans délai et la simulation donne de super résultats. L'énergie du système varie toujours un peu, mais elle ne fait qu'osciller et revient à la même valeur de manière cyclique, donc c'est vraiment moins problématique. Merci à vous pour votre aide et vos conseils.