Calcul de coordonnées force gravitationelle

[Azghar]

Je ne comprend pas pourquoi ça ne marche pas bien :
#programme de calcul de coordonnees du point C (point B fixe sur l'origine) Point C demarre sur axe des abscisses

from math import *
G=6.6738*10**(-11)
m=10000000
          #Masse de l'objet B               
Cxt=0            #Coordonnee x initiale du point C
Cyt=50              #Coordonnee y initiale du point C
d=50 
V0=1  
cosinusbeta=1                                      #beta : orientation du vecteur vitesse à t=0   (horizontal)
sinusbeta=0                                       #alpha : orientation du vecteur acceleration a t=0  (ici vertical) 
Vx=V0
Vy=0                           
for t in range (1,2000):
    cosinusalpha=-Cxt/d
    sinusalpha=-Cyt/d
    V0=sqrt((cosinusalpha*G*t*m/(d**2)+cosinusbeta*V0)**2+(sinusalpha*G*m*t/d**2+sinusbeta*V0)**2)
    Vx=cosinusalpha*G*t*m/(d**2)+cosinusbeta*V0
    Vy=sinusalpha*G*m*t/d**2+sinusbeta*V0
    sinusbeta=Vy/V0
    cosinusbeta=Vx/V0
    d=sqrt((Cxt)**2+(Cyt)**2)
    Cxt=cosinusalpha*G*m*(t**2)/(2*d**2)+cosinusbeta*V0*t+0     #Les termes additionnes a la fin viennent des coordonnes initiales.
    Cyt=sinusalpha*G*m*(t**2)/(2*d**2)+sinusbeta*V0*t+50
    print(Cxt), print(Cyt)
n=input('pressez entrer pour fermer')

Une idée ?
-----

[invite9dc7b526]

bonjour,

je ne connais pas ce langage mais je vois une instruction "for t in ..." et manifestement il y a plusieurs choses à faire dans cette boucle et je ne vois pas de marque de début et fin de boucle (l'équivalent des {} de C ou R)

[invitef35ebd48]

ca ressemble a du python sans indentation....bref c'est illisible

[invite1c6b0acc]

Bonjour, Azghar

C'est difficile de vous aider : vous ne comprenez pas ce que fait votre programme, mais au moins vous savez en quel langage il est écrit, ce que vous avez voulu faire et ce qui se passe quand vous le lancez.
Nous on ne sait rien de tout ça ...
Un problème informatique n'est pas forcément facile à résoudre, alors rajouter des énigmes en plus, ça ne va pas aider ...

[A voir en vidéo sur Futura]

[invite43901482]

Équation horaire du mouvement plan ?

[invited9b9018b]

Bonsoir,

1) langage non précisé
2) code dégueu balancé comme ça sans balise [code]
3) un simple "ça marche pas bien" lâché comme seule indication du problème ce qui n'aide pas

ça ne donne pas envie d'aider c'est sur.

A+

[Azghar]

Bonsoir,
dsl pour la lisibilité, c'était pas voulu mes espaces ne sont pas passés.
En fait c'est un programme qui calcul la position au cours du temps de deux corps, un fixe à l'origine de masse m et l'autre en mouvement.
à terme je veux simuler des système à n corps mais je commence par ça.
j'ai des difficultés à exprimer l'angle du vecteur vitesse au cours du temps
C'est en python.
quand je fais mes tests l'objet en mouvement ne semble pas vraiment être influencé par la masse, pourtant grande. (j'ai fais plusieurs test)
je suppose que l'erreur vient de l'expression de ma vitesse. Pouvez vous m'aider ?

[JPL]

Utilise la balise Code : regarde le résultat dans ton premier message.

[invite1c6b0acc]

Il me semble que tu devrais revoir l'analyse de ton problème : je ne vois vraiment pas ce qui motive l'utilisation de sinus et de cosinus ...
Est-ce que tu as fait un dessin ?
Est-ce que tu as écrit clairement les formules que tu veux implémenter ?

Question subsidiaire : on pourrait voir l'un et les autres ?

Parce que soit tu t'es trompé dans ton analyse, soit tu t'es trompé dans l'implémentation, soit tu te trompes sur le résultat que ton programme devrait donner. Et on peut difficilement le savoir avec uniquement le programme.

[Azghar]

OK merci JPL. Je ne connaissais pas cette balise.


Je donne toutes les infos ce soir alors, les dessins et les formules qui m'ont fait aboutir a ça.
Peut on intégrer des photos sur le forum? Car je ne sais pas encore écrire en latex.
Bonne journée

[invite1c6b0acc]

Oui, on peut intégrer des photos. Il faut les mettre en pièce jointe. Attention, il est interdit d'utiliser un serveur externe (sinon les images risquent de disparaître au bout d'un certain temps et rendre le message incompréhensible).
Pour Latex, tu peux trouver un message dont les formules ont vaguement la même forme que les tiennes, le copier-coller (tu cliques d'abord sur "citer" comme si tu voulais répondre, et tu annules ensuite), et bricoler pour y mettre tes propres formules. C'est un peu lourdingue comme méthode, mais ça marche.

[invite4bf147f6]

Essai avec les deux print sur 2 lignes et enleve la virgule

[Azghar]

Un début (j'ai pas tout mis)

[invite1c6b0acc]

Tu te trompes dans tes formules :

Tu écris :
et tu en déduis :

Mais c'est faux. On n'a pas , sauf si est constante, ce qui n'est clairement pas le cas.

[Azghar]

Miiince exact, j'ai confondu avec mes anciens exos de terminal ou g était constant
Bon je repart de 0

[Azghar]

Je reviens quelques mois plus tard avec une proposition.
J'ai appris un peu de python entre temps. Pour expliquer rapidement ce que fais le programme, il calcul à chaque instant t le vecteur accélération, puis il en déduit le vecteur vitesse, puis le vecteur position et renvoit une liste des coordonnés du point m
C'est bien sûr une approximation, car l'objet en mouvement est "laissé libre" entre chaque intervalle de temps t. Le programme "recadre " et recalcul les variables toutes les secondes dans mon exemple.
J'ai regardé sur géogebra, les coordonnés renvoyées me semblent plus ou moins correctes, sauf que le mouvement de m est "décentré" par rapport à M, ce que je ne comprend pas.
J'aimerai lié ce programme (je sais pas comment on dit, ni comment on fait) à une interface de type géogebra, qui me placerait les coordonnées renvoyées sur un graphique. Savez vous comment je peux faire ?
Merci d'avance

def simulation(mx,my,vxm,vym,M):     #mx , my coordonnees initiales du point m, vx,vy vitesse initiale du point m,  M = masse de l astre attracteur.
    h=1                              #precision (intervalle de temps entre deux calculs) 
    C=[]                             #liste qui donne les coordonnees de M a chaque instant t
    t=0  
    G=6.67384*10**(-11)
    while t<100000:
        d=sqrt(mx**2+my**2)
        if d>50:                             #distance minimale entre m et M (collision)
            alpha=asin(mx/d) 
            axm=sin(alpha)*G*M/(d**2)        #calcul du vecteur acceleration a l instant t
            aym=cos(alpha)*G*M/(d**2)
            
            vxm+=-cos(alpha)*G*t*M/(d**2)    #calcul du vecteur vitesse  a l instant t
            vym+=sin(alpha)*G*t*M/(d**2)
                
            mx+=-sin(alpha)*G*t**2*M/(2*(d**2))+vxm*t         #calcul du vecteur position a  l instant t
            my+=-cos(alpha)*G*t**2*M/(2*(d**2))+vym*t
            C.append([mx,my,t])
            t+=h
        else :
            break
    return C
    
    
    
print(simulation(-2000,1000,0.1,0,1000000))

[invite43901482]

Tu ne peux pas, dans ton cas, regarde du côté de matplotlib...

[Azghar]

voilà un essai avec matplotlib, j'ai mis des données qui sont proches du système terre lune, donc il y'a surement une erreur dans mon code :
(durée = 24h)

[Azghar]

Je reviens vers vous, avec quelques connaissances en plus ;
J'aimerai appliquer la méthode d'Euler pour une approximation de mes équa diffs, mais je ne sais pas comment faire. (2 ED à coefficients non constants et dont l'une dépend de l'autre)
Je vous donne l'état actuel des calculs :

[Azghar]

Si j'ai bien compris, je dois séparer chacune des ED d'ordre 2 en 2 ED d'ordre 1.
C'est bien ça ?

[invitef35ebd48]

Bonjour,


Oui c'est l'idée et derrière tu résous chaque équations avec une méthode d'Euler ( ou éventuellement RK4)

[invite1cb4eccb]

Bonjour
Euler me semble une bonne approche car peu compliqué a coder.
après si tu veux une solution plus précise, utilise Euler implicite ou rk2/rk4

[Azghar]

ça y est ! j'y suis arrivé pour un système à n=2 corps ! (pour le moment en 2D, la 3D arrive bientôt)
J'ai toujours un problème : quand je prends des valeurs initiales trop grandes (système terre lune par exemple) les calculs partent dans les choux.
Avec des valeurs plus raisonnables on observe des courbes convaincantes, dites ce que vous en pensez :

from math import*
import matplotlib.pyplot as plt

def orbite(M0,M1,x0,y0,vx0,vy0,x1,y1,vx1,vy1,ti,tf,N):
    Lx0=[x0]
    Ly0=[y0]
    Lx1=[x1]
    Ly1=[y1]
    t=ti
    G=6.67384*10**-11
    x0old,y0old,vx0old,vy0old,x1old,y1old,vx1old,vy1old=x0,y0,vx0,vy0,x1,y1,vx1,vy1
    h=(tf-ti)/N
    while t<tf:
        #Mouvement de l'objet 0:
        vx0new=G*M1*(x1old-x0old)/((x1old-x0old)**2+(y1old-y0old)**2)*h+vx0old
        vy0new=G*M1*(y1old-y0old)/((x1old-x0old)**2+(y1old-y0old)**2)*h+vy0old
        x0new=vx0old*h+x0old
        y0new=vy0old*h+y0old
        Lx0.append(x0new)
        Ly0.append(y0new)
        
        #Mouvement de l'objet 1:
        vx1new=-G*M0*(x1old-x0old)/((x1old-x0old)**2+(y1old-y0old)**2)*h+vx1old
        vy1new=-G*M0*(y1old-y0old)/((x1old-x0old)**2+(y1old-y0old)**2)*h+vy1old
        x1new=vx1old*h+x1old
        y1new=vy1old*h+y1old
        Lx1.append(x1new)
        Ly1.append(y1new)
        
        #changement des variables:
        x0old,y0old,vx0old,vy0old,x1old,y1old,vx1old,vy1old=x0new,y0new,vx0new,vy0new,x1new,y1new,vx1new,vy1new
        t+=h
        
    plt.plot(Lx0,Ly0)
    plt.plot(Lx1,Ly1)
    plt.show
    
orbite(6*10**10,10**10,0,0,0,1,1000,0,0,3,0,60*60,1000000)

[Azghar]

Il y avait une petite erreur dans le code, voilà la nouvelle version qui marche parfaitement avec les paramètres du système terre/lune.
Maintenant je passe en 3D et à 3 corps.

from math import*
import matplotlib.pyplot as plt

def orbite(M0,M1,x0,y0,vx0,vy0,x1,y1,vx1,vy1,ti,tf,N):
    Lx0=[x0]
    Ly0=[y0]
    Lx1=[x1]
    Ly1=[y1]
    t=ti
    G=6.67384*10**-11
    x0old,y0old,vx0old,vy0old,x1old,y1old,vx1old,vy1old=x0,y0,vx0,vy0,x1,y1,vx1,vy1
    h=(tf-ti)/N
    while t<tf:
        #Mouvement de l'objet 0:
        vx0new=G*M1*(x1old-x0old)/(((x1old-x0old)**2+(y1old-y0old)**2)**1.5)*h+vx0old
        vy0new=G*M1*(y1old-y0old)/(((x1old-x0old)**2+(y1old-y0old)**2)**1.5)*h+vy0old
        x0new=vx0old*h+x0old
        y0new=vy0old*h+y0old
        Lx0.append(x0new)
        Ly0.append(y0new)
        
        #Mouvement de l'objet 1:
        vx1new=-G*M0*(x1old-x0old)/(((x1old-x0old)**2+(y1old-y0old)**2)**1.5)*h+vx1old
        vy1new=-G*M0*(y1old-y0old)/(((x1old-x0old)**2+(y1old-y0old)**2)**1.5)*h+vy1old
        x1new=vx1old*h+x1old
        y1new=vy1old*h+y1old
        Lx1.append(x1new)
        Ly1.append(y1new)
        
        #changement des variables:
        x0old,y0old,vx0old,vy0old,x1old,y1old,vx1old,vy1old=x0new,y0new,vx0new,vy0new,x1new,y1new,vx1new,vy1new
        t+=h
    plt.plot(Lx0,Ly0)
    plt.plot(Lx1,Ly1)
    plt.show
    
orbite(6*10**24,7*10**22,0,0,0,0,385000000,0,0,1023,0,60*60*24*27,10000000)

[Azghar]

J'aimerais faire une animation du mouvement, mais j'ai pas bien compris comment on peut faire.
Si quelqu'un peut m'aider, je lui en serais très reconnaissant

[invite1cb4eccb]

Pour commencer choisi ta méthode de rendu
Alors la, tu as deux solutions : soit tu fais un truc par toi même (renseigne toi sur les tk et canvas du module tkinter (de base dans python)) soit tu utilises un module externe (comme pygame par exemple (je connais pas très bien les modules de python))
Et après tu va vouloir mettre ton système dans chaque position une par une.

[Azghar]

Voilà ça marche en 3D
Prochaine étape, animation avec Tkinter (comme laentheorist me l'a proposé), et mettre tout ça en forme avec des fenêtres.
exemple avec le système Jupiter-Soleil (données approximatives)
code :

###Programme de simulation d'attraction gravitationelle entre deux corps
#Les masses sont considérées ponctuelles, 

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,x0,y0,z0,vx0,vy0,vz0,x1,y1,z1,vx1,vy1,vz1,ti,tf,N):
    Lx0=[x0]
    Ly0=[y0]
    Lz0=[z0]
    Lx1=[x1]
    Ly1=[y1]
    Lz1=[z1]
    t=ti
    G=6.67384*10**-11
    x0old,y0old,z0old,vx0old,vy0old,vz0old,x1old,y1old,z1old,vx1old,vy1old,vz1old=x0,y0,z0,vx0,vy0,vz0,x1,y1,z1,vx1,vy1,vz1
    h=(tf-ti)/N
    
    while t<tf:
        #Mouvement de l'objet 0:
        vx0new=G*M1*(x1old-x0old)/(((x1old-x0old)**2+(y1old-y0old)**2+(z1old-z0old)**2)**1.5)*h+vx0old
        vy0new=G*M1*(y1old-y0old)/(((x1old-x0old)**2+(y1old-y0old)**2+(z1old-z0old)**2)**1.5)*h+vy0old
        vz0new=G*M1*(z1old-z0old)/(((x1old-x0old)**2+(y1old-y0old)**2+(z1old-z0old)**2)**1.5)*h+vz0old
        
        x0new=vx0old*h+x0old
        y0new=vy0old*h+y0old
        z0new=vz0old*h+z0old
        
        Lx0.append(x0new)
        Ly0.append(y0new)
        Lz0.append(z0new)
        
        
        
        #Mouvement de l'objet 1:
        vx1new=-G*M0*(x1old-x0old)/(((x1old-x0old)**2+(y1old-y0old)**2+(z1old-z0old)**2)**1.5)*h+vx1old
        vy1new=-G*M0*(y1old-y0old)/(((x1old-x0old)**2+(y1old-y0old)**2+(z1old-z0old)**2)**1.5)*h+vy1old
        vz1new=-G*M0*(z1old-z0old)/(((x1old-x0old)**2+(y1old-y0old)**2+(z1old-z0old)**2)**1.5)*h+vz1old
        
        x1new=vx1old*h+x1old
        y1new=vy1old*h+y1old
        z1new=vz1old*h+z1old
        
        Lx1.append(x1new)
        Ly1.append(y1new)
        Lz1.append(z1new)
        
        #changement des variables:
        x0old,y0old,z0old,vx0old,vy0old,vz0old,x1old,y1old,z1old,vx1old,vy1old,vz1old=x0new,y0new,z0new,vx0new,vy0new,vz0new,x1new,y1new,z1new,vx1new,vy1new,vz1new
        t+=h
        
###Mesures (ne fontionne qu'en 2D pour le moment)
    '''#Mesures
    #periode:
    for i in range(len(Ly1)):
        if Ly1[i]-Ly0[i]<0:
            T=2*i*h
            if T<=60:
                print('T=',T,'s')
            elif 60<T<=3600:
                print('T=',T/60,'min')
            elif 3600<T<=86400:
                print('T=',T/(60**2),'heures')
            elif 86400<T<=31536000:
                print('T=',T/(60**2*24),'jours')
            else :
                print('T=',T/(60**2*24*365),'annees')
                
            break
            
    #distance minimale:
    m=2*sqrt((Lx1[0]-Lx0[0])**2+(Ly1[0]-Ly0[0])**2)
    for i in range(len(Lx0)):
        u=sqrt((Lx1[i]-Lx0[i])**2+(Ly1[i]-Ly0[i])**2)
        if u<m:
            k=i
            m=u
    print('rmin=',m/1000,'km')

    #distance maximale:
    M=0
    for i in range(len(Lx0)):
        u=sqrt((Lx1[i]-Lx0[i])**2+(Ly1[i]-Ly0[i])**2)
        if u>M:
            j=i
            M=u
    print('rmax=',M/1000,'km')        
    
    #vitesse maximale:
    if k==0:
        vmax=sqrt(((Lx1[k+1]-Lx1[k])/h)**2+((Ly1[k+1]-Ly1[k])/h)**2)
    elif k==len(Lx0):
        vmax=sqrt(((Lx1[k]-Lx1[k-1])/h)**2+((Ly1[k]-Ly1[k-1])/h)**2)
    else :
        vmax=sqrt(((Lx1[k+1]-Lx1[k-1])/(2*h))**2+((Ly1[k+1]-Ly1[k-1])/(2*h))**2)
    print('vmax=',vmax,'m/s')
    #vitesse minimale:
    if j==0:
        vmin=sqrt(((Lx1[j+1]-Lx1[j])/h)**2+((Ly1[j+1]-Ly1[j])/h)**2)
    elif j==len(Lx0):
        vmin=sqrt(((Lx1[j]-Lx1[j-1])/h)**2+((Ly1[j]-Ly1[j-1])/h)**2)
    else :
        vmin=sqrt(((Lx1[j+1]-Lx1[j-1])/(2*h))**2+((Ly1[j+1]-Ly1[j-1])/(2*h))**2)
    print('vmin=',vmin,'m/s')
    #demi grand axe:
    a=(m+M)/2
    print('a=',a/1000,'km')
    
    #excentricite:
    e=((M-a)/a+(a-m)/a)/2
    print('e=',e)'''
    
    
    #affichage
    ax.plot(Lx0,Ly0,Lz0)
    ax.plot(Lx1,Ly1,Lz1)
    plt.show() 
    
orbite(1.989*10**30,1.898*10**27,0,0,0,0,0,0,800000000000,0,0,0,12000,0,0,3600*24*365*12,1000000)
'''orbite(M0,M1,x0,y0,z0,vx0,vy0,vz0,x1,y1,z1,vx1,vy1,vz1,ti,tf,N)'''

###Exemples types :
'''Données système :
                    Terre-Lune :
                    Soleil-Jupiter:

[invitef35ebd48]

Bon amusement !

Tu pourras rajouter un 3ème corps pour plus de fun.

[JPL]

Là ça devient sadique !

[Azghar]

Oui c'est ce que je voulais faire. Mais cela ne fait que reporter le problème car ensuite je voudrais faire à 4 corps, puis à 5...
J'aimerai faire un algorithme qui fait du N corps (objectif initial de mon programme à terme)

C'est pour ça que pour le moment je m’entraîne avec des solutions simples pour voir comment ça fonctionne et j'apprend pleins de trucs.