питона: прерывание odeint при выполнении условия

Question

Я использую odeint функцию из scipy.integrate пакета:

r0 = np.array([1,2,3,4]) 
t=np.linspace(0,1,20) 
def drdt(r,t): 
    return r # or whatever else 
r = odeint(drdt,r0,t) 

r0 является NumPy массива, который содержит начальные позиции определенного количества очков. В конце скрипта, как и ожидалось, я получаю позиции точек по 20 меток времени.

Теперь я хотел бы остановить решатель odeint, когда выполняется условие на r. В частности, я хотел бы остановить odeint, когда 2 из этих точек ближе определенного порога, внести некоторые изменения в вектор r и продолжить решение odeint с новыми начальными позициями. Есть ли способ реализовать это?

Возможное решение, о котором я думал, заключается в том, чтобы запустить odeint до конца, а позже проверить, выполнено ли условие, но это, конечно, не очень эффективно.

Спасибо всем за помощь, Nicola

Answer 1

У меня есть ответ на C++. Возможно, это не лучшее место для публикации, но это может быть интересно. (Я не нашел лучшего места, чтобы выразить это, и именно здесь я приземлился при поиске решения на C++).

Ниже приведен пример C++, который останавливает интеграцию, когда переменная становится равной или ниже нуля.

#include <iostream> 
#include <boost/range/algorithm.hpp> 
#include <boost/numeric/odeint.hpp> 

using namespace std; 
using namespace boost::numeric::odeint; 


typedef double state_type; 

typedef runge_kutta_cash_karp54<state_type> error_stepper_type; 


class Myode{ 
public: 
    void operator()(const state_type& x, state_type& dxdt, const double t){dxdt=-1-x;} 
    static bool done(const state_type &x){return x<=0.0;} 
}; 

int main(int argc, char* argv[]){ 
    Myode ode; 
    state_type x=10.0; 
    auto stepper = make_controlled<error_stepper_type>(1.0e-10 , 1.0e-6); 
    auto iter= boost::find_if(make_adaptive_range(stepper,ode,x, 0.0 , 20.0 , 0.01), 
       Myode::done); 

    cout<<"integration stopped at"<<x<<endl; 
    return 1; 
} 

Интеграция останавливается в первый раз, когда достигает значения x, меньшего или равного нулю (см. Выполненную функцию). Поэтому, в зависимости от вашего текущего размера шага, он может быть намного ниже нуля.

Обратите внимание, что это использует конструкции C++ 11, поэтому вам нужно включить это в свой компилятор. В моем случае (gcc 4.4) это было достигнуто добавлением -std = gnu ++ 0x в команду компиляции.