限號:: 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()); 

我們也正在開發一種特殊的正是這個任務整合程序。但在完成之前還需要幾周的時間。此外，我們開發了也可以使用並且很快就會準備就緒的迭代器（我希望在下週的下一週）。