On 14.01.2016 22:48, Jules Tamagnan wrote:

Hi Karsten,

Thank you for your answer.

...

How do you use the solver? Runge Kutta 78 is an adaptive stepper, but you need to put it into the controlled_runge_kutta stepper.

This is how I used the runge_kutta_fehlberg78 stepper initally.

...

void Printer_Observer( const model_type &x , double t ){ printf("%g\n", t); } typedef runge_kutta_fehlberg78< model_type , double , model_type , double , vector_space_algebra > stepper; double abs_err = .1; double rel_err = .00001; double step_size = .02; int steps = integrate_const( make_controlled<stepper>( abs_err , rel_err ) , model2, x, 0.0 , 10.0 , step_size, Printer_Observer);

Running it this way works but it takes a long time (which is fine and is to be expected). The results though seem to be depending on the abs_err and rel_err and stepsize very heavily. How do I know at what point the stepsize, abs_err and rel_err are small enough that I'll get the "right" answer? Also it gets stuck at t=0 if I set step_size=1.0.

In what way depend the results on the step size? In you example, the step size that you pass to the integrate function is in principle the observation time step. But it is also the initial time step which might be in appropriate for your use case. Can you try to solve the above example with integrate times: vector< double > obs_times = { 0.0 , 0.1 , 0.2 , 0.3 , 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 1.0 }; double step_size = 0.02; // now it is really the initial step size integrate_times( make_controlled< stepper >( abs_err , rel_err ) , model2 , x , step_size , times.begin() , times.end() , Print_observer );

...

Furthermore, it is not a dense output solvers. This means that the solver need to reach to the time points where you want to observe your solution. Therefore, it needs to change the step size from the adaptive scheme just to reach the appropriate observation time.

This makes sense but I don't feel like that should cause the solver to return an entirely different answer. I guess that is just the way of numeric solutions to ODEs though.

Yes, this is in principle correct. The problem with the above controlled stepper and integrate_const that they need to reach exactly the observation points. Hence, if you want to observe at say t = 0.5 the controlled stepper will break the time step adaption scheme to reach exactly at t = 0.5. This is not the case for the dense output steppers. Here intermediate values of your ode are interpolated (with the same accuracy as the solver).

...

Maybe should try the dopri5 or bulirsh_stoer stepper, since they can use dense output and they should not be sensitive on the observation time and the time step.

I've since tried running the bulirsch_stoer stepper like this:

...

void Printer_Observer( const model_type &x , double t ){ printf("%g\n", t); } double abs_err = .1; double rel_err = .00001; double step_size = .02; bulirsch_stoer_dense_out< model_type , double , model_type , double , vector_space_algebra > stepper {abs_err, rel_err}; int steps = integrate_const(stepper, model2, x, 0.0 , 10.0 , step_size, Printer_Observer);

The program compiles and runs but it seems like it gets stuck at t=0.02 (10+ minutes without moving anywhere).

Hmm, this really looks like you have a problem with your ODE. You might be right, it could be too stiff. For analysis purposes you could also try to unroll the inner integration look to see what is happening here in detail: while( t < end_time ) { observ( x , t ); controlled_step_result res; size_t trials = 0; do { res = stepper.try_step( system , x , t , dt ); ++trials; // add some analysis to see the step size, e.g. // cout << t << " " << dt << endl; } while( ( res == fail ) && ( trials < max_attempts ) ); if( trials == max_attempts ) throw std::runtime_error( "Too many steps" ); } If you system is stiff you should also look at the Rosenbrock method. It is much better suited for stiff problems. But its needs to calculate the Jacobian.

I've also tried with the runge_kutta_dopri5 with dense output:

...

void Printer_Observer( const model_type &x , double t ){ printf("%g\n", t); } double abs_err = .1; double rel_err = .00001; double step_size = .02; int steps = integrate_const( make_dense_output( abs_err , rel_err , runge_kutta_dopri5< model_type , double , model_type, double , vector_space_algebra >()) , model2, x, 0.0 , 10.0 , step_size, Printer_Observer);

Here it also gets stuck at t=.02 and if I try with a step_size=1.0 it will get stuck at t=0.

It seems like these steppers are also affected by the step_size and that they fail in the same way that the runge_kutta_fehlberg78 stepper did. Did these fail because they were improperly instantiated? Is the problem somehow too large or too stiff?

Many thanks for your continued help, Jules _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users

...

vector< double > obs_times = { 0.0 , 0.1 , 0.2 , 0.3 , 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 1.0 }; double step_size = 0.02; // now it is really the initial step size integrate_times( make_controlled< stepper >( abs_err , rel_err ) , model2 , x , step_size , times.begin() , times.end() , Print_observer );

I used integrate_times as you described above but it also froze after the .1 step. After having waited for a very long time I eventually get this error message """ terminate called after throwing an instance of 'boost::exception_detail::clone_implstd::overflow_error >' what(): Integrate adaptive : Maximal number of iterations reached. A step size could not be found. """

Might it be possible that you reached something like a singular point? If the stepper does not find a step size it usually hints to some problems with the ODE. Can you show me roughly your equations?

...

If you system is stiff you should also look at the Rosenbrock method. It is much better suited for stiff problems. But its needs to calculate the Jacobian.

The biggest issue is that the Jacobian is too complicated to calculate, it depends on some population dynamics and higher order interpolation.

Is there no other way to solve stiff systems. I've read about other methods such as CVODE and LSODA, but couldn't figure out how to get those to work.

These methods exists but not within odeint and there are currently not reachable from odeint. You might have a look at here: https://computation.llnl.gov/casc/sundials/description/description.html