Runge-Kutta方法无法给您完美闭合的轨迹,因为存在能量漂移(能量逐渐减小)。
我试图创建围绕太阳的行星的投影。我正在尝试使用RungeKutta来投影和创建matplotlib图。但是,椭圆向外没有闭合。您能帮助我,我到底在哪里出错?
使用龙格·库塔(Runge Kutta)查找位置为时间函数的向量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)
输出:
您犯了一个不太常见的复制粘贴错误。在[>]中有两个错误的除法
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
。
您可以通过使用类似in中的元组分配来缩短这些计算,避免冗余或重复的计算,并将错误的位置减少一半。
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
获得闭环,对于细分的选择,通过来获得逐渐闭合的轨道]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)
Runge-Kutta方法无法给您完美闭合的轨迹,因为存在能量漂移(能量逐渐减小)。
如果需要能量来保持接近初始值,则可以使用辛积分器。您可以在Chris Rackauckas's answer on "What does 'symplectic' mean[...]"中阅读有关此内容的更多信息。
[然而
,您的时间步长足够小,所以这不是为什么您的轨迹没有明显闭合的原因。该信息可能会让其他人知道,所以,我在这里留下这个答案。Runge-Kutta方法无法给您完美闭合的轨迹,因为存在能量漂移(能量逐渐减小)。