0

我一直在尝试在 python 中使用 rk4 解决两个身体问题。但是，我的情节有些问题，似乎无法获得轨道，而只是一条略微弯曲的线。

我试图改变步长和初始值，但它没有让我到任何地方。这是我的太阳和地球系统代码。

from __future__ import division 
import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt
#from VPython import * 
#from VPython.graph import *


AU=1.5e11
a=AU  #semi major axis
e=0.1 #eccentricity
ms = 2E30
me = 5.98E24
h=10**(-6) 
G=6.67e-11
step=600 #timestep

#sun=sphere(pos=vec(0,0,0),radius=0.1*AU,color=color.yellow)
#earth=sphere(pos=vec(1*AU,0,0),radius=0.1*AU)

sunpos=np.array([0,0,0])
earthpos=np.array([a,0,0])

#scene.range=1.3*AU
#scene.autoscale=0

earthv=np.array([0,35000,0])
sunv=np.array([0,-35000,0])

norme=np.sqrt(abs(earthpos[0]**2+earthpos[1]**2+earthpos[2]**2-sunpos[0]**2+sunpos[1]**2+sunpos[2]**2))
norma=np.sqrt(abs(-(earthpos[0]**2+earthpos[1]**2+earthpos[2]**2)+(sunpos[0]**2+sunpos[1]**2+sunpos[2]**2)))

eartha=G*ms*(earthpos-sunpos)/norme**3
suna=G*me*(sunpos-earthpos)/norma**3

xarray=[]
yarray=[]
zarray=[]
xarray.append(earthpos[0])
yarray.append(earthpos[1])
zarray.append(earthpos[2])

#trail = curve(color = earth.color)

#k1 = [vector(0,0,0),vector(0,0,0), vector(0,0,0), vector(0,0,0)] 
#k2 = [vector(0,0,0),vector(0,0,0), vector(0,0,0), vector(0,0,0)] 
#k3 = [vector(0,0,0),vector(0,0,0), vector(0,0,0), vector(0,0,0)] 
#k4 = [vector(0,0,0),vector(0,0,0), vector(0,0,0), vector(0,0,0)]
t=0
T=10**9
while t<T:
    k1v1=h*earthv
    k1v2=h*sunv
    k1a1=eartha
    k1a2=suna
    earthpos=earthpos+.5*k1v1
    sunpos=sunpos+.5*k1v2
    earthv=earthv+.5*k1a1
    sunv=sunv+.5*k1a2
    norme=np.sqrt(abs(earthpos[0]**2+earthpos[1]**2+earthpos[2]**2-sunpos[0]**2+sunpos[1]**2+sunpos[2]**2))
    norma=np.sqrt(abs(-(earthpos[0]**2+earthpos[1]**2+earthpos[2]**2)+(sunpos[0]**2+sunpos[1]**2+sunpos[2]**2)))
    eartha=G*ms*(earthpos-sunpos)/norme**3
    suna=G*me*(sunpos-earthpos)/norma**3


    k2v1=h*earthv
    k2v2=h*sunv
    k2a1=eartha
    k2a2=suna
    earthpos=earthpos+.5*k2v1
    sunpos=sunpos+.5*k2v2
    earthv=earthv+.5*k2a1
    sunv=sunv+.5*k2a2
    norme=np.sqrt(abs(earthpos[0]**2+earthpos[1]**2+earthpos[2]**2-sunpos[0]**2+sunpos[1]**2+sunpos[2]**2))
    norma=np.sqrt(abs(-(earthpos[0]**2+earthpos[1]**2+earthpos[2]**2)+(sunpos[0]**2+sunpos[1]**2+sunpos[2]**2)))
    eartha=G*ms*(earthpos-sunpos)/norme**3
    suna=G*me*(sunpos-earthpos)/norma**3

    k3v1=h*earthv
    k3v2=h*sunv
    k3a1=eartha
    k3a2=suna
    earthpos=earthpos+k3v1
    sunpos=sunpos+k3v2
    earthv=earthv+k3a1
    sunv=sunv+k3a2
    norme=np.sqrt(abs(earthpos[0]**2+earthpos[1]**2+earthpos[2]**2-sunpos[0]**2+sunpos[1]**2+sunpos[2]**2))
    norma=np.sqrt(abs(-(earthpos[0]**2+earthpos[1]**2+earthpos[2]**2)+(sunpos[0]**2+sunpos[1]**2+sunpos[2]**2)))
    eartha=G*ms*(earthpos-sunpos)/norme**3
    suna=G*me*(sunpos-earthpos)/norma**3

    k4v1=h*earthv
    k4v2=h*sunv
    k4a1=eartha
    k4a2=suna
    earthpos=earthpos+1/6*(k1v1+2*k2v1+2*k3v1+k4v1)
    sunpos=sunpos+1/6*(k1v2+2*k2v2+2*k3v2+k4v2)
    earthv=earthv+1/6*(k1a1+2*k2a1+2*k3a1+k4a1)
    sunv=sunv+1/6*(k1a2+2*k2a2+2*k3a2+k4a2)
    norme=np.sqrt(abs(earthpos[0]**2+earthpos[1]**2+earthpos[2]**2-sunpos[0]**2+sunpos[1]**2+sunpos[2]**2))
    norma=np.sqrt(abs(-(earthpos[0]**2+earthpos[1]**2+earthpos[2]**2)+(sunpos[0]**2+sunpos[1]**2+sunpos[2]**2)))
    eartha=G*ms*(earthpos-sunpos)/norme**3
    suna=G*me*(sunpos-earthpos)/norma**3
    xarray.append(earthpos[0])
    yarray.append(earthpos[1])
    zarray.append(earthpos[2])

    t=t+step

plt.plot(xarray,yarray)

这是我得到的情节。

地球“轨道”图。

4

1 回答 1

0

  • 请再次阅读欧几里得距离是什么。在最简单的python实现中

    norme=sum( (earthpos-sunpos)**2 )**0.5

    也就是说，首先计算坐标差，然后将它们平方，然后将它们相加，然后取和的平方根。

  • 此外，您似乎永远不会将加速度与步长相乘。

  • 力的方向似乎是向外的，而不是像吸引力那样向内的。

  • 在 RK4 步骤期间，基本状态不会改变，所有的修改都是相对于未改变的基本状态。

为了更正所有这些并使代码更紧凑，您可以使用以下想法：

将重力/加速度计算定义为

def accelerations(earthpos, sunpos):
    norme=sum( (earthpos-sunpos)**2 )**0.5
    gravit = G*(earthpos-sunpos)/norme**3
    suna = me*gravit
    eartha = -ms*gravit
    return eartha, suna

然后可以将 RK4 步骤缩短为

eartha, suna = accelerations(earthpos, sunpos)
k1v1 = h*earthv
k1v2=h*sunv
k1a1=h*eartha
k1a2=h*suna

eartha, suna = accelerations(earthpos+0.5*k1v1, sunpos+0.5*k1v2)
k2v1 = h*(earthv+0.5*k1a1)
k2v2=h*(sunv+0.5*k1a2)
k2a1=h*eartha
k2a2=h*suna

eartha, suna = accelerations(earthpos+0.5*k2v1, sunpos+0.5*k2v2)
k3v1 = h*(earthv+0.5*k2a1)
k3v2=h*(sunv+0.5*k2a2)
k3a1=h*eartha
k3a2=h*suna

等等

于 2019-07-29T12:56:35.733 回答