Forum |  HardWare.fr | News | Articles | PC | S'identifier | S'inscrire | Shop Recherche
1612 connectés 

  FORUM HardWare.fr
  Programmation
  Python

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

 
 


 Mot :   Pseudo :  
 
Bas de page
Auteur Sujet :

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

n°2301577
maths1
Posté le 04-06-2017 à 19:44:08  profilanswer
 

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.

mood
Publicité
Posté le 04-06-2017 à 19:44:08  profilanswer
 

n°2301589
maths1
Posté le 05-06-2017 à 10:46:03  profilanswer
 

Personne ne voit l'erreur ?

n°2301594
MaybeEijOr​Not
but someone at least
Posté le 05-06-2017 à 12:08:35  profilanswer
 

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.
n°2301600
maths1
Posté le 05-06-2017 à 16:28:59  profilanswer
 

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)

n°2301601
rat de com​bat
attention rongeur méchant!
Posté le 05-06-2017 à 16:36:21  profilanswer
 

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

n°2301602
maths1
Posté le 05-06-2017 à 17:01:11  profilanswer
 

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
n°2301603
MaybeEijOr​Not
but someone at least
Posté le 05-06-2017 à 17:37:22  profilanswer
 

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.
n°2301606
maths1
Posté le 05-06-2017 à 18:22:16  profilanswer
 

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.

n°2301609
MaybeEijOr​Not
but someone at least
Posté le 05-06-2017 à 18:36:06  profilanswer
 

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.
n°2301611
maths1
Posté le 05-06-2017 à 19:57:30  profilanswer
 

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
mood
Publicité
Posté le 05-06-2017 à 19:57:30  profilanswer
 

n°2301615
MaybeEijOr​Not
but someone at least
Posté le 05-06-2017 à 20:29:22  profilanswer
 

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.
n°2301618
maths1
Posté le 05-06-2017 à 21:10:54  profilanswer
 

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.

n°2301619
MaybeEijOr​Not
but someone at least
Posté le 05-06-2017 à 21:32:26  profilanswer
 

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.
n°2301675
MaybeEijOr​Not
but someone at least
Posté le 06-06-2017 à 16:43:10  profilanswer
 

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.
n°2301708
maths1
Posté le 06-06-2017 à 22:57:53  profilanswer
 

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

n°2301711
MaybeEijOr​Not
but someone at least
Posté le 06-06-2017 à 23:11:45  profilanswer
 

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.
n°2301712
maths1
Posté le 06-06-2017 à 23:21:00  profilanswer
 

Ca fait quoi comme trajectoire

n°2301728
maths1
Posté le 07-06-2017 à 12:18:00  profilanswer
 

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

n°2301733
MaybeEijOr​Not
but someone at least
Posté le 07-06-2017 à 13:26:19  profilanswer
 

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.

Aller à :
Ajouter une réponse
  FORUM HardWare.fr
  Programmation
  Python

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

 

Sujets relatifs
[Python] Questions pratiques installation et utilisation d'OpenCVExporter données Python vers Gnuplot
Aidez moi svp urgent isn Pythonpb avec python
Probleme avec un programe python[Python] Parser un CSV vers un format custom
[Python]Utiliser Socket pour app de gestion réseauAide Programmation python
boucle while avec affectation en python[Python] Quel cours choisir ?
Plus de sujets relatifs à : résolution d'équation différentielle avec python

Forum MesDiscussions.Net, Version 2010.2
(c) 2000-2011 Doctissimo
Page générée en 0.084 secondes

Copyright © 1997-2022 Hardware.fr SARL (Signaler un contenu illicite / Données personnelles) / Groupe LDLC / Shop HFR