Root function where g(x) = g(x + d) = 0

[guibou]

Hello.

General problem

In the sundial documentation, at least for IDA, it is noted that the root finding logic will generate an error if after finding a root such as g(x) = 0, sundials will try g(x + d) (where d is small) and if it is equal to 0 too, it will fail to ensure that at any point in time we know from which side of the root we are.

The documentation for IDASolve, https://sundials.readthedocs.io/en/latest/ida/Usage/index.html#c.IDASolve, states that it will return IDA_ILL_INPUT in this context.

What is the recommended way to handle such a case?

My context

I'm solving arbitrary problems: users provides a set of ode / algebraic rules, and my program use sundial to solve that.

The ode provided by the user may contain discontinuities (e.g. dx / dt = if g > 20 then y else z). In order to not miss some possibly important discontinuities, we use root finding in order to detect the discontinuity, stop the solver, change the relevant values, restart. Without this, we can miss important / short phenomenums. For example if we have dx/dt = if t > 20 and t < 20.1 then 10000 else 0 we expect an important change of x, however a step larger than 0.1 may miss it.

The root function are not provided by the user, we extract them from the different equations of the problem. For the two previous examples, we introduce roots functions r such as r_0 = g - 20 = 0, r_1 = t - 20 = 0 and r_2 = t - 20.1 = 0 as root functions and stop the solver accordingly.

Unfortunately, we can generate functions for which r(x) = r(x + d) = 0, this is a bit out of our hand (e.g. because it depends on user input) and until now I had not been able to find a robust solution.

Robust solution

The only idea I do have right now is to store in user data a context about the root value and sign and when the root function is called and a 0 is found, lookup in the context to see if a 0 was found previously at a close distance, and if yes, instead of returning a 0, return a small value of a different sign than the value observed in previous calls on the left of the current timepoint. This is impressively ugly.

Not robust solution, detecting IDASolve fail and ignoring the problem

So right now my only solution is to handle the problem when it arise, e.g. when IDASolve fails with IDA_ILL_INPUT.

I hence do have multiples questions:

a) When this happen, how do I detect that I'm in this context? I can get the error message and hope it matchs Root found at and very near, but that's not robust.
b) Am I guaranteed to be stopped at the root location? The documentation states that the value of tret is the farthest point reached during the integration.
c) How do I detect which condition is 0? IDAGetRootInfo tells me which function had changed sign or are ==0, but how can I detect which one have triggered the error? Is the result of IDAGetRootInfo even guranteed to be relevant?
d) Once I detected that I may be in a weird state, I may decide to restart the solver. If I restart it at the current time, I suppose it will fail (because the root function is =0, so it does not know how to initialize it). A hack could be to advance by a small amount and hope for the best. I may also disable the root function which are problematic, but I would want to reenable them as soon as possible.

Of course I can experiment with sundials to get most of the answer to these questions, but it will be implementation based and may not be robust, hence why I asking theses questions here.

In summary, do you have any suggestion to handle this very specific case?

Thank you.

[guibou]

Ping @eltix (because you work in the same company, so you may be interested ;)

[guibou]

For the record, I think I found a way to avoid the 0 value while keeping the same root properties.

The idea is that we could introduce:

g'(x) = if
   g(x) != 0 then g(x)
   else small_epsilon # See note about the sign

Details need to be checked for numerical accuracy, but the idea is that we map g(x) to a function which is exactly the same BUT for the 0 value which is remove and replaced by the smallest non zero value.

Note that it means that the function is never 0, which does seem to contradicte the idea of "root finding". However, for sundials, we only care about a change of sign and this function behaves exactly like the previous one with respect to the sign except on 0 where anyway, the notion of sign is open to debate.

Negative/positive zero

A small digression about the sign. In IEEE float, the zero do have a sign and actually I'm wondering if this property could not be used to remove the special case described here. If we take the sign of the 0 into account, our equation always have an explicit sign change and hence being stuck in a context where we don't know the sign of the previous timepoint is now impossible. However it would render the root detection of function which does not change sign (but do have a perfect 0) a bit more unreliable (e.g. a positive function with a 0 value won't be considered as sign changing, when a negative function where the computed 0 is not strictly a negative 0 would see its sign change).

Small epsilon

The sign of the small epsilon is important to ensure that the root finding algorithm actually triggers immediately and not after a small delay.

For example, if you want to trigger an event when x > 10. The solution is to look for the root of x - 10. This function is negative when x < 10, positive when x > 10 and ==0 when x==10.

Let's assume that the epsilon is positive:

If we start with x < 10, the function is hence negative. When x == 10, x - 10 gives 0 + epsilon which is >0 hence the root is detected exactly at the right time where x == 10, however we were looking for x strictly above 10.

Let's assume that the epsilon is negative:

If we start with x < 10, the function is hence negative. When x == 10, x-10 gives 0 - epsilon which is <0, hence the root is not detected when x == 10. The solver may go a step further when x = 10.000000001, where x - 10 is now positive, detecting the root.

On most context, the different is irrelevant. But it can be of importance if the root function stays a long time on 0. Actually, it also completly change how the solver will behave when trying to locate the root of function which have root, but do not change sign. For example, if you track the root of x^2, the function never change sign anyway. However, on x==0, if the small epsilon is positive, we will never detect the root. If the small epsilon is negative, it could be detected (providing that x stays ==0 for long enough so a solver step won't jump above it)