How to get a more accurate cancelation

I shall try to get to the point, so let me know if there is something left and you need more details.

I am solving a couple of equations that are not coupled explicitly, but their corresponding unknown variables, say $$x$$ and $$y$$ must satisfy a differential equation:

$$\dot x = x + y,$$

where dots denote derivative with respect to an independent variable, say $$t$$.

The equation for $$x$$ is of second-order, so one gets $$x$$ and $$\dot x$$ from it, and one can check if the equation above is consistently satisfied. However (see the plot attached), it turns out that, no matter which integration method I use from SciPy (the ones that are already implemented), the equality above stops being satisfied at some point. This is due to the fact that $$x$$ and $$y$$ cancel each other up to a very high precision, which seems not to be attained by none of the methods provided by SciPy (I have checked this by taking every method and by lowering the absolute and relative tolerance as much as possible. In the plot attached, the method employed is DOP853, which is supposed to be very useful when very low tolerances are required).

My question is if you know of any way to improve the accuracy so the cancelation gets more precise (I would like that the equation will be conveniently satisfied all the way through the entire calculation). The only parameters I have changed so far were the relative and absolute tolerances (and of course the different methods at our disposal). Is there any parameter that I am missing and that might be useful for that?

• Unless you use something like GNU Multiple Precision library and increase the precision or rewrite the formula so it avoids subtraction, there are no parameters of ODE solvers that you can tune in this case. This looks like a textbook catastrophic cancellation to me, if it is not, it is another floating point arithmetic issue. Aug 20 '20 at 19:39
• To me, it looks like you've reached machine accuracy: about 1e-16 times the initial value. There is nothing you can do about that if you use double precision floating point accuracy, other than using something like multiple precision libraries (see above). The question that would be interesting to investigate, however, is why you care about such extreme accuracy? Aug 20 '20 at 23:31
• @WolfgangBangerth, if $x$ and $y$ are FP numbers of opposite sign and of order $10^{-12}$ and their mantissas agree except for the last place then it can be also catastrophic cancellation too, can't it? In that case, you can not really except to gain any more accuracy because the result of each subtraction will be garbage except for the first significant digit (e.g. up to %50 relatively). That could occur before machine accuracy is reached. Nevertheless, the solution is pretty much the same. I guess I am asking this out of curiosity, to get your opinion. Aug 21 '20 at 2:07
• If you get $x$ and $\dot{x}$ from the equation for $x$, why to solve the equation for $y$? Just use $y = \dot{x} -x$, that would satisfy the coupling constraint identically. Aug 21 '20 at 3:49

You can see this in action in the Feagin convergence plots which demonstrate accuracy of a 14th order method to $$10^{-50}$$ via BigFloats. In Julia with BigFloats or ArbFloats, you can do setprecision to then change the precision of the number types to get the accuracy you need.