In [94]:
import numpy
from math import log
from matplotlib import pyplot
%matplotlib inline
In [123]:
T = 3600*24*27.3
dt = 100
mT = 5.972*10**24
mL = 7.348*10**22
G = 6.67408*10**-11
dataL = numpy.empty((N,4))
dataT = numpy.empty((N,4))
dataL0 = [-362600000,0,0,-1083.4] # conditions initiales
dataT0 = [0,0,0,0]
dataL[0] = dataL0
dataT[0] = dataT0
Codons la méthode "Runge-Kutta 4"
La fonction ci-dessous ne sera pas généralisable à plus de deux corps, ce qui évite de lourds calculs. La première ligne de out est la trajectoire de la Terre, la seconde celle de la lune. Nous avons pris u' = f(u) avec u = [x , y , vx , vy] Nous traiterons simultanément les trajectoires de la lune et de la Terre.
In [40]:
def f(uT,uL):
out = numpy.ones((2,4))
out[0] = numpy.array([uT[2],uT[3],-G*mL*(uT[0]-uL[0])/((uT[0]-uL[0])**2+(uT[1]-uL[1])**2)**(3/2),-G*mL*(uT[1]-uL[1])/((uT[0]-uL[0])**2+(uT[1]-uL[1])**2)**(3/2)])
out[1] = numpy.array([uL[2],uL[3],-G*mT*(uL[0]-uT[0])/((uT[0]-uL[0])**2+(uT[1]-uL[1])**2)**(3/2),-G*mT*(uL[1]-uT[1])/((uT[0]-uL[0])**2+(uT[1]-uL[1])**2)**(3/2)])
return out
Ci-dessous la méthode "Runge-Kutta 4" appliquée directement aux trajectoires des deux astres.
In [41]:
def rk4(f,u,v,dt):
k1 = f(u,v)
k2 = f(u+dt/2*k1[0,:],v+dt/2*k1[1,:])
k3 = f(u+dt/2*k2[0,:],v+dt/2*k2[1,:])
k4 = f(u+dt*k3[0,:],v+dt*k3[1,:])
out = numpy.zeros((2,4))
out[0] = u+dt/6*(k1[0,:]+2*k2[0,:]+2*k3[0,:]+k4[0,:])
out[1] = v+dt/6*(k1[1,:]+2*k2[1,:]+2*k3[1,:]+k4[1,:])
return out
On calcule maintenant la trajectoire. Notons que dataL[t] est le vecteur uL qui est l'input de f (idem pour dataT[t]).
In [42]:
N = 1+int(T/dt)
for t in range(1,N):
W = rk4(f,dataT[t-1],dataL[t-1],dt)
dataT[t] = W[0]
dataL[t] = W[1]
posLT = dataL[:,0:2] - dataT[:,0:2] #position de la lune dans le repère de la Terre
In [43]:
# Plot the position of the moon and the earth in the initial frame
pyplot.figure(figsize=(10,6))
pyplot.grid(True)
pyplot.xlabel('$x$')
pyplot.ylabel('$y$')
pyplot.plot(dataL[:,0],dataL[:,1], 'b-', label='RK4')
pyplot.plot(dataT[:,0],dataT[:,1], 'b-', label='RK4')
pyplot.title('Position of the moon and the earth in the initial frame')
pyplot.legend();
pyplot.show();
In [44]:
# Plot the position of the moon in the rest frame of the earth
pyplot.figure(figsize=(10,6))
pyplot.grid(True)
pyplot.xlabel('$x$')
pyplot.ylabel('$y$')
pyplot.plot(posLT[:,0],posLT[:,1], 'b-', label='RK4')
pyplot.title('Position of the moon in the rest frame of the earth')
pyplot.legend();
pyplot.show();
Problèmes éventuels
In [45]:
# Compute the parameters a,b,e
minX = 0
for i in range(1,N-1):
if (posLT[i+1,0]>=posLT[i,0] and posLT[i,0]<posLT[i-1,0]):
minX = posLT[i,0]
break
maxX = 0
for i in range(1,N-1):
if (posLT[i+1,0]<=posLT[i,0] and posLT[i,0]>posLT[i-1,0]):
maxX = posLT[i,0]
break
minY = 0
for i in range(1,N-1):
if (posLT[i+1,1]>=posLT[i,1] and posLT[i,1]<posLT[i-1,1]):
minY = posLT[i,1]
break
maxY = 0
for i in range(1,N-1):
if (posLT[i+1,1]<=posLT[i,1] and posLT[i,1]>posLT[i-1,1]):
maxY = posLT[i,1]
break
print("minX =", minX)
print("maxX =", maxX)
print("minY =", minY)
print("maxY =", maxY)
a = abs(maxX-minX)/2
b = abs(maxY-minY)/2
e = numpy.sqrt(1-b**2/a**2)
print("a =", a)
print("b =", b)
print("e =", e)
In [46]:
# Compute the period
for i in range(2,N-1):
if (posLT[i+1,0]>=posLT[i,0] and posLT[i,0]<posLT[i-1,0]):
period = i*dt
break
print("period in seconds =", period)
print("period in 24h-days =", period/24/3600)
In [62]:
# plagiat ?
def get_diffgrid(data_current, data_fine, dt):
N_current = len(data_current[:,0])
N_fine = len(data_fine[:,0])
grid_size_ratio = int(numpy.ceil(N_fine/N_current))
diffgrid = dt * numpy.sum( numpy.abs(\
data_current[:,0]- data_fine[::grid_size_ratio,0]))
return diffgrid
In [139]:
# Show the order of the method
r = 2
h = 10
dt_values = numpy.array([h, r*h, r**2*h, r**3*h,r**4*h,r**5*h,r**6*h,r**7*h,r**8*h])
#dt_values = numpy.array([1000,300,100,30,10]) #,3,1
dataL_values = numpy.empty_like(dt_values, dtype=numpy.ndarray)
for i, dt in enumerate(dt_values):
N = int(T/dt)+1
dataL = numpy.empty((N, 4))
dataT = numpy.empty((N, 4))
dataL[0] = dataL0
dataT[0] = dataT0
for t in range(1,N):
W = rk4(f, dataT[t-1], dataL[t-1], dt)
dataT[t] = W[0]
dataL[t] = W[1]
dataL_values[i] = dataL
In [140]:
# Compute the differences
diffgridL = numpy.empty_like(dt_values, dtype=numpy.ndarray)
for i, dt in enumerate(dt_values):
diffgridL[i] = get_diffgrid(dataL_values[i], dataL_values[0], dt)
In [141]:
print(diffgridL)
In [146]:
alpha = numpy.empty_like(dt_values, dtype=numpy.ndarray)
for i in range(0, len(dt_values)-2):
alpha[i] = (log(get_diffgrid(dataL_values[i+2], dataL_values[i+1], dt_values[i+2]))
- log(get_diffgrid(dataL_values[i+1], dataL_values[i], dt_values[i+1]))) / log(r)
print(alpha)
In [143]:
pyplot.figure(figsize=(6,6))
pyplot.grid(True)
pyplot.xlabel(r'$\Delta t$', fontsize=18)
pyplot.ylabel(r'$L_1$-norm of the grid differences', fontsize=18)
pyplot.xlim(1e-4,1)
pyplot.ylim(1e-4,1)
pyplot.axis('equal')
pyplot.loglog(dt_values[1:], diffgridL[1:], color='k', ls='--', lw=2, marker='o');
In [ ]: