Решение векторного дифференциального уравнения второго порядка при индексации в массив - PullRequest
Новая
       12

Решение векторного дифференциального уравнения второго порядка при индексации в массив

0 голосов
/ 05 ноября 2019

Я пытаюсь решить дифференциальное уравнение:

м (т) = M ( x ) x '' + C ( x , x ') + B x'

где x и x ' - векторы с 2 записями, представляющими углы и угловую скорость вдинамическая система. M ( x ) - матрица 2x2, которая является функцией компонентов тета, C - вектор 2x1, который является функцией тета и тета ', и B - матрица 2x2. постоянных. m (t) - это массив 2 * 1001, содержащий крутящие моменты, приложенные к каждому из двух соединений с шагом 1001 времени, и я хотел бы рассчитать эволюцию углов в зависимости от этих шагов 1001.

Я преобразовал его в стандартную форму так, что:

x '' = M ( x )^ -1 (м (т) - C ( x , x ') - B x')

Затем подставим y_1 = x и y_2 = x ' дает линейную систему уравнений первого порядка:

y_2 = y_1 '

y_2' = M ( y_1 ) ^ - 1 (м (т) - C ( y_1 , y_2 ) - B y_2 )

(я использовал тета и фи в своем коде для x и 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, что имеет смысл, однако я не уверен, как заставить его использовать текущее значение крутящих моментов на каждом временном шаге.

Я также пытался поместить команду odeint в цикл for и оценивать ее только по одному временному шагу за раз, используя решение функции в качестве начальных условий для следующего цикла, однако функция просто возвращала исходнуюусловия, означающие, что каждый цикл был идентичен. Это заставляет меня подозревать, что я допустил ошибку в своей реализации стандартной формы, но я не могу понять, что это такое. Однако было бы предпочтительнее не каждый раз вызывать решатель odeint в цикле for, а делать все это как единое целое.

Если полезно, мои начальные условия и постоянные значения: 

    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

Спасибо за помощь!

1 Ответ

1 голос
/ 06 ноября 2019

Решатель использует внутренний степпинг, адаптированный к проблемам. Данный временной список представляет собой список точек, в которых внутреннее решение интерполируется для выходных выборок. Внутренние и внешние временные списки никоим образом не связаны, внутренний список зависит только от заданных допусков.

Нет фактической естественной связи между индексами массива и временем выборки.

Преобразование заданного времени в индекс и построение выборочного значения из записей окружающей таблицы называется интерполяцией (кусочно-полиномиальной функцией).

Крутящий момент как физическое явление, по крайней мере, непрерывный, кусочно-линейная интерполяция - это самый простой способ преобразовать данную таблицу значений функции в фактическую непрерывную функцию. Конечно, также нужен массив времени.

Поэтому используйте numpy.interp1d или более продвинутые подпрограммы scipy.interpolate, чтобы определить функцию крутящего момента, которая может быть оценена в произвольные моменты времени, как того требует решатель и его метод интегрирования. .

...