假设我有以下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(),但这可能意味着集成将持续到我感兴趣的结束时间之后。

最佳答案

目前没有用于此任务的集成例程。但是,您有几种选择:

首先,在Integrated()中使用观察者,如果超过最大步数,则会在其中引发异常。当然,这不是很优雅:

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迭代器,该迭代器也可以使用,并且很快就会准备就绪(我希望在下周的下一个星期)。

10-07 19:17
查看更多