Accurate computation of the current time in time integrator

Question

I implement Runge--Kutta method for time integration of the system of nonlinear conservation laws

$$ u_t + f(u)_x = 0. $$

As the system is nonlinear, we have to recompute time step dt on each time step. Then we have to find new current time by using code like this:

cur_time = cur_time + dt

However, if we want to integrate to a very large time, then cur_time will be large while time step dt will be small, therefore, the error due to the floating-point arithmetics can accumulate significantly.

Is there a way to avoid this? If time step were constant, I could handle it, but what to do with the nonconstant time step?

UPDATE 2016-04-17. I add an example to demonstrate that the floating-point arithmetics introduce significant error even for a small number of step and small ratio cur_time/dt. For simplicity, I use constant time step dt here, nevertheless, the issue with nonconstant time step should be similar. Code is in Python programming language:

dt = 0.2
n = 100
cur_time = [0]
for i in range(n):
    cur_time.append(cur_time[-1] + dt)
cur_time[-5:]  # Shows last five time points.

which leads to

19.199999999999964,
19.399999999999963,
19.599999999999962,
19.79999999999996,
19.99999999999996

Therefore, even with small time of integration and small number of time steps, the time points are slightly corrupted by accumulation of error due to the floating-point arithmetics.

Answer 1

This shouldn't happen. You get floating point errors if dt is, say, $10^{-13}$ times cur_time. But this would mean that you have already done around $10^{13}$ time steps, which I think you are unlikely to have done, just because that would have taken a humongously long time.

In practice, the likelihood that you do so many time steps that round-off is an error is pretty small.

Answer 2

Your updated question implies that you need to know the current time very precisely -- so precisely that a difference on the order of $10^{-14}$ is unacceptable. If you need to track something to an accuracy better than a few times double-precision roundoff error, the only alternative is to use something other than double-precision numbers. You could use higher-precision floating-point numbers (see mpmath) or rational arithmetic (see sympy).

(If possible, you may want to cast $\Delta t$ to a double-precision float when using it elsewhere in order to avoid slowing down the whole code. But then when your time variable reaches $T$, your solution may as well be considered as a solution at $\hat{T} = T + \epsilon$ where $\epsilon$ is related to double-precision roundoff anyway).

Answer 3

Also don't get misleaded by the display of your computer regarding precision and remember that error does not stack the same way depending on your calculations for the same values.

Consider the following:

import numpy as np # Personnal comfort, not relevant here.
dt = 0.2
n  = 100
add_time, mul_time = np.zeros(n), np.zeros(n)
for i in range(n):
    if i == 0:
        add_time[i], mul_time[i] = 0.2
    else:
        add_time[i] = add_time[i-1] + dt
        mul_time[i] = (i-1)*dt
add_time[-5:], mul_time[-5]

Which leads to the following:

(array([ 19.2,  19.4,  19.6,  19.8,  20. ]), array([ 19.2,  19.4,  19.6,  19.8,  20. ]))

But for demonstration purpose I had the following print:

print "add_time: {0:18.18f}\n {1:18.18f}\n {2:18.18f}\n {3:18.18f}\n {4:18.18f}\n mul_time: {5:18.18f}\n {6:18.18f}\n {7:18.18f}\n {8:18.18f}\n {9:18.18f}\n".format(add_time[-5], add_time[-4], add_time[-3], add_time[-2], add_time[-1], mul_time[-5], mul_time[-4], mul_time[-3], mul_time[-2], mul_time[-1])
add_time: 19.199999999999963762
19.399999999999963052
19.599999999999962341
19.799999999999961631
19.999999999999960920
mul_time: 19.200000000000002842
19.400000000000002132
19.600000000000001421
19.800000000000000711
20.000000000000000000

My initial guess is also that this works well because dt is set at the value of 0.2 yet you'll always end up with a round-up error due to floats precision. But still, if you can store your values in advance with a minimum set of operation, as for time-step in this short example, do it. The lesser number of operations, the smaller the error.

EDIT: One can note the following on simple precision when substracting the expected value from the arrays for the printing:

add_time: -0.000000000000035527
0.000000000000035527
0.000000000000039080
0.000000000000039080
0.000000000000039080
mul_time: 0.000000000000003553
0.000000000000003553
0.000000000000000000
0.000000000000000000
0.0000000000000000000000000

We are at a precision setup where the errors are not even relevant to be printed. Hence, I actually advise to set a flag at each time step change and to calculate the new time step with something like, say: $$t_n = t_o + i*dt_n$$ where $t_o$ is the time at which time-step was changed, $dt_n$ the new time step and $i$ the number of iteration with this new value.