bugGNU Scientific Library - Bugs: bug #21933, adaptive step size control in ODE...

 
 

You are not allowed to post comments on this tracker with your current authentication level.

bug #21933: adaptive step size control in ODE solving

Submitted by:  -Deleted Account- <bjg>
Submitted on:  Wed 02 Jan 2008 07:00:15 PM UTC  
 
Category: Runtime errorSeverity: 3 - Normal
Operating System: Status: Fixed
Assigned to: -Deleted Account- <bjg>Open/Closed: Closed
Release: 1.10

Fri 04 Jan 2008 10:57:38 AM UTC, comment #1:

added test cases and patch (with modifications)

-Deleted Account- <bjg>In charge of this item.
Wed 02 Jan 2008 07:00:15 PM UTC, original submission:

From: Frank Reininghaus <frank78ac@googlemail.com>
To: -unavailable-
Subject: Bug report and fix for adaptive step size control in ODE solving
Date: Wed, 26 Dec 2007 23:14:17 +0100

[1 <text/plain; ISO-8859-1 (7bit)>]
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

-Deleted Account- <bjg>In charge of this item.

 

Attached Files
file #14739:  testcase1.c added by bjg (1KiB - application/octet-stream)
file #14740:  testcase2.c added by bjg (1KiB - application/octet-stream)
file #14741:  testcase3.c added by bjg (1KiB - application/octet-stream)
file #14742:  evolve.c.patch added by bjg (1KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

CC list is empty

 

Please enter the title of George Orwell's famous dystopian book (it's a date):

 

 

Follow 6 latest changes.

Date Changed By Updated Field Previous Value => Replaced By
Fri 04 Jan 2008 10:57:58 AM UTCbjgOpen/ClosedOpen=>Closed
Fri 04 Jan 2008 10:57:38 AM UTCbjgStatusConfirmed=>Fixed
Wed 02 Jan 2008 07:00:15 PM UTCbjgAttached File-=>Added testcase1.c, #14739
  Attached File-=>Added testcase2.c, #14740
  Attached File-=>Added testcase3.c, #14741
  Attached File-=>Added evolve.c.patch, #14742

Back to the top


Powered by Savane 3.1-cleanup1