在索引到数组时求解向量二阶微分方程 [英] Solving vector second order differential equation while indexing into an array

问题描述

我正在尝试解微分方程:

I'm attempting to solve the differential equation:

m(t) = M(x)x'' + C(x, x') + B x'

m(t) = M(x)x'' + C(x, x') + B x'

其中 x 和 x' 是向量，其中 2 个条目表示动力学中的角度和角速度系统.M(x) 是一个 2x2 矩阵，它是 theta 分量的函数，C 是一个 2x1 向量，它是 theta 和 theta' 的函数，B 是一个2x2 常量矩阵.m(t) 是一个 2*1001 数组，包含在 1001 个时间步长处施加到两个关节中的每一个的扭矩，我想计算角度的演变，作为这 1001 个时间步长的函数.

where x and x' are vectors with 2 entries representing the angles and angular velocity in a dynamical system. M(x) is a 2x2 matrix that is a function of the components of theta, C is a 2x1 vector that is a function of theta and theta' and B is a 2x2 matrix of constants. m(t) is a 2*1001 array containing the torques applied to each of the two joints at the 1001 time steps and I would like to calculate the evolution of the angles as a function of those 1001 time steps.

我已将其转换为标准形式:

I've transformed it to standard form such that :

x'' = M(x)^-1 (m(t) - C(x、x') - B x'>)

x'' = M(x)^-1 (m(t) - C(x, x') - B x')

然后替换 y_1 = x 和 y_2 = x' 给出一阶线性方程组:

Then substituting y_1 = x and y_2 = x' gives the first order linear system of equations:

y_2 = y_1'

y_2' = M(y_1)^-1 (m(t) - C(y_1、y_2) - B y_2)

(我在代码中为 x 和 y 使用了 theta 和 phi)

(I've used theta and phi in my code for x and y)

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 数组索引扭矩，这很有意义，但是我不确定如何让它在每个时间步使用扭矩的当前值.

I get the error that I cannot index into torques using the t array, which makes perfect sense, however I am not sure how to have it use the current value of the torques at each time step.

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

I also tried putting odeint command in a for loop and only evaluating it at one time step at a time, using the solution of the function as the initial conditions for the next loop, however the function simply returned the initial conditions, meaning every loop was identical. This leads me to suspect I've made a mistake in my implementation of the standard form but I can't work out what it is. It would be preferable however to not have to call the odeint solver in a for loop every time, and rather do it all as one.

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

If helpful, my initial conditions and constant values are:

    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

感谢您的帮助！

推荐答案

求解器使用适合问题的内部步进.给定的时间列表是一个点列表，在这些点中，内部解决方案为输出样本进行了插值.内部和外部时间列表没有任何关系，内部列表只取决于给定的容差.

The solver uses an internal stepping that is problem adapted. The given time list is a list of points where the internal solution gets interpolated for output samples. The internal and external time lists are in no way related, the internal list only depends on the given tolerances.

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

There is no actual natural relation between array indices and sample times.

将给定时间转换为索引并从周围表条目构建样本值称为插值(通过分段多项式函数).

The translation of a given time into an index and construction of a sample value from the surrounding table entries is called interpolation (by a piecewise polynomial function).

<小时>

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

---

Torque as a physical phenomenon is at least continuous, a piecewise linear interpolation is the easiest way to transform the given function value table into an actual continuous function. Of course one also needs the time array.

因此使用 numpy.interp1d 或更高级的 scipy.interpolate 例程来定义扭矩函数，该函数可以根据求解器及其要求在任意时间进行评估集成方法.

So use numpy.interp1d or the more advanced routines of scipy.interpolate to define the torque function that can be evaluated at arbitrary times as demanded by the solver and its integration method.

这篇关于在索引到数组时求解向量二阶微分方程的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！