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


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.

  • $\begingroup$ Depending on what language you're working in, it may not be too difficult to use higher precision for the time value. $\endgroup$ – David Ketcheson Apr 15 '16 at 5:48
  • $\begingroup$ You could always try to keep track of time using a fraction representation using unsigned integers. I have this in my simulation setup, and it's useful for accurate time tracking. $\endgroup$ – spektr Apr 17 '16 at 0:51
  • $\begingroup$ Can you explain why an error of a few times machine epsilon is problematic for your computations? $\endgroup$ – David Ketcheson Apr 17 '16 at 7:01
  • $\begingroup$ @DavidKetcheson, here I asking about this mainly because of the curiosity. Previously, I was checking linear solver (with constant time step) for self-convergence (by taking twice smaller time steps and comparing norm of error between solutions on corresponding time layers). $\endgroup$ – Dmitry Kabanov Apr 17 '16 at 20:16
  • $\begingroup$ Continuation: There, it was crucial to compute current time like this: cur_time = time_step * dt. This way, time layers match each other: that is, each time layer of the solution with time step 2dt matches every second time layer of the solution with time step dt to 16 significant digits. Then, I could be sure that I compute error correctly. That way I became aware of the accumulation error in the computation of current time (before I didn't think about it). Now, I'm curious how it affects the computation in the case of nonconstant time step. $\endgroup$ – Dmitry Kabanov Apr 17 '16 at 20:16

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.

| cite | improve this answer | |
  • 4
    $\begingroup$ I suppose if he's already run into problems with this, he could try something like Kahan summation to accumulate a more accurate final time. However, I agree with you that he should be more concerned about why his time step was so small. $\endgroup$ – Tyler Olsen Apr 13 '16 at 19:25
  • $\begingroup$ But to do so, you have to store all time steps and that would be a very large number. $\endgroup$ – Wolfgang Bangerth Apr 13 '16 at 21:30
  • 4
    $\begingroup$ You don't need to store all time steps to use Kahan summation, which is a form of running compensation. For something like pairwise summation you would. $\endgroup$ – horchler Apr 13 '16 at 22:19
  • 1
    $\begingroup$ I seem to recall running into this problem in a significant way when doing ~$10^9$ time steps. I guess it depends on how accurate you need your current time to be. $\endgroup$ – David Ketcheson Apr 14 '16 at 8:58
  • 1
    $\begingroup$ @DavidKetcheson, and how big a difference between your largest and smallest so far has been. If you have a tiny time step to start, an adaptive time stepper, and a section of large steps, you may get the current total very far from the next step size in a small $O(1)$ number of steps. $\endgroup$ – Bill Barth Apr 14 '16 at 16:29

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

| cite | improve this answer | |
  • $\begingroup$ I prefer to stay with double-precision arithmetics. Thanks for the mpmath link. I didn't know about it. $\endgroup$ – Dmitry Kabanov Apr 18 '16 at 9:42

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
        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
mul_time: 19.200000000000002842

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
mul_time: 0.000000000000003553

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.

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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