尝试正则化二体系统的运动方程（稍后更新为三体）。我已经编写了一个 RK4 循环，并首先使用牛顿方程对其进行了尝试，并且效果很好。我现在正在尝试使用正则化方法来获得位置的二阶导数，即加速度。但是我只得到一条平坦的曲线。
from __future__ import division
import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
AU=1.5e11
a=AU
e=0
ms = 2E30
me = 5.98E24
mv=4.867E24
yr=3.15e7
h=100
mu1=ms*me/(ms+me)
mu2=ms*me/(ms+me)
G=6.67E11
#initial velocity
vi=np.sqrt(G*ms*(2/(a*(1-e))-1/a))
#initial position&velocity
sunpos=np.array([0,0,0])
earthpos=np.array([a*(1-e),0,0])
earthv=np.array([0,vi,0])
sunv=np.array([0,0,0])
#magnitude of position vector
norme=sum( (earthpos-sunpos)**2 )**0.5
norma=sum( (sunpos-earthpos)**2 )**0.5
#System of non-linear equations to get regularized position
def equations(p,qf):
q1, q2, q3, q4 = p
r1=np.sqrt((qf[0]-mu2)**2+qf[1]**2+qf[2]**2)
return q1**2-q2**2-q3**2+q4**2-qf[0]+mu2, 2*q1*q2-2*q3*q4-qf[1], 2*q1*q3+2*q2*q4-qf[2], q1**2+q2**2+q3**2+q4**2-r1
#system of non-linear equations to get regularized velocity
def equationsp(pp,qfp,qs):
dq1, dq2, dq3, dq4 = pp
D=qs[0]**2+qs[1]**2+qs[2]**2+qs[3]**2
return 2*(qs[0]*dq1-qs[1]*dq2-qs[2]*dq3+qs[3]*dq4)-D*qfp[0],2*(qs[1]*dq1+qs[0]*dq2-qs[3]*dq3-qs[2]*dq4)-D*qfp[1], 2*(qs[2]*dq1+qs[0]*dq3+qs[3]*dq2+qs[1]*dq4)-D*qfp[2],2*(qs[3]*dq1-qs[2]*dq2+qs[1]*dq3-qs[0]*dq4)
#Actual regularization method to get acceleration
def DotR(qf,dqf):
r1=np.sqrt((qf[0]-mu2)**2+qf[1]**2+qf[2]**2)
r2=np.sqrt((qf[0]+mu1)**2+qf[1]**2+qf[2]**2)
Q1,Q2,Q3,Q4=fsolve(equations, (1, 1, 1, 1),qf)
qf[0]=Q1**2-Q2**2-Q3**2+Q4**2+mu2
qf[1]=2*Q1*Q2-2*Q3*Q4
qf[2]=2*Q1*Q3+2*Q2*Q4
r1=Q1**2+Q2**2+Q3**2+Q4**2
D=4*(Q1**2+Q2**2+Q3**2+Q4**2)
dQ1,dQ2,dQ3,dQ4=fsolve(equationsp, (1, 1, 1, 1),(dqf,(Q1, Q2, Q3, Q4)))
dq0=2*(Q1*dQ1-Q2*dQ2-Q3*dQ3+Q4*dQ4)
dq1=2*(Q2*dQ1+Q1*dQ2-Q4*dQ3-Q3*dQ4)
dD=8*(Q1*dQ1+Q2*dQ2+Q3*dQ3+Q4*dQ4)
ddq0=(D**3*(qf[0]-mu1/r1**3*(qf[0]-mu2)-mu2/r2**3*(qf[0]+mu1))+dD*dq0+2*D**2*dq1)/D
ddq1=(D**3*(qf[1]-mu1*qf[1]/r1**3-mu2*qf[1]/r2**2)+dD*dq1+2*D**2*dq0)/D
ddq2=(D**3*(-mu1*qf[2]/r1**3-mu2*qf[2]/r2**3))/D
return np.array([ddq0,ddq1,ddq2])
#initial acceleration
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])
t=0
T=10**3
#Runge-kutta while loop
while t<T:
k1v1=h*earthv
k1v2=h*sunv
k1a1=h*eartha
k1a2=h*suna
earthpos=earthpos+.5*k1v1
sunpos=sunpos+.5*k1v2
earthv=earthv+.5*k1a1
sunv=sunv+.5*k1a2
norme=sum( (earthpos-sunpos)**2 )**0.5
norma=sum( (sunpos-earthpos)**2 )**0.5
print(DotR(earthpos,earthv))
eartha=DotR(earthpos,earthv)
suna=DotR(sunpos,sunv)
k2v1=h*earthv
k2v2=h*sunv
k2a1=h*eartha
k2a2=h*suna
earthpos=earthpos+.5*k2v1
sunpos=sunpos+.5*k2v2
earthv=earthv+.5*k2a1
sunv=sunv+.5*k2a2
norme=sum( (earthpos-sunpos)**2 )**0.5
norma=sum( (sunpos-earthpos)**2 )**0.5
eartha=DotR(earthpos,earthv)
suna=DotR(sunpos,sunv)
k3v1=h*earthv
k3v2=h*sunv
k3a1=h*eartha
k3a2=h*suna
earthpos=earthpos+k3v1
sunpos=sunpos+k3v2
earthv=earthv+k3a1
sunv=sunv+k3a2
norme=sum( (earthpos-sunpos)**2 )**0.5
norma=sum( (sunpos-earthpos)**2 )**0.5
eartha=DotR(earthpos,earthv)
suna=DotR(sunpos,sunv)
k4v1=h*earthv
k4v2=h*sunv
k4a1=h*eartha
k4a2=h*suna
earthpos=earthpos+h/6.0*(k1v1+2*k2v1+2*k3v1+k4v1)
sunpos=sunpos+h/6.0*(k1v2+2*k2v2+2*k3v2+k4v2)
earthv=earthv+h/6.0*(k1a1+2*k2a1+2*k3a1+k4a1)
sunv=sunv+h/6.0*(k1a2+2*k2a2+2*k3a2+k4a2)
norme=sum( (earthpos-sunpos)**2 )**0.5
norma=sum( (sunpos-earthpos)**2 )**0.5
eartha=DotR(earthpos,earthv)
suna=DotR(sunpos,sunv)
xarray.append(earthpos[0])
yarray.append(earthpos[1])
# print(earthpos,sunpos)
# zarray.append(earthpos[2])
t=t+h
plt.plot(xarray,yarray)
plt.savefig('orbit.png')
问题似乎是打印时 eartha 和 suna 没有“更新”我每次迭代都获得相同的加速度，但我不确定是什么原因造成的。