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

Subject: Re: [Boost-users] [Odeint] Stiff System of differential equations
From: Karsten Ahnert (karsten.ahnert_at_[hidden])
Date: 2016-01-15 02:00:39

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_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/boost-users
>