1

我正在使用 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 求解器。特别是当其中 2 个点比某个阈值更接近时，我想停止 odeint，对 r 向量进行一些更改，并使用新的初始位置继续 odeint 求解器。有没有办法实现这个？

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

谢谢大家的帮助，尼古拉

4

1 回答 1

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

于 2015-03-23T16:42:11.590 回答