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]

[lg_53]

Vous avez des problèmes également sur vos paramètres.
v0 à 20% de la vitesse de la lumière c'est possible ca ?

Si je prends ca :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
t0 = 0
#tmax = 1*60*60
#tmax = 60
tmax=1.0e-11
npoint = 1000
t = np.linspace(t0, tmax, npoint) 
 
# Conditions initiales et résolution
x0= 1.496e11
#v0= 0.2*c
v0 = 1.0
theta = np.pi
#omega = 30000
omega = 3

voyez que pour un temps très petit on fait déjà plus d'une révolution...



Et sinon faite votre plot avec '+' plutot qu'avec des lignes. Car tant que vos points sont aussi éloignés les uns des autres vous n'allez rien voir.

[quentinserveille]

J'étudie le projet breakthrough starshot, qui prévoit de propulser des voiles solaires très légères à très grande vitesse initiale pour qu'elles puissent atteindre alpha du centaure. La vitesse 0.2c est celle qu'ils souhaiteraient atteindre. https://breakthroughinitiatives.org/initiative/3

Je n'avais pas remarqué que je traçais avec v au lieu de theta, merci de m'avoir indiqué cette erreur !
Vous avez aussi raison concernant le F1/r² que j'avais oublié de recopier

Malheureusement en effectuant toutes les corrections, j'obtiens une trajectoire rectiligne pour des temps très courts comme très longs

[lg_53]

Ok mais meme si ca va très très vite je doute que 10^-11 secondes soit suffisante pour faire plus d'une orbite ...
Donc il y a des problèmes dans les ordres de grandeur des paramètres.

Essayer de mettre F2 à 0, et de prendre les paramètres de la Terre, pour voir déjà si dans ce cas classique si vous arrivez à reproduire correctement la dynamique de la rotation de la terre autour du soleil (et donc là en 365 jours)

[quentinserveille]

J'ai essayé mais j'ai toujours le même probleme, lorsque je regarde les valeurs de theta trouvées par odeint, theta est toujours constant ..

[lg_53]

Désolé mais ca ne suffit pas de dire ca marche pas.

Il faut présenter le(s) code(s) tester, commenter ce qu'on a fait par rapport au précédent, etc ...

Car moi je viens de faire en virant les forces de radiation, et en prenant les constantes terre soleil et ca marche bien. J'obtiens une belle ellipse que dessine la terre autour du soleil. Vous devriez être capable d'en faire autant.

Si votre code ne fonctionne pas dans ce cas simple, alors normal qu'il ne fonctionne pas non plus dans un cas plus compliqué ...

[quentinserveille]

J'ai trouvé mon erreur dans le cas de la trajectoire de la terre, elle résidait dans les conditions initiales, sur la vitesse angulaire initiale. J'obtiens moi aussi une ellipse.

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
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
 
from scipy.integrate import odeint
 
G = 6.67*(10**-11)
M = 1.989*(10**30)
m= 5.972*(10**24)
 
def fonction(S,t):
    r = S[0]
    theta = S[1]
    v = S[2]
    omega = S[3]
 
    dr = v
    dtheta = omega
 
    dv = -G*M*m/(r**2) /m + r*(omega**2)
    dw = - 2*v*omega/r
 
    return (dr,dtheta,dv,dw)
 
t0 = 0
tmax = 365*24*60*60
npoint = 1000
t = np.linspace(t0, tmax, npoint) 
 
# Conditions initiales et résolution
x0= 149600000000
v0= np.sqrt(G*Mt/x0)
theta = np.pi
omega = 2*np.pi/(365*24*60*60)
 
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[:, 1]
 
 
plt.polar(u,x)
plt.show()

Il me semble donc que je dois revoir la détermination des conditions initiales dans mon projet originel. merci pour votre aide

[lg_53]

Voilà c'était pas si dur ! Et voyez que j'avais moi même tester un code très semblable :

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
import numpy as np
import matplotlib.pyplot as plt
 
from scipy.integrate import odeint
 
 
c = 3e8       ### vitesse de la lumière (m/s)
G = 6.67e-11  ### constante gravitationnelle
M = 1.989e30  ### masse du soleil (Kg)
m = 5.972e24  ### masse de la terre (Kg)
 
def derivee_systeme(S,t):
    r, theta, v, w = S
 
    dv = - G*M / r**2 + r*w**2
    dw = - 2*v*w/r
 
    return [v,w,dv,dw]
 
# Discrétisation temporelle 
t0 = 0
tmax=365*24*3600 ## 1 an en secondes
npoint = 1000
t = np.linspace(t0, tmax, npoint) 
 
# Conditions initiales
r0= 1.496e11         ## distance terre soleil (en mètres)
theta0 = 0           ## l'angle de la position départ importe peu
omega0 = 1.99e-7     ## vitesse angulaire terrestre (en radian par seconde)
v0 = 0               ## vitesse radiale
syst_CI=np.array([r0,theta0,v0,omega0])    # Vecteur des CI
 
# Résolution du système différentiel
Sols=odeint(derivee_systeme,syst_CI,t)   # Résolution numérique des équations différentielles
 
# Récupération des solutions
r = Sols[:, 0]
theta = Sols[:, 1]
 
# Visualisation du réultat
plt.plot(r*np.cos(theta),r*np.sin(theta), '.')
plt.show()

Donc effectivement comme je le supputtais, les paramètres doivent obligatoirement avoir les bons ordre de grandeur sinon ca donne effectivement n'importe quoi. Là vous avez la vitesse de lancement de votre satellite. C'est une vitesse radiale ? une vitesse angulaire ? Simplement la norme et c'est à vous de la décomposer ?

[quentinserveille]

Cette subtilité n'a pas été abordée par les personnes gérant le projet, du moins lors de mes recherches je n'ai pas trouvé de précisions à ce sujet. Pour moi, c'est une vitesse radiale, qui pourrait être tellement forte que la voile solaire ne ressentirait plus visiblement les effets de la gravitation et obtiendrait ainsi une trajectoire rectiligne. Je ne sais pas si cette hypothèse est fondée scientifiquement, mais elle expliquerait pourquoi je n'obtient jamais de trajectoire en forme de spirale autour du soleil

[lg_53]

Essaie d'obtenir une trajectoire de ce type pour un satelitte "normal" déjà. C'est l'étape intremédiaire entre juste la terre autour du soleil, et ton satelite nouvelle génération. Une fois que tu auras cette intermède tu pourras bouger ses paramètres pour voir coment le système répond et éventuellement voir aussi si une augmentation de la vitesse initiale peut compenser par qqch (la diminuation du poids, mais peut etre autre chose) pour conserver ce profil de trajectoire.