Méthode d'Euler ordre 2 trajectoire satellite

Un problème, une question, un nouveau théorème ?

Messages : 2

Inscription : 26 sept. 2020 17:23

Profil de l'utilisateur : Élève de lycée

Méthode d'Euler ordre 2 trajectoire satellite

Message par AlexQui » 23 sept. 2021 19:27

Bonjour,

Je souhaite programmer la représentation de la trajectoire d'un satellite autour de la Terre en Python, je me concentre pour le moment sur la trajectoire du satellite et non sur l'aspect graphique que je programmerai avec Tkinter par la suite. J'ai obtenu deux équations différentielles d'ordre 2, une selon x et une autre selon y en coordonnées cartésiennes à l'aide d'un PFD, qui sont les suivantes :

d²x/dt²+GM_T*x/((x²+y²)^3/2)=0
d²y/dt²+GM_T*y/((x²+y²)^3/2)=0

(j'ai simplifié par la masse m du satellite et j'ai négligé les frottements pour le moment) avec G la constante de gravitation et M_T la masse de la Terre.
J'ai appliqué la méthode d'Euler et j'ai le programme suivant :

Code : Tout sélectionner

from scipy import*
from pylab import*
 
def euler(x0,y0,n,tf):
    T=[0]
    X=[0]
    Y=[0]
    h=tf/float(n)
    x,y,G,M_T=x0,y0,6.67*10**(-11),5.972*10**24
    t=0
    while t<tf:
        x+=h*(-G*M_T*x)/((x**2+y**2)**3/2)
        y+=h*(-G*M_T*y)/((x**2+y**2)**3/2)
        t+=h
        X.append(x)
        Y.append(y)
        T.append(t)
    plot(Y,X)
    show()
Mon programme ne fonctionne pas (il renvoi une droite) et je comprends bien sûr que mon problème est dû au fait que je ne prends pas en compte le temps t dont dépendent x et y mais je ne vois pas comment faire. Pourriez-vous m'éclairer ?
2020-2021 : MPSI Chrestien de Troyes
2021-2022 : MP Chrestien de Troyes

Messages : 3

Inscription : 06 sept. 2021 19:51

Profil de l'utilisateur : Professionnel

Re: Méthode d'Euler ordre 2 trajectoire satellite

Message par JohanB » 23 sept. 2021 23:14

Il n'y a pas de vitesse initiale dans tes données. Pas étonnant que le satellite s'écrase sur Terre.

En terme informatique, valorise à tout prix la lisibilité de ton code, quitte à être un peu redondant. Il y avait dans ton code originel, une erreur de priorité entre l'opérateur puissance ( ** en Python ) et la division. Il fallait parenthéser le 3/2 pour que celui ci soit bien pris en compte. On cubait la distance avec ta version.

Code : Tout sélectionner

import matplotlib.pyplot
from math import sqrt

G, M_T, R_T = 6.67 * 10 ** (-11), 5.972 * 10 ** 24, 6.4 * 10 **6


def sat_in_geocentric_ref(x0, y0, vx0, vy0, n, tf):
    """ Le calcul est mené dans le cadre de la mécanique Newtonienne dans le référientiel géocentrique.
    x0,y0  correspondent aux coordonnées de départ
    vx0, vy0 correspondent aux cordonnées du vecteur vitesse du satellite
    """
    T = [0]
    X = [0]
    Y = [0]
    VX = [0]
    VY = [0]
    V2X = [0]
    h = tf / float(n)

    x, y, vx, vy = x0, y0, vx0, vy0
    t = 0
    while t < tf:
        d2 = x ** 2 + y ** 2
        x_n, y_n = x/sqrt(d2), y/sqrt(d2)  # vecteur normalisé
        vx += h * (-G * M_T * x_n) / d2
        vy += h * (-G * M_T * y_n) / d2
        x += h * vx
        y += h * vy
        t += h

        X.append(x)
        Y.append(y)
        VX.append(vx)
        VY.append(vy)
        T.append(t)
        if d2 < R_T ** 2:
            print("CRASH")
            t = tf
    return X, Y, T


if __name__ == "__main__":
    #Valeurs pour l'orbite géostationnaire
    X,Y,T= sat_in_geocentric_ref(36000* 10 ** 3 + R_T, 0, 0, 3*10**3, 86400, 86400)
    matplotlib.pyplot.plot(Y,X)
    matplotlib.pyplot.show()

Messages : 2

Inscription : 26 sept. 2020 17:23

Profil de l'utilisateur : Élève de lycée

Re: Méthode d'Euler ordre 2 trajectoire satellite

Message par AlexQui » 24 sept. 2021 19:01

Je n'imaginais pas avoir une telle réponse avec un code qui fonctionne très bien, je vais essayer de le comprendre. Merci !
2020-2021 : MPSI Chrestien de Troyes
2021-2022 : MP Chrestien de Troyes

Messages : 4

Inscription : 13 oct. 2019 16:35

Profil de l'utilisateur : Élève de lycée

Re: Méthode d'Euler ordre 2 trajectoire satellite

Message par autobox » 24 sept. 2021 20:42

Il y a aussi quelques améliorations qu'on peut y apporter

* initialiser directement les listes avec r+1 zéros et les remplir dans la boucle sera plus efficace que de faire à chaque fois un append et réallouer toute la mémoire.

* précalculer G*M_T et R_T**2

* écrire plutôt x, y = x + vx*h, y + vy*h est susceptible de faire comprendre à l'interpréteur qu'il peut paralléliser ces instructions. Efficacité marginale ici, mais ça peut servir ailleurs, je précise au cas où ;)

* remplacer le t = tf par un break évite une assignation, mais en fait c'est inutile si tu remplaces ton while par un for avec le bon nombre d'itérations, puisque t = i*h où i est l'itération.
Tu auras automatiquement t>=tf dès que i*h >= tf ie i >= tf/h, ie i >= ceil(tf/h). Donc la bonne borne est effectivement r = ceil(tf/h) (for i in range(r): ...)

* aussi, on remplit VX et VY mais on ne les renvoie pas et on n'en fait rien, ce qui est du gâchis de mémoire. V2X n'apparaît même pas dans la suite du code, il faut le dégager aussi :wink:

* remplacer les *10**... par une notation scientifique qui est faite pour ça et donne directement la mantisse et l'exposant à Python. L'avantage d'utiliser 10** est que la masse et le rayon restent des entiers ce qui peut être utile si tu fais du calcul exact (fractions) mais ici on se contente d'un float

Code : Tout sélectionner

import numpy as np
import math

G,M_T,R_T = 6.67e-11, 5.972e24, 6.4e6
K, RT2 = G*M_T, R_T**2


def sat_in_geocentric_ref(x0, y0, vx0, vy0, n, tf):
	""" Le calcul est mené dans le cadre de la mécanique Newtonienne dans le référientiel géocentrique.
	x0,y0  correspondent aux coordonnées de départ
	vx0, vy0 correspondent aux cordonnées du vecteur vitesse du satellite
	"""
	h = tf / float(n)
	r = int(math.ceil(tf/h))
	T = [i*h for i in range(r)]

	X = np.zeros(r+1) # remplaçable par [0]*(r+1) ici parce que 0 est une constante. Sinon gare au piège de la référence
	Y = np.zeros(r+1)
	#VX = np.zeros(r+1)
	#VY = np.zeros(r+1)
	#V2X = np.zeros(r+1)

	x, y, vx, vy = x0, y0, vx0, vy0	
	hK = h*K
	
	for i in range(r): # xrange si tu es encore sur Python 2
		d2 = x**2+y**2
		if d2 < RT2:
			print("CRASH")
			break # ou mieux : return X[:i],Y[:i],T[:i]
		
		d = d2**.5
		x_n, y_n = x/d, y/d  # vecteur normalisé
		vx, vy = vx - (hK*x_n)/d2, vy - (hK*y_n)/d2
		x, y = x + h*vx, y + h*vy
		
		X[i], Y[i] = x, y
		#VX[i], VY[i] = vx, vy
		
		
	return X, Y, T
	
	
	
X, Y, T = sat_in_geocentric_ref(36000e3 + R_T, 0, 0, 3e3, 86400, 86400)

for u in (T,X,Y): print(u)
* enfin dernière remarque, le break dans le for peut te faire retourner un tuple de listes ayant de nombreux zéros. Tu peux le remplacer par return X[:i],Y[:i],T[:i]

Répondre