Моделирование траектории движения pl anet вокруг Солнца. Эллипс не закрывается - PullRequest
Новая
       116

Моделирование траектории движения pl anet вокруг Солнца. Эллипс не закрывается

2 голосов
/ 27 января 2020

Я пытался создать проекцию планет вокруг Солнца. Используя RungeKutta, я пытаюсь спроектировать и создать граф matplotlib. Тем не менее, внешний эллипс не закрывается. Не могли бы вы помочь мне, где именно я делаю ошибку?

Использовал Рунге Кутту, чтобы найти вектор R, имеющий положение как функцию времени. И когда я использую график, чтобы нарисовать траекторию, он не образует эллипс, который я должен найти. 

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)

Вывод:

enter image description here

Ответы [ 2 ]

4 голосов
/ 27 января 2020

Вы сделали не очень редкую ошибку копирования-вставки. В 

k4=F(t+h,r+h*k3,theta+h*l3/2)[0]
l4=F(t+h,r+h*k3,theta+h*l3/2)[1]

имеется ошибочное деление на два. Обратите внимание также на отсутствующее F в строке l2.

Вы можете сократить эти вычисления и избежать избыточных или дублированных вычислений. и уменьшить количество мест для ошибок наполовину, используя назначения кортежей, как в 

k4,l4 = F(t+h,r+h*k3,theta+h*l3)

и более поздних 

ri, thetai = RK4(F1,h,ri,thetai,ti)

Изменение вычисления размера шага на h=(tb-ta)/N, как, вероятно, первоначально предполагалось, и используя tb=150 для получения замкнутого l oop, для выбора подразделений можно получить все более замкнутые орбиты с помощью 

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)

enter image description here

2 голосов
/ 27 января 2020

Метод Рунге-Кутты не может дать вам идеально замкнутую траекторию, потому что есть дрейф энергии (энергия постепенно уменьшается).

Если вам нужна энергия, чтобы оставаться близкой к начальному значению, вы можете используйте симплектический c интегратор. Вы можете прочитать больше об этом в ответе Криса Ракауцкаса "Что означает" symplecti c "[...]" .

Однако , как отметил Лутц Леманн, ваш временной шаг достаточно мал, поэтому не поэтому ваша траектория не видна замкнутой. Эта информация может быть интересна для других, поэтому я оставляю этот ответ здесь.

...