I would like to set the current_state
of a stepper in boost::odeint.
I have the following MWE:
#include <queue>
#include <boost/numeric/odeint.hpp>
#include <array>
typedef std::array< double , 2 > state_type;
void eqsys(const state_type &x, state_type &dxdt, double t) {
dxdt[0] = 0.3*x[1];
dxdt[1] = 4*x[0];
}
int main() {
int steps=0;
auto stepper=make_dense_output(1.0e-6,1.0e-6,boost::numeric::odeint::runge_kutta_dopri5< state_type >() );
state_type x={2,2};
double stoptime=10;
stepper.initialize(x,0,0.01);
while(true){
//Run stepper
while( stepper.current_time()<stoptime && steps<10000 ) {
stepper.do_step(eqsys);
++steps;
}
x=stepper.current_state();
x[0]=2;
stepper.set_current_state(x);
stoptime+=10;
}
}
The line stepper.set_current_state(x);
is the one that doesn't work. Is there a command I can use to accomplish this?
I realize I could probably just use:
stepper.initialize(x, stepper.current_time(), 0.01);
but's it's unclear if this is really the best strategy, or if I'm losing some internal state of the stepper (such as the current stepsize or error estimate) which it would be better to keep.