我正在尝试求解微分方程：

m（t）= M（ x ） x'' + C（ x ， x' ）和B x'

其中 x 和 x' 是具有2个条目的矢量，分别表示动力学中的角度和角速度系统。 M（ x ）是2x2矩阵，它是theta分量的函数，C是2x1向量，它是theta和theta'的函数，而B是a 2x2常数矩阵。 m（t）是一个2 * 1001数组，其中包含在1001个时间步长处施加到两个关节中的每个关节的扭矩，我想根据这些1001个时间步长来计算角度的变化。

我已经将其转换为标准格式，例如：

x'' = M（ x ）^-1（m（t）-C（ x ， x' ）-B x' ）

然后替换 y_1 = x 和 y_2 = x' 给出方程的一阶线性系统：

y_2 = y_1'

y_2' = M（ y_1 ）^-1（m（t）-C（ y_1 ， y_2 ）-B y_2 ）

（我在代码中为x和y使用了theta和phi）

def joint_angles(theta_array,t,torques,B):

    phi_1 = np.array([theta_array[0],theta_array[1]])
    phi_2 = np.array([theta_array[2],theta_array[3]])

    def M_func(phi):
        M = np.array([[a_1+2.*a_2*np.cos(phi[1]),a_3+a_2*np.cos(phi[1])],[a_3+a_2*np.cos(phi[1]),a_3]])
        return np.linalg.inv(M)

    def C_func(phi,phi_dot):
        return a_2 * np.sin(phi[1]) * np.array([-phi_dot[1] * (2. * phi_dot[0] + phi_dot[1]),phi_dot[0]**2])

    dphi_2dt = M_func(phi_1) @ (torques[:,t] - C_func(phi_1,phi_2) - B @ phi_2)

    return dphi_2dt,phi_2

t = np.linspace(0,1,1001)
initial = theta_init[0],theta_init[1],dtheta_init[0],dtheta_init[1]

x = odeint(joint_angles,initial,args = (torque_array,B))

我得到了一个错误，我无法使用t数组来索引扭矩，这是很有意义的，但是我不确定如何在每个时间步中使用扭矩的当前值。

我还尝试将odeint命令放在for循环中，并且一次只能一次评估它，使用函数的解作为下一个循环的初始条件，但是该函数只是返回了初始条件，这意味着每个循环都是相同的。这使我怀疑我在执行标准表单时犯了一个错误，但是我无法弄清楚它是什么。但是，最好不必每次都在for循环中调用odeint求解器，而将它们作为一个整体来完成。

如果有帮助，我的初始条件和常数值是：

    theta_init = np.array([10*np.pi/180,143.54*np.pi/180])
    dtheta_init = np.array([0,0])

    L_1 = 0.3
    L_2 = 0.33
    I_1 = 0.025
    I_2 = 0.045
    M_1 = 1.4
    M_2 = 1.0
    D_2 = 0.16

    a_1 = I_1+I_2+M_2*(L_1**2)
    a_2 = M_2*L_1*D_2
    a_3 = I_2

感谢您的帮助！