Интегрируйте жесткие ODE с Python

Question

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

EDIT: Конкретная проблема, с которой я столкнулась с PyGSL, заключается в выборе правильной функции шага. Есть несколько из них, но нет прямых аналогов ode15s или ode23s (формула bdf и модифицированный Rosenbrock, если это имеет смысл). Итак, что такое хорошая функция для выбора жесткой системы? Я должен решить эту систему на очень долгое время, чтобы убедиться, что она достигнет устойчивого состояния, и решатели GSL либо выбирают минимальный временной шаг, либо слишком большой.

Answer 1

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

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

Добавлен здесь, потому что не будет вписываться в комментариях: Это из док заголовка DLSODE:

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 2

Я в настоящее время изучает немного ОДУ и ее решателей, поэтому ваш вопрос очень интересен для меня ...

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

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

Answer 3

Если вы можете решить вашу проблему с помощью 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 4

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