Mouvement de particules

[invite77a86db2]

Bonjour,
Sous un titre générique (j'ai mis 5 minutes à me décider ), j'ai une série de questions...
Je cherche à créer un modèle assez simple de particules (pas de modélisation trop proche de la réalité). En l'occurrence, elle sont chargées, ont une masse, mais c'est à peu près tout.
J'ai donc une matrice NxN de particules élémentaires, et je cherche à appliquer les lois gravitationelle et électromagnétique.
Là, je bloque complètement, car comment calculer la résultante, pour chaque particule, des forces appliquées, de manière à connaître la position future ?
Par un vecteur (ou une matrice) particulier(e) ?

Merci de m'aider !
-----

[yat]

Je ne suis pas sur de bien comprendre la question, mais la méthode bourrin devrait bien marcher :
Chaque particule est associée à un vecteur position, un vecteur vitesse et un vecteur accélération/force.
A chaque itération tu remets les vecteurs acceleration/force de toutes tes particules à zero.
Pour chaque particule i, tu parcours le tableau des particules de 0 à i-1, pour chacune tu calcules la distance entre les deux particules, et les vecteurs force réciproques opposés que tu ajoutes dans le vecteur force des deux particules.
Une fois que c'est fait, les vecteurs force contiennent pour chaque particule le bilan des forces extérieures, tu transformes ça en accélération, tu modifies la vitesse en fonction puis tu modifies la position en fonction, tu dessines et tu recommences.

J'applique cette méthode quand il n'y a pas trop de particules (pas plus de quelques milliers), et le temps de calcul reste correct.

Si tu veux accélérer les choses ou taper dans les centaines de milliers de particules, tu découpes ton espace en cubes, chaque cube contient une liste de particules, et tu ne testes pour chaque particule que les interactions avec les particules des cubes voisins. C'est plus rapide mais plus compliqué et c'est une approximation.

Si quelqu'un connait des solutions plus efficaces, je suis preneur aussi

[Coincoin]

Salut,
Effectivement, la méthode bourrine marche. La complexité doit être en n², mais pour n raisonnable, ça passe très bien.
Tu peux raffiner un peu en ne calculant qu'une fois chaque force (quand tu considères la particule i et que tu regardes la force exercée par la particule j, tu peux en profiter pour dire que la force que j subit de i, c'est l'opposé). Du coup au final, t'as une boucle du genre for i from 1 to N-1 do for j from i+1 to N do ... end end. C'est toujours ça de gagné.

Sinon, il faut voir la concentration de ton système par rapport à la portée des forces, mais se limiter aux plus proches voisins peut énormément alléger le calcul.

Une fois que tu as les accélérations, pour calculer la nouvelle position, essaye peut-être de faire un peu plus poussé que la méthode d'Euler (ordre 1) : Runge-Kutta 2, ...


Sinon, c'est quoi exactement ta matrice NxN ? Tu fais ça en combien de dimensions ?

Et pou l'affichage, tu sais comment tu vas faire ?

[invite77a86db2]

Mmh...
Je comprends, en gros !
As-tu un exemple de code (matlab ou autre) ?
Je code en Python pour ma part.
Les vecteurs position, vitesse et accélération/force sont des scalaires associés aux particules ?
Tu parles d'un tableau de particules, c'est quoi ?
Un vecteur contenant chaque particule ?

Voilà, pour résumer, un bout de code, je suis preneur .
Merci de ton aide

[A voir en vidéo sur Futura]

[invite77a86db2]

Sinon, il faut voir la concentration de ton système par rapport à la portée des forces, mais se limiter aux plus proches voisins peut énormément alléger le calcul.

J'y ai pensé, mais quelles forces doivent être prises en compte ?
Quelle est la force prédominante à une distance faible ?

Une fois que tu as les accélérations, pour calculer la nouvelle position, essaye peut-être de faire un peu plus poussé que la méthode d'Euler (ordre 1) : Runge-Kutta 2, ...

Là, j'ai rien compris .

Pour ma matrice NxN, c'est une représentation du plan des particules, en considérant qu'une particule occupe une case.
Ce qui me limite rapidement pour l'affichage, car du coup, au mieux, une case=1pixel, donc pas plus de 600x600 pour rester assez visualisable (et encore).

[yat]

Tu peux raffiner un peu en ne calculant qu'une fois chaque force (quand tu considères la particule i et que tu regardes la force exercée par la particule j, tu peux en profiter pour dire que la force que j subit de i, c'est l'opposé). Du coup au final, t'as une boucle du genre for i from 1 to N-1 do for j from i+1 to N do ... end end. C'est toujours ça de gagné.

C'est bien ce que j'ai fait

Une fois que tu as les accélérations, pour calculer la nouvelle position, essaye peut-être de faire un peu plus poussé que la méthode d'Euler (ordre 1) : Runge-Kutta 2, ...

Ca ça m'intéresse... je me débrouille avec google ou tu peux me résumer ce que c'est ?

[Coincoin]

Je suis loin d'être un spécialiste...
Personnellement, dans le cadre de mes études, je fais un petit stage dans un labo où je fais une modélisation en Fortran de particules colloïdales. Elles sont soumises à la force magnétique et à la force électrostatique et y a du mouvement brownien, mais au fond, ça correspond exactement à ce que veut faire Progfou.

Le truc le plus simple, c'est de faire un tableau x de dimension N, regroupant les abscisses des particules, un autre y pour les ordonnées (et éventuellement un troisième pour z). Tu fais aussi des tableaux pour la vitesse : vx et vy. Pour les forces, tu fais un tableau Fx (taille N) des forces selon x, et un Fy. Et si les masses et les charges ne sont pas les mêmes, tu peux faire un tableau m et un tableau q...
Ensuite, à chaque étape, tu calcules les forces (tu calcules la distance au carré d2=(x(i)-x(j))²+(y(i)+y(j))², tu calcules la force correspondante, tu calcules la projection selon l'axe qui t'intéresses), faut juste faire un peu de maths avant pour trouver les bons calculs.

Pour ce qui est de l'intégration des équations du mouvement, tu as . La méthode la plus simple (et la plus intuitive) consiste à dire à dire que et . C'est ce qu'on appelle la méthode d'Euler. Mathématiquement, c'est juste si dt est infiniment petit. En pratique dt n'est pas nul, c'est juste une approximation au premier ordre. Y a moyen de faire un peu plus juste en compliquant le calcul. Par exemple au deuxième ordre, on a : .
Pour ce qui est des méthodes de Runge-Kutta, elles utilisent en fait plusieurs points précédents au lieu d'un seul (Runge-Kutta 2 utilisent deux points, ...). Là où Euler suppose que la pente est constante entre deux points, RK fait une moyenne (pondérée) des dernières pentes.
C'est un peu du raffinement, je vous l'accorde, mais pour donner un exemple concret, si on modélise un oscillateur harmonique (un ressort) avec Euler, les oscillations augmentent avec le temps. Alors que cet artefact numérique n'existe pas avec RK.
Mais pour un modèle ne cherchant pas à coller à tout prix à la réalité, on peut s'en passer...

[invite77a86db2]

Donc, je calcule a à partir de la somme des forces appliquées.
Puis, v, puis enfin, x.
Mais c'est uniquement l'abscisse ?
J'ai de très vagues souvenirs de méca du point, mais où retrouver ça...

[Coincoin]

J'ai séparé les composantes. Donc tu calcules Fx et Fy, puis tu intégres pour obtenir vx et vy, et enfin x et y.

Pour ce qui est de l'aspect mathématique, essaye de trouver des cours de physique niveau 1ère / Tale S par exemple.

[invite77a86db2]

En fait, je m'en suis souvenu 5 minutes après avoir tappé le message !
Mais j'ai retrouvé mon cours de méca du point .
Vais creuser !

[yat]

C'est un peu du raffinement, je vous l'accorde, mais pour donner un exemple concret, si on modélise un oscillateur harmonique (un ressort) avec Euler, les oscillations augmentent avec le temps. Alors que cet artefact numérique n'existe pas avec RK.
Mais pour un modèle ne cherchant pas à coller à tout prix à la réalité, on peut s'en passer...

Coller à la réalité est une chose... le problème avec ma méthode d'Euler (tu viens de m'apprendre que c'était la méthode d'Euler d'ordre 1 ), c'est qu'à long terme mes planètes finissent par se barrer et échapper à l'attraction de leur étoile

Je vais m'empresser de tester cette méthode de Runge-Kutta. Merci !

[invite77a86db2]

Oki, donc, si je comprends bien, dans la méthode expliquée ici :
http://en.wikipedia.org/wiki/Runge-Kutta_methods

Le h c'est ton dt.
Je vais voir ce que ça donne en éclatant le calcul, mais ça me parle un peu plus .

[Coincoin]

à long terme mes planètes finissent par se barrer et échapper à l'attraction de leur étoile

Je pense que RK devrait régler ce problème (cf ce que je dis sur l'oscillateur harmonique).
J'avais essayé de simuler le système Soleil-Terre-Lune, mais la Terre se précipitait sur le Soleil, suivie par la Lune. La fin du monde avait lieu dans quelques semaines à peine... Faudra que je réessaye à l'occasion.

Le h c'est ton dt.

Oui, c'est le pas de temps.

[invite77a86db2]

Donc, il faut faire dérouler un temps...
Sinon, c'est un temps discret qu'il faut faire intervenir, avc dt=Te la période d'échantillonnage ?
Si tel était le cas, je pourrais exploiter la transformée en Z, et sortir l'expression de mon abscisse (et ordonnée) en fonction des temps précédents et de l'accélération ?

[Coincoin]

Il faut effectivement un temps discret. Un temps continu demanderait une résolutionde s équations, pas une simulation.
Pour ce qui est de la tranformée en z, j'en ai déjà entendu parler, mais je maîtrise pas vraiment. Mais je ne vois vraiment pas ce que ça pourrait t'apporter de plus que les équations de base.

[invite77a86db2]

Si tu le fais par TZ, tu as ceci :
x(n)=int(int(a(n))) avec int l'intégrale.
=>X(z)=(1/(1-z^-1)^2)A(z)

Donc, après développement :
X(z)=A(z)-2*z^-1X(z)+z^-2X(z)

Donc, en repassant en n:
x(n)=a(n)-2x(n-1)+x(n-2)

Or, x(n-1) et x(n-2) sont justes les positions antérieures.
Donc, ils sont stockables, et a(n) est l'accélération.

[Coincoin]

Euh... ça me paraît louche comme résultat : le pas de temps n'intervient pas, si tu prends a=0 et x(n-1)=x(n-2) (différents de 0), tu as quand même un mouvement, ...

Tu es peut-être plus à l'aise avec la transformée en z qu'avec les dérivées, mais je pense sincérement que ça ne fait que compliquer le problème.

[invite77a86db2]

Non, n=nTe où Te est la période d'échantillonage.
Par convention, on écrit n, c'est plus rapide .

[yat]

Je me suis un peu documenté sur la méthode RK d'ordre 2, plus précisément la méthode de Heun (elle m'a l'air relativement simple, et les résultats obtenus sont hallucinants par rapport à mon vieil Euler)

A ce que j'en ai compris, on a une fonction y(t) qu'on ne peut pas exprimer directement et une fonction f telle que y'(t)=f(y(t),t), et on cherche à approcher numériquement la valeur de y(t+h) à partir de y(t) et de f.

Le principe est apparemment de calculer une valeur k1 par la méthode d'Euler (c'est à dire k1=y(t)+hf(y(t),t) ), et une valeur k2 par une méthode qui ressemble, mais en utilisant la valeur de f en k1. (k2=y(t)+hf(k1,t+h)). On garde ensuite la moyenne de ces deux valeurs. J'ai donc trois questions :

-est-ce que je suis complêtement à coté de la plaque ?

-il n'est apparemment pas question ici de tenir compte des deux valeurs précédentes. Mais en regardant un peu ce qui se passe, et avec un pas deux fois plus grand, il me semble que ça revient à peu près à chercher y(t+2h), de prendre k1=y(t+h)+(y(t+h)-y(t)), calculer la pente en k1 et faire directement y(t+2h)=y(t+h)+hf(k1,t+2h). J'ai bon ?

-on parle ici de dérivées, moi j'ai des dérivées secondes. Est-ce qu'il faut appliquer une fois la méthode pour approximer la dérivée, puis une deuxième fois pour approximer la valeur à partir de cette dérivée, ou bien directement calculer mes k1 et k2 en appliquant deux fois la méthode d'Euler, puis faire la moyenne ?

[invite77a86db2]

Comment implémenter ces fonctions ?
Je suis complètement bloqué !
Matlab ou autre chose, si vous avez des exemples de code, ça m'aiderait à y voir plus clair...

[hterrolle]

tu utilise quel outil pour réaliser ta programation 3D ?

[Coincoin]

Yat, pour ce qui est de Runge-Kutta, je m'en servi quelques fois mais je ne m'en souviens pas bien. Je ne peux pas te dire si ce que tu racontes est juste ou pas.

-on parle ici de dérivées, moi j'ai des dérivées secondes. Est-ce qu'il faut appliquer une fois la méthode pour approximer la dérivée, puis une deuxième fois pour approximer la valeur à partir de cette dérivée, ou bien directement calculer mes k1 et k2 en appliquant deux fois la méthode d'Euler, puis faire la moyenne ?

La première solution, c'est intégrer deux fois avec RK. La deuxième solution, c'est un Euler puis un RK. Ce que j'en pense, c'est que c'est dommage de gâcher avec un Euler. Mieux vaut faire 2 RK, ça sera toujours plus précis...

[invite77a86db2]

Je me permets de relancer ce sujet .
Si quelqu'un a une implémentation, ne serait-ce que d'Euler déjà, voir de RK, je prend, on simplement un algorithme en pseudo-langage ?
J'essaie de le faire, mais je n'y arrive pas, c'est compliqué de voir comment le faire...

Merci d'avance !

[Chup]

Bonjour,
vous trouverez ça dans les Numerical Recipes au lien suivant
http://www.library.cornell.edu/nr/

[zoup1]

Bonjour,

voici un exemple de code qui utilise Soit Euler, Runge et Kutta d'ordre 4 pour calculer le mouvement d'un oscillateur harmonique amorti en fonction du temps.

Il y a les sources en C
http://zoup1.free.fr/oscillateur.c
Et le fichier des paramètres.
http://zoup1.free.fr/in

C'est destiné à tourner sous un environnement avec ligne de commande genre Linux...
Cela devrait pouvoir aussi tourner sous DOS même si je n'ai pas testé cela.

[invite77a86db2]

Merci à vous deux, zoup1: je suis sous linux, donc pas d'inquiétude .

[zoup1]

cool...
A ce moment, les instructions pour la compilation, l'execution et l'affichage sont sommairement décrites dans les premières lignes du fichier.

Par ailleurs tu peux aller voir ce lien si tu as besion d'un peu de documentation sur les outils...
http://www-ph332.script.univ-paris7.fr/index.html

[invite77a86db2]

Je connais bien les outils .
Par contre, il y a un truc bizarre :
Tu dis que ton équation est de la forme :
m d²x/dt²= -k x - a dx/dt

Et dans la fonction accel (qui au passage m'a fait tout comprendre !!) tu as :
return -k/m*x - a*v;

Il ne te manque pas un m dans a*v => a*v/m ?

Et puis, tu ne travailles qu'avec x, dans ton implémentation. Ca se fait de la même manière pour connaitre les nouvelles coordonnées y et z ?
Je dois seulement trouver les équa diff qui régissent mon système pour chaque coordonnée ?

[zoup1]

Je connais bien les outils .
Par contre, il y a un truc bizarre :
Tu dis que ton équation est de la forme :
m d²x/dt²= -k x - a dx/dt

Et dans la fonction accel (qui au passage m'a fait tout comprendre !!) tu as :
return -k/m*x - a*v;

Il ne te manque pas un m dans a*v => a*v/m ?

Oui, tu as raison, c'est une erreur...

Et puis, tu ne travailles qu'avec x, dans ton implémentation. Ca se fait de la même manière pour connaitre les nouvelles coordonnées y et z ?
Je dois seulement trouver les équa diff qui régissent mon système pour chaque coordonnée ?

Oui, cela se fait de la même façon avec les autres coordonnées... Simplement lorsque l'on commence a avoir plusieurs équations différentielles à résoudre, on a interêt mettre les différentes variables dans un tableau et à faire des boucles. Pour cela tu peux regarder ce qui est fait dans Numerical Recipes (le lien qu'a donné Chup).

[yat]

Zoup1, j'ai des problèmes avec ton algo...

Je m'intéresse particulièrement à la méthode de Heun, qui est apparemment celle que tu appelles méthode du point millieu.

Mon souçi est le passage à la dérivée seconde d'une méthode qui est présentée partout avec la dérivée première. J'ai donc fait un test sous excel, une particule en orbite autour de l'origine, dans un espace 2D avec une froce centripète en 1/d².

A chaque étape, je calcule donc la nouvelle vitesse en lui appliquant linéairement l'accélération, puis j'applique linéairement à la position du point la moyenne de l'ancienne et de la nouvelle vitesse.

Au bout d'une dizaine de milliers d'opérations, la distance entre la particule et l'origine a été multipliée par 4 alors qu'avec une méthode d'Euler brutale elle a été multipliée par 5...

Cette amélioration ne me semble vraiment pas trancendante, d'autant plus qu'en faisant des essais en travaillant sur la dérivée première, j'obtiens un joli cercle là ou Euler me donne une lamentable spirale.