Why the values are not correctly iterated using range arrays?

Question

I came across this problem

julia> Tinit = 0.0
    0.0  

julia> for ii in 1:10
           timeinterval = Tinit:0.05:(Tinit + 0.05)
           println(Tinit + 0.05)
           for item in timeinterval
                   println("+++++++++ ", item)
           end
           Tinit = timeinterval[end]
       end

And for some strange reason, I cannot obtain values bigger than 0.15

Answer 1

You're encountering a classic problem in floating-point arithmetic:

julia> 0.05 + 0.05
0.1

julia> 0.05 + 0.05 + 0.05
0.15000000000000002

julia> 0.05 + 0.05 + 0.05 + 0.05
0.2

So your first problem is that 0.1 + 0.5 == 0.15000000000000002 > 0.15. See Is floating point math broken?. So you're not getting the exact real number you were hoping to get. The second problem you're encountering is that you're constructing a range from that:

julia> r = 0.15000000000000002:0.05:0.2
0.15000000000000002:0.05:0.15000000000000002

julia> length(r)
1

julia> r[end]
0.15000000000000002

What?!? Why? This happens because your starting point is larger than you expect (your end point is also larger than you expect, but not significantly so: 0.2 == 2.00000000000000011102230246251565404236316680908203125). In particular, note that

julia> (0.2-0.15000000000000002)/0.05
0.9999999999999998

In other words, the range you're constructing isn't quite big enough to have two points. However, this is the case as well:

julia> 0.15000000000000002 + 0.05
0.2

Arguably, if you can add the step n times and hit the end point exactly, then the range should include that end point, so this could be considered a bug (issue filed: #20373). I'm fairly certain, however, that we could find other cases where your algorithm will accumulate enough floating-point error to fail.

There are a couple of possible solutions:

1) Don't do absolute floating-point calculations instead of iterative ones. For example, you could express your range as (0:0.05:0.05) + Tinit + (ii-1)*0.05:

julia> Tinit = 0.0; for ii in 1:10
           println(Tinit + (ii-1)*0.05 + (0:0.05:0.05))
       end
0.0:0.05:0.05
0.05:0.05:0.1
0.1:0.05:0.15
0.15000000000000002:0.05:0.2
0.2:0.05:0.25
0.25:0.05:0.3
0.30000000000000004:0.05:0.35000000000000003
0.35000000000000003:0.05:0.4
0.4:0.05:0.45
0.45:0.05:0.5

This doesn't entirely address the floating-point error problem, as you can see, but the error doesn't accumulate from one iteration to the next.

2) Use integer arithmetic, which is exact, and divide by 100 or whatever divisor to lower that computation down to the closest representable floating-point values. Since floating-point would only be used for presentation, not the core computation, the core computation would be correct and never accumulate error.

julia> Tinit = 0.0; for ii in 1:10
           println((Tinit+(ii-1)*5)/100:0.5:(Tinit+ii*5)/100)
       end
0.0:0.5:0.0
0.05:0.5:0.05
0.1:0.5:0.1
0.15:0.5:0.15
0.2:0.5:0.2
0.25:0.5:0.25
0.3:0.5:0.3
0.35:0.5:0.35
0.4:0.5:0.4
0.45:0.5:0.45

As you can see, this actually eliminates the inexact problem entirely.