URL: <https://savannah.gnu.org/bugs/?66849>
Summary: gsl_odeiv2_evolve_apply() may exceed final time
Group: GNU Scientific Library
Submitter: mangrove84
Submitted: Wed 26 Feb 2025 08:57:01 AM CET
Category: None
Severity: 3 - Normal
Operating System:
Status: None
Assigned to: None
Open/Closed: Open
Discussion Lock: Any
Release: v2.8
_______________________________________________________
Follow-up Comments:
-------------------------------------------------------
Date: Wed 26 Feb 2025 08:57:01 AM CET By: Christian Krueger <mangrove84>
The function gsl_odeiv2_evolve_apply() may exceed final time in rare cases
(loss of accuracy) as described below.
When evolving an ODE system using gsl_odeiv2_evolve_apply(), the user provides
the current time t, the maximum time t1, and a suggested time step h. At some
point, the final step is reached when the suggested time step h is larger than
the remaining time dt = t1 - t. Now gsl_odeiv2_step_apply() is called
providing the current time t and the requested time step dt. Then, the stepper
evaluates the right-hand side of the ODE at the times t + a[i]*dt, where a[i]
are the intermediate RK steps. Importantly, the final evaluation time is t +
dt.
Consider now this problem: In my case, I integrate an ODE system from a "time"
t = 1 in small steps toward zero (i.e. backwards). However, my final time t1
is not exactly zero but something like t1 = 1e-12. Let's assume a fixed step
size of 1e-3, which results in 1000 steps. In the last step,
gsl_odeiv2_evolve_apply() will calculate dt = t1 - t = 1e-12 - 0.001 (loss of
accuracy!), and then pass t = 0.001 and dt to gsl_odeiv2_step_apply().
Here is the issue: The time of the last RK substep is t + dt = 0.001 + (1e-12
- 0.001). This calculation suffers from loss of accuracy and *is not
guaranteed to not exceed the final time t1*.
I will try to create an MWE with a simple ODE soon. Second, given the current
implementation of the functions in GSL, I do not see how this could easily be
fixed; also, the problem should appear only in such cases where the final time
is very close to but not exactly zero.
For the moment, my solution is to pass the final time t1 also as a parameter
to my RHS function and check if the time is out of bounds and return GSL_EDOM
in such case. This means that the adaptive step sizing will create a few more
steps of smaller size instead of one single last step.
_______________________________________________________
Reply to this item at:
<https://savannah.gnu.org/bugs/?66849>
_______________________________________________
Message sent via Savannah
https://savannah.gnu.org/
signature.asc
Description: PGP signature
