Использование scipy odeint для уравнений с переменной фазой

Question

В основном ... Мне нужен способ включить фазовый сдвиг в мои дифференциальные уравнения. То есть, я имею в определении моей системной функции, которая возвращает dY/dt что-то вроде Y (t-3). Как это дифференциальное уравнение:

dY/dt = a*Y(t) + b*Y(t-tau) 

Теперь, если я пытаюсь писать это как функции определения системы для перехода к scipy.odeint, я потерял:

def eqtnSystem(A,t): 
    Y = A 
    a = 1 
    b = 5 
    tau = 3 
    return a*Y + b*???  # how do I Y(t-tau) ? 

Это в основном это. Я действительно надеюсь, что есть простой ответ, но я не мог отследить его.

Конкретно ... Я пытаюсь численно вычислить решение для системы, определяемой следующей функцией:

def etaFunc(A,t): 
    #...definition of all those constants is here... 
    return array([(gamma[0,0]*xi(t-theta[0])[0] - eta[0] + zeta[0])/tau[0],\ 
      (gamma[1,1]*xi(t-theta[1])[1] - eta[1] + zeta[1])/tau[1],\ 
      (gamma[2,2]*xi(t-theta[2])[2] - eta[2] + zeta[2])/tau[2],\ 
      ( beta[3,0]*pastEta(t-theta[3])[0] \ 
      + beta[3,1]*pastEta(t-theta[4])[1] \ 
      + beta[3,2]*pastEta(t-theta[5])[2] -eta[3]+ zeta[3])/tau[3],\ 
      ( beta[4,3]*pastEta(t-theta[6])[3] \ 
      + beta[4,2]*pastEta(t-theta[7])[2] - eta[4] + zeta[4])/tau[4]]) 

Эта функция затем позже в данной odeint так:

ETA = integrate.odeint(etaFunc,initCond,time) 

, а затем я могу получить каждый отдельный компонент ETA (например, eta_0) следующим образом: ETA[:,0].

Проблема, которая возникает у меня здесь, связана с pastEta(t-theta[?]). На данный момент это функция, которая пытается найти уже рассчитанные значения eta (когда start_time < t-theta[?] < t и theta[?] > 0.Это не очень хорошо работает.

Я вижу, что в этом случае я мог бы найти каждый компонент eta индивидуально и затем вычислить вычисленные значения для ранее вычисленных этанных компонентов (eta_0, eta_1, eta_2) для вычисления eta_3 и аналогично для eta_4, но это не идеально, так как оно отнимает у меня способность «подключать и играть» любые общие формулы.

Answer 1

Для этого существует ряд существующих библиотек и примеров.

http://www.google.fi/search?q=python+delay+differential+equation дает мне:

• http://zulko.wordpress.com/2013/03/01/delay-differential-equations-easy-with-python/ (использует scipy.integrate.ode)
• https://pypi.python.org/pypi/PyDDE
• http://pydelay.sourceforge.net/

Answer 2

Один из способов сделать это с помощью integrate.odeint() - это запустить integrate.odeint() за много коротких промежутков времени между начальным временем и временем окончания, сохраняя значение времени и выходное значение Y после каждого короткого интервала l в списках. Это позволит вам интерполировать значение Y в списках, используя, например, scipy.interpolate.interp1d(), каждый раз, когда вам нужно Y(t-3).

Вы получите только приблизительное значение для Y(t-3), если вы сделаете это, конечно, но если временные интервалы будут достаточно близко друг к другу, этот подход может быть для вас удовлетворительным. В конце концов, значения Y(t), рассчитанные численными решателями ОДУ, также являются приблизительными.

Answer 3

Задержки не совсем линейные функции. Обычная задержка шага представлена ​​в области Лапласа как e**(a*s)/s, где a - это задержка.

Это означает, что «нормальные» решатели ODE не будут работать, если у вас нет обходного пути. Обычно это обходное решение не очень легко сделать, так как для жестких задач вы обычно не можете интерполировать с достаточным приближением.

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

Другое решение делает это симметрично (если можно, вы можете попробовать SymPy).

Третье решение - это хранение прошлых результатов и интерполяция, чтобы найти точное прошлое, которое вам нужно (может быть, недостаточно).

Четвертое решение может быть рекомендовано некоторыми документами для симуляторов: используйте c2d() и моделируйте всю модель в дискретном времени и сохраняйте прошлые переменные в списке/массиве (без интерполяции, но вам может потребоваться использование небольших шагов для улучшения точность).

Пятое решение использует Padé approximation для представления задержки вашей модели (может работать в зависимости от вашего случая). В python-контроле есть pade() function, чтобы приблизиться именно к этому.