bugGNU Scientific Library - Bugs: bug #52336, odeiv2/msbdf crashes with method...

 
 

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

bug #52336: odeiv2/msbdf crashes with method order 3 and maximum step length increase

Submitted by:  Tuomo Keskitalo <tkeskita>
Submitted on:  Sat 04 Nov 2017 07:44:57 AM UTC  
 
Category: Runtime errorSeverity: 3 - Normal
Operating System: Status: Fixed
Assigned to: Tuomo Keskitalo <tkeskita>Open/Closed: Closed
Release: 1.16

Sat 04 Nov 2017 08:24:25 AM UTC, comment #2:

As a fix, changed max step scaling to 4.9 in cscal.c. Implemented in commit 9cd8a8d58c499232689d6ee849e24724c3a181cc

Tuomo Keskitalo <tkeskita>
Project MemberIn charge of this item.
Sat 04 Nov 2017 07:50:14 AM UTC, comment #1:

Added debug output from Andrew to bug report.

It looks like you have found a rare collision. The errors of your function just before crashing are suddely all smaller continuously, so controller is accelerating and finally reaches maximum acceleration, which hits the singularity of the solver at ord=3 as you described.

Also you are right about the easy cure: r could be limited in cscal.c to something else than 5.0. I think that 4.9 is better than 5.1, because it eliminates possibility for hitting 5.0 accidentally. The maximum value 5.0 is anyway more or less arbitrary, so I'm not worried about changing it in any case.

(file #42326)

Tuomo Keskitalo <tkeskita>
Project MemberIn charge of this item.
Sat 04 Nov 2017 07:44:57 AM UTC, original submission:

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

Tuomo Keskitalo <tkeskita>
Project MemberIn charge of this item.

 

Attached Files
file #42326:  benson_debug.txt added by tkeskita (172KiB - text/plain - Debug of case from Andrew Benson)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -unavailable- added by tkeskita (Submitted the item)
  •  

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

     

     

    Follow 3 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Sat 04 Nov 2017 08:24:25 AM UTCtkeskitaStatusNone=>Fixed
      Open/ClosedOpen=>Closed
    Sat 04 Nov 2017 07:50:14 AM UTCtkeskitaAttached File-=>Added benson_debug.txt, #42326

    Back to the top


    Powered by Savane 3.1-cleanup1