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

[Frank Reininghaus <frank78ac@googlemail.com>]

[-- Attachment #1: Type: text/plain, Size: 3309 bytes --]

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

[-- Attachment #2: testcase1.c --]
[-- Type: text/x-csrc, Size: 1067 bytes --]

/* compile and link with
gcc -o testcase1 testcase1.c -lgsl -lgslcblas -lm
*/

#include <stdio.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_errno.h>

int rhs (double t, const double * y, double * dydt, void * params) {
  if (t >= 1.0)
    dydt [0] = 1;
  else
    dydt [0] = 0;

  return 0;
}

int main (int argc, char ** argv) {
  /* infinite loop for epsabs = 1e-18, but not for 1e-17 */
  double epsabs = 1e-18;
  double epsrel = 1e-6;

  const gsl_odeiv_step_type * T = gsl_odeiv_step_rk2;
  gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 1);
  gsl_odeiv_control * c = gsl_odeiv_control_y_new (epsabs, epsrel);
  gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (1);
  gsl_odeiv_system sys = {rhs, 0, 1, 0};
     
  double t = 0.0;
  double h = 1e-6;
  double y = 0.0;
     
  while (t < 2.0) {
    int status = gsl_odeiv_evolve_apply (e, c, s,
					 &sys, 
					 &t, 2,
					 &h, &y);
    if (status != GSL_SUCCESS)
      break;
  }

  printf ("%g\n", y);
     
  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);
  return 0;
}


[-- Attachment #3: testcase2.c --]
[-- Type: text/x-csrc, Size: 1074 bytes --]

/* compile and link with
gcc -o testcase2 testcase2.c -lgsl -lgslcblas -lm
*/

#include <stdio.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_errno.h>

int rhs (double t, const double * y, double * dydt, void * params) {
  if (t >= 0.0)
    dydt [0] = 1e300;
  else
    dydt [0] = 0;

  return 0;
}

int main (int argc, char ** argv) {
  /* infinite loop for epsabs = 1e-25, but not for 1e-24 */
  double epsabs = 1e-25;
  double epsrel = 1e-6;

  const gsl_odeiv_step_type * T = gsl_odeiv_step_rk2;
  gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 1);
  gsl_odeiv_control * c = gsl_odeiv_control_y_new (epsabs, epsrel);
  gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (1);
  gsl_odeiv_system sys = {rhs, 0, 1, 0};
     
  double t = -1.0;
  double h = 1e-6;
  double y = 0.0;
     
  while (t < 1.0) {
    int status = gsl_odeiv_evolve_apply (e, c, s,
					 &sys, 
					 &t, 1.0,
					 &h, &y);
    if (status != GSL_SUCCESS)
      break;
  }

  printf ("%g\n", y);
     
  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);
  return 0;
}


[-- Attachment #4: testcase3.c --]
[-- Type: text/x-csrc, Size: 1076 bytes --]

/* compile and link with
gcc -o testcase3 testcase3.c -lgsl -lgslcblas -lm
*/

#include <stdio.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_errno.h>

int rhs (double t, const double * y, double * dydt, void * params) {
  if (t >= 0.0)
    dydt [0] = 1e300;
  else
    dydt [0] = 0;

  return 0;
}

int main (int argc, char ** argv) {
  /* infinite loop for epsabs = 1e-26, but not for 1e-25 */
  double epsabs = 1e-26;
  double epsrel = 1e-6;

  const gsl_odeiv_step_type * T = gsl_odeiv_step_rkf45;
  gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 1);
  gsl_odeiv_control * c = gsl_odeiv_control_y_new (epsabs, epsrel);
  gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (1);
  gsl_odeiv_system sys = {rhs, 0, 1, 0};
     
  double t = -1.0;
  double h = 1e-6;
  double y = 0.0;
     
  while (t < 1.0) {
    int status = gsl_odeiv_evolve_apply (e, c, s,
					 &sys, 
					 &t, 1.0,
					 &h, &y);
    if (status != GSL_SUCCESS)
      break;
  }

  printf ("%g\n", y);
     
  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);
  return 0;
}


[-- Attachment #5: evolve.c.patch --]
[-- Type: text/x-diff, Size: 1236 bytes --]

--- cvs/gsl/ode-initval/evolve.c	2007-12-14 23:53:49.000000000 +0100
+++ cvs-no-ode-bug/gsl/ode-initval/evolve.c	2007-12-26 22:51:34.000000000 +0100
@@ -199,15 +199,28 @@
   if (con != NULL)
     {
       /* Check error and attempt to adjust the step. */
+      const double h0_old = h0;
       const int hadjust_status 
         = gsl_odeiv_control_hadjust (con, step, y, e->yerr, e->dydt_out, &h0);
 
       if (hadjust_status == GSL_ODEIV_HADJ_DEC)
         {
-          /* Step was decreased. Undo and go back to try again. */
-          DBL_MEMCPY (y, e->y0, dydt->dimension);
-          e->failed_steps++;
-          goto try_step;
+	  /* Step size h0 was decreased. Make sure it does not get too small:
+	     It should be positive (i.e., h0 >= GSL_DBL_MIN) and not 
+	     smaller than GSL_DBL_EPSILON * t0 (because t0 + h0 has to be 
+	     distinguishable from t0) */
+
+	  h0 = GSL_MAX_DBL (h0, GSL_MAX_DBL (GSL_DBL_MIN, GSL_DBL_EPSILON * t0));
+
+          /* If the step size h0 is still smaller than in the previous step, 
+	     undo and go back to try again. */
+
+	  if (h0 < h0_old) 
+	    {
+	      DBL_MEMCPY (y, e->y0, dydt->dimension);
+	      e->failed_steps++;
+	      goto try_step;
+	    }
         }
     }