Double variable as loop counter

Question

I often have to code explicit schemes which means that I have to look at the evolution of a function by incrementing the time variable t <- t+ dt. It is therefore only natural to have my loops increment on dt:

int N = 7;
double t=0., T = 1., dt=T/N; // or I could have dt=0.03 for example
for(; t<T; t+= dt){
   if(T - t < dt){
      dt = T-t;
   }
   //some functions that use dt, t, T etc
} 

The rationale behind this is that I'm incrementing t by a constant dt at each step, except at the last iteration, where if my current time t is such that T- dt < t < T then I modify my time increment by dt <- T-t.

What are the possible pitfalls of such a procedure or ways I could improve it? I do realise that I might get a very small time increment.

Are there any floating problems that might appear (should I stick to incrementing on integers)?

In terms of optimisation, I assume that this technique is not costly at all, since a basic branch prediction would almost always skip the if block.

EDIT

I realise my question wasn't really good. Usually the dt is given by a CFL condition i.e. it is given so that it is small enough compared to some other parameters.

So from a logical point of view, dt is first given and afterwards we can define an integer N=floor(T/dt), loop with integers up to N, then deal with the leftover time interval N*dt --- T.

The code would be:

double dt = //given by some function;
double t=0., T = 1.;
for(; t<T; t+= dt){
   if(T - t < dt){
      dt = T-t;
   }
   //some functions that use dt, t, T etc
} 

Answer 1

First at all the compensation if (T - t < dt) is not needed, because it's only purpose appears to set the last value to t == T, which won't be processed due to inequality ...;t < T;... in the for loop condition.

That being said, finite difference method doesn't work that well with floats unless N is a power of two. If e.g. one wishes to evaluate a function at steps of 0.1f, one will most likely miss a few integral points.

Branch prediction may skip the condition evaluation, but there may also be a penalty/latency in mixing floating point operations with flow control operations.

Due to the cumulating rounding errors, it's possible that the iteration count can't be easily determined by the optimizer, disallowing some optimizations (loop unrolling or even vectorizing).

The inaccuracy can be mitigated simply by linear interpolation: t = c * dt;, but not perfectly, because not for all cases (dbl / N) * N == dbl. In practice the error should be in the epsilon magnitude. To get exact ending value, one has to calculate t = (range * N) / N; this time making sure range * N doesn't drop the least significant bits.

Answer 2

With the new information that dt must be set to a fixed, predetermined value (at least for all but the last step), here is my recommendation:

double T0 = 0.0;
double T  = 1.0;
int    N  = floor((T - T0)/dt);
double t  = T0;
for (int step_number = 0; step_number < N; ++step_number, t += dt)
{
  t = T0 + step_number * dt;
  do_one_step(t, T, dt);
}
if (t < T)
{
  do_one_step(t, T, T - t);
}

The function do_one_step performs the necessary calculations using t, T, and dt for each iteration. The data that must be updated by the function can either be made member variables of the class or can be included in the function parameter list as non-const references.

By the way, I have the last call to the function outside the loop not in order to save the possible cost of the branch condition, but because I find the code to be better organized and easier to understand that way.

---

The old answer:

You can easily get a very small time increment at the end, as you say, because the result T/N is typically not exact (and 0.03 certainly is not exact).

I would prefer to develop t and dt like this:

int    N  = 7;
double t  = 0.0;
double T  = 1.0;
double dt = (T - t)/N;
for (int step_number = 0; step_number < N; ++step_number, t += dt)
{
  // ... calculations with t, T, dt, etc.
}

(Notice that this says dt = (T - t)/N just in case you ever decide to start the iterations at a non-zero value of t.)

Alternatively, and potentially slightly more accurate if N is very large (because t += dt effectively has to round off dt as soon as t gets much larger):

int    N  = 7;
double T  = 1.0;
double dt = T/N;
for (int step_number = 0; step_number < N; ++step_number)
{
  double t = T0 + step_number * dt;
  // ... calculations with t, T, dt, etc.
}