您的状态有 4 个组件，因此您在每个阶段都需要 4 个斜坡.很明显，z 的斜率/更新不能来自 k1_zdot，它必须是 k1_z 之前要计算的

Your state has 4 components, thus you need 4 slopes in each stage. It should be obvious that the slope/update for z can not come from k1_zdot, it has to be k1_z which is to be computed previously as

k1_z = zdot

进入下一阶段

k2_z = zdot + dt/2*k1_zdot

等

但***使用矢量化界面

def derivs(t,u):
    z, theta, dz, dtheta = u
    ddz = -omega_z * z - epsilon / 2 / m * theta
    ddtheta = -omega_theta * theta - epsilon / 2 / I * z 
    return np.array([dz, dtheta, ddz, ddtheta]);

然后使用 RK4 的标准公式

and then use the standard formulas for RK4

i = 0
while True:
    # runge_kutta
    k1 = derivs(t[i], u[i])
    k2 = derivs(t[i] + dt/2, u[i] + dt/2 * k1)
    k3 = derivs(t[i] + dt/2, u[i] + dt/2 * k2)
    k4 = derivs(t[i] + dt, u[i] + dt * k3)

    u.append (u[i] + (k1 + 2*k2 + 2*k3 + k4) * dt / 6)
    i += 1

然后解压为

z, theta, dz, dtheta = np.asarray(u).T