基于Runge-Kutta法的单摆周期仿真异常：代码错误排查

问题根源 #

你的代码存在两个核心错误，导致仿真结果不符合预期：

1. 摆长计算错误 #

单摆小角度周期公式为：
T = 2π√(l/g)
若预期周期T=2秒，代入求解摆长l：
l = g*(T/(2π))² = g*(1/π)² = g/π²
你当前代码中l = g/pi，摆长是正确值的π倍，直接导致周期偏离预期。

2. 未正确实现Runge-Kutta方法 #

你写的迭代逻辑不是标准的Runge-Kutta（如RK4）步骤，自定义的更新规则（包括t==0的特殊分支）完全错误，导致数值积分精度极低，甚至会让系统能量不守恒，结果发散。

修正后的代码 #

以下是实现标准RK4方法的修正代码，同时修正摆长计算、优化变量命名：

import math
import matplotlib.pyplot as plt

g = 9.81
pi = math.pi
# 修正摆长计算，确保周期接近2秒
l = g / (pi ** 2)

t0 = 0.0
t = t0
t_list = [t0]

# v是角速度，theta是角位移
v0 = 0.0
v = v0
v_list = [v0]

theta0 = 0.1  # 小角度，符合单摆近似条件
theta = theta0
theta_list = [theta0]

dt = 0.02
tfinal = 20.0

# 定义微分方程：d(theta)/dt = v，d(v)/dt = -(g/l)*sin(theta)
def d_theta_dt(t, theta, v):
    return v

def dv_dt(t, theta, v):
    return -(g / l) * math.sin(theta)

while t <= tfinal:
    # RK4四个斜率计算
    k1_theta = d_theta_dt(t, theta, v)
    k1_v = dv_dt(t, theta, v)
    
    k2_theta = d_theta_dt(t + dt/2, theta + k1_theta*dt/2, v + k1_v*dt/2)
    k2_v = dv_dt(t + dt/2, theta + k1_theta*dt/2, v + k1_v*dt/2)
    
    k3_theta = d_theta_dt(t + dt/2, theta + k2_theta*dt/2, v + k2_v*dt/2)
    k3_v = dv_dt(t + dt/2, theta + k2_theta*dt/2, v + k2_v*dt/2)
    
    k4_theta = d_theta_dt(t + dt, theta + k3_theta*dt, v + k3_v*dt)
    k4_v = dv_dt(t + dt, theta + k3_theta*dt, v + k3_v*dt)
    
    # 用加权平均更新状态
    theta_new = theta + (dt/6) * (k1_theta + 2*k2_theta + 2*k3_theta + k4_theta)
    v_new = v + (dt/6) * (k1_v + 2*k2_v + 2*k3_v + k4_v)
    
    # 更新变量
    theta = theta_new
    v = v_new
    t += dt
    
    t_list.append(t)
    theta_list.append(theta)
    v_list.append(v)

# 绘制结果
plt.plot(t_list, theta_list)
plt.xlabel('时间 (s)')
plt.ylabel('角位移 (rad)')
plt.title('单摆角位移随时间变化（周期≈2秒）')
plt.show()

效果说明 #

修正后，仿真得到的单摆周期会接近预期的2秒，角位移曲线为稳定的简谐运动，不会出现原代码中的发散或周期偏移问题。

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