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

[invite88ef51f0]

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]

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).

[A voir en vidéo sur Futura]

[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 ?

[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

[invite88ef51f0]

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...

[invite88ef51f0]

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 .

[invite88ef51f0]

à 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 ?

[invite88ef51f0]

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.

[invite88ef51f0]

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 .