Euler/Runge-Kunta: Comment les implémenter

Bonjour tout le monde !
je cherche à implémenter, afin de modéliser le comportement de particules chargées, les méthodes d'Euler dans un premier temps, puis de Runge-Kunta. Bien sûr, si je pouvais implémenter la seconde directement, ce serait mieux...
J'ai trouvé un tas de documents qui expliquent Runge-Kunta (RK), mais je n'arrive pas à en déduire une implémentation :x .
Mon but est d'utiliser l'équation fondamentale: Somme des forces=m.accélération .
Donc, en calculant F (sur x, y et z), je peux en déduire a. Ca, ça peut encore aller. Par contre, pour connaître la position de la particule à l'instant t+dt (connaissant sa position à l'instant t), il me faudrait utiliser RK, et là, je bloque :? .

Si vous pouviez méclairer !
Merci d'avance.

Les "Numerical Recipes" sont ma bible dans ce domaine.

Outre leurs très bonnes explications (en anglais, par contre), elles donnent des implémentations qui m'ont parfois permis de comprendre ou appliquer des méthodes que mon (faible) baggage mathématique ne me permettait pas d'appréhender formellement.

L'ouvrage est disponible en librairie spécialisée, mais aussi en ligne, ici pour la version avec sources en C (lire les mentions légales, d'usage et de licence, notamment).

Voir plus précisément les détails de Runge-Kutta qui m'ont servi dans une autre vie d'informaticien scientifique ;)

Bonjour,
je ne sais pas si le langage est imposé, mais sur Maple (langage de calcul formel) les calculs et affichages se font tres facilement.
Par ailleur, tu doit meme pouvoir trouver les sources de ce que tu cherches.
Bon courage.

La méthode de Runge-Kutta (d'ordre 4 par exemple, qui est très facile à implémenter) permet de résoudre des équations différentielles ordinaires d'ordre 1. Le problème est que l'équation fondamentale de la mécanique F = m.accel est d'ordre 2 (à cause de l'acceleration qui est une dérivée seconde). Il faut donc commencer par transformer l'équation en équation d'ordre 1, et pour cela il suffit de remplacer la fonction inconnue (disons f) par le couple (f,f'). On réduit ainsi le degré à 1, mais en contre partie les solutions sont à valeurs dans un espace de dimension double (6 si f prend ses valeurs dans R^3).

Tout comme une réponse précedente, le conseille vivement les
Numerical Recipes et surtout le livre qui l'accompagne.

Merci pour vos réponses !
DrTopos: pourrais-tu m'expliquer un peu plus ?
Si j'applique RK4 sur mes coordonnées x, puis y (puis z), je n'aurais pas le bon résultat ?

tout comme pour le méthode d'Euler, ces techniques résolvent un système du 1er ordre
on résoud y' = f( x,y)
avec y(t+dt) = y(t) + dt* f(t)+ O(dt^2) [ pour Euler ]
ou une expression en 4 membres pour Runge-Kutta
( y(t+dt) = y(t) +dt * [ k1/6+ k2/3+ k3/3+ k4/6 ] + O(dt^5) )
si l'équation est d'1 ordre + élevé, la technique standard est la suivante
exemple ordre 2
on a
Y''. A + Y' .B + C.Y + D =0

on note U=Y et V=Y'

alors on a

U' = V
et
V' = -1/A * ( V*B + C*U+D)

on baisse d'un ordre en doublant le nombre d'équations
Ce système de 2 équations ( ou plus ) se résouds alors avec un méthode standard . on a les solutions U(t) et V(t) soit Y(t) et Y'(t) en //.

Il arrive que des changements de variables puissent serieusement simplifier la solution.

Les détails viennent d'être donnés par j.p.mignot. Pour résumer la méthode en deux mots: il faut prendre la position (x,y,z) et la vitesse (x',y',z') comme inconnues (qu'on renomme par exemple (u,v,w)) et non pas la position seulement, pour abaisser le degré de l'équation. Il faut donc appliquer RK4 à 6 coordonnées et non pas 3.

Sur ma page enseignement, il y a un document qur les équations différentielles. Je n'y parle pas de RK4, mais il contient des informations d'intérêt général sur le sujet.

Bon, bon, bon...
C'est bien ce que je pensais, mon problème, c'est de trouver les équa diff qui régissent le système...
J'ai un système avec deux particules chargées (charge 1.6e-19C), de charges opposées (elles devraient s'attirer) et j'ai un peu n'importe quoi.
Pouvez-vous m'aider ?
Je cherche l'équation de (x,y,z) pour chaque particule (bien sûr !).
Je pense que si j'ai x j'ai les autres, mais là...

Ah Ah ! Il ne s'agit pas d'un point matériel, mais de deux points matériels. Il faut donc trois coordonnées (x1,y1,z1) pour le premier point, et trois autres (x2,y2,z2) pour le deuxième point. Il nous faut aussi trois coordonnées (u1,v1,w1) poir le vecteur vitesse diu premier point et (u2,v2,w2) pour le vecteur vitesse du deuxième point. Au total on a donc 12 fonctions inconnues du temps.

On a déjà bien sûr les 6 équations:

Code:

---

1 2 3 4 5 6 7

x1'(t) = u1(t) y1'(t) = v1(t) z1'(t) = w1(t) x2'(t) = u2(t) y2'(t) = v2(t) z2'(t) = w2(t)

---

Les autres équations vont être données par la loi de la mécanique, en faisant attention d'écrire les accélerations comme des dérivées premières de u1,v1 etc... (sans jamais écrire de dérivées secondes):

Code:

---

1 2 3

u1'(t) = .... etc....

---

On se retrouve alors avec un système différentiel linéaire de 12 équations qui est du premier ordre et qu'on peut traiter par RK4.

pour compléter la réponse ci-dessus, en plus des 12 équations (Xi(t),Yi(t),Zi(t), X'i(t),Y'i(t) et Z'i(t) i=1 et 2 il faut de plus impérativement définir les conditions initiales du problème tant sur les Xi,Yi,Zi que sur leurs dérivées.
Mathématiquement toute intégration fait intervenir des "constantes" dont les valeurs peuvent changer drastiquement les résultats. Leurs affectations sont liées aux conditions initiales et le jeu d'équations différentielles à lui seul n'est pas suffisant pour résoudre le problème.
La force électrique d'attraction
(- 1/(4.pi.epsilon0.epsilon r) )* ( (q1*q2)/r^2 )
va définir en chaque instant l'accélération - donc la correction de vitesse - dans la direction du vecteur P1P2. les composantes des vitesses dans le plan orthogonal à ce vecteur restent elles inchangées. A chaque étape de l'intégration, la direction P1P2 change.
pour débuter le calcul il faut connaître P1(0), P1(0) ainsi que V1(0) et V2(0)
ce qui redonne les 12 conditions initiales aux 12 équations du 1er ordre précédemment citées.

Le fait de considérer les vitesses comme des fonction inconnues du temps a deux avantages:

(1) Cela permet d'appliquer RK4 en abaissant le degré de l'équation différentielle.

(2) RK4 va donner à chaque étape non seulement les positions des deux particules, mais aussi les vecteurs vitesse. Si tu veux faire une animation 3D par exemple, tu pourras très facilement représenter les particules ainsi que leurs vecteurs vitesse. Cela peut être intéressant pédagogiquement.

Merci à vous deux (on sent les profs de maths ;) ) !
Juste une dernière chose, et je me remets dans mon implémentation, la direction du vecteur P1P2 s'obtient bien par (X2-X1, Y2-Y1, Z2-Z1) (et inversement pour P2P1) ?
J'entends par là que la force F12(Fx1, Fy1, Fz1) de P1 sur P2 s'exprime :

Citation:

---

Avec k=1/(4*PI*Eps0)*(q1*q2)/(r^2),
Fx1=k*(X2-X1)
Fy1=k*(Y2-Y1)
Fz1=k*(Z2-Z1)

---

On a alors:

Citation:

---

u'1(t)=Fx1/m
...

---

Non ?

Sujet interessant étant donné que je souhaiterais utilisé également Somme des forces=m.accélération pour calculer mon deplacement, mon problème se passe dans un univers 2D, donc il sera moins complexe mais je suis content d'avoir trouvé autant d'explication à la fois sur ce thread et par l'intermediaires des liens/ressources liées.

voud avancez

Citation:

---

Avec k=1/(4*PI*Eps0)*(q1*q2)/(r^2),
Fx1=k*(X2-X1)
Fy1=k*(Y2-Y1)
Fz1=k*(Z2-Z1)

---

NON Fx = k * (X2-X1) / |P1P2|
ou alors

Fx=k * (X2-X1) avec k = k=1/(4*PI*Eps0)*(q1*q2)/(r^3)

!! idem pour les autres composantes !!

Oki doki !
J'ai réussi à avoir mes 12 équations (simplissime en fin de compte), j'ai maintenant à peu près compris le principe, et je pourrais (sous peu) implémenter tout ça.
Quelles méthodes permettraient de résoudre directement des systèmes d'ordre 2 (ou plus), si elles existent ?

Citation:

---

Envoyé par progman

Quelles méthodes permettraient de résoudre directement des systèmes d'ordre 2 (ou plus), si elles existent ?

---

Autant que je sache, on est obligé de se ramener au premier ordre en considérant la vitesse comme une fonction inconne. Si l'équation est d'ordre 3, il faut faire de même avec l'accelération etc... Je ne connais pas d'autre méthode. Mais c'est facile à faire donc pas un problème.

Oui, c'est facile en effet, sans optimiser, j'ai un programme qui me résoud une équation du second degré en passant par deux rk4 (un par équation du premier degré) et je retrouve la théorie (avec de minimes erreurs). J'ai essayé avec l'équation :
d²x/dt²=1 (ce qui conduit, je crois, si les CI sont nulles, à x=t²/2). Eh bien j'ai bien t²/2 en sortie du calcul :wink: .
J'adore :D .

Je vous conseille d'aller sur:
http://library.lanl.gov/numerical/bookcpdf.html
et d'en downloader les pdf du chapitre 16 qui traitent de la résolution des ODE
Une partie importante traite de RK4 et des techniques de pas adaptatifs d'intégration, d'autres approches sont aussi proposées.

Oui, j'ai lu cela hier soir, merci !
Il existe beaucoup de méthodes, mais, si j'ai bien compris, on se ramène toujours à des équations du 1er ordre.
Maintenant, repartons dans l'implémentation ;) .

Merci à vous deux.

le site
http://www.sciences.univ-nantes.fr/physique/perso/aloui/m_numeri/61eqdiff/61eqdiff.htm
propose au chapitre 2.6 une méthode directe pour résoudre
A(x) y'' + B(x) y' + C(x) + D(x) = 0
je n'ai pas testé mais probablement en ordre.

Pfiou, j'ai mis une journée à comprendre (à peu près) le principe...
Je ne suis pas prêt d'implémenter ça, déjà que RK4 dans mon problème s'avère complexe à vérifier...