Discussion :

[quentinserveille]

Bonjour à tous,
J’aimerais modéliser le trajet d’un satellite dans l’espace, celui-ci étant seulement soumis aux forces de gravitation et à la pression de radiation. J’obtiens ces équations de trajectoires couplées. m(d²r/dt² - r(dθ/dt)²) = F1/r - GMm/r² et m(2*dr/dt*dθ/dt +F2/r²)

J’ai essayé de résoudre ce système et d’afficher la trajectoire grâce aux modules pyplot et odeint de python. Je ne pense pas avoir fait d’erreur lors de la détermination des équations de trajectoire. Pourtant les résultats que j’obtiens sont totalement incohérents. Sauriez-vous si j’ai fait une erreur dans mon code ?

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
 
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
 
from scipy import *
from scipy.integrate import odeint
 
c= 3*(10**8)
F1 = 2*(10**8)
F2 = 0.6
G = 6.67*(10**-11)
M = 1.989*(10**30)
m=1
 
def fonction(S,t):
    r = S[0]
    theta = S[1]
    v = S[2]
    omega = S[3]
 
    dr = v
    dtheta = omega
 
    dv = (F1 - G*M*m) /(m*(r**2)) + r*(dtheta**2)
    dw = F2/(m*(r**3)) - 2*dr*dtheta/r
 
     return [dr,dtheta,dv,dw]
 
t0 = 0
tmax = 1*60*60
npoint = 100000
t = np.linspace(t0, tmax, npoint) 
 
# Conditions initiales et résolution
x0= 149600000000
v0= 0.2*c
theta = np.pi
omega = 30000
 
syst_CI=np.array([x0,theta,v0,omega])                  # Tableau des CI
Sols=odeint(fonction,syst_CI,t)            # Résolution numérique des équations différentielles
 
# Récupération des solutions
x = Sols[:, 0]
u = Sols[:, 2]
 
 
plt.plot(x*np.cos(u),x*np.sin(u), linewidth=0.25)
plt.show()

[lg_53]

Lorsque vous poster un message, mettez les prérequis, on est pas censé les devinez ...

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
import numpy as np
from matplotlib import pyplot as plt
from scipy.integrate import odeint

Quand bien même on rajouterais le code ci dessus, votre code ne fonctionne toujours pas : il nous manque la valeur de c...
Et ensuite dans votre fonction (qui ne devrait pas avoir ce nom là d'ailleurs), il va nous manquer F1 et F2 ...


Fournissez un code qui tourne si vous voulez qu'on vous aide. Même si le résultat produit n'est pas celui escompté, c'est tout de même la 1ere des choses à faire.

[quentinserveille]

Lorsque vous poster un message, mettez les prérequis, on est pas censé les devinez ...

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
import numpy as np
from matplotlib import pyplot as plt
from scipy.integrate import odeint

Quand bien même on rajouterais le code ci dessus, votre code ne fonctionne toujours pas : il nous manque la valeur de c...
Et ensuite dans votre fonction (qui ne devrait pas avoir ce nom là d'ailleurs), il va nous manquer F1 et F2 ...


Fournissez un code qui tourne si vous voulez qu'on vous aide. Même si le résultat produit n'est pas celui escompté, c'est tout de même la 1ere des choses à faire.

Bonjour, je ne m'étais pas aperçu que le code ne fonctionnait pas, je l'ai modifié dans le post original il devrait être complet maintenant. merci

[lg_53]

Ok et quelle est le problème ? Que ca diverge ?

Si vous vous arrêtez à un t plus petit (genre 2), en mettant juste 1000 points, le résultat n'est pas celui attendu ?

Quel serait le résultat attendu ? Une ellipse ?

Lorsque vous présenter vos équations dans votre 1er message, la 2eme est égale à rien du tout ...

[quentinserveille]

Lorsque je mets t=2 je m'obtiens pas le résultat attendu. Et justement j'essaie de modéliser dans l'idéal le trajet d'un satellite sur plusieurs années et le temps étant en seconde il faut que le programme fonctionne pour des valeurs de t très grande. Le résultat que je souhaite obtenir ressemblerait à une spirale autour du centre, partant de la position initiale et s'éloignant du centre petit à petit.

Les équations exactes sont : m(d²r/dt² - r(dθ/dt)²) = F1/r - GMm/r² qui représente la projection de la seconde loi de newton sur Ur (radiale) en coordonnées polaires. la deuxième est m(2*dr/dt*dθ/dt + r*d²θ/dt²) = F2/r² qui est la projection sur Uθ.

[lg_53]

Si moi je traduis tes equations je trouve cela :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
    dv = (F1 - G*M*m/r) /(m*r) + r*(dtheta**2)
    dw = F2/(m*(r**2)) - 2*dr*dtheta/r

Donc déjà doit avoir un souci là. Mais ca reste bizarre d'avoir d'un coté F1/r et de l'autre F2/r² ...

Ensuite c'est pas parce que tu veux pour un t grand que ca doit t'empêcher de regarder pour un t petit ! D'ailleurs tu devrais commencer avec un t petit ! Car si ca ne marche pas pour un t petit, aucune chance que ca marche pour un t grand.

Enfin pourquoi tracer r en fonction de v ?
Car tu as écris :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
x = Sols[:, 0]
u = Sols[:, 2] ### 3eme composante !

Ca ne serait pas plutot r et theta que tu veux ?

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
x = Sols[:, 0]
u = Sols[:, 1]