Array elements not equal despite only difference being their respective end points

2014 7 17
2 个回答
回答已采纳
更新时间：2014 7 17
3 次查看（30 天）

0 个投票

I have two time arrays t1400 and t1500 defined as follows: The only difference between these two arrays is that t1500 has 100000 more cycles than t1400.
t1400=0:1e-3:1400;
t1500=0:1e-3:1500;
I would have assumed that the first 1400001 elements of both arrays would be absolutely equal but alas this is not the case as can be seen from the output of the following code:
equal=ones(length(t1400),1)*36;
for n=1:length(t1400)
equal(n)=(t1400(n)==t1500(n));
end
for n=1:length(equal)
if equal(n)==0
break;
end
end
n
n =
700003
The 700003rd being the first of many positions of the two arrays that are ever so slightly unequal and it's very slight:
>> format long
>> t1400(n)
ans =
7.000020000000000e+02
>> t1500(n)
ans =
7.000020000000001e+02
Unfortunately the system I am modelling is incredibly sensitive and by changing the end time of the time array can lead to very different results. I have accepted this as the case but I was just wondering whether someone could explain to me as to why Matlab gives different values depending upon the end time? Thanks in advance :)

请先登录，再回答此问题。

 采纳的回答

Peter
Peter 2014-7-17

1 个投票

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)
% COLONOP Demonstrate how the built-in a:d:b is constructed.
%
% v = colonop(a,b) constructs v = a:1:b.
% v = colonop(a,d,b) constructs v = a:d:b.
%
% v = a:d:b is not constructed using repeated addition. If the
% textual representation of d in the source code cannot be
% exactly represented in binary floating point, then repeated
% addition will appear to have accumlated roundoff error. In
% some cases, d may be so small that the floating point number
% nearest a+d is actually a. Here are two imporant examples.
%
% v = 1-eps : eps/4 : 1+eps is the nine floating point numbers
% closest to v = 1 + (-4:1:4)*eps/4. Since the spacing of the
% floating point numbers between 1-eps and 1 is eps/2 and the
% spacing between 1 and 1+eps is eps,
% v = [1-eps 1-eps 1-eps/2 1 1 1 1 1+eps 1+eps].
%
% 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.
%
% Ideally, in exact arithmetic, for b > a and d > 0,
% v = a:d:b should be the vector of length n+1 generated by
% v = a + (0:n)*d where n = floor((b-a)/d).
% In floating point arithmetic, the delicate computatations
% are the value of n, the value of the right hand end point,
% c = a+n*d, and symmetry about the mid-point.
if nargin < 3
b = d;
d = 1;
end
tol = 2.0*eps*max(abs(a),abs(b));
sig = sign(d);
% Exceptional cases.
if ~isfinite(a) || ~isfinite(d) || ~isfinite(b)
v = NaN;
return
elseif d == 0 || a < b && d < 0 || b < a && d > 0
% Result is empty.
v = zeros(1,0);
return
end
% n = number of intervals = length(v) - 1.
if a == floor(a) && d == 1
% Consecutive integers.
n = floor(b) - a;
elseif a == floor(a) && d == floor(d)
% Integers with spacing > 1.
q = floor(a/d);
r = a - q*d;
n = floor((b-r)/d) - q;
else
% General case.
n = round((b-a)/d);
if sig*(a+n*d - b) > tol
n = n - 1;
end
end
% c = right hand end point.
c = a + n*d;
if sig*(c-b) > -tol
c = b;
end
% v should be symmetric about the mid-point.
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

更多回答（1 个）

dpb
dpb 2014-7-17

1 个投票

The "why" of the difference here I can't explain, either. It seems the internals for colon are somehow being generated somewhat like the definition of linspace such that the actual delta step is being computed and that generates a rounding error that eventually accumulates.
Alternatively, it would seem the "deadahead" implementation would interpret the delta from the character constant value in the expression and that should, it would seem as you note, be invariant to the upper bound.
I'd submit this query to TMW Tech Support for an interpretation; it may not be considered a bug but it surely is peculiar. (Not that there's an eventual roundoff error accumulated; only that the two have different as your point is made).

4 个评论
显示 2更早的评论 隐藏 2更早的评论

Peter
Peter 2014-7-17
Thanks for looking at this for me. What do you mean by the "deadahead" implementation?
I think I will submit my query to Tech Support in a bid to get to the bottom of this.
Thanks again!
dpb
dpb 2014-7-17
What do you mean by the "deadahead" implementation?
The for loop with the delta computed from the literal constant converted to internal representation. What it seems must be happening is akin to precomputing the end length and then linspace(first,end,nPoints) wherein the delta would be computed and one could see a rounding difference there in that value whereas it's hard to grok why the literal constant difference doesn't lead to the same result up to the matching length as you noted.
Peter
Peter 2014-7-17
Thank you again for clarifying. Tech support got back to me with a very good and insightful answer which I have posted. Thanks for suggesting to get in touch with them.
Peter
dpb
dpb 2014-7-17
It helps when you have access to the code... :)
Makes sense when one knows how colon is actually implemented, indeed.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Logical 的更多信息

产品

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by