MATLAB 的冒号运算符如何工作? [英] How does MATLAB's colon operator work?
问题描述
如
<小时>冒号:
function out = Colonop(start,step,stop)% COLONOP 演示如何构建内置 a:d:b.%% v = Colonop(a,b) 构造 v = a:1:b.% v = Colonop(a,d,b) 构造 v = a:d:b.%% v = a:d:b 不是使用重复加法构造的.如果% 源代码中 d 的文本表示不能% 以二进制浮点精确表示，然后重复% 加法似乎具有累积的舍入误差.在% 某些情况下，d 可能小到浮点数% 最接近 a+d 实际上是 a.这里有两个重要的例子.%% v = 1-eps : eps/4 : 1+eps 是九个浮点数% 最接近 v = 1 + (-4:1:4)*eps/4.由于间距1-eps 和 1 之间的 % 浮点数是 eps/2 并且1 和 1+eps 之间的 % 间距是 eps，% v = [1-eps 1-eps 1-eps/2 1 1 1 1 1+eps 1+eps].%% 即使 0.01 没有完全用二进制表示，% v = -1 : 0.01 : 1 由 201 个浮点数组成% 以零对称居中.%% 理想情况下，在精确算术中，对于 b >a 和 d >0,% v = a:d:b 应该是长度为 n+1 的向量，由% v = a + (0:n)*d 其中 n = floor((b-a)/d).% 在浮点运算中，精细的计算% 是 n 的值，右手端点的值，% c = a+n*d，并且关于中点对称.如果 nargin <3停止=步骤；步 = 1;结尾tol = 2.0*eps*max(abs(start),abs(stop));sig = 符号(步骤)；% 特殊情况.if ~isfinite(start) ||~isfinite(step) ||~无限(停止)出 = NaN;返回elseif 步骤 == 0 ||开始<停止&&步<0 ||停止 <开始&&步骤 >0% 结果为空.出 = 零(1,0)；返回结尾% n = 间隔数 = length(v) - 1.if start == floor(start) &&步骤 == 1% 连续整数.n = 楼层(停止)- 开始；elseif 开始 == 地板(开始)&&步骤==地板(步骤)% 带间距的整数 >1.q = 楼层(开始/步骤)；r = 开始 - q*步；n = floor((stop-r)/step) - q;别的% 一般情况.n = 轮((停止-开始)/步)；如果 sig*(start+n*step - stop) >托尔n = n - 1;结尾结尾% last = 右手边的终点.最后 = 开始 + n*步；如果 sig*(last-stop) >-tol最后=停止；结尾% out 应该关于中点对称.出 = 零(1，n+1)；k = 0:地板(n/2);out(1+k) = start + k*step;out(n+1-k) = last - k*step;如果 mod(n,2) == 0出(n/2+1)=(开始+最后)/2；结尾As noted in this answer by Sam Roberts and this other answer by gnovice, MATLAB's colon operator (start:step:stop) creates a vector of values in a different way that linspace does. In particular, Sam Roberts states:
The colon operator adds increments to the starting point, and subtracts decrements from the end point to reach a middle point. In this way, it ensures that the output vector is as symmetric as possible.
However, offical documentation about this from The MathWorks has been deleted from their site.
If Sam's description is correct, wouldn't the errors in the step sizes be symmetric?
>> step = 1/3;
>> C = 0:step:5;
>> diff(C) - step
ans =
1.0e-15 *
Columns 1 through 10
0 0 0.0555 -0.0555 -0.0555 0.1665 -0.2776 0.6106 -0.2776 0.1665
Columns 11 through 15
0.1665 -0.2776 -0.2776 0.6106 -0.2776
Interesting things to note about the colon operator:
Its values depend on its length:
>> step = 1/3; >> C = 0:step:5; >> X = 0:step:3; >> C(1:10) - X ans = 1.0e-15 * 0 0 0 0 0 -0.2220 0 -0.4441 0.4441 0It can generate repeated values if they are rounded:
>> E = 1-eps : eps/4 : 1+eps; >> E-1 ans = 1.0e-15 * -0.2220 -0.2220 -0.1110 0 0 0 0 0.2220 0.2220There is a tolerance for the last value, if the step size creates a value just above the end, this end value is still used:
>> A = 0 : step : 5-2*eps(5) A = Columns 1 through 10 0 0.3333 0.6667 1.0000 1.3333 1.6667 2.0000 2.3333 2.6667 3.0000 Columns 11 through 16 3.3333 3.6667 4.0000 4.3333 4.6667 5.0000 >> A(end) == 5 - 2*eps(5) ans = logical 1 >> step*15 - 5 ans = 0
The deleted page referred to by Sam's answer is still archived by the Way Back Machine. Luckily, even the attached M-file colonop is there too. And it seems that this function still matches what MATLAB does (I'm on R2017a):
>> all(0:step:5 == colonop(0,step,5))
ans =
logical
1
>> all(-pi:pi/21:pi == colonop(-pi,pi/21,pi))
ans =
logical
1
I'll replicate here what the function does for the general case (there are some shortcuts for generating integer vectors and handling special cases). I'm replacing the function's variable names with more meaningful ones. The inputs are start, step and stop.
First it computes how many steps there are in between start and stop. If the last step exceeds stop by more than a tolerance, it is not taken:
n = round((stop-start)/step);
tol = 2.0*eps*max(abs(start),abs(stop));
sig = sign(step);
if sig*(start+n*step - stop) > tol
n = n - 1;
end
This explains the last observation mentioned in the question.
Next, it computes the value of the last element, and makes sure that it does not exceed the stop value, even if it allowed to go past it in the previous computation.
last = start + n*step;
if sig*(last-stop) > -tol
last = stop;
end
This is why the lasat value in the vector A in the question actually has the stop value as the last value.
Next, it computes the output array in two parts, as advertised: the left and right halves of the array are filled independently:
out = zeros(1,n+1);
k = 0:floor(n/2);
out(1+k) = start + k*step;
out(n+1-k) = last - k*step;
Note that they are not filled by incrementing, but by computing an integer array and multiplying it by the step size, just like linspace does. This exaplains the observation about array E in the question. The difference is that the right half of the array is filled by subtracting those values from the last value.
As a final step, for odd-sized arrays, the middle value is computed separately to ensure it lies exactly half-way the two end points:
if mod(n,2) == 0
out(n/2+1) = (start+last)/2;
end
The full function colonop is copied at the bottom.
Note that filling the left and right side of the array separately does not mean that the errors in step sizes should be perfectly symmetric. These errors are given by roundoff errors. But it does make a difference where the stop point is not reached exactly by the step size, as in the case of array A in the question. In this case, the slightly shorter step size is taken in the middle of the array, rather than at the end:
>> step=1/3;
>> A = 0 : step : 5-2*eps(5);
>> A/step-(0:15)
ans =
1.0e-14 *
Columns 1 through 10
0 0 0 0 0 0 0 -0.0888 -0.4441 -0.5329
Columns 11 through 16
-0.3553 -0.3553 -0.5329 -0.5329 -0.3553 -0.5329
But even in the case where the stop point is reached exactly, some additional error accumulates in the middle. Take for example the array C in the question. This error accumulation does not happen with linspace:
C = 0:1/3:5;
lims = eps(C);
subplot(2,1,1)
plot(diff(C)-1/3,'o-')
hold on
plot(lims,'k:')
plot(-lims,'k:')
plot([1,15],[0,0],'k:')
ylabel('error')
title('0:1/3:5')
L = linspace(0,5,16);
subplot(2,1,2)
plot(diff(L)-1/3,'x-')
hold on
plot(lims,'k:')
plot(-lims,'k:')
plot([1,15],[0,0],'k:')
title('linspace(0,5,16)')
ylabel('error')
colonop:
function out = colonop(start,step,stop)
% 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
stop = step;
step = 1;
end
tol = 2.0*eps*max(abs(start),abs(stop));
sig = sign(step);
% Exceptional cases.
if ~isfinite(start) || ~isfinite(step) || ~isfinite(stop)
out = NaN;
return
elseif step == 0 || start < stop && step < 0 || stop < start && step > 0
% Result is empty.
out = zeros(1,0);
return
end
% n = number of intervals = length(v) - 1.
if start == floor(start) && step == 1
% Consecutive integers.
n = floor(stop) - start;
elseif start == floor(start) && step == floor(step)
% Integers with spacing > 1.
q = floor(start/step);
r = start - q*step;
n = floor((stop-r)/step) - q;
else
% General case.
n = round((stop-start)/step);
if sig*(start+n*step - stop) > tol
n = n - 1;
end
end
% last = right hand end point.
last = start + n*step;
if sig*(last-stop) > -tol
last = stop;
end
% out should be symmetric about the mid-point.
out = zeros(1,n+1);
k = 0:floor(n/2);
out(1+k) = start + k*step;
out(n+1-k) = last - k*step;
if mod(n,2) == 0
out(n/2+1) = (start+last)/2;
end
这篇关于MATLAB 的冒号运算符如何工作?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
