python:满足条件时中断odeint

Question

我使用的是scipy.integrate包中的odeint函数：

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个时间步长上的点的位置。

现在，我想在满足r上的条件时停止odeint求解器。特别是，当其中两个点接近某个阈值时，我希望停止odeint，对r向量进行一些更改，并使用新的初始位置继续odeint求解器。有没有办法实现这一点？

我一直在考虑的一种可能的解决方案是运行odeint到最后，然后再检查条件是否满足，但这当然不是很有效。

谢谢你们的帮助，尼古拉

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来实现的。