使用四阶Runge-Kutta法求解非线性摆微分方程遇异常结果求助

一、异常图像问题排查 #

大角度非线性摆的RK4求解出现“突然停止”，大概率是数值不稳定或RK4迭代逻辑错误，常见问题及修复方案：

常见错误点 #

1. 步长过大：大角度摆的周期比小角度摆更长（比如160°初始角的周期约是小角度的1.8倍），步长太大时RK4无法准确捕捉振荡细节，甚至导致数值发散。建议将步长h设置为0.01~0.05之间。
2. 一阶方程组导数计算错误：非线性摆的二阶方程需转化为一阶方程组：
   • θ' = ω
   • ω' = -g/L * sinθ
     若符号错误或变量映射错误，会导致运动方向异常。
3. 初始条件处理失误：确保初始角度已正确转为弧度（160°=8π/9≈2.7925弧度）。

正确实现代码 #

import numpy as np
import matplotlib.pyplot as plt

# 定义非线性摆的一阶导数方程组
def pendulum_dydt(t, y, g, L):
    theta, omega = y
    dtheta_dt = omega
    domega_dt = -(g / L) * np.sin(theta)
    return np.array([dtheta_dt, domega_dt])

# RK4单步迭代函数
def rk4_update(f, t, y, h, *args):
    k1 = h * f(t, y, *args)
    k2 = h * f(t + h/2, y + k1/2, *args)
    k3 = h * f(t + h/2, y + k2/2, *args)
    k4 = h * f(t + h, y + k3, *args)
    return y + (k1 + 2*k2 + 2*k3 + k4) / 6

# 参数与初始条件
g = 9.81
L = 1.0
theta0 = np.deg2rad(160)  # 初始角度转弧度
omega0 = 0.0
t_start = 0.0
t_end = 50.0
h = 0.01  # 足够小的步长保证数值稳定

# 初始化数组
t = np.arange(t_start, t_end, h)
y = np.zeros((len(t), 2))
y[0] = [theta0, omega0]

# 执行RK4迭代
for i in range(1, len(t)):
    y[i] = rk4_update(pendulum_dydt, t[i-1], y[i-1], h, g, L)

theta = y[:, 0]

# 绘制曲线
plt.figure(figsize=(10, 6))
plt.plot(t, np.rad2deg(theta))
plt.xlabel('时间 (s)')
plt.ylabel('摆角 (°)')
plt.title('非线性摆时间-摆角曲线')
plt.grid(True)
plt.show()

运行上述代码应该能得到预期的非均匀周期类正弦曲线。

二、输出时间-角度表格数据 #

保存为CSV文件（方便查看和后续处理） #

# 组合时间与摆角（转成角度）
data = np.column_stack((t, np.rad2deg(theta)))
# 保存到本地文件
np.savetxt(
    'pendulum_data.csv',
    data,
    header='时间(s),摆角(°)',
    delimiter=',',
    fmt='%.4f'
)

直接打印前N行数据 #

print("时间(s)\t摆角(°)")
for i in range(10):  # 打印前10行
    print(f"{t[i]:.4f}\t{np.rad2deg(theta[i]):.4f}")

三、查询特定时间对应的摆角 #

方法1：线性插值（精准） #

借助scipy的插值函数，可获取任意时间点的摆角：

from scipy.interpolate import interp1d

# 创建插值函数
theta_interp = interp1d(t, theta, kind='linear')

# 查询目标时间的摆角
target_time = 10.5
target_theta_deg = np.rad2deg(theta_interp(target_time))
print(f"时间{target_time:.2f}s对应的摆角：{target_theta_deg:.2f}°")

方法2：查找最近时间点（无需额外库） #

如果不想安装scipy，可以直接找最接近目标时间的计算点：

target_time = 10.5
# 找到与目标时间最接近的索引
closest_idx = np.argmin(np.abs(t - target_time))
closest_theta_deg = np.rad2deg(theta[closest_idx])
print(f"最接近{target_time:.2f}s的时间点{t[closest_idx]:.4f}s，摆角：{closest_theta_deg:.2f}°")

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