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