résolution d'équation différentielle avec python

résolution d'équation différentielle avec python - Python - Programmation

Marsh Posté le 04-06-2017 à 19:44:08    

Bonjour, j'essaye de coder un algorithme pour résoudre les équations du problème à N corps mais quand je trace la trajectoire avec la méthode d'Euler j'obtiens des droites ce qui est manifestement faux. Si vous pouviez trouver l'erreur ça serait très sympathique car je la détecte pas.
 

Code :
  1. def mullist(coeff,L):
  2.     for k in range(0,len(L)):
  3.         L[k]=coeff*L[k]
  4.     return(L)
  5.  
  6. def somlist(L,M):
  7.     N=[]
  8.     for k in range(0,len(L)):
  9.         N.append(L[k]+M[k])
  10.     return(N)
  11.  
  12. def diflist(L,M):
  13.     N=[]
  14.     for k in range(0,len(L)):
  15.         N.append(L[k]-M[k])
  16.     return(N)
  17.  
  18. def normevecteur(V):
  19.     normecarre=0
  20.     for k in range(0,len(V)):
  21.         normecarre+=V[k]**2
  22.     return(sqrt(normecarre))
  23.  
  24.  
  25.  
  26. def Force_entre_deux_corps(m1,p1,m2,p2):#force de 2 sur 1 à un instant fixé
  27.     G=39.4
  28.     Force=[0,0,0]
  29.     for k in range(0,3):
  30.         Force[k]=(-G*m1*m2*(p2[k]-p1[k]))/(normevecteur(diflist(p2,p1))**3)
  31.     return(Force)
  32.  
  33. def Force_entre_N_corps(j,m,position): #Résultante des forces exercées sur j à un instant fixé
  34.     Resultante=[0,0,0]
  35.     for k in range(0,len(position)):
  36.         if k!=j:
  37.             Force=Force_entre_deux_corps(m[j],position[j],m[k],position[k])
  38.             Resultante[0]=Resultante[0]+Force[0]
  39.             Resultante[1]=Resultante[1]+Force[1]
  40.             Resultante[2]=Resultante[2]+Force[2]
  41.     return((Resultante))
  42.  
  43. def Acceleration(m,position):
  44.     L=[]
  45.     for k in range(0,len(position)):
  46.         L.append(mullist(1/m[k],Force_entre_N_corps(k,m,position))) #Calcul du vecteur acceleration pour chaque corps"
  47.     return(L)
  48.  
  49. def fonction_de_passage(m,Y):
  50.     f=[Y[1],Acceleration(m,Y[0])]
  51.     return(np.array(f))
  52.  
  53. def suite_Euler_pas_fixe(m,t_min,t_max,n,position_initiale,vitesse_initiale):
  54.     t=t_min
  55.     temps=[t_min]
  56.     h=(t_max-t_min)/n  #pas de discrétisation
  57.     Y_corps=np.array([position_initiale,vitesse_initiale])
  58.     Y_temps=[list(Y_corps)]
  59.     while t<=t_max:
  60.         Y_corps=Y_corps+h*fonction_de_passage(m,Y_corps)
  61.         Y_temps.append(list(Y_corps))
  62.         t+=h
  63.         temps.append(t)
  64.     return(temps,np.array(Y_temps))


 
Merci d'avance(j'ai même codé les méthodes de Runge-Kutta et Verlet mais ça ne fonctionne pas)
 
Pour les unités j'ai tous mis en unité astronomique.

Reply

Marsh Posté le 04-06-2017 à 19:44:08   

Reply

Marsh Posté le 05-06-2017 à 10:46:03    

Personne ne voit l'erreur ?

Reply

Marsh Posté le 05-06-2017 à 12:08:35    

Bonjour,
 
Cela aiderait déjà de savoir si le problème est lié à l'algo (résultat non attendu) ou au langage (erreur Python).


---------------
C'est en écrivant n'importe quoi qu'on devient n'importe qui.
Reply

Marsh Posté le 05-06-2017 à 16:28:59    

L'algo ne donne pas le résultat attendu(je l'ai dit dans le premier message puisque ça donne des droites comme trajectoires, les coordonnées sont toujours croissantes en valeur absolue)

Reply

Marsh Posté le 05-06-2017 à 16:36:21    

Ca pourrait être utile aussi de donner un lien (Wikipédia?) vers cet algo. "méthode d'Euler" ça me dit quelque chose mais j'ai pas ça en tête. :o

Reply

Marsh Posté le 05-06-2017 à 17:01:11    

http://desaintar.free.fr/python/tp/tp_euler.pdf
 
Dans mon cas c'est le cas d'une EDL d'ordre 2 vectorielle.


Message édité par maths1 le 05-06-2017 à 17:01:42
Reply

Marsh Posté le 05-06-2017 à 17:37:22    

Oui désolé j'ai oublié la première ligne de ton message pendant que je commençais à regarder.
 
Je vais peut-être dire de la merde mais pourquoi c'est au cube dans le calcul de la force de gravité entre deux corps :
 

Code :
  1. Force[k]=(-G*m1*m2*(p2[k]-p1[k]))/(normevecteur(diflist(p2,p1))**3)


---------------
C'est en écrivant n'importe quoi qu'on devient n'importe qui.
Reply

Marsh Posté le 05-06-2017 à 18:22:16    

Normalement c'est au carré et on multiplie par le vecteur directeur de ri-rj.
Or le vecteur directeur est égale à (ri-rj)/(norme(ri-rj))
Ainsi la norme est au cube.

Reply

Marsh Posté le 05-06-2017 à 18:36:06    

EDIT : ok je viens de comprendre les unités utilisées pour G, je retrouve pareil, une erreur sur ce dernier ou dans la conversion des unités peut conduire à des erreurs de trajectoires donc vérifie bien ça. Je ne vois pas plus d'erreur sinon.


Message édité par MaybeEijOrNot le 05-06-2017 à 19:13:35

---------------
C'est en écrivant n'importe quoi qu'on devient n'importe qui.
Reply

Marsh Posté le 05-06-2017 à 19:57:30    

m=[1.0,2.988621344*10**(-6)]
vitesse_initiale=[[0,0,0],[0,6.27,0]]
position_initiale=[[0,0,0],[0.997332336,0,0]]
G=39.34
 
Il y avait quelques légères erreurs de conversion mais même en les corrigeant j'obtiens encore des droites.


Message édité par maths1 le 05-06-2017 à 19:58:13
Reply

Marsh Posté le 05-06-2017 à 19:57:30   

Reply

Marsh Posté le 05-06-2017 à 20:29:22    

Oui ben a+ le soleil, on part en voyage. :D  
 
Je ne vois pas non plus, j'avais pensé à un problème de référentiel mais ça ne semble pas non plus être le cas. Essaye toujours de récupérer la norme de ta vitesse voir si elle reste constante.


---------------
C'est en écrivant n'importe quoi qu'on devient n'importe qui.
Reply

Marsh Posté le 05-06-2017 à 21:10:54    

Merci quand même car je n avais pas pensé aux unités. J ai l impression que pour toutes les conditions initiales ca donne des droites.

Reply

Marsh Posté le 05-06-2017 à 21:32:26    

Je pense qu'il va falloir débuguer en mettant un n de 2 et un delta t de 0.2 puis ressortir tous les résultats intermédiaires pour voir là où ça déconne.


---------------
C'est en écrivant n'importe quoi qu'on devient n'importe qui.
Reply

Marsh Posté le 06-06-2017 à 16:43:10    

J'ai lancé une petite simu de ton algo et du coup j'ai vite trouvé l'erreur, dans ton calcul de la force ton vecteur unitaire est dans le sens contraire. En fait à cette étape tu calcules la force exercée par l'autre corps sur ton corps, il faut donc faire p1[k]-p2[k].
 
Et après ça fonctionne.


---------------
C'est en écrivant n'importe quoi qu'on devient n'importe qui.
Reply

Marsh Posté le 06-06-2017 à 22:57:53    

J ai fais la modification mais ça fait encore des droites.Sur ton ordi ça marche ?

Reply

Marsh Posté le 06-06-2017 à 23:11:45    

Points exportés sur Excel.
 

Code :
  1. import os
  2. import numpy as np
  3. from math import sqrt
  5. def mullist(coeff,L):
  6.    for k in range(0,len(L)):
  7.       L[k]=coeff*L[k]
  8.    return(L)
  10. def somlist(L,M):
  11.    N=[]
  12.    for k in range(0,len(L)):
  13.       N.append(L[k]+M[k])
  14.    return(N)
  16. def diflist(L,M):
  17.    N=[]
  18.    for k in range(0,len(L)):
  19.       N.append(L[k]-M[k])
  20.    return(N)
  22. def normevecteur(V):
  23.    normecarre=0
  24.    for k in range(0,len(V)):
  25.       normecarre+=V[k]**2
  26.    return(sqrt(normecarre))
  28. def Force_entre_deux_corps(m1,p1,m2,p2):#force de 2 sur 1 à un instant fixé
  29.    G=39.49070145
  30.    Force=[0,0,0]
  31.    for k in range(0,3):
  32.       Force[k]=(-G*m1*m2*(p1[k]-p2[k]))/(normevecteur(diflist(p2,p1))**3)
  33.    return(Force)
  35. def Force_entre_N_corps(j,m,position): #Résultante des forces exercées sur j à un instant fixé
  36.    Resultante=[0,0,0]
  37.    for k in range(0,len(position)):
  38.       if k!=j:
  39.          Force=Force_entre_deux_corps(m[j],position[j],m[k],position[k])
  40.          Resultante[0]=Resultante[0]+Force[0]
  41.          Resultante[1]=Resultante[1]+Force[1]
  42.          Resultante[2]=Resultante[2]+Force[2]
  43.    return((Resultante))
  45. def Acceleration(m,position):
  46.    L=[]
  47.    for k in range(0,len(position)):
  48.       L.append(mullist(1/m[k],Force_entre_N_corps(k,m,position))) #Calcul du vecteur acceleration pour chaque corps"
  49.    return(L)
  51. def fonction_de_passage(m,Y):
  52.    f=[Y[1],Acceleration(m,Y[0])]
  53.    return(np.array(f))
  55. def suite_Euler_pas_fixe(m,t_min,t_max,n,position_initiale,vitesse_initiale):
  56.    t=t_min
  57.    temps=[t_min]
  58.    h=(t_max-t_min)/n  #pas de discrétisation
  59.    Y_corps=np.array([position_initiale,vitesse_initiale])
  60.    Y_temps=[list(Y_corps)]
  61.    while t<=t_max:
  62.       Y_corps=Y_corps+h*fonction_de_passage(m,Y_corps)
  63.       Y_temps.append(list(Y_corps))
  64.       t+=h
  65.       temps.append(t)
  66.    return(np.array(Y_temps))
  68. m = [1, 3.00251*10**(-6)]
  69. position_initiale = [[0,0,0], [1,0,0]]
  70. vitesse_initiale = [[0,0,0], [0,6.27,0]]
  71. t_min = 0
  72. t_max = 1
  73. n = 2999
  75. a = suite_Euler_pas_fixe(m,t_min,t_max,n,position_initiale,vitesse_initiale)
  77. b = "xs;ys;zs;xt;yt;zt" + "\n"
  78. for i in range(0, len(a)):
  79.    b += str(a[i][0][0][0]) + ";" + str(a[i][0][0][1]) + ";" + str(a[i][0][0][2]) + ";" + str(a[i][0][1][0]) + ";" + str(a[i][0][1][1]) + ";" + str(a[i][0][1][2]) + "\n"
  81. os.chdir("C:/Users/***/Desktop" )
  82. txt = open("datas.txt", "w" )
  83. txt.write(b)
  84. txt.close()


 
Tout tourne.


---------------
C'est en écrivant n'importe quoi qu'on devient n'importe qui.
Reply

Marsh Posté le 06-06-2017 à 23:21:00    

Ca fait quoi comme trajectoire

Reply

Marsh Posté le 07-06-2017 à 12:18:00    

Non c est bon c est juste que j avais modifié un truc. Merci beaucoup

Reply

Marsh Posté le 07-06-2017 à 13:26:19    

Du coup si ton problème ne vennait pas de là, heureusement qu'il y avait une erreur dans le script sinon on aurait cherché longtemps...
 
Sans ce changement de signe tu as dans le cas Soleil - Terre, les deux corps qui se fuient donc ça finit en effet en droites mais ça commence quand même par une trajectoire arrondie.


---------------
C'est en écrivant n'importe quoi qu'on devient n'importe qui.
Reply

Sujets relatifs:

Leave a Replay

Make sure you enter the(*)required information where indicate.HTML code is not allowed