[Bug-gsl] [bug #52336] odeiv2/msbdf crashes with method order 3 and maximum step length increase

[Tuomo Keskitalo]

URL:
  <http://savannah.gnu.org/bugs/?52336>

                 Summary: odeiv2/msbdf crashes with method order 3 and maximum
step length increase
                 Project: GNU Scientific Library
            Submitted by: tkeskita
            Submitted on: Sat 04 Nov 2017 07:44:57 AM UTC
                Category: Runtime error
                Severity: 3 - Normal
        Operating System: 
                  Status: None
             Assigned to: tkeskita
             Open/Closed: Open
                 Release: 1.16
         Discussion Lock: Any

    _______________________________________________________

Details:

>From Andrew Benson on 26 Oct 2017:

Hi Tuomo,

I've been making good use of ode-initval2 for the past seven years(!), but
recently ran into a problem. Using the msbdf algorithm I sometimes get
floating point exceptions. I tracked this down to the following lines in
msbdf.c (I'm currently using GSL v1.16 for this):

*ordm1coeff = fabs (c1 / (x / l[ord]));

*ordp1coeff = fabs (a2 / (l[ord] * (h / hsum) / x));

The problem is that x and l[ord] are both equal to zero, leading to a floating
exception. This seems to be caused by the fact that h/hprev[0]=5 - for that
exact ratio, and for ord=3, the calculation leads to l[1]=1+5/6=1.833333 and
coeff1=-1-5/6=-1.83333 - as a result x=0.

The reason that h/hprev[0]=5 precisely is that in cscal.c sc_control_hadjust()
limits step size increases to a factor of 5.

I've worked around this problem by changing that factor in
sc_control_hadjust() from 5 to something slightly different (e.g. 5.1). But,
this seems like an ugly solution!

Is there a more intelligent solution that could be taken in msbdf.c in the
case thath/hprev[0]=5?

Thanks,

Andrew




    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?52336>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/