Интеграция жестких ODE с Python

Question

Я ищу хорошую библиотеку, которая будет интегрировать жесткие ODE в Python.Проблема в том, что одисциплина Сципи дает мне хорошие решения иногда , но малейшее изменение начальных условий заставляет его падать и сдаваться.Та же проблема решается довольно успешно с помощью жестких решателей MATLAB (ode15s и ode23s), но я не могу использовать это (даже из Python, потому что ни одна из привязок Python для MATLAB C API не реализует обратные вызовы, и мне нужно передать функциюрешателю ODE).Я пытаюсь PyGSL, но это ужасно сложно.Любые предложения будут с благодарностью.

РЕДАКТИРОВАТЬ: Конкретная проблема, с которой я столкнулся с PyGSL, это выбор правильной функции шага.Их несколько, но прямых аналогов к ode15 или ode23 нет (формула bdf и модифицированный Розенброк, если это имеет смысл).Так какую же хорошую шаговую функцию выбрать для жесткой системы?Я должен решить эту систему в течение очень долгого времени, чтобы гарантировать, что она достигнет стационарного состояния, и решатели GSL либо выбирают минимальный временной шаг, либо слишком большой.

Answer 1

Если вы можете решить вашу проблему с помощью Matlab's ode15s, вы сможете решить ее с помощью vode решателя scipy Для имитации ode15s я использую следующие настройки:

ode15s = scipy.integrate.ode(f)
ode15s.set_integrator('vode', method='bdf', order=15, nsteps=3000)
ode15s.set_initial_value(u0, t0)

и тогда вы можете счастливо решить свою проблему с помощью ode15s.integrate(t_final). Это должно очень хорошо сработать для решения сложной проблемы.

(см. Также http://www.scipy.org/NumPy_for_Matlab_Users)

Answer 2

Python может вызывать C. Промышленный стандарт - LSODE в ODEPACK. Это общественное достояние. Вы можете скачать C версию . Эти решатели чрезвычайно сложны, поэтому лучше использовать хорошо проверенный код.

Добавлено: Убедитесь, что у вас действительно жесткая система, т. Е. Если значения (собственные значения) отличаются более чем на 2 или 3 порядка. Кроме того, если система жесткая, но вы ищете только стационарное решение, эти решатели дают вам возможность решить некоторые уравнения алгебраически. В противном случае хороший решатель Рунге-Кутты, такой как DVERK , будет хорошим и намного более простым решением.

Добавлен здесь, потому что он не помещается в комментарии: это из документа DLSODE header:

C     T     :INOUT  Value of the independent variable.  On return it
C                   will be the current value of t (normally TOUT).
C
C     TOUT  :IN     Next point where output is desired (.NE. T).

Кроме того, да, кинетика Михаэлиса-Ментена является нелинейной. Ускорение Aitken работает с этим, все же. (Если вам нужно краткое объяснение, сначала рассмотрите простой случай, когда Y является скаляром. Вы запускаете систему, чтобы получить 3 точки Y (T). Проведите через них экспоненциальную кривую (простая алгебра). Затем установите Y в асимптоту и повторите. Теперь просто обобщите, чтобы Y был вектором. Предположим, что 3 точки находятся на плоскости - это нормально, если их нет.) Кроме того, если у вас нет функции форсирования (например, постоянная капля IV), исключение ММ затухнет прочь и система приблизится к линейности. Надеюсь, это поможет.

Answer 3

PyDSTool упаковывает решатель Radau, который является превосходным неявным жестким интегратором.Это имеет больше настроек, чем odeint, но намного меньше, чем PyGSL.Наибольшим преимуществом является то, что ваша функция RHS указана в виде строки (как правило, хотя вы можете построить систему с помощью символических манипуляций) и конвертируется в C, поэтому нет медленных обратных вызовов Python и все будет очень быстро.*

Answer 4

Я сейчас немного изучаю ODE и его решатели, поэтому ваш вопрос мне очень интересен ...

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

Итак, короче говоря, похоже, что метод bsimp стоит того, чтобы сделать снимок, хотя он требует функции якобиана. Если вы не можете вычислить якобиан, я попытаюсь использовать rk2imp, rk4imp или любой из методов экипировки.