求解向量二阶微分方程,同时索引到数组中

问题描述 投票:0回答:1

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

m(t)= M(x x'' + C(xx'] >)+ B x'] >>

其中x []]]和x'

是具有2个条目的矢量,代表动态系统中的角度和角速度。 M([[[x)是一个2x2矩阵,它是theta分量的函数,C是一个2x1向量,它是theta和theta'的函数,而B是一个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 ['y_2'

= M(y [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, t, 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

感谢您的帮助!

我正在尝试求解微分方程:m(t)= M(x)x''+ C(x,x')+ B x'其中x和x'是具有2个项的矢量,分别表示角度和动力系统中的角速度。 M(x)... 

求解器使用适合问题的内部步进。给定的时间列表是内部解决方案针对输出样本进行插值的点的列表。内部和外部时间列表没有任何关系,内部列表仅取决于给定的公差。

数组索引和采样时间之间没有实际的自然关系。


将给定时间转换为索引并从周围的表条目中构建样本值称为

interpolation(通过分段多项式函数)。

作为物理现象的扭矩至少是连续的,分段线性插值是将给定函数值表转换为实际连续函数的最简单方法。当然也需要时间数组。


因此,请使用numpy.interp1d或更高级的例程scipy.interpolate定义可根据求解器及其积分方法的要求在任意时间评估的转矩函数。
python-3.x ode
1个回答
0
投票

数组索引和采样时间之间没有实际的自然关系。


将给定时间转换为索引并从周围的表条目中构建样本值称为

interpolation(通过分段多项式函数)。

作为物理现象的扭矩至少是连续的,分段线性插值是将给定函数值表转换为实际连续函数的最简单方法。当然也需要时间数组。


因此,请使用numpy.interp1d或更高级的例程scipy.interpolate定义可根据求解器及其积分方法的要求在任意时间评估的转矩函数。
© www.soinside.com 2019 - 2024. All rights reserved.