Méthode d'Euler trajectoire satellite

Messages : 2

Inscription : 26 sept. 2020 17:23

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

Méthode d'Euler trajectoire satellite

Message par AlexQui » 02 nov. 2021 13:24

Bonjour,

J'ai programmé la méthode d'Euler sur Python pour un vecteur avec une fonction f qui renvoie la représentation d'un satellite autour de la Terre en prenant en compte les frottements. J'ai le code suivant :

Code : Tout sélectionner

import numpy as np
from pylab import plot,show
from math import sqrt

G, M_T, R_T,Cx = 6.67*10**(-11), 5.972e24, 6.4e6,0.47

def P(h):
    """
    argument : altitude en km : renvoie la masse volumique de l'air correspondant à l'altitude donnée en argument
    """
    val=np.array([[0,50,100,150,200,250,300,350,400,450,500,550,600,650,700,750,800,850,900,950,1000],[1.29e-3,8.261e-7,5.173e-10,1.561e-12,1.583e-13,2.644e-14,5.63e-15,1.382e-15,3.868e-16,1.308e-16,5.704e-17,3.163e-17,2.045e-17,1.428e-17,1.034e-17,7.625e-18,5.691e-18,4.289e-18,3.263e-18,2.507e-18,1.945e-18]])
    m,M=[],0
    for k in range(0,21):
        m.append(abs(h-val[0][k])) #compare l'altitude donnée en argument avec la valeur la plus proche
    M=np.min(m)
    for k in range(0,21):
        if M==m[k]:
            return val[1][k]

def f(vecteur,fonction2):
    """
    argument : une matrice colonne qui contient les CI de la position et de la vitesse en x et y et la fonction P : renvoie la dérivée et l'accélération selon x et y
    """
    Fx=-(G*M_T*vecteur[0])/((vecteur[0]**2+vecteur[1]**2)**(3/2)) #force de gravitation en x
    Fy=-(G*M_T*vecteur[1])/((vecteur[0]**2+vecteur[1]**2)**(3/2)) #force de gravitation en y
    Frx=((fonction2(vecteur[0]))*0.58*Cx*sqrt(vecteur[2]**2+vecteur[3]**2)*vecteur[2])/2 #force de trainée en x
    Fry=((fonction2(vecteur[1]))*0.58*Cx*sqrt(vecteur[2]**2+vecteur[3]**2)*vecteur[3])/2 #force de trainée en y
    return([vecteur[2],vecteur[3],Fx+Frx/(43.6),Fy+Fry/(43.6)])

def euler(x0,y0,vx0,vy0,n,tf,fonction2):
    h=tf/float(n)
    v=[np.array(([x0,y0,vx0,vy0]))]
    t=0
    while t<tf:
        new=np.dot(h,f(v[-1],fonction2(x0)))+v[-1]
        v.append(new)
        t=t+h
    return v

def extract(fonction1,x0,y0,vx0,vy0,n,tf,fonction2):
    """
    argument : une fonction (ici euler), les CI de la position et de la vitesse en x et y, la fonction P : renvoie la trajectoire du débris
    """
    X=[0]
    Y=[0]
    E=fonction1(x0,y0,vx0,vy0,n,tf,P)
    for k in range(0,len(E)):
        X.append(E[k][0])
        Y.append(E[k][1])
    plot(X,Y)
    show()

"""extract(euler,400e3 + R_T, 0, 0, 7.6e3, 86400, 86400,P)"""
Ce qu'il y'a entre guillemets à la fin est ce que j'exécute. L'erreur envoyée est 'numpy.float64' object is not callable.

Mon programme fonctionne sans la fonction P() et cela ressemble à ça :

Code : Tout sélectionner

from numpy import*
from pylab import*

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

def f(vecteur):
    """
    argument : une matrice colonne qui contient les CI de la position et de la vitesse en x et y : renvoie la dérivée et l'accélération selon x et y
    """
    Fx=-(G*M_T*vecteur[0])/((vecteur[0]**2+vecteur[1]**2)**(3/2)) #force de gravitation en x
    Fy=-(G*M_T*vecteur[1])/((vecteur[0]**2+vecteur[1]**2)**(3/2)) #force de gravitation en y
    Frx=((1.29e-20)*0.58*Cx*sqrt(vecteur[2]**2+vecteur[3]**2)*vecteur[2])/2 #force de trainée en x
    Fry=((1.29e-20)*0.58*Cx*sqrt(vecteur[2]**2+vecteur[3]**2)*vecteur[3])/2 #force de trainée en y
    return([vecteur[2],vecteur[3],Fx+Frx,Fy+Fry])

def euler(x0,y0,vx0,vy0,n,tf):
    h=tf/float(n)
    v=[array(([x0,y0,vx0,vy0]))]
    t=0
    while t<tf:
        new=dot(h,f(v[-1]))+v[-1]
        v.append(new)
        t=t+h
    return v

def extract(euler,x0,y0,vx0,vy0,n,tf):
    """
    argument : une fonction (ici euler) et les CI de la position et de la vitesse en x et y : renvoie la trajectoire du débris
    """
    X=[0]
    Y=[0]
    E=euler(x0,y0,vx0,vy0,n,tf)
    for k in range(0,len(E)):
        X.append(E[k][0])
        Y.append(E[k][1])
    plot(X,Y)
    show()
La valeur (1.29e-20) représente le coefficient de trainée et la fonction P() dans le programme plus haut modifie cette valeur selon l'altitude du satellite pour rendre cela plus complexe. La faible valeur de ce coefficient est due au fait que pour observer quelque chose il faut avoir une faible valeur de ce coefficient sinon le satellite tombe.

Pourriez-vous me dire pourquoi le premier programme ne fonctionne pas ?
2020-2021 : MPSI Chrestien de Troyes
2021-2022 : MP Chrestien de Troyes

Messages : 2

Inscription : 26 sept. 2020 17:23

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

Re: Méthode d'Euler trajectoire satellite

Message par AlexQui » 02 nov. 2021 16:47

J'ai trouvé mon erreur, il ne faut pas mettre d'argument dans la fonction f dans la fonction euler, il faut simplement mettre fonction2.
2020-2021 : MPSI Chrestien de Troyes
2021-2022 : MP Chrestien de Troyes

Messages : 9679

Inscription : 30 juil. 2008 16:59

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

Re: Méthode d'Euler trajectoire satellite

Message par fakbill » 04 nov. 2021 11:38

Même si ton code tourne c'est voué à l'échec :(
Tu vas vite te rendre compte que, même sans frottement, la méthode d'Euler ne va pas te donner des trajectoires fermées. C'est du au fait que cette méthode n'est pas adaptée pour résoudre ce genre de problème. Utiliser un algo de type RK serait déjà plus pertinant mais pas encore vraiment adapté. Il existe des intégrateurs qui sont fait pour ce genre de problème : regarde un peu la biblio sur le sujet.
https://en.wikipedia.org/wiki/Verlet_integration
https://en.wikipedia.org/wiki/Leapfrog_integration

Quel est ton but? Faire de la physique? des maths appliquées? Pourquoi recoder des choses qui existent, bien faites, dans scipy (les méthodes de résolutions d'équa diff)?
Pas prof.
Prépa, école, M2, thèse (optique/images) ->ingé dans le privé.

Messages : 41

Inscription : 08 mars 2020 19:07

Profil de l'utilisateur : Élève de CPGE

Re: Méthode d'Euler trajectoire satellite

Message par Sherlock12 » 05 nov. 2021 11:50

fakbill a écrit :
04 nov. 2021 11:38
Quel est ton but? Faire de la physique? des maths appliquées?
Ou reproduire la prouesse évoquée dans un excellent film ?

(cf. https://youtu.be/v-pbGAts_Fg )
2020/2021 : PCSI 834 Masséna
2021/2022 : PC Masséna
2022/2023 : 5/2
2023/2026 : Télécom SudParis

Répondre