Scipy odeint неотрицательный раствор

Question

По-видимому, getting a non-negative solution from an ODE solver is non-trivial. В Matlab для некоторых решателей есть NonNegative option, чтобы получить неотрицательное решение. Есть ли аналогичный вариант в 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

Проблема заключается не только в том, что вам нужно избегать квадратного укоренения отрицательного x. Проблема в том, что «лучший» способ наложения ограничения по-прежнему зависит от того, что приложение вашей системы и какое поведение вы считаете «разумным». Если ваша система не имеет равновесия в 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, о котором я хотел бы вам посоветовать, если вы публикуете его на форумах помощи.