This is the mail archive of the
gsl-discuss@sourceware.org
mailing list for the GSL project.
Re: Bug report and fix for adaptive step size control in ODE solving
- From: Brian Gough <bjg at network-theory dot co dot uk>
- To: Frank Reininghaus <frank78ac at googlemail dot com>
- Cc: gsl-discuss at sourceware dot org
- Date: Fri, 04 Jan 2008 09:46:23 +0000
- Subject: Re: Bug report and fix for adaptive step size control in ODE solving
- References: <4772D239.9070206@googlemail.com>
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 */