Machine precision and local error

I'm working with an RKF45 integrator that I have programmed using CUDA C++ on my GPU and am pondering a few questions as I'm trying to track down some issues with my code.

1. I'm using double precision values in the calculation. Is it feasible of me to demand of my integrator that the difference between the 4th and 5th order estimates be < 1 x 10 ^-16 if machine precision is that of a double?

2. If I do impose local error to be 10E-16 and I integrate over, for example, 20,000 time steps, should I then be expecting a maximal global error of the product of these two numbers : 2E-12?

Indeed if that is the case, that's not too bad, however I'm trying to track the source of a leak where my numbers become unrealistic and I'm trying to rule out machine precision as the source of that error.

EDIT

My RHS of my equations look like this.

double frag_term = 0;

double flux = 0;
if (r == ((maxlength)-1))
flux = -km*(r)*conc[r]+2*(ka)*conc[r-1]*conc[0];

else if ( (r > ((nc)-1)) && (r != ((maxlength)-1)) )
{
frag_term = conc2[maxlength-1] - conc2[r];

flux = -(km)*(r)*conc[r] + 2*(km)*frag_term - 2*(ka)*conc[r]*conc[0] + 2*(ka)*conc[r-1]*conc[0];
}
else if (r == ((nc)-1))
{
frag_term = conc2[maxlength-1] - conc2[r];

flux = (kn)*pow(conc[0],(nc)) + 2*(km)*frag_term - 2*(ka)*conc[r]*conc[0];
}
else if (r < ((nc)-1))
flux = 0;


where conc is an array of concentrations whose values I am integrating over time, conc2 is an array of partial sums such that conc2[r] = sum of all conc array items less than and including the rth term.

• Have you compared your results to matlab's ODE45 or other established implementations?
– Jan
Commented Aug 22, 2014 at 8:13
• I have not, yet. I've had a hard enough time getting to know C++, haha. I've been trying to devote more time to physics rather than all the computational aspects as my advisor isn't too keen on the how of producing the simulation, only the results. Commented Aug 22, 2014 at 18:51
• I see... Nevertheless, you may want to check or use existing implementations. See Ernest Hairar's website for a list of (probably) high performant code. Most is in Fortran but there are also C++ versions.
– Jan
Commented Aug 23, 2014 at 14:31

I'm using double precision values in the calculation. Is it feasible of me to demand of my integrator that the difference between the 4th and 5th order estimates be < 1 x 10 ^-16 if machine precision is that of a double?

Assuming you're talking about absolute error tolerances, in general, yes, it makes sense.. A unit in the last place (ulp) would be the 16th decimal place in scientific notation, which does not necessarily correspond to $10^{-16}$. Suppose you have a scalar ODE, and for a given time step, the 5th order method calculates a value of $1.01 \times 10^{-15}$, while the 4th order method calculates a value of $1 \times 10^{-15}$. The estimated local truncation error would be $1 \times 10^{-17}$, without loss of precision due to floating point arithmetic, so you could use an absolute tolerance of $5 \times 10^{-17}$ for this example time step without having to resort to reducing the size of the time step.

If I do impose local error to be 10E-16 and I integrate over, for example, 20,000 time steps, should I then be expecting a maximal global error of the product of these two numbers : 2E-12?

That all depends on how the summation of time steps is implemented. For naive summation, I'd expect the rounding errors to grow roughly linearly with the number of time steps. Compensated summation methods can be used to control the rounding error well when a large number of time steps are taken. (The link references Accuracy and Stability of Numerical Algorithms, 2nd edition by Higham, page 87; see that book for more details regarding how rounding errors propagate for a variety of different numerical methods.)

I wouldn't necessarily expect the global error to grow linearly, however, purely on the basis of local truncation error.

• OK, I think I see what you're saying. My IC's for my initial concentrations look like [5E-6, 0, 0...0] and the RKF45 integrator subtracts the 4th and 5th order calculations for all the elements, picks the max and uses that as a worst case error. If I'm understanding you, then I wouldn't have to worry about floating point arithmetic errors unless I set my local error tolerance to ~ 1E-22 or thereabouts since the first five zeroes are not contributing to the significance of the calculation? Commented Aug 22, 2014 at 18:44
• For time steps near or shortly after your initial condition, yes. Roughly speaking, given a time step, the maximum value of the state variables in the 5th order predictor times $10^{-16}$ is going to be the minimum realizable value of the error greater than zero. Smaller values of the error can't occur, because you couldn't represent them using floating-point arithmetic. So if your state variables decay with time to smaller values, an absolute error tolerance of $10^{-25}$ might be reasonable; if they increase with time, larger error tolerances might be appropriate. Commented Aug 22, 2014 at 19:24
• I think I am experiencing this as I try to test my integrator. I changed all my RHS to be conc[r], with the expectation that my calculations should return e^t. Stepping through it with the debugger, I was finding that while my concentration was 5E-6, I was calculating "zero" local error, which I assume comes from the fact that e^t is relatively flat and unchanging at that point. After letting it run up to a concentration in the realm of ~1, I noticed that the aggregated time was off and that there was difficulty in the stepper permitting a legal step of <10^-16 error. Commented Aug 22, 2014 at 19:43
• Generally, you want to set the absolute tolerance to be the largest value you'd be willing to call "zero" when comparing solutions. You might also consider using a relative error tolerance instead of an absolute error tolerance; relative error tolerances tend to work better with values close to zero. Commented Aug 22, 2014 at 20:19
1. As I see it, it makes sense to compare on the level of machine precision. In the worst case you will not be able to satisfy your criterion because of rounding errors, but this is a Type I error (the time step is rejected although it meets the requirement) which is 'safe'.
2. No, this depends on how much an error is amplified by your incremental function. This quantity depends mainly on the Lipshitz constant of your right hand side.

Finally, since you are going for high precision, keep in mind, that every single-step method suffers from an inherent instability which is $\mathcal O(\frac{\epsilon}{\tau})$, i.e., an error that scales with machine precision divided by the time-step size. This prevents you from reaching very high precision and, for very small time steps, may destroy any convergence.