Python蛙跳法计算地球轨道能量溢出及结果异常问题

编辑补充：目前我已得到地球轨道仿真的相关绘图，但总能量结果不符合预期，请问总能量是否应当在0值附近振荡？
感谢解答！

我目前正在完成数值方法课程作业的b部分题目：地球轨道仿真题目说明
目标是修改a部分已完成的代码，使其能够输出地球轨道的动能与势能，但运行时出现如下报错：

(8.12).py:31: RuntimeWarning: overflow encountered in multiply
  z_vec[i+1,:]=z_vec[i,:] + h*f(z_half_step, t_vec[i])

目前编写的动能、势能求解代码存在明显错误，核心问题是不知道如何正确修改a部分的代码来实现势能与动能的计算。

a部分可正常运行的代码 #

以下a部分代码运行后可得到5年周期内略偏离正圆的地球运行轨迹：

import matplotlib.pyplot as plt
import numpy as np
import scipy.integrate as spi


"""a"""


G = 6.6738*10**-11
M = 1.9891*10**30
h = 3600
y = 1.4710*10**11
vx = 3.0287*10**4

def LeapFrog(f, t_start, t_stop, z0, h):

    t_vec = np.arange(t_start, t_stop, h) 
    n = len(t_vec)
    d = len(z0)  
    z_vec = np.zeros((n, d))
    z_vec[0,:] = z0
    z_half_step=z_vec[0 , :] + (1/2)*h*f(z0,t_vec[0]) 
    
    
    for i in range(0, n - 1):
        z_vec[i+1,:]=z_vec[i,:] + h*f(z_half_step, t_vec[i])
        z_half_step += h*f(z_vec[i+1,:], t_vec[i])   


    return t_vec, z_vec, 


def f(z,t):   #dz/dt
    
    x=z[0]
    y=z[1] 
    vx=z[2] 
    vy=z[3] 
    r=np.sqrt(x**2+y**2)

    dz=np.zeros(4)
    
    dz[0]=vx 
    dz[1]=vy
    dz[2]=-G*M*x/r**3
    dz[3]=-G*M*y/r**3

    return dz

t_start = 0
t_stop = 24*365*h*5 #5 years
z0 = np.array([0,y,vx,0])
t_vec, z_vec = LeapFrog(f, t_start, t_stop, z0, h)



figure, ax = plt.subplots()

plt.title("Orbit of the Earth")
plt.plot(0,0,'yo', label = 'Sun positon')  
plt.plot(z_vec[:,0],z_vec[:,1], 'g', markersize=1, label='Earth trajectory')
ax.set_aspect('equal')

ax.set(xlim=(-1.55*10**11, 1.55*10**11), ylim = (-1.55*10**11, 1.55*10**11))
a_circle = plt.Circle((0, 0), 1.4710*10**11, fill=False,color='r')
ax.add_artist(a_circle)
ax.set_aspect('equal')

legend = plt.legend(['Sun position','Earths 5-year trajectory','Perfect circle'],loc='center left', bbox_to_anchor=(1, 0.5))

plt.show()



"""
We can see that the trajectory of the earth for 5 years is very slightly non-circular.
"""

b部分错误代码 #

完成b部分时错误修改了dz[2]和dz[3]项，导致代码无法正常运行，错误代码如下：

"""b"""

m = 5.9722*10**24
def f(z,t):   #dz/dt

    x=z[0]
    y=z[1] 
    vx=z[2] 
    vy=z[3] 
    r=np.sqrt(x**2+y**2)

    dz=np.zeros(4)

    dz[0]=vx 
    dz[1]=vy
    dz[2]=-G*M*m/r
    dz[3]=0.5*m*y**2

    return dz

t_start = 0
t_stop = 24*365*h*5 #5 years
z0 = np.array([0,y,vx,0])
LeapFrog(f, t_start, t_stop, z0, h)

解答 #

总能量取值说明 #

地球绕太阳公转属于束缚轨道，总能量（动能+引力势能）恒为负值，不会在0附近振荡：

• 动能计算公式为 $E_k=\frac{1}{2}m(v_x^2+v_y2)$，恒为正
• 引力势能计算公式为 $E_p=-\frac{GMm}{r}$（以无穷远为势能零点），恒为负
  束缚轨道下引力势能的绝对值始终大于动能，因此总能量为负。你用的蛙跳法（Leapfrog）是辛积分器，总能量只会在真实值附近做小幅度的有界振荡，不会出现无限制漂移。如果总能量在0附近甚至为正，要么是能量计算逻辑错误，要么是数值仿真发散，轨道已经变成非束缚的抛物线/双曲线轨道。

报错与代码错误原因 #

核心错误是误将能量计算写进了运动方程的导数项。
你a部分写的运动方程f(z,t)是完全正确的，描述的是万有引力产生的加速度，不需要做任何修改。动能和势能是仿真得到位置、速度结果后，通过后处理计算的导出量，不是用来更新轨道运动状态的变量。你错误替换了加速度项dz[2]和dz[3]，相当于完全改乱了地球的受力规则，导致数值计算快速发散，才会触发数值溢出的RuntimeWarning。

正确实现方法 #

不需要改动a部分的Leapfrog积分函数、运动方程f(z,t)，只需要在跑完轨道仿真、得到t_vec和z_vec之后，追加以下后处理代码即可计算能量：

# 地球质量
m = 5.9722*10**24
# 从仿真结果中提取位置、速度序列
x = z_vec[:, 0]
y = z_vec[:, 1]
vx = z_vec[:, 2]
vy = z_vec[:, 3]
r = np.sqrt(x**2 + y**2)
# 分别计算动能、势能、总能量
Ek = 0.5 * m * (vx**2 + vy**2)
Ep = -G * M * m / r
E_total = Ek + Ep

# 可选：绘制能量随时间变化曲线
plt.figure(figsize=(10,6))
plt.plot(t_vec, Ek, label='动能')
plt.plot(t_vec, Ep, label='引力势能')
plt.plot(t_vec, E_total, label='总能量')
plt.xlabel('仿真时间 (s)')
plt.ylabel('能量 (J)')
plt.legend()
plt.grid(alpha=0.3)
plt.show()

注：蛙跳法的速度定义在半时间步上，如果追求更高精度，可以用半时间步的速度值计算动能，不过以1小时为步长仿真地球轨道时，直接用整步输出的速度计算的能量误差完全可以满足课程作业要求。

内容的提问来源于stack exchange，提问作者Grayham_T