I was trying to create the projection of planets around the sun. Using the RungeKutta I'm trying to project and create the matplotlib graph. However, the out ellipse is not closing. Could you please help me, where exactly I'm doing the mistake?
Used Runge Kutta to find a vector R having the position as a function of time. And when I use the plot to draw the trajectory it doesn't make an ellipse which is what I am supposed to find.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
#unités de normalisation
UA=149.59787e6 #distance moyenne Terre-Soleil
MASSE=6.0e24 #masse Terre
JOUR=24*3600
#données :
k=0.01720209895
G=(k**2) # constante de gravitation en km^3/kg/s²
r0= 1 # distance initiale terre soleil en m
Ms = 2.0e30/MASSE #masse Soleil
Mt = 1 #masse Terre
w0 = 30/(UA/JOUR) #vitesse angulaire en Km/s
C = r0*(w0**2)
m = (Ms*Mt/Ms+Mt) #masse réduite
K = G*m #on choisit K > 0
b = 2 #beta
def RK4(F, h, r, theta, t, *args):
k1=F(t,r,theta,)[0]
l1=F(t,r,theta,)[1]
k2=F(t+h/2,r+h*k1/2,theta+h*l1/2)[0]
l2=(t+h/2,r+h*k1/2,theta+h*l1/2)[1]
k3=F(t+h/2,r+h*k2/2,theta+h*l2/2)[0]
l3=F(t+h/2,r+h*k2/2,theta+h*l2/2)[1]
k4=F(t+h,r+h*k3,theta+h*l3/2)[0]
l4=F(t+h,r+h*k3,theta+h*l3/2)[1]
return [r +(h/6)*(k1+2*k2+2*k3+k4),theta +(h/6)*(l1+2*l2+2*l3+l4)]
def F1(t,r,theta):
return np.array([r[1],r[0]*theta[1]**2-K/r[0]**b]),np.array([theta[1],-2*r[1]*theta[1]/r[0]])
def RK4N(F1,N,r0,r1,theta0,theta1,ta,tb):
h=0.05
ri=np.array([r0,r1])
thetai=np.array([theta0,theta1])
ti=ta
R=np.zeros((N,2))
THETA=np.zeros((N,2))
lt=np.zeros(N)
lt[0], R[0],THETA[0] = ta , ri ,thetai
for i in range (1, N):
a = ri
ri = RK4(F1,h,ri,thetai,ti)[0]
thetai=RK4(F1,h,a,thetai,ti)[1]
ti=ti+h
R[i]=ri
THETA[i]=thetai
lt[i]=ti
return R,THETA,lt
def trace_position(F1,N,r0,r1,theta0,theta1,ta,tb):
R,THETA,lt=RK4N(F1,N,r0,r1,theta0,theta1,ta,tb)
lx,ly=[],[]
for i in range(N):
lx.append(R[i][0]*np.cos(THETA[i][0]))
ly.append(R[i][0]*np.sin(THETA[i][0]))
plt.plot(lx,ly)
plt.plot(0,0)
plt.show()
# rapport a/b
max_x,min_x,max_y,min_y= max(lx),min(lx),max(ly),min(ly)
rapport = (max_x-min_x)/(max_y-min_y)
print("rapport a/b = ",rapport)
if abs(rapport-1)< 10e-2:
print("le mouvement est presque circulaire")
else:
print("le mouvement est elliptique")
def trace_Ep(F1,N,r0,r1,theta0,theta1,ta,tb):
R,THETA,lt=RK4N(F1,N,r0,r1,theta0,theta1,ta,tb)
lEp = []
for i in range(N):
lEp.append(-K/R[i][0]**(b-1))
#fig, (ax1, ax2) = plt.subplots(1, 2)
#ax1.plot(lt,lEp)
#ax2.plot(R[:,0],lEp)
plt.plot(lt,lEp)
plt.show()
trace_position(F1,380,r0,0,0,w0,0,1)
Output:
You made a not so uncommon copy-paste error. There is an erroneous division by two left over in
k4=F(t+h,r+h*k3,theta+h*l3/2)[0]
l4=F(t+h,r+h*k3,theta+h*l3/2)[1]
Note also the missing F
in the l2
line.
You can shorten these computation and avoid redundant or duplicated computations, and reduce places for errors by one half, by using tuple assignments like in
k4,l4 = F(t+h,r+h*k3,theta+h*l3)
and later
ri, thetai = RK4(F1,h,ri,thetai,ti)
Changing the step size computation to h=(tb-ta)/N
as probably initially intended, and using tb=150
to get a closed loop, for a selection of subdivisions one gets the increasingly closed orbits via
for k in range(4):
N = [16,19,25,120][k]
plt.subplot(2,2,k+1,aspect='equal')
trace_position(F1,N,r0,0,0,w0,0,150)