I have started using boost::odeint in my C++ code, and I think I'm missing a simple feature available in other integrators, namely Scipy's odeint.
scipy.odeint lets the user specify the times at which the state must be added to the output state history.scipy.odeint is a variable timestep integrator whose one-liner call looks like this (the state is integrated from the initial condition X0 and interpolated/stored at the times specified in times)
X = scipy.odeint(dxdt,X0,times,atol = 1e-13,rtol = 1e-13)
where X is a matrix that has as many rows as there are elements in times
Basically, I am looking for a similar feature in boost::odeint in order to do two things:
t0 to tf, but only retrieve the final value of the state. I think I could write an observer only storing the state if the internal time satisfies t == tf, but this looks like a rather ugly hack. If I want to let the integrator choose the proper internal timestep to meet the tolerance values, storing intermediate states is a unnecessary burden.t0 to tf but store the state at times specified in advance, that are not necessarily evenly distributed, in a similar fashion to the call to scipy.odeint hereabove, while also letting the integrator choose the proper internal timestep.The closest I've been to achieving that is with the following
size_t steps = integrate_adaptive( make_controlled< error_stepper_type >( 1.0e-10 , 1.0e-16 ) ,
dynamics , x , 0.0 , 10.0 , 1. , push_back_state_and_time( x_vec , times ));
The tolerances are met, but all the states are stored into x_vec by the observer, without letting me specify what the storage times should be.
How should I proceed?
It seems you are looking for the integrate_times function:
It lets you specify a range of exact times for which the observer will be invoked, adjusting the step size to reach every time step exactly, if necessary.
Especially for adaptive methods, this is quite useful as it computes the solution at the exact times you specified while still controlling the time step size to not exceed the error bounds.
So your current call could be modified to something like
auto stepper = make_controlled<error_stepper_type>( 1.0e-10 , 1.0e-16 );
// std::vector<time> times;
// std::vector<state> x_vec;
// time tstep;
auto tbegin = times.begin();
auto tend = times.end();
integrate_times(stepper, dynamics, x, tbegin, tend, tstep, push_back_state(x_vec));