Re: [Boost-users] [Odeint] Stiff System of differential equations

Subject: Re: [Boost-users] [Odeint] Stiff System of differential equations
From: Jules Tamagnan (jtamagnan_at_[hidden])
Date: 2016-01-14 16:48:34

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.

> 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.

> 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).

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