From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 21615 invoked by alias); 3 Jan 2008 05:15:32 -0000 Received: (qmail 21598 invoked by uid 22791); 3 Jan 2008 05:15:28 -0000 X-Spam-Check-By: sourceware.org Received: from vms044pub.verizon.net (HELO vms044pub.verizon.net) (206.46.252.44) by sourceware.org (qpsmtpd/0.31) with ESMTP; Thu, 03 Jan 2008 05:14:56 +0000 Received: from PageXP ([71.175.123.164]) by vms044.mailsrvcs.net (Sun Java System Messaging Server 6.2-6.01 (built Apr 3 2006)) with ESMTPA id <0JU100DRSYKL1IE3@vms044.mailsrvcs.net> for gsl-discuss@sourceware.org; Wed, 02 Jan 2008 23:14:46 -0600 (CST) Date: Thu, 03 Jan 2008 05:15:00 -0000 From: "Michael Stauffer" Subject: RE: Bug report and fix for adaptive step size control in ODE solving In-reply-to: <4772D239.9070206@googlemail.com> To: "'Frank Reininghaus'" , Message-id: <007801c84dc7$8efd3a50$0202fea9@PageXP> MIME-version: 1.0 X-Mailer: Microsoft Outlook, Build 10.0.6822 Content-type: text/plain; charset=us-ascii Content-transfer-encoding: 7bit Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org X-SW-Source: 2008-q1/txt/msg00003.txt.bz2 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 >