When does an ODE integrator reach the best numerical accuracy?

I'm using an ODE Matlab integrator (Bulirsch-Stoer) that is (supposedly, I never used it before) quite good for solving ODEs with high accuracy. I'm using it to solve highly non-linear diff. eqs. of which one does not know the true solution and which are highly senstivie to initial conditions and numerical accuracy. To be sure that it does the right thing I integrate the curves back and forth and see whether the back-integration returns to the intial conditions. It does. Of course, if the time interval (the ODE's time-parameter "t") is large enough one sees it diverging (i.e., it doesn't return to the intital conditions.) This is normal, and expected. However, since the integrator's accuracy also depends on three independent parameters (the error tolerance, and two internal loops, midpoint and nr. of segmentation points) and it is difficult to find the optimal setting in a large 3D parameter space, I'm wondering whether the divergence from a true solution at some time, say t_max, beyond which the integrator no longer funrishes the correct values, is due to a non-optimal setting of the accuracy parameters (that is, I could still increase the accuracy), or wheteher it has reached the precision that is limited by the internal 15 digits number representation (that is, I reached the max. accuracy, and can't do anything about it as long I work with double precision.)

So, my question is: Is there a method to know when one has effectively reached the best numerical evaluation in solving an ODE inside the limitation of a double precision integratioon? A method that essentially tells me: "Yes, that's the best integration, beyond which you can't go with a 15 digits internal representation, no matter how you fine-tune your solver."

I hope that I expressed clearly what my issue is. Feel free to ask for more information.

[Steven Lord]

What is the problem you're trying to solve that requires a down-to-the-last bit solution? Are you solving an orbital mechanics problem over the course of years or decades? If we know the specific type of problem you're trying to solve, I'm wondering if there's a domain-specific tool or algorithm that can solve that problem more easily or effectively than a general purpose ODE solver. In a pinch you could use a wrench to hammer in a nail, but you're probably going to be more effective if you use a hammer instead.

If you require an exact solution you could always go to a symbolic solution (depending on whether or not your system of ODEs has an analytical solution) but that likely will take more time if it is possible.

Are you looking for a theoretical best or an actual best approach (based on the specific problem or type of problem you're trying to solve)?

In a comment on Torsten's answer you wrote "As to the irreversibility, I'm dealing with ODEs (y'=f(y,t)) having a unique solution y(t). Integrating it from y(0) to y(T) with T>0 must coincide with the same integration from y(T) to y(0)." Is that "must" in theory, in practice, or both? As the old saying says, "In theory there is no difference between theory and practice. In practice there is."

[Marco]

I'm trying to solve the equations of a damped forced pendulum. It is a non-linear chaotic system highly sensitive to intial conditions (which means, also to numerical accuracy). I'm looking for the exact solution on the longest possible time interval. As far as I undestand it, a very accurate integrator is the Bulrisch-Stoer one, and that I found here: https://www.mathworks.com/matlabcentral/fileexchange/55528-bulirsch-stoer Indeed, it works. However, it doesn't look like to be very accurate. I expected a better accuracy (it integrates for about 12-14 seconds and then diverges, slightly better than the Runge-Kutta, but not that much.) But, as I said, that might be a matter of internal precision. That's why I'm asking.

[Steven Lord]

Over what time scale are you hoping to integrate? If you're trying to model that system over the span of years or even days, my guess is you're going to be out of luck if the system is that sensitive.

[Marco]

In the literature I found people integrating up to 200 seconds.Would like to get there to.

Answer 1

[Torsten]

Solve the problem for different error tolerances and see when the solution stabilizes.

If available, also test different solvers for the problem.

That's all you can do.

Some problems are irreversable. So you shouldn't use the back-integration method you suggested as a test for a successful integration.

[Marco]

Ok, I should have clarified this better. Let me rephrase this in another way.

The solutions converge, this isn't the real question here. But, even if a solution converges, it always converges inside a certain time interval to the true solution. Beyond a certain integration time (number of iterations) the numerical error (due to the limited internal numerical precision) tends to propagate exponentially, and the numerical solution diverges from the real one. This will always happen, no matter how accurate one's integrator is.

As to the irreversibility, I'm dealing with ODEs (y'=f(y,t)) having a unique solution y(t). Integrating it from y(0) to y(T) with T>0 must coincide with the same integration from y(T) to y(0).

[Torsten]

Beyond a certain integration time (number of iterations) the numerical error (due to the limited internal numerical precision) tends to propagate exponentially, and the numerical solution diverges from the real one. This will always happen, no matter how accurate one's integrator is.

Your prescribed tolerances and the stepsizes computed by the solver from these tolerances should prevent this from happening - independent of the time interval over which you integrate. The error control of the solver is constructed to guarantee that the difference between true solution and numerical solution remains in prescribed bounds. If the solver is no longer able to satisfy these bounds, it will usually stop with the message that the minimum stepsize has been reached. If it doesn't do this, it's a bad solver.

[Marco]

Yes, but in double precision, the tolerance is a 15 digits number, it can't be much less than 1e-15/6, or alike, and all the estimations are limited to that. This limits the tolerance due to inherent internal precision and a numerical error grows exponentially, no matter how minute one sets it. Which means that you wil see the curve converge for some time t<T, but from t>T onward it will no longer converge to the true solution. This will always be the case for a sufficiently large T.

[Torsten]

What is the aim of your investigation in applying numerical methods to chaotic systems ?

[Marco]

To check how noise scales up in time.

[Torsten]

Hard to check if the noise is introduced by the solution method itself and you don't have a reference solution.