Система дифференциальных уравнений с постоянными во времени в массивах с использованием odeint - PullRequest
Новая
       39

Система дифференциальных уравнений с постоянными во времени в массивах с использованием odeint

0 голосов
/ 21 февраля 2019

Допустим, у меня есть система дифференциальных уравнений, и я хочу решить ее с помощью odeint.Некоторые из системных констант зависят от времени, и их значения хранятся в массивах (a, b, c и d с формой (8000,)).Я хочу, чтобы система использовала разные значения этих констант на каждом временном шаге.См. Пример упрощенного кода: 

t = np.linspace(0,100,8000)

a = np.array([0, 5, 6, 12, 1.254, ..., 0.145])     # shape (8000, )
b = np.array([1.45, 5.9125, 1.367, ..., 3.1458])
c = np.array([0.124, 0.258, 0.369, ..., 0.147])
d = np.array([7.145, 5.123, 6.321, ..., 0.125])

def system(k,t):
    vcx_i = k[0]
    vcy_i = k[1]
    psi_i = k[2]
    wz_i = k[3]

    vcx_i_dot = a*np.cos(psi_i)-b*np.sin(psi_i)
    vcy_i_dot = b*np.cos(psi_i)+a*np.sin(psi_i)
    psi_i_dot = wz_i
    wz_i_dot = c*vcx_i-a*np.sin(psi_i)-d*np.sin(psi_i)-b*np.cos(psi_i)

    return [vcx_i_dot, vcy_i_dot, psi_i_dot wz_i_dot]

k0 = [0.1257, 0, 0, 0]

k = odeint(system, k0, t)

vcx_i = k[:,0]
vcy_i = k[:,1]
psi_i = k[:,2]
wz_i = k[:,3]

psi_i = [system(t_i, k_i)[2] for k_i, t_i in zip(k,t)]
wz_i = [system(t_i, k_i)[3] for k_i, t_i in zip(k,t)]

Наиболее подходящее решение, которое я нашел до сих пор, это: Решение системы од (с изменением константы!) С использованием scipy.integrate.odeint? Но так как у меня есть только значения моих переменных в массивах, а не уравнения переменных, которые зависят от времени (например, a = f (t)), я попытался применить интерполяцию между значениями в моих массивах, как показано здесь: ODEINT с несколькими параметрами (зависит от времени) Мне удалось заставить код работать без ошибок, но общее время резко возросло, и результаты решенной системы совершенно неверны.Я попробовал любой возможный тип интерполяции, который я нашел здесь: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html, но все еще неверный результат.Это означает, что моя интерполяция не самая лучшая, или мои точки в массивах (8000 значений) слишком велики, чтобы интерполировать между ними и правильно решать систему.Как правильно решить такую ​​систему?Я новичок в Python и использую Python 2.7.12 на Ubuntu 16.04 LTS.Заранее спасибо.

1 Ответ

0 голосов
/ 21 февраля 2019

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

В последнее время ode и odeint были заменены другими, более гибкими функциями (см. здесь )

Iбудет начинаться с и явного метода вместо неявного (по умолчанию для solve_ivp это runge kutta, а по умолчанию для odeint - LSODA): 

interp = scipy.interpolate.interp1d(t,(a,b,c,d))

def system(t,k):
    vcx_i = k[0]
    vcy_i = k[1]
    psi_i = k[2]
    wz_i = k[3]
    a, b, c, d = interp(t)

    vcx_i_dot = a*np.cos(psi_i)-b*np.sin(psi_i)
    vcy_i_dot = b*np.cos(psi_i)+a*np.sin(psi_i)
    psi_i_dot = wz_i
    wz_i_dot = c*vcx_i-a*np.sin(psi_i)-d*np.sin(psi_i)-b*np.cos(psi_i)
    return [vcx_i_dot, vcy_i_dot, psi_i_dot, wz_i_dot]

k0 = [0.1257, 0, 0, 0]
steps = 1
method = 'RK23'
atol = 1e-3
s = solve_ivp(dydt, (0, 100), k0, method=method, t_eval=t, atol=atol, vectorized=True)

вы можете увеличить / уменьшить atol или изменитьметод.

...