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

[Brian Gough <bjg@network-theory.co.uk>]

At Wed, 26 Dec 2007 23:14:17 +0100,
Frank Reininghaus wrote:
> 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.

I've implemented it as follows, using GSL_COERCE_DBL to avoid any
optimisation or excess precision problems.  This is checked in along
with the test cases.

-- 
Brian Gough

 if (con != NULL)
    {
      /* Check error and attempt to adjust the step. */

      double h_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)
        {
          /* Check that the reported status is correct (i.e. an actual
             decrease in h0 occured) and the suggested h0 will change
             the time by at least 1 ulp */

          double t_curr = GSL_COERCE_DBL(*t);
          double t_next = GSL_COERCE_DBL((*t) + h0);

          if (fabs(h0) < fabs(h_old) && t_next != t_curr) 
            {
              /* Step was decreased. Undo step, and try again with new h0. */
              DBL_MEMCPY (y, e->y0, dydt->dimension);
              e->failed_steps++;
              goto try_step;
            }
          else
            {
              h0 = h_old; /* keep current step size */
            }
        }
    }

  *h = h0;  /* suggest step size for next time-step */