Об обыкновенных дифференциальных уравнений (ОДУ) и оптимизации, в Python

Question

Я хочу, чтобы решить такую ​​проблему:

dy/dt = 0.01*y*(1-y), find t when y = 0.8 (0<t<3000) 

Я попытался функцию ода ​​в Python, но он может рассчитывать только тогда, когда у Т данный.

Итак, есть ли какие-либо простые способы решения этой проблемы в Python?

---

PS: Эта функция является просто примером. Моя реальная проблема настолько сложна, что ее невозможно решить аналитически. Поэтому я хочу знать, как его решить численно. И я думаю, что эта проблема больше похожа на проблему оптимизации:

Objective function y(t) = 0.8, Subject to dy/dt = 0.01*y*(1-y), and 0<t<3000 

---

PPS: Моя реальная проблема заключается в:

objective function: F(t) = 0.85, 
subject to: F(t) = sqrt(x(t)^2+y(t)^2+z(t)^2), 
      x''(t) = (1/F(t)-1)*250*x(t), 
      y''(t) = (1/F(t)-1)*250*y(t), 
      z''(t) = (1/F(t)-1)*250*z(t)-10, 
      x(0) = 0, y(0) = 0, z(0) = 0.7, 
      x'(0) = 0.1, y'(0) = 1.5, z'(0) = 0, 
      0<t<5 

Answer 1

Использование scipy.integrate.odeint для обработки вашей интеграции, а также анализировать результаты позже.

import numpy as np 
from scipy.integrate import odeint 

ts = np.arange(0,3000,1) # time series - start, stop, step 

def rhs(y,t): 
    return 0.01*y*(1-y) 

y0 = np.array([1]) # initial value 

ys = odeint(rhs,y0,ts) 

Затем анализируют Numpy массив ys, чтобы найти ответ на свой вопрос (размеры массива ts спичек ys). (Это может не работать в первый раз, потому что я строю из памяти).

Это может быть связано с использованием функции scipy interpolate для массива ys, так что вы получите результат в t.

EDIT: Я вижу, что вы хотите решить весну в 3D. Это должно быть хорошо с вышеуказанным методом; Odeint на сайте scipy есть примеры для таких систем, как связанные пружины, которые можно решить, и они могут быть расширены.

Answer 2

Это дифференциальное уравнение может быть решено аналитический довольно легко:

д/дт = 0,01 * у * (1-у)

Перестановки собирать у и т условия на противоположных сторонах

100 dt = 1/(y * (1-y)) dy

Lhs интегрируется тривиально в 100 * t, rhs немного сложнее. Мы всегда можем написать произведение двух дробей в виде суммы двух дробей * некоторые константы:

1/(у * (1-у)) = А/г + В/(1-у)

Значения для A и B могут быть разработаны путем помещения rhs на тот же знаменатель и сравнения констант и первого порядка y с обеих сторон. В этом случае это просто, A = B = 1. Таким образом, мы должны интегрировать

1/у + 1/(1-у) ау

Первый член интегрируется с Ln (у), второй член может быть интегрирован с заменой переменных и = 1- y до -ln (1-y).Наше интегральное уравнение для этого выглядит следующим образом:

100 * T + C = Ln (у) - п (1-у)

не забывая при этом постоянную интегрирования (удобно записать его на LHS здесь) , Мы можем объединить два логарифмов термины:

100 * T + C = п (г/(1-у))

Для того, чтобы решить т для точного значения у, в первую очередь необходимо разработать значение C. Мы делаем это, используя начальные условия. Ясно, что если y начинается с 1, dy/dt = 0, а значение y никогда не изменяется. Таким образом подключить значения для у и т в начале

100 * 0 + С = п (у (0)/(1 - у (0))

Это даст значение для C (при условии, y не 0 или 1), а затем используйте y = 0,8 для получения значения для t. Обратите внимание, что из-за логарифма и коэффициента 100 умножение ty достигнет 0,8 в относительно коротком диапазоне значений t, если только начальное значение y . невероятно мал это, конечно, также просто переставить уравнение выше, чтобы выразить у через т, то вы можете построить график функции, а также

Edit:. Численное интегрирование

Для более сложной ODE, которая не может быть решена аналитически, вам придется попробовать численно. Первоначально мы знаем только значение функции в нулевое время y (0) (мы должны знать, по крайней мере, то, чтобы однозначно определить траекторию функции) и как оценить градиент. Идея численного интегрирования состоит в том, что мы можем использовать наши знания градиента (который говорит нам, как изменяется функция), чтобы определить, какое значение функции будет находиться в непосредственной близости от нашей начальной точки. Самый простой способ сделать это состоит в Эйлера интеграции:

у (дт) = у (0) + Dy/дт * дт

Эйлера интеграция предполагает, что градиент является постоянным между Т = 0 и Т = DT. Как только y (dt) известно, градиент также может быть рассчитан и в свою очередь используется для вычисления y (2 * dt) и т. Д., Постепенно наращивая полную траекторию функции. Если вы ищете конкретное целевое значение, просто подождите, пока траектория не пройдет мимо этого значения, затем интерполируйте между двумя последними позициями, чтобы получить точный t.

Проблема с интеграцией Эйлера (и со всеми другими методами численного интегрирования) заключается в том, что ее результаты являются точными только тогда, когда его допущения действительны. Поскольку градиент не является постоянным между парами временных точек, для каждого шага интеграции возникает определенное количество ошибок, которое со временем будет нарастать, пока ответ не будет полностью неточным. Чтобы улучшить качество интеграции, необходимо использовать более сложные приближения к градиенту. Ознакомьтесь, например, с методами Runge-Kutta, которые представляют собой семейство интеграторов, которые удаляют прогрессивные порядки ошибок за счет увеличения времени вычислений. Если ваша функция дифференцируема, знание второй или даже третьей производных может также использоваться для уменьшения ошибки интеграции.

К счастью, кто-то еще проделал здесь тяжелую работу, и вам не нужно слишком беспокоиться о решении таких проблем, как численная стабильность или глубокое понимание всех деталей (хотя понимание примерно того, что происходит на помогает много). Обратитесь к http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html#scipy.integrate.ode за примером класса интегратора, который вы сможете использовать сразу. Например,

from scipy.integrate import ode 
def deriv(t, y): 
    return 0.01 * y * (1 - y) 
my_integrator = ode(deriv) 
my_integrator.set_initial_value(0.5) 
t = 0.1 # start with a small value of time 
while t < 3000: 
    y = my_integrator.integrate(t) 
    if y > 0.8: 
     print "y(%f) = %f" % (t, y) 
     break 
    t += 0.1 

Этот код напечатает первое значение t, когда y пройдет 0.8 (или ничего, если оно никогда не достигает 0,8). Если вы хотите получить более точное значение t, сохраните y предыдущего t и интерполируйте между ними.

Answer 3

Что вы просите, является интегратором ODE с возможностями поиска root. Они существуют и низкоуровневый код для таких интеграторов снабжен scipy, но они еще не были обернуты в привязки python.

Для получения дополнительной информации см этого списка рассылки пост, который обеспечивает несколько вариантов: http://mail.scipy.org/pipermail/scipy-user/2010-March/024890.html

Вы можете использовать следующий пример реализацию, который использует возвраты (следовательно, он не является оптимальным, так как это на болты дополнения к интегратору что не имеет корней вывод сам по себе): https://github.com/scipy/scipy/pull/4904/files

Answer 4

в дополнение к Krastanov`s ответ:

в стороне PyDSTool есть и другие пакеты, такие как Pysundials и Assimulo, которые обеспечивают привязку к й e решатель IDA от Sundials. Этот решатель имеет возможности поиска корней.