耦合摆微分方程求解代码报错：IndexError索引越界问题求助

问题背景 #

编程经验尚浅，课程中需实现耦合摆模型（二阶微分方程组），已将其简化为4个一阶微分方程，尝试用scipy.integrate.solve_ivp求解时遇到错误。

原代码 #

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

%matplotlib inline

g = 9.8
L = 2 #length of the pendulum in meters
k = 3
m1 = 2
m2 = 2

s = (20,20)
state = np.zeros(s)

def der_state(t, state):
        """compute the derivative of the given state"""
        der = np.zeros_like(state)
        der[0] = state[1]
        der[1] = (-(g/L)*np.sin(state[0]))-((k/m1)*((np.sin(state[0])-(np.sin(state[2])))))
        der[2] = state[3]
        der[3] = (-(g/L)*np.sin(state[0]))-((k/m2)*((np.sin(state[2])-(np.sin(state[1])))))
        return der

tf = 25 #simulation for tf seconds
n = 1000 #number of evaluation points
dt = tf/n
T = np.linspace(0.0, tf, n+1)

state0 = ([np.pi/4, 0.2]) #this is the initial state

sol = integrate.solve_ivp(der_state, (0, tf), state0, t_eval=T)
ang_pos = sol.y[0]

错误信息 #

3         der = np.zeros_like(state)
      4         der[0] = state[1]
----> 5         der[1] = (-(g/L)*np.sin(state[0]))-((k/m1)*((np.sin(state[0])-(np.sin(state[2])))))
      6         der[2] = state[3]
      7         der[3] = (-(g/L)*np.sin(state[0]))-((k/m2)*((np.sin(state[2])-(np.sin(state[1])))))

IndexError: index 2 is out of bounds for axis 0 with size 2

错误原因及修正方案 #

1. 初始状态维度不匹配
   你将二阶方程组拆分为4个一阶方程，意味着状态向量state需要包含4个元素：

   • state[0]: 第一个摆的角度θ₁
   • state[1]: 第一个摆的角速度dθ₁/dt
   • state[2]: 第二个摆的角度θ₂
   • state[3]: 第二个摆的角速度dθ₂/dt
     但你定义的state0 = ([np.pi/4, 0.2])只有2个元素，导致solve_ivp运行时state的长度为2，访问state[2]自然会触发越界错误。
     修正：给state0补充完整的初始值，比如：

   state0 = [np.pi/4, 0.2, 0.0, 0.0]  # θ₁初始值、dθ₁/dt初始值、θ₂初始值、dθ₂/dt初始值
   

2. 微分方程公式错误
   原代码中der[3]的重力项写错了，应该对应第二个摆的角度state[2]，而非第一个摆的state[0]；同时耦合项里错误地使用了角速度state[1]，应该用第一个摆的角度state[0]。正确的der_state函数应该是：

   def der_state(t, state):
       """compute the derivative of the given state"""
       der = np.zeros_like(state)
       # 第一个摆的角度导数=角速度
       der[0] = state[1]
       # 第一个摆的角速度导数（加速度）
       der[1] = -(g/L)*np.sin(state[0]) - (k/m1)*(np.sin(state[0]) - np.sin(state[2]))
       # 第二个摆的角度导数=角速度
       der[2] = state[3]
       # 第二个摆的角速度导数（加速度）
       der[3] = -(g/L)*np.sin(state[2]) - (k/m2)*(np.sin(state[2]) - np.sin(state[0]))
       return der
   

3. 多余代码清理
   代码中s = (20,20)和state = np.zeros(s)完全没有作用，solve_ivp会根据传入的state0自动确定状态向量的维度，建议删除这两行避免混淆。

修正后的完整代码 #

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

%matplotlib inline

g = 9.8
L = 2 #length of the pendulum in meters
k = 3
m1 = 2
m2 = 2

def der_state(t, state):
        """compute the derivative of the given state"""
        der = np.zeros_like(state)
        der[0] = state[1]
        der[1] = -(g/L)*np.sin(state[0]) - (k/m1)*(np.sin(state[0]) - np.sin(state[2]))
        der[2] = state[3]
        der[3] = -(g/L)*np.sin(state[2]) - (k/m2)*(np.sin(state[2]) - np.sin(state[0]))
        return der

tf = 25 #simulation for tf seconds
n = 1000 #number of evaluation points
T = np.linspace(0.0, tf, n+1)

state0 = [np.pi/4, 0.2, 0.0, 0.0] # 完整初始状态：θ1, dθ1/dt, θ2, dθ2/dt

sol = integrate.solve_ivp(der_state, (0, tf), state0, t_eval=T)
# 可以分别提取两个摆的角度和角速度
ang_pos1 = sol.y[0]
ang_pos2 = sol.y[2]
ang_vel1 = sol.y[1]
ang_vel2 = sol.y[3]

# 简单绘制角度变化曲线验证
plt.figure(figsize=(10,6))
plt.plot(T, ang_pos1, label='摆1角度')
plt.plot(T, ang_pos2, label='摆2角度')
plt.xlabel('时间(s)')
plt.ylabel('角度(rad)')
plt.legend()
plt.show()

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