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

["Michael Stauffer" <mstauff at verizon dot net>]

Hi,

Thanks for the patch. I've had a problem with an infinite, or at least
near-infinite loop, using the ODE solver, so I'm hoping this will help me.
I'm new to GSL and ODE solvers, so I appreciate someone who understands the
complexities taking a whack at it.

Cheers,
Michael

> -----Original Message-----
> From: Frank Reininghaus [mailto:frank78ac@googlemail.com] 
> Sent: Wednesday, December 26, 2007 5:14 PM
> To: gsl-discuss@sourceware.org
> Subject: Bug report and fix for adaptive step size control in 
> ODE solving
> 
> 
> Hi,
> 
> I've found some problems with adaptive step size control 
> which can lead to an infinite loop under special 
> circumstances. I've attached three test case programs where 
> the rhs of the ODE to be solved is a step function. All 
> programs lead to infinite loops either in the the loop in my 
> code (which should be equivalent to the loop in the example 
> program
> in the documentation) or in gsl_odeiv_evolve_apply () (in the file
> ode-initval/evolve.c) for different reasons.
> 
> The bug ocurred on an AMD Athlon 64 in 32-bit mode with the 
> current CVS version of GSL. 'gcc --version' yields gcc (GCC) 
> 4.1.3 20070929 (prerelease) (Ubuntu 4.1.2-16ubuntu2). I've 
> also checked the testcases on a Pentium III with gcc (GCC) 
> 4.1.2 20061115 (prerelease) (SUSE Linux).
> 
> In testcase1.c, the function rhs (t) yields 1 for t>=1.0 and 
> 0 otherwise. Start and end points for the ODE solver are 
> t=0.0 and t=2.0, respectively, and the initial value is 
> y=0.0. The problem is that for very small epsabs, the step 
> size h0 in gsl_odeiv_evolve_apply
> () is reduced so much [by calls to gsl_odeiv_control_hadjust 
> ()] that t0 + h0 seems to be equal  to t0 due to the limited 
> precision (h0/t0 is smaller than GSL_DBL_EPSILON, defined in 
> gsl_machine.h). Therefore, the ODE solving process gets stuck 
> at the point of discontinuity (t=1.0).
> 
> This behaviour occurs not only with gsl_odeiv_step_rk2, but 
> also with other solvers (the value of epsabs necessary to get 
> into the infinite loop might vary), except for the implicit 
> solvers, but some of these yield very inaccurate results.
> 
> In testcase2.c, the rhs is 1e300 for t>=0 and 0 otherwise. 
> Integrating from t=-1.0 to t=1.0 causes step size control to 
> reduce h0 successively at the point of discontinuity until it 
> becomes zero, and this makes the ODE solver get into the same 
> infinite loop.
> 
> In testcase3.c, the problem is different. The only difference 
> from testcase2.c is the choice of a different solver (rkf45 
> instead of rk2). Now the step size h0 is not reduced to zero, 
> but only to 4.94066e-324, which seems to be the smallest 
> possible double on my machine. The problem here is that the 
> function std_control_hadjust () multiplies the step size in 
> line 113 in ode-initval/cstd.c with r=0.535431. This is 
> rounded to 4.94066e-324 by the CPU, i.e., the step size is 
> unchanged. Nevertheless, the function returns 
> GSL_ODEIV_HADJ_DEC, and this makes gsl_odeiv_evolve_apply () 
> believe that the step size was reduced, and therefore it 
> jumps back to the label 'try_step' and does exactly the same 
> thing as before with the unchanged step size, so we end up in 
> an infinite loop in gsl_odeiv_evolve_apply ().
> 
> I've attached a suggestion for a fix. In the current CVS 
> version, gsl_odeiv_evolve_apply () jumps back to the label 
> 'try_step' if gsl_odeiv_control_hadjust () returns 
> GSL_ODEIV_HADJ_DEC, but in my version, it first makes sure 
> that h0 does not get smaller than the maximum of GSL_DBL_MIN 
> and GSL_DBL_EPSILON * t0 to avoid the problems described 
> above. Before it jumps back to 'try_step', it checks if h0 is 
> really smaller than in the previous step (we would get into 
> an infinite loop otherwise).
> 
> Note that std_control_hadjust () still returns 
> GSL_ODEIV_HADJ_DEC although the step size is unchanged in 
> testcase3.c, but the new version of gsl_odeiv_evolve_apply () 
> detects this.
> 
> Regards,
> Frank
>