I received an answer from Tech Support which clarifies this issue:
This is a very interesting question. This is an expected behaviour, to understand why this happens we need to understand how the colon operator works.
v = a:d:b is not constructed using repeated addition (as one would think) as that might accumulate round-off errors.
In fact, this kind of round-off error is inherent to software or hardware systems that perform calculations using floating point numbers. Since floating point numbers have to be represented within a limited number of bits, the result of some calculations is rounded off in order to fit the result into the finite representation (sign, exponent, mantissa) for floating-point numbers.
An example of a number that cannot be represented exactly in binary floating point is the number 0.1, which, while finite in decimal representation, would require an infinite number of bits to be exactly represented in floating-point notation. For instance, executing the following in MATLAB:
0.1+0.1+0.1 == 0.3
will not return true due to such round-off errors.
To counteract such error accumulation, the algorithm of the COLON operator dictates that:
1. The first half of the output vector (in this example ‘v’) is calculated by adding integer multiples of the step ‘d’ to the left-hand endpoint ‘a’. 2. The second half is calculated by subtracting multiples of the step ‘d’ from the right-hand endpoint.
The COLON operator ensures symmetry of the output vector about its midpoint and also that the round off error is minimized at both endpoints. For example, even though 0.01 is not exactly represented in binary,
v = -1 : 0.01 : 1
consists of 201 floating points numbers centered symmetrically about zero.
For more detailed information on how the COLON operator works, refer to the attached function (COLONOP).
function v = colonop(a,d,b)
if nargin < 3
b = d;
d = 1;
end
tol = 2.0*eps*max(abs(a),abs(b));
sig = sign(d);
if ~isfinite(a) || ~isfinite(d) || ~isfinite(b)
v = NaN;
return
elseif d == 0 || a < b && d < 0 || b < a && d > 0
v = zeros(1,0);
return
end
if a == floor(a) && d == 1
n = floor(b) - a;
elseif a == floor(a) && d == floor(d)
q = floor(a/d);
r = a - q*d;
n = floor((b-r)/d) - q;
else
n = round((b-a)/d);
if sig*(a+n*d - b) > tol
n = n - 1;
end
end
c = a + n*d;
if sig*(c-b) > -tol
c = b;
end
v = zeros(1,n+1);
k = 0:floor(n/2);
v(1+k) = a + k*d;
v(n+1-k) = c - k*d;
if mod(n,2) == 0
v(n/2+1) = (a+c)/2;
end
Best Answer