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)"""
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()
Pourriez-vous me dire pourquoi le premier programme ne fonctionne pas ?