3
$\begingroup$

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).

enter image description here

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?

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ 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. $\endgroup$
    – Abdullah Ali Sivas
    Aug 20, 2020 at 19:39
  • $\begingroup$ 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? $\endgroup$
    – Wolfgang Bangerth
    Aug 20, 2020 at 23:31
  • $\begingroup$ @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. $\endgroup$
    – Abdullah Ali Sivas
    Aug 21, 2020 at 2:07
  • 1
    $\begingroup$ 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. $\endgroup$
    – Maxim Umansky
    Aug 21, 2020 at 3:49

1 Answer 1

Reset to default
4
$\begingroup$

I am not sure this is possible with the Python libraries since they are using Fortran under the hood and that can't be easily recompiled, but the Julia DifferentialEquations.jl JIT compile specializes the solvers based on the number types that you give it. Here's a demonstration of some weird types like rational numbers, MPFR BigFloats, and ArbFloats (based on the Arb library).

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.

While the Julia methods are very fast in comparison to SciPy (order and a half of magnitude), and even though they compile to specialize on input types so a special optimize code is generated for the high precision case, high precision arithmetic is still fairly expensive and you should keep this in mind. Specializing the integrator for this high accuracy range will be fairly important. Note that if you want to do this, I would probably recommend Vern9 or one of the multithreaded extrapolation methods like ExtrapolationMidpointDeuflhard (they will multithread between f calls which will be more important as tolerance decreases), or maybe the new 16th order symplectic integrator IRKGL16.

Also, if you need validated arithmetic, you can use TaylorIntegration.jl for high order Taylor methods with floating point accuracy bounds on the solution.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.