如何用Scipy的odeint求解二阶微分方程组？代码结果不符求助

老哥，当然可以用odeint求解二阶微分方程组啦！这正是它的拿手好戏之一，不过这里有个关键前提——你得先把二阶微分方程转化为一阶微分方程组，这大概率就是你代码出问题的核心原因。

核心思路：二阶转一阶 #

对于任何二阶微分方程 $\ddot{\theta} = f(\theta, \dot{\theta}, t)$，我们可以引入新变量把它拆成两个一阶方程：

• 令 $y_1 = \theta$（原变量）
• 令 $y_2 = \dot{\theta}$（原变量的一阶导数）
  这样原方程就转化为：
  $\dot{y_1} = y_2$
  $\dot{y_2} = f(y_1, y_2, t)$

如果是多个二阶方程，就以此类推，把每个变量的导数都作为新的状态元素加入状态向量——比如两个二阶方程的话，状态向量应该有4个元素（每个变量+它的导数）。

你的代码疑点分析 #

从你给出的片段来看，有几个明显的可能出错的地方：

• 你的初始状态initial_state = [0.1, 0, 0]只有3个元素，但如果是二阶方程组，状态向量长度应该对应「变量数 + 导数数」，比如一个二阶方程需要2个元素，两个二阶方程需要4个，这里的长度不太对。
• 你在my_system里把current_state解包成theta1, theta2, theta3，但没看到对应的导数方程定义——这部分是核心，一旦方程写错，结果肯定和预期不符。
• 你定义了参数a=4，但没看到你把它传入odeint的过程，如果直接在函数里用全局变量，容易出现作用域问题。

示例：用odeint解双摆（两个二阶方程） #

给你一个完整的可运行示例，帮你理解正确的写法：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def double_pendulum(state, t, g, L1, L2, m1, m2):
    # 解包状态变量：theta1, theta1的导数, theta2, theta2的导数
    theta1, dtheta1_dt, theta2, dtheta2_dt = state
    
    # 计算二阶导数（方程组的右端项）
    c = np.cos(theta1 - theta2)
    s = np.sin(theta1 - theta2)
    denominator = m1 + m2 * s**2
    
    ddtheta1_dt = (m2 * g * np.sin(theta2) * c - m2 * s * (L1 * dtheta1_dt**2 * c + L2 * dtheta2_dt**2) - (m1 + m2) * g * np.sin(theta1)) / (L1 * denominator)
    ddtheta2_dt = ((m1 + m2) * (L1 * dtheta1_dt**2 * s - g * np.sin(theta2) + g * np.sin(theta1) * c) + m2 * L2 * dtheta2_dt**2 * s * c) / (L2 * denominator)
    
    # 返回每个状态元素的一阶导数，顺序要和state一致
    return [dtheta1_dt, ddtheta1_dt, dtheta2_dt, ddtheta2_dt]

# 物理参数
g = 9.81
L1 = 1.0
L2 = 1.0
m1 = 1.0
m2 = 1.0

# 初始状态：theta1初始角, theta1初始角速度, theta2初始角, theta2初始角速度
initial_state = [np.pi/4, 0, np.pi/2, 0]
time_points = np.linspace(0, 20, 1000)

# 求解：注意用args传递额外参数
solution = odeint(double_pendulum, initial_state, time_points, args=(g, L1, L2, m1, m2))

# 提取结果
theta1 = solution[:, 0]
theta2 = solution[:, 2]

# 绘图展示
plt.figure(figsize=(10, 6))
plt.plot(time_points, theta1, label='摆角1')
plt.plot(time_points, theta2, label='摆角2')
plt.xlabel('时间 (s)')
plt.ylabel('角度 (rad)')
plt.legend()
plt.show()

排查你的代码的建议 #

• 先确认所有二阶方程都完成了「一阶化」，状态向量的长度符合要求。
• 检查my_system的返回值：必须是一个列表，每个元素对应状态向量中对应元素的一阶导数，顺序完全一致。
• 用odeint的args参数传递额外参数（比如你的a=4），避免用全局变量。
• 核对初始状态：确保包含了所有变量的初始值和初始导数，顺序和函数解包的顺序一致。

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