Calcul de coordonnées force gravitationelle

[JPL]

Sûr qu'il vaut mieux commencer simple parce qu'à partir de trois corps tu risques de rencontrer rapidement le chaos, ce qui n'est pas favorable pour suivre le comportement du programme !
-----

[Azghar]

J'ai avancé l'algorithme, mais il y a encore un problème
J'aimerais sauvegarder une image de l'état du système tous les 10 itérations (ou plus) par exemple. Mais j'ai un message d'erreur : 'float' object has no attribute 'scatter'

from random import random
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

        
def mvt(n,m,L,N,tf):        #n nombre de particules, m leurs masses, L côté du cube qui contient les particules, N nombre d'itérations par particules (précision) tf durée totale simulation
    x=[]
    y=[]
    z=[]
    G=6.67384*10**-11
    #coordonnées aléatoires du nuage de particule initial
    for i in range(n):
        x.append(random()*L)
        y.append(random()*L)
        z.append(random()*L)
    h=tf/N
    t=0
    while t<tf:
        Lx=[]
        Ly=[]
        Lz=[]
        for i in range(n):
            xold,yold,zold,vxold,vyold,vzold=x[i],y[i],z[i],0,0,0
            for k in range(n):
                d=sqrt((x[k]-x[i])**2+(y[k]-y[i])**2+(z[k]-z[i])**2)
                if d>0:
                    ax=G*m*(x[k]-x[i])/d**1.5
                    ay=G*m*(y[k]-y[i])/d**1.5
                    az=G*m*(z[k]-z[i])/d**1.5
    
            #Création des listes de coordonnées temporaires
            Lx=[]
            Ly=[]
            Lz=[]
            #Mouvement de l'objet i:
            vxnew=G*m*ax*h+vxold
            vynew=G*m*ay*h+vyold
            vznew=G*m*az*h+vzold
        
            xnew=vxold*h+xold
            ynew=vyold*h+yold
            znew=vzold*h+zold
        
            Lx.append(xnew)
            Ly.append(ynew)
            Lx.append(znew)
                
            #changement de variable
            xold,yold,zold,vxold,vyold,vzold=xnew,ynew,znew,vxnew,vynew,vznew
        t+=h    
    #enregistrement image  
    ax.scatter(Lx, Ly, Lz, c=c, marker=m)
    plt.savefig("test.png")

Avez vous une idée ?

[Arzhur]

Bonjour,


Le message d'erreur parle de lui meme : tu tentes d'appeler la méthode/attribut "scatter" sur un objet de type float....qui n'a pas ce genre de methode/attribut.

[CM63]

Bonjour,

Effectivement, tu écrases l'objet ax à la ligne 30. Choisis un autre nom pour cet objet.

A plus.

[Azghar]

ah oui effectivement, ax était le nom d'une de mes variables.

Par contre ça ne marche toujours pas pour l'image : Il me dis qu'il ne connait pas c (dans scatter), il ne comprend pas que c'est la couleur.
Je ne sais pas utilisé ce module d'ailleurs j'ai utilisé du code que j'ai trouvé sur un tutoriel :

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

def randrange(n, vmin, vmax):
    return (vmax-vmin)*np.random.rand(n) + vmin

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
n = 100
for c, m, zl, zh in [('r', 'o', -50, -25), ('b', '^', -30, -5)]:
    xs = randrange(n, 23, 32)
    ys = randrange(n, 0, 100)
    zs = randrange(n, zl, zh)
    ax.scatter(xs, ys, zs, c=c, marker=m)

ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')

plt.show()

Bizarrement, là il comprend ce qu'est c, le code fonctionne bien.
Et encore pire : j'ai enlevé c dans le scatter (j'imagine qu'il va mettre une couleur par défaut) et maintenant il me dit que xs et ys ne sont pas de même taille , variables que je n'utilise pas dans mon programme, mais qui sont utilisées dans le code trouvé sur internet.
Je ne comprend vraiment pas. Moi je veux juste afficher un nuage de point en 3D, et j'ai les coordonnées de chaque point sous forme de 3 grandes listes.

[Arzhur]

ax.scatter(Lx, Ly, Lz, c=c, marker=m)

c et marker sont optionnels, t'as essayé de les enlever ?

Il me dis qu'il ne connait pas c (dans scatter), il ne comprend pas que c'est la couleur.
Je ne sais pas utilisé ce module d'ailleurs j'ai utilisé du code que j'ai trouvé sur un tutoriel :

oui mais ici c'est le second c qu'il ne connait pas : logique, tu ne l'as pas déclaré.
Attention avec tes variables marker=m --> m c'est ta masse !

Utiliser du code sans comprendre ce qu'on fait...ça donne souvent un résultat aléatoire. Essaye de jeter un coup d'oeil a la documentation de matplotlib ou de trouver un tuto simple sur le net.

[Azghar]

J'ai mis à jour le programme qui prend maintenant 3 corps en compte.
J'ai mis plusieurs système déjà prédéfini si vous voulez vous amuser chez vous. Attention cependant, c'est gourmand en calcul.

###Programme de simulation d'attraction gravitationelle entre trois corps
#Les masses sont considérées ponctuelles
from pylab import *
from math import*
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')


def orbite(M0,M1,M2,x0,y0,z0,vx0,vy0,vz0,x1,y1,z1,vx1,vy1,vz1,x2,y2,z2,vx2,vy2,vz2,ti,tf,N):
    Lx0=[x0]
    Ly0=[y0]
    Lz0=[z0]
    Lx1=[x1]
    Ly1=[y1]
    Lz1=[z1]
    Lx2=[x2]
    Ly2=[y2]
    Lz2=[z2]
    t=ti
    i=int(N/100000)
    k=0 #compteur d'opérations
    G=6.67384*10**-11
    x0old,y0old,z0old,vx0old,vy0old,vz0old,x1old,y1old,z1old,vx1old,vy1old,vz1old,x2old,y2old,z2old,vx2old,vy2old,vz2old=x0,y0,z0,vx0,vy0,vz0,x1,y1,z1,vx1,vy1,vz1,x2,y2,z2,vx2,vy2,vz2
    h=(tf-ti)/N
    
    while t<tf:
        d01=sqrt((x1old-x0old)**2+(y1old-y0old)**2+(z1old-z0old)**2)
        d02=sqrt((x0old-x2old)**2+(y0old-y2old)**2+(z0old-z2old)**2)
        d12=sqrt((x1old-x2old)**2+(y1old-y2old)**2+(z1old-z2old)**2)
        
        #Mouvement de l'objet 0:
        ax0=-G*M1*(x0old-x1old)/d01**3-G*M2*(x0old-x2old)/d02**3
        ay0=-G*M1*(y0old-y1old)/d01**3-G*M2*(y0old-y2old)/d02**3
        az0=-G*M1*(z0old-z1old)/d01**3-G*M2*(z0old-z2old)/d02**3
        
        vx0new=ax0*h+vx0old
        vy0new=ay0*h+vy0old
        vz0new=az0*h+vz0old
        
        x0new=vx0old*h+x0old
        y0new=vy0old*h+y0old
        z0new=vz0old*h+z0old
        
        if k%i==0: 
            Lx0.append(x0new)
            Ly0.append(y0new)
            Lz0.append(z0new)
        
        
        
        #Mouvement de l'objet 1:
        ax1=G*M0*(x0old-x1old)/d01**3-G*M2*(x1old-x2old)/d12**3
        ay1=G*M0*(y0old-y1old)/d01**3-G*M2*(y1old-y2old)/d12**3
        az1=G*M0*(z0old-z1old)/d01**3-G*M2*(z1old-z2old)/d12**3
        
        vx1new=ax1*h+vx1old
        vy1new=ay1*h+vy1old
        vz1new=az1*h+vz1old
        
        x1new=vx1old*h+x1old
        y1new=vy1old*h+y1old
        z1new=vz1old*h+z1old
        
        if k%i==0: 
            Lx1.append(x1new)
            Ly1.append(y1new)
            Lz1.append(z1new)
        
        
        #Mouvement de l'objet 2:
        ax2=-G*M1*(x2old-x1old)/d12**3+G*M0*(x0old-x2old)/d02**3
        ay2=-G*M1*(y2old-y1old)/d12**3+G*M0*(y0old-y2old)/d02**3
        az2=-G*M1*(z2old-z1old)/d12**3+G*M0*(z0old-z2old)/d02**3
        
        vx2new=ax2*h+vx2old
        vy2new=ay2*h+vy2old
        vz2new=az2*h+vz2old
        
        x2new=vx2old*h+x2old
        y2new=vy2old*h+y2old
        z2new=vz2old*h+z2old
        
        if k%i==0: 
            Lx2.append(x2new)
            Ly2.append(y2new)
            Lz2.append(z2new)
        
        
        #changement des variables:
        x0old,y0old,z0old,vx0old,vy0old,vz0old,x1old,y1old,z1old,vx1old,vy1old,vz1old,x2old,y2old,z2old,vx2old,vy2old,vz2old=x0new,y0new,z0new,vx0new,vy0new,vz0new,x1new,y1new,z1new,vx1new,vy1new,vz1new,x2new,y2new,z2new,vx2new,vy2new,vz2new
        t+=h
        k+=1
        
 
        
    #affichage
    ax.plot(Lx0,Ly0,Lz0)
    ax.plot(Lx1,Ly1,Lz1)
    ax.plot(Lx2,Ly2,Lz2)
    plt.show() 
    
    
        
orbite(1.989*10**30,5.9736*10**24,1.898*10**27,0,0,0,0,0,0,0,147098074*10**3,0,30.287*10**3,0,0,0,740520000*10**3,0,13.72*10**3,0,0,0,3600*24*365*12,100000000)
'''orbite(M0,M1,M2,x0,y0,z0,vx0,vy0,vz0,x1,y1,z1,vx1,vy1,vz1,x2,y2,z2,vx2,vy2,vz2,ti,tf,N)'''

###Exemples types :
'''Terre-soleil-Jupiter:1.989*10**30,5.9736*10**24,1.898*10**27,0,0,0,0,0,0,0,147098074*10**3,0,30.287*10**3,0,0,0,740520000*10**3,0,13.72*10**3,0,0,0,3600*24*365*12,1000000
Terre Lune Soleil :1.989*10**30,5.9736*10**24,7.3477*10**22,0,0,0,0,0,0,0,147098074*10**3,0,30.287*10**3,0,0,0,147098074*10**3+363104*10**3,0,1.052*10**3+30.287*10**3,0,0,0,3600*24*365,1000000'''
'''3 soleils:5.972*10**24,5.972*10**24,5.972*10**24,0,149597887.5*10**3,0,0,0,0,149597.5*10**3,149597887.5*10**3,0,0,0,0,0,149597887.5*10**3+363104*10**3,149597887.5*10**3,0,0,0,0,3600*24*365,1000000'''
'''Pluton Terre Soleil 1.989*10**30,5.9736*10**24,1.314*10**22,0,0,0,0,0,0,0,147098074*10**3,0,30.287*10**3,0,0,0,4436824613*10**3,0,6116.7,0,0,0,3600*24*365*247,10000000'''

Et voilà un exemple sur le système {Soleil,Terre,Jupiter}

[Chanur]

Joli ...
Et est-ce que ça marche si tu prends tes 3 corps dans le cas où les trajectoires sont chaotiques, c'est à dire des corps de masse proches et "également éloignés" (je veux dire que ce n'est pas chaotique si deux corps sont très proches et le troisième très loin) ?

[Azghar]

Oui biensûr,
voici par exemples 3 corps de masse solaire, avec des vitesses et une répartition choisies arbitrairement :
J'ai fais une simulation sur 10 ans et voici ce qu'on obtient:

[Azghar]

même simulation (même conditions initiales) sur 1 an maintenant

[Azghar]

Et voilà, après plusieurs heures à trouver et apprendre un bon module, j'ai finalement choisi Visual Python.
Et ça marche en 3D ! Simulation à n corps pour très bientôt

[Azghar]

Bonjour/bonsoir,
voici la version n corps 'figée' sans animation. Je travail sur l'animation avec matplotlib mais j'ai du mal (je suis entrain d'essayer de faire une classe)
Si vous voulez tester chez vous, attention la complexité est quadratique, le temps de calcul s'envole très vite.

# -*- coding: utf-8 -*-
"""
Created on Wed Feb 10 23:25:25 2016

@author: Jean 
jpl96@orange.fr
"""

from math import*
import numpy as np
import mpl_toolkits.mplot3d.axes3d as p3
import matplotlib.pyplot as plt
import random 


fig = plt.figure()
ax = p3.Axes3D(fig)


###Nbody sim using Eulerian method

def nbody(n,m,L,N,tf):
    x=[]    #listes contenant les
    y=[]    #coordonnees des points
    z=[]    #elles se mettent a jour a chaque boucle
    vx=[]   #idem pour les vitesses
    vy=[]
    vz=[]
    G=1
    #coordonnées aléatoires du nuage de particule initial
    for i in range(n):
        x.append(random.random()*L)
        y.append(random.random()*L)
        z.append(random.random()*L)
        vx.append(0)
        vy.append(0)
        vz.append(0)
    dt=tf/N
    t=0
    while t<tf:
        #listes temporaires qui figent les coordonnées de tous les points pour calculer les mvts à partir du même instant
        Lx=x
        Ly=y
        Lz=z
        for i in range(n):
            xold,yold,zold,vxold,vyold,vzold=x[i],y[i],z[i],vx[i],vy[i],vz[i]
            Ax=Ay=Az=0
            for k in range(n):
                d=sqrt((Lx[k]-Lx[i])**2+(Ly[k]-Ly[i])**2+(Lz[k]-Lz[i])**2)
                if d>0:
                    Ax+=G*m*(Lx[k]-Lx[i])/d**3
                    Ay+=G*m*(Ly[k]-Ly[i])/d**3
                    Az+=G*m*(Lz[k]-Lz[i])/d**3
            #Mouvement de l'objet i
            vxnew=Ax*dt+vxold
            vynew=Ay*dt+vyold
            vznew=Az*dt+vzold
            vx[i]=vxnew
            vy[i]=vynew
            vz[i]=vznew
        
            xnew=vxold*dt+xold
            ynew=vyold*dt+yold
            znew=vzold*dt+zold
        
            x[i]=xnew
            y[i]=ynew
            z[i]=znew
                
            #changement de variable
            xold,yold,zold,vxold,vyold,vzold=xnew,ynew,znew,vxnew,vynew,vznew
        t+=dt    
        
    #affichage
    for i in range(n):
        xi=x[i]
        yi=y[i]
        zi=z[i]
        ax.scatter(xi, yi, zi, c='r')
        
nbody(30,1,100,10000,60)

[Chanur]

le temps de calcul s'envole très vite.

Quelques pistes si tu cherches à accélérer les choses :
- utiliser un langage compilé (par ex. C++) au lieu d'un langage interprété
- optimiser l'algorithme : par exemple, tu calcules deux fois la distance entre deux corps. c'est sans doute ce qui prends le plus de temps à cause de sqrt()
- tu utilises l'opérateur puissance (**), mais x² = x*x qui est sûrement plus rapide à calculer que x**2
- de même, tu calcules d**3, mais en fait tu commence par calculer d², et par en prendre la racine carrée pour avoir d, alors que d³ = d*d² : il suffirait d'avoir mémorisé d²
- je ne vois pas l'intérêt, à la fin de la boucle de la ligne "changement de variables", puisque tu va les réinitialiser au début de la boucle suivante.
- je ne vois pas l'intérêt des variables vxold, vyold et vzold : il suffirait d'écrire vx[i] = vx[i] + Ax*dt

[Azghar]

Bonjour Chanur et merci pour ton aide. Je suis encore novice en programmation , le passage en C se fera plus tard quand j'aurais plus de temps cad après mes concours. Pour le calcul des distances (le seul truc dont j'avais conscience que je pouvais l'optimiser) je vais regarder.
Je n'étais pas au courant pour le x**2 et le x*x, bon à savoir.
Pour la distance au cube je corrige ça en effet ça économise des calculs.
Pour les deux derniers points ça vient de mes anciens algorithmes un peu laborieux, d'ailleurs je comptais faire de même pour x[i] car je nai plus besoin forcément de garder la trajectoire en mémoire avec l'animation.
Sinon j'avais d'autres idées. Vu que j'ai deux cartes graphiques(nvidia) puissantes un processeur également très bon. J'avais pensé faire fonctionner le tout au lieu de faire fonctionner qu'un seul coeur de mon CPU(cas actuel). J'ai vu qu'on pouvait faire marcher cuda avec python.

Sinon j'ai vu l'algorithme de Barnes-Hut qui semble rependu et permet de passer a une complexité en (nlog(n)).
Qu'en penses-tu?

[Chanur]

Sinon j'avais d'autres idées. Vu que j'ai deux cartes graphiques(nvidia) puissantes un processeur également très bon. J'avais pensé faire fonctionner le tout au lieu de faire fonctionner qu'un seul coeur de mon CPU(cas actuel). J'ai vu qu'on pouvait faire marcher cuda avec python.

Je n'y connais pas grand'chose, mais je sais que les cartes graphiques sont très douées pour faire du calcul vectoriel en 3 dimensions (en parallélisant le calcul des différentes coordonnées). En allant plus loin, on peut peut-être (sûrement, même) paralléliser le calcul des différents corps mais ça doit être plus dur ...

Sinon j'ai vu l'algorithme de Barnes-Hut qui semble rependu et permet de passer a une complexité en (nlog(n)).
Qu'en penses-tu?

Je ne connaissais pas. Ça a l'air intéressant, mais comme il faut calculer l'arbre, le bénéfice ne se manifestera qu'à partir d'un certain nombre de corps. Combien ? A part en essayant, je ne vois pas trop de moyen de le savoir.

Une chose que j'oubliais : tes corps peuvent s'approcher l'un de l'autre autant qu'il veulent, alors qu'en fait il y a une distance minimale en dessous de laquelle l'intervalle dt devient trop grand (et, par exemple, des corps se font éjecter de façon absolument pas réaliste).
Tu peux contrôler si c'est le cas en calculant l'énergie totale (somme des énergies cinétiques + somme des énergies potentielles des corps 2 à 2) qui doit rester constante : si elle varie c'est que ton algo est pris en défaut. Ce n'est pas difficile à calculer, puisque tu connais les masses, les vitesses et les distances.
On peut pallier à ça :
- en diminuant dt quand ça se produit. (mais c'est plus dur d'avoir une animation régulière)
- en faisant fusionner deux corps s'ils s'approchent trop l'un de l'autre : somme des masses, barycentre des positions et, pour la vitesse, conservation de la quantité de mouvement. (mais c'est moyennement réaliste : par exemple quand la Terre a percuté un astre de la masse de Mars, il n'en est pas résulté un astre unique, mais deux : la Terre et la Lune, et ça, pour le coup, c'est un peu plus dur à simuler )

[Azghar]

Haha tu lis dans mes pensées. Mon objectif à terme est justement de prendre en compte les collisions pour arriver à simuler la formation d'une planète à partir d'une ceinture d’astéroïde. Mais on en est encore loin.
Effectivement je voulais mettre un compteur d'énergie du système dans sa globalité(qui ne sera pas constant vu la méthode de résolution de mes eq diff). Simplement je ne savais pas vraiment comment faire; donc l'e totale si je reprend ce que tu dis, c'est la somme de l'énergie cinétique de chaque corps(ça je veux bien le comprendre) + la somme des ep des corps 2à2 (ça c'est beaucoup moins évident pour moi). A moins que tu parlais seulement du cas n=2(là ok) Comment le justifier ?

Pour revenir a des choses plus réalistes, dans mon programme j'ai effectivement toujours un corps qui se promène super loin du groupe, qui a été éjecté, tu as vu juste. Je vais intégrer un système de dt adaptatif car la deuxième solution parait plus compliqué effectivement mais aussi plus fun, donc je l'essayerais quand même. Mais parcontre est-ce judicieux d'être a quantité de mvt constante avant après choc? Car dans ce genre d'évènement, une grande partie de l'énergie doit être dissipée sous forme d'énergie thermique, mais pour savoir combien... ? ça doit dépendre de pas mal de choses

[Chanur]

donc l'e totale si je reprend ce que tu dis, c'est la somme de l'énergie cinétique de chaque corps(ça je veux bien le comprendre) + la somme des ep des corps 2à2 (ça c'est beaucoup moins évident pour moi). A moins que tu parlais seulement du cas n=2(là ok) Comment le justifier ?

Non, je parlais des n corps.
Ça se justifie par le principe de superposition : le champ en un point est égal à la somme des champs dus aux différents corps, et le champ est le gradient du potentiel. Donc le potentiel en un point est égal à la somme des potentiels dus aux différents corps, sous réserve d'avoir une référence commune (potentiel nul), ce qui se fait facilement en prenant comme référence le potentiel à l'infini. Dans ce cas, le potentiel d'une masse m dans le champ d'une masse M est égal à l'intégrale de l'infini à d de MmG/r²dr, soit -MmG/r (ne pas oublier le signe moins : par rapport à l'infini, toutes les énergies potentielles sont négatives).
Note que l'énergie potentielle de m par rapport à M est la même que l'énergie de M par rapport à m : inutile de la calculer 2 fois

Mais parcontre est-ce judicieux d'être a quantité de mvt constante avant après choc? Car dans ce genre d'évènement, une grande partie de l'énergie doit être dissipée sous forme d'énergie thermique, mais pour savoir combien... ? ça doit dépendre de pas mal de choses

C'est justement pour ça qu'on ne peut pas se baser sur l'énergie. Mais la quantité de mouvement, elle, reste constante (en pratique, il faudrait évidemment tenir compte de la quantité de mouvement des gaz et poussières éjectés). Du coup, bien sûr, dans le cas d'une collision l'énergie totale diminue.
Ça revient à supposer un choc parfaitement mou : l'intégralité de l'énergie cinétique relative entre les deux corps est transformée en chaleur.

[Azghar]

rien qu'en ajustant les deux premiers points que tu m'as conseillé je passe de 17,17sec d’exécution à 9.8.

Par contre j'essaye de régler le problème des distances (jusqu'ici calculées deux fois) mais là j'ai quelques soucis. Ma nouvelle méthode de calcul fait passer le temps à 28sec d'execution et en plus donne un résultat aberrant.
J'ai raisonné sur la symétrie de la distance donc pour tout corps i,k appartenant à la boite , d(i,k)=d(k,i).
Donc j'ai pensé que ça revient a remplir la moitié d'une matrice, et remplir l'autre moitié par symétrie. Donc en fait je rempli un triangle, donc le k commence à i à chaque fois
mais c'est manifestement faux

ça donne ça :
[...]

 while t<tf:
        #listes temporaires qui figent les coordonnées de tous les points pour calculer les mvts à partir du même instant
        Lx=x
        Ly=y
        Lz=z
        #Calcul des distances
        for i in range(n):
            k=i
            while k<=n-1:
                D[i][k]=((Lx[k]-Lx[i])*(Lx[k]-Lx[i])+(Ly[k]-Ly[i])*(Ly[k]-Ly[i])+(Lz[k]-Lz[i])*(Lz[k]-Lz[i]))**1.5
                if k!=i:
                    D[k][i]=D[i][k]
                k+=1

            #D est une matrice symétrique.
            #en effet par symétrie d(k,i)==d(i,k)

[...]

[Chanur]

Je ne vois pas ce qui cloche, mais je ne connais rien à python ...
Je me suis contenté de lire attentivement ton code, ce qui est notoirement insuffisant pour trouver un bug.
Comment est déclarée la matrice D ?

Sinon, pour les performances, je crains que le **1.5 ne soit particulièrement gourmand,
tu fais deux fois chaque soustraction,
il est inutile de calculer D[i][i] : tu ne t'en serviras pas (et de toutes façon ça vaut 0),
et k<=n-1 fait une soustraction inutile à chaque boucle.

J'aurais écrit :

for i in range(n-1):
    k=i+1
    while k<n:
        dx = Lx[k]-Lx[i];
        dy = Ly[k]-Ly[i];
        dz = Lz[k]-Lz[i];
        d2 = dx*dx + dy*dy + dz*dz
        D[k][i]=D[i][k]=d2*sqrt(d2)
        k+=1

(je ne suis pas sûr de la syntaxe D[k][i]=D[i][k]=...)

[Azghar]

Bonjour,
Voila tout fonctionne bien: n corps en 3d https://youtu.be/sQsFcBHf57I
Prochaine étape: prendre en compte les collisions.

[Azghar]

Un petit test avec une simulation de saturne, j'ai mis 10 000 corps.
Pour celle ci j'ai supprimé l'interaction entre chaque corps, cela aurait été trop long sinon. (et c'est largement légitime de dire que les corps constituants les anneaux interagissent peu entre eux comparé à l'interaction avec Saturne)

https://youtu.be/Xlcf5XY5aZo

[Arzhur]

En tout cas, ça a de la gueule

[Azghar]

Oui en effet c'est pas mal
Néanmoins j'ai quelques difficultés à intégrer un bon système de collision au programme.

En effet, vu qu'il y pleins de forces appliquées aux corps, je ne peux pas utiliser la conservation de la quantité de mouvement avant et après le choc.
Après tout je ne sais pas si c'est une bonne chose les collisions car ça va tout ralentir encore plus. Et j'ai trouvé une alternative sympa qui évite que les corps se catapultent lorsqu'ils sont trop proches :
La KS regularization : http://kiaa.pku.edu.cn/~kouwenhoven/...seth_talk2.pdf

Mon but à terme est de former des planétésimaux à partir d'un nuage de poussière (il faut que j'utilise Barnes Hut, impossible en brute force)
Donc c'est pas pour tout de suite...

[Azghar]

Voici une des rares configurations stables en 3 corps.
https://youtu.be/mJkD1W0pbmY

[cxFreeze]

Du très bon, on voit très bien la montée en puissance depuis la première itération du programme...

De plus, l'interface est particulièrement esthétique et bien pensée, un grand bravo à vous.

[Azghar]

merci ^^ et c'est pas fini

[Azghar]

Un compteur d'énergie mécanique et un compteur de temps a été ajouté. L'interface est faite par un pote
voici la simulation du célèbre problème pythagoricien avec la méthode d'euler en prenant dt=0.2*10**-4 s
Je suis entrain d'implémenter une autre méthode d'intégration : Verlet
Affaire à suivre

[Azghar]

Simulation du système solaire, dommage le compteur de temps ne s'affiche pas quand je film.
La simulation s'étend sur 30 ans :