Скажи, что у меня есть следующее 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()
, но это может означать, что интеграция продолжается после того времени, которое меня интересует.
в настоящий момент для этой задачи нет процедуры интеграции. Тем не менее, у вас есть несколько вариантов:
Первый, используйте наблюдателя в 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()
вместо вызова наблюдателя. Кроме того, если ваш компилятор не поддерживает auto, то есть у вас есть компилятор 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, которые также могут быть использованы и которые будут готовы в ближайшее время (надеюсь на следующей неделе).
Других решений пока нет …