2

I have an RKF7(8) integrator whose output I've verified with several simple test functions. However, when I use it on the equation of motion I'm interested in, the local truncation errors are suddenly very small. For a timestep of around 1e-1, my errors are all around 1e-18 or 1e-19. For the simple test functions (so far, sine and an exponential), the errors are always reasonable, ie 1e-7 or so for the same timestep.

The only difference between the simple test functions and the problem function is that it's a huge chunk of code, like maybe 1000 different terms, some with relatively large exponents (like 9 or 10). Could that possibly affect the precision? Should I change my code to use long doubles?

14
  • 1
    How do the errors behave when you change the step size by a factor of 2, h=0.05, h=0.2 etc? May 27, 2020 at 16:56
  • Did you have also some vector cases among the test problems, for instance the simple circular motion or some orbit simulation for circular orbits? May 27, 2020 at 17:11
  • Interesting. When I set the step size to 0.5, 0.25, and 0.2, the error quickly goes to 1e-11, which is still low, but the step size rapidly becomes >1. The sin test problem was technically a vector case, as I did it as two coupled 1st-order equations. Should I do a more complex case, like an orbit?
    – David
    May 27, 2020 at 17:27
  • Ok, that seems right then for the method implementation. Looking for existing implementations I found one where the indentation in the loop over the stages was wrong, which might lead to errors that can only be seen in ODE systems. // What are the scales of the components that show this error, that is, what is the relative error? 1e-18 makes only sense in floating point if the component is of size 1e-4 or so. May 27, 2020 at 17:36
  • With a method of order 7, the optimal step size is at Lh=(2e-16)^(1/8)~0.01. If the Lipschitz constant L of the system is small, then h=0.1 can be in the region where the accumulation of floating point errors dominates the step error. In this range the error will be rather random, so large jumps like this would still be exceptional, but not unexpected. May 27, 2020 at 17:36

1 Answer 1

0

Very interesting question. The problem you are facing might be related to issues (or limitations) of the floating point arithmetic. Since your function contains coefficients in a wide numerical interval, it is likely that you have some loss of precision in your calculations. In general, these problems can come in the form of:

  • Overflow
  • Underflow
  • Multiplication and division
  • Adding numbers of very different magnitudes
  • Subtracting numbers of similar magnitudes

Overflow and underflow occur when the numbers you are dealing with are too large or too small with respect to the machine precision, and it would be my bet that this is not what happens in your system. Nevertheless, one must take into account that multiplication and division operations can lead to overflow and underflow. On the other hand, adding numbers of very different magnitudes (or subtracting numbers of similar magnitudes) can lead to severe loss of precision due to the roundoff errors. From my experience in optimization problems that involve large and small numbers, I would say this could be a reasonable explanation of the behavior of your integrator.

I have two suggestions for you. The fist one is of course increasing the precision of your numbers to the maximum available. This might help or not depending on how ill conditioned your problem is. The second one is to use a better algorithm to perform the sums in your numerical method. In contrast to the naive addition of all number sequentially, you could use a more elaborated strategy by dividing your sums into sub-sums, effectively reducing roundoff errors. Notable examples of these algorithms are the pairwise summation and the Kahan summation.

I hope this answer offers you some clues. Good luck!

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.