I am on the fence on this. Looking at your proposal, it seems you want to accept a root when in fact it really does not belong to the interval. The fact the exception message shows 0 (or -0) does not mean the value is exactly zero, but rather a small number not displayed (like +1.0e-15 or something like that). The tests on the boundary are written in such a way they match “perfect root” (this is what is written in the comment), and by “perfect” we mean exact 0, not close to zero within some tolerance.
any reason to not do that ?
Yes, there are reasons not to do that, and these reasons are in fact related to event detection. In many cases, events are set up to separate processing of the “before” and “after” phases. This sometimes implies that a computation is changed and restarted or reset at the detected event, which implies that one of the boundary for the next phase (the after phase if propagating forward, the before phase if propagating backward) becomes by construction one of the detected zero. So basically, just when we restart, we are already at zero. If we consider that when we are close to zero we ignore the sign and consider we have found a root, then we will stall at the zero and never cross it. This was exactly what caused years ago the unstability on events detection we corrected since then. Now, we take great care to make sure we cross the zero and once we have crossed it we don’t want to detect the same zero again. If for example we know we have already considered a transition from negative to positive, if we start the next phase with a value y(t1) = +1.0e-300 and y(t2) = +0.2, then we consider this interval does not bracket a root. The root was in the previous interval and we consider a root at the boundary between two intervals is in one of the intervals only.
the epsilon is arbitrary set to 1e-12 , is there a value that is already used somewhere which would be more significant or configurable ?
There is already a configurable value that has this meaning. You can get it by calling getFunctionValueAccuracy()
, the default value is 1.0e-14 and it can be set to another value by using one of the constructors.
can this behaviour be included in a future release of hipparchus ? should I contact directly hipparchus ?
This has to be discussed, because as I explained before, there is a reason for this very strict handling of bracketing intervals. For discussing this, yes the Hipparchus project is the place to go. You should open an issue here https://github.com/Hipparchus-Math/hipparchus/issues. You will stay in well-known community since the two Hipparchus developers that are more familiar with events detection and root solving are @evan.ward and myself (and I designed the algorithm behind BracketingNthOrderBrentSolver
).
I would like to come back to your original problem. You say you encounter this exception quite often so there is definitely a problem there. Normally, as long as the function is really continuous in the neighbourhood of roots, the logic between EventState and the solver is very stable by now. With functions continuous around zero crossings, the solver should not be called with an interval where the zero is just before or after the boundaries, the solver should either have been called for the interval before or it will be called for the interval after, where the event really lies. We even handle properly functions that just touch zero without crossing it (and which don’t trigger events as there is no sign change). In many cases, problems in events detection were in fact problems in the function definition and coupling between the function definition and the dynamics changes triggered by the event. In most of these cases, the function was not really continuous and a zero crossing triggered by the before event logic changed the behaviour in the after event logic (say you switch off a device, or a force, or change an attitude) and as a side effect of this change the event function definition changed it sign again, hence breaking continuity. Are you sure this is not what happens in your case?