我是一个编程新手，必须在我的一个课程中编写一些模型。在这种情况下，我必须对耦合摆建模，这是一个二阶微分方程系统。

正如在课堂上看到的，我将方程简化为4个一阶微分方程的系统，并尝试在我的代码中实现它们并使用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

我似乎不明白是什么导致了这个问题。