Re: Bug report and fix for adaptive step size control in ODE solving

[Jochen Küpper <jochen at fhi-berlin dot mpg dot de>]

Dear Frank, All,


On 26.12.2007, at 23:14, Frank Reininghaus wrote:

> I've found some problems with adaptive step size control which can  
> lead to an infinite loop under special circumstances.


Nice write-up of this problem!

Although I have not looked through your cases in detail, I can  
conform that we have had similar infinite-loop problems, for example,  
on AMD Opteron systems running Linux and operated in 64bit mode  
(gcc-4.x with -m64).

Because I did not have time to look into it in detail, and because  
the q&d-hack below works for us, so far we use the following code  
snippet to work around it:

                while(t < tend) {
                    if(GSL_SUCCESS != gsl_odeiv_evolve_apply(evolve,  
control, step,
                                                              
&system, &t, tend, &h, y)) {
                        throw Error("Flight::operator(): integration  
problem in GSL");
                    }
                    if(h < hmin) {
                        h = hinit;
                    }
                }

with (for example)
  hinit = tstep * 0.1
  hmin  = tstep * 0.01

We use this with gsl_odeiv_step_rkf45, gsl_odeiv_step_rkck, or  
gsl_odeiv_step_rk8pd.


We would be quite interested in seeing this fixed in GSL, obviously.


I'll try to find some time to test your patch. However, why don't you  
test what you really describe your comment? This would look like the  
following (untested):


  /* Step size h0 was decreased. Make sure it does not get too small:
     It should be positive (i.e., h0 >= GSL_DBL_MIN) and large enough
     to make t0 + h0 distinguishable from t0, i.e., at least of size
     GSL_DBL_EPSILON*t0 */

  h0 = GSL_MAX_DBL(GSL_DBL_MIN, (t0 != t0+h0) ? h0 :  
GSL_DBL_EPSILON*t0);


This also avoids the multiplication in most cases.

Greetings,
Jochen
--
Einigkeit und Recht und Freiheit                http://www.Jochen- 
Kuepper.de
    Liberté, Égalité, Fraternité                GnuPG key: CC1B0B4D
        Sex, drugs and rock-n-roll

Attachment: PGP.sig
Description: This is a digitally signed message part