4

次の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()たが、これは、関心のある終了時間を過ぎて統合が進行することを意味する可能性があります。

4

1 に答える 1

5

現時点では、このタスクの統合ルーチンはありません。それにもかかわらず、いくつかのオプションがあります。

まず、 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()andを直接操作します。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 イテレータを開発し、すぐに準備が整う予定です (来週の来週になることを願っています)。

于 2013-01-22T20:10:30.187 に答える