Python中RK4法求解复杂系统运动微分方程的代码修复求助

核心问题分析 #

1. 循环未执行：main函数中先用np.linspace(t0, T, n)生成了完整的时间序列，导致t_n[-1]直接等于T，后续的while循环根本不会运行，所以y_n只保留了初始条件。
2. 微分方程计算顺序错误：equations函数中计算f4时用到了还未定义的f3和f5，导致依赖关系混乱，计算结果异常。
3. 微分方程表达式不完整：f5的分子最后一项缺少对应的导数项，会导致方程描述的物理关系不成立。

修复后的完整代码 #

import numpy as np

# 系统参数
p = 7.850
A = (0.01) ** 2
m1 = 1
m3 = 1
g = 9.81
r2 = 1
r3 = 1
r1 = 1
rcm = r1
r = 4 * r2
m2 = 1


def equations(t, y):
    theta, phi1, phi2, x1, w1, p1 = y

    # 定义各阶导数的对应关系
    f0 = x1  # d(theta)/dt = x1
    f1 = w1  # d(phi1)/dt = w1
    f2 = p1  # d(phi2)/dt = p1

    # 先计算f3（d(x1)/dt），消除依赖顺序问题
    f3 = ((m1 * r2 * rcm * w1**2 * np.sin(theta + phi1)) + 
          (g * m1 * rcm * np.sin(theta)) - 
          (m1 * r2 * rcm * np.cos(theta + phi1) * w1)) / (m1 * r2**2)
    # 计算f4（d(w1)/dt），此时f3已定义
    f4 = ((m1 * r2 * rcm * x1**2 * np.sin(theta + phi1)) - 
          (m3 * r3 * (r - r2) * p1**2 * np.sin(phi2 - phi1)) - 
          (g * m1 * r2 * np.sin(phi1)) + 
          (g * m2 * ((r/2) - r2) * np.sin(phi1)) + 
          (g * m3 * (r - r2) * np.sin(phi1)) - 
          (m1 * r2 * rcm * np.cos(theta + phi1) * f3) + 
          (m3 * r3 * (r - r2) * np.cos(phi2 - phi1) * p1)) / (m1 * r2**2)
    # 修复f5的表达式，补充缺失的导数项（需根据实际方程调整）
    f5 = ((m3 * r3 * (r - r2) * np.sin(phi2 - phi1) * w1**2) + 
          (g * m3 * r3 * np.sin(phi2)) + 
          (m3 * r3 * (r - r2) * np.cos(phi2 - phi1) * w1)) / (m3 * r3**2)

    return np.array([f0, f1, f2, f3, f4, f5])


def RK4_step(t, y, h, f):
    # 标准RK4四阶龙格库塔计算
    k1 = f(t, y)
    k2 = f(t + h/2, y + h/2 * k1)
    k3 = f(t + h/2, y + h/2 * k2)
    k4 = f(t + h, y + h * k3)
    return y + h/6 * (k1 + 2*k2 + 2*k3 + k4)


def main():
    t0 = 0
    T = 100
    n = 10000
    # 初始条件
    y0 = np.array([np.pi - np.arccos((2 - 1.5)/r2), 
                   np.arccos((2 - 1.5)/r2), 
                   np.arcsin(1.5/(r - r2)), 
                   0, 0, 0])
    
    h = (T - t0) / n
    t_n = np.zeros(n)
    y_n = np.zeros((n, 6))
    t_n[0] = t0
    y_n[0] = y0

    # 逐步迭代计算RK4
    for i in range(1, n):
        t_n[i] = t_n[i-1] + h
        y_n[i] = RK4_step(t_n[i-1], y_n[i-1], h, equations)

    # 输出最后几个结果验证
    print("最后5个时间步的系统状态：")
    print(y_n[-5:])


if __name__ == "__main__":
    main()

关键修复点说明 #

• 循环逻辑修正：改用for循环直接迭代n次，避免np.linspace导致的循环失效问题，同时用数组存储结果更高效。
• 微分方程计算顺序调整：先计算f3，再计算f4，确保变量依赖关系正确，不会出现未定义变量的问题。
• RK4函数优化：将RK4函数改为直接返回更新后的状态，简化主循环的计算逻辑，去掉不必要的y解包操作。
• f5表达式修复：补充了缺失的导数项（这里假设是乘以w1，如果实际方程不同，需要根据你的微分方程附图调整这部分）。

系统模拟方法说明 #

1. 微分方程建模：确保所有导数项的依赖关系正确，高阶导数的计算必须基于已定义的低阶导数或已知变量。
2. 数值方法选择：RK4是四阶精度的数值积分方法，适合求解非线性微分方程组，每一步需要计算四个斜率值并加权平均。
3. 结果验证：可以先取较小的时间范围（比如T=10）和较少的步数，观察结果是否符合物理直觉，也可以对比解析解（如果存在）或简化模型的结果。
4. 参数调整：如果结果出现发散，可能需要减小步长h，或者检查微分方程的表达式是否正确。

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