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

Я использую функцию 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 до конца, а потом только проверить, было ли выполнено условие, но это, конечно, не очень эффективно.

Спасибо всем за помощь, Николай


odeint
person Nicola Gritti    schedule 03.12.2014    source источник

Ответы (1)


arrow_upward
1
arrow_downward

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

Ниже приведен пример 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, меньшего или равного нулю (см. функцию done). Таким образом, в зависимости от вашего текущего размера шага, он может быть значительно ниже нуля.

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

person citronrose    schedule 23.03.2015