Bug report and fix for adaptive step size control in ODE solving
Michael Stauffer
mstauff@verizon.net
Thu Jan 3 05:15:00 GMT 2008
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
>
More information about the Gsl-discuss
mailing list