From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 29763 invoked by alias); 26 Dec 2007 22:14:29 -0000 Received: (qmail 29754 invoked by uid 22791); 26 Dec 2007 22:14:28 -0000 X-Spam-Check-By: sourceware.org Received: from fk-out-0910.google.com (HELO fk-out-0910.google.com) (209.85.128.187) by sourceware.org (qpsmtpd/0.31) with ESMTP; Wed, 26 Dec 2007 22:14:22 +0000 Received: by fk-out-0910.google.com with SMTP id 26so3392240fkx.8 for ; Wed, 26 Dec 2007 14:14:19 -0800 (PST) Received: by 10.82.158.12 with SMTP id g12mr13262884bue.23.1198707259389; Wed, 26 Dec 2007 14:14:19 -0800 (PST) Received: from ?192.168.178.21? ( [87.189.25.32]) by mx.google.com with ESMTPS id j2sm21256824mue.3.2007.12.26.14.14.17 (version=TLSv1/SSLv3 cipher=RC4-MD5); Wed, 26 Dec 2007 14:14:18 -0800 (PST) Message-ID: <4772D239.9070206@googlemail.com> Date: Wed, 26 Dec 2007 22:14:00 -0000 From: Frank Reininghaus User-Agent: Thunderbird 2.0.0.6 (X11/20071022) MIME-Version: 1.0 To: gsl-discuss@sourceware.org Subject: Bug report and fix for adaptive step size control in ODE solving Content-Type: multipart/mixed; boundary="------------020808050805000801090908" 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: 2007-q4/txt/msg00047.txt.bz2 This is a multi-part message in MIME format. --------------020808050805000801090908 Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit Content-length: 3309 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 --------------020808050805000801090908 Content-Type: text/x-csrc; name="testcase1.c" Content-Transfer-Encoding: 7bit Content-Disposition: inline; filename="testcase1.c" Content-length: 1067 /* compile and link with gcc -o testcase1 testcase1.c -lgsl -lgslcblas -lm */ #include #include #include 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; } --------------020808050805000801090908 Content-Type: text/x-csrc; name="testcase2.c" Content-Transfer-Encoding: 7bit Content-Disposition: inline; filename="testcase2.c" Content-length: 1074 /* compile and link with gcc -o testcase2 testcase2.c -lgsl -lgslcblas -lm */ #include #include #include 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; } --------------020808050805000801090908 Content-Type: text/x-csrc; name="testcase3.c" Content-Transfer-Encoding: 7bit Content-Disposition: inline; filename="testcase3.c" Content-length: 1076 /* compile and link with gcc -o testcase3 testcase3.c -lgsl -lgslcblas -lm */ #include #include #include 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; } --------------020808050805000801090908 Content-Type: text/x-diff; name="evolve.c.patch" Content-Transfer-Encoding: 7bit Content-Disposition: inline; filename="evolve.c.patch" Content-length: 1236 --- 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; + } } } --------------020808050805000801090908--