Scipy odeint Неотрицательное решение

Question

По-видимому, получение неотрицательного решения из решателя ODE нетривиально. В Matlab есть опция NonNegative для некоторых решателей, чтобы получить неотрицательное решение. Есть ли подобная опция в scipy?

Если нет, каков "лучший" способ наложения ограничения неотрицательности? На данный момент у меня есть что-то вроде следующего:

def f(x, t, params):
     ... ... ...
     ... ... ...
     x_dot[(x <= 0) * (x_dot <= 0)] = 0.0
     return x_dot
... ... ...
x = odeint(f, x0, t, args=params)

Однако это приводит к численной нестабильности. Мне нужно было установить mxstep на 1e8 и hmin=1e-15.

Answer 1

Проблема не только в том, что вы должны избегать квадратного корня отрицательного х. Проблема в том, что "наилучший" способ наложения ограничения по-прежнему зависит от того, каково приложение вашей системы и какое поведение вы считаете "разумным". Если ваша система не имеет равновесия в 0, то ваша проблема может быть некорректной. Что может означать, что он движется с ненулевой скоростью в область минус-х? Если интерпретация заключается в том, что решение должно оставаться на нуле, то фактически у вас больше нет системы ODE в качестве вашей предполагаемой модели: у вас есть гибридная динамическая система с негладким компонентом, т.е. когда траектория x(t) достигает 0 при t = t_1, он должен оставаться в точке x(t) для всех t > t_1. Это может быть легко достигнуто с помощью соответствующего пакета динамических систем, такого как PyDSTool.

Альтернативно, x=0 является устойчивым равновесием, и вам просто нужно предотвратить оценку f для x<0. Это также можно взломать с помощью обнаружения событий.

В любом случае, обнаружение события при x = 0 сложно, когда ваш f не определен для x<0. Существует несколько стандартных решателей ODE, которые могут быть буквально вынуждены избегать оценки в подобласти при любых обстоятельствах, и большая часть обнаружения событий будет включать в себя оценки по обе стороны границы. Прагматичное решение состоит в том, чтобы выбрать небольшое число для x, ниже которого безопасно (в контексте вашего приложения) объявить x = 0. Затем заставить событие обнаруживать, когда x достигает этого (что, учитывая, что вы можете контролировать размер шага оставаться достаточно маленьким) должно препятствовать тому, чтобы x когда-либо оценивался в отрицательном значении. Затем вы бы поставили условие, чтобы x = 0 после этой точки, если это поведение, которое вы хотите. Опять же, это немного шумно в scipy/python, но вы можете сделать это. Также довольно легко настроить желаемое поведение в PyDSTool, о котором я хотел бы посоветовать, если вы будете публиковать сообщения на его справочных форумах.