限制 boost::odeint 集成的步数

Question

假设我有以下boost::odeint 代码：

#include <iostream>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>
using namespace std;
using namespace boost::numeric::odeint;

const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;

typedef boost::array< double , 3 > state_type;

void lorenz( const state_type &x , state_type &dxdt , double t ){
    dxdt[0] = sigma * ( x[1] - x[0] );
    dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
    dxdt[2] = -b * x[2] + x[0] * x[1];
}

void write_lorenz( const state_type &x , const double t ){
    cout << t << '\t' << x[0] << '\t' << x[1] << '\t' << x[2] << endl;
}

int main(int argc, char **argv){
    state_type x = { 10.0 , 1.0 , 1.0 }; // initial conditions
    cout<<"Steps: "<<integrate( lorenz , x , 0.0 , 25.0 , 0.1 , write_lorenz )<<endl;
}

如何修改代码以使集成在一定数量的步骤后中断？我正在运行大量集成，并希望避免在集成任何特定系统上花费太多时间。

我曾考虑过使用integrate_n_steps()，但这可能意味着集成会持续到我感兴趣的结束时间。

Answer 1

目前没有此任务的集成例程。不过，您有多种选择：

首先，在integrate() 中使用观察者，如果超过最大步数，则在此处抛出异常。当然，这不是很优雅：

struct write_lorenz_and_check_steps
{
    size_t m_steps;
    write_lorenz_and_check_steps( void ) : m_steps( 0 ) { }
    void operator()( const state_type &x , const double t ) const {
       cout << t << '\t' << x[0] << '\t' << x[1] << '\t' << x[2] << endl;
       ++m_steps;
       if( m_steps > max_steps ) throw runtime_error( "Too much steps" );
    }
};

// ...

size_t steps = 0;
try {
    steps = integrate( lorenz , x , 0.0 , 25.0 , 0.1 , write_lorenz );
} catch( ... ) { steps = max_steps; }
cout << steps << endl;

第二，可以自己写步进循环：

// Attention: the code has not been check to compile
double tmax = 25.0;
size_t imax = 1000;
size_t i = 0;
auto stepper = make_dense_output( 1.0e-6 , 1.0e-6 , runge_kutta_dopri5< state_type >() );
stepper.initialize( x , t , dt );
while ( ( stepper.current_time() < tmax ) && ( i < imax ) )
{
    observer( stepper.current_state() , stepper.current_time() );
    stepper.do_step( lorenz() );
    ++i;
}
x = stepper.current_state();

在此示例中，您还可以直接使用stepper.current_state() 和stepper.current_time()，而不是调用观察者。此外，如果您的编译器不支持自动，即您有一个 C++03 编译器，只需使用

typedef runge_kutta_dopri5< state_type > stepper_type;
result_of::make_dense_output< stepper_type >::type stepper =
    make_dense_output( 1.0e-6 , 1.0e-6 , stepper_type() );

我们还专门为此任务开发了一个特殊的集成例程。但它仍然需要几个星期才能完成。此外，我们还开发了可以使用的 ode 迭代器，并且很快就会准备就绪（我希望在下周的下周）。