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
|