为什么在Matlab中sin(pi)不是精确的，但sin(pi/2)是精确的？

Question

我在使用 matlab 进行计算时遇到了问题。 我知道"pi"是一个浮点数，而不是精确的数字。因此，在matlab中，sin(pi)并不完全等于零。 我的问题是，如果"pi"不是精确的数字，那么为什么sin(pi/2)恰好等于1。

sin(pi)--> 由于pi是浮点数，所以不精确。 但是， sin(pi/2) 恰好等于1。

我感到很惊讶和困惑！

- Nafiseh Eftekharian

---

它不完全等于1，但误差实际上太小以至于不能被表示。 - Paul Stelian

但是Matlab显示确切的1。 - Nafiseh Eftekharian

尝试运行sin(sym('pi'))并访问http://www2.math.umd.edu/~jmr/241/introsession.html。 - Malick

坦白说，这可能只是运气！只要结果在eps范围内，你应该会很高兴 :) - gregswiss

嗯，运气？我想知道是否有人能提供一个好的答案，其中使用了pi / 2的双精度算术表示和计算sine所使用的精确算法，以证明这种情况在浮点精度范围内收敛到1.0。 - Łukasz Rogalski

也许它是1.00000000000000000000000000000143..（而浮点格式只能容纳前10-15个小数位，类似于）。这真的是关键。 - Paul Stelian

Answer 1

我不知道Matlab计算sin(x)的确切方法，但你可以通过使用幂级数来进行探究，即

sin x = x - (x^3)/3! + (x^5)/5! - (x^7)/7! + (x^9)/9! ...

将此转换为Matlab代码，我们表示为：

clc
x = pi;     %  or x = pi/2
res = x;
factor = -1;
for ii=3:2:19
  res = res + factor*power(x,ii)/factorial(ii);
  factor = factor*-1;  
  fprintf ( 'iteration %2i  sin(x)=%1.16f\n', (ii-1)/2, res );
end
res

运行这段代码，对于x=pi和x=pi/2两种情况，你会发现x=pi/2很快（9次迭代）就收敛到了正确的结果（误差在eps之内），而x=pi的情况则没有在同一时间范围内收敛。

需要注意的是，在第9次迭代时，正在计算的是factorial(19)。在这个序列中下一个要计算的阶乘是21。由于双精度的限制（参见help factorial），这是最后一个可以用100%准确表示的阶乘。

因此，我认为发生的情况是，对于pi/2，数学解法比pi更快地收敛于1，以双精度结果存储的数学和精度限制导致pi无法完全收敛。

尽管如此，sin(pi)的误差在eps之内，所以你应该利用这个事实。

以下是我得到的结果（R2015b）：

Results for PI/2
iteration  1  sin(x)=0.9248322292886504
iteration  2  sin(x)=1.0045248555348174
iteration  3  sin(x)=0.9998431013994987
iteration  4  sin(x)=1.0000035425842861
iteration  5  sin(x)=0.9999999437410510
iteration  6  sin(x)=1.0000000006627803
iteration  7  sin(x)=0.9999999999939768
iteration  8  sin(x)=1.0000000000000437
iteration  9  sin(x)=1.0000000000000000
Final Result: 1.0000000000000000


Results for PI
iteration  1  sin(x)=-2.0261201264601763
iteration  2  sin(x)=0.5240439134171688
iteration  3  sin(x)=-0.0752206159036231
iteration  4  sin(x)=0.0069252707075051
iteration  5  sin(x)=-0.0004451602382092
iteration  6  sin(x)=0.0000211425675584
iteration  7  sin(x)=-0.0000007727858894
iteration  8  sin(x)=0.0000000224195107
iteration  9  sin(x)=-0.0000000005289183
Final Result: -0.0000000005289183

- matlabgui

---

包括使用 pi 和 pi/2 的标准输出会很好。MATLAB 内部使用哪个算法来计算正弦函数是否有文档记录？ - Łukasz Rogalski

完成。我不知道Matlab内部使用哪种算法。它可能会根据所请求的值和平台等使用不同的算法组合...（请参见c++中类似问题的此答案）。 - matlabgui

Answer 2

原因是 sin(pi)=0.0，所以即使是微小的误差也与0相比非常大，因此是可见的。而对于 sin(pi/2)=1，如果算法产生的误差小于eps（约为2.220446e-16），你将看不到这个误差，因为1+eps=1。

误差部分是由于不精确的输入（pi值不精确）和计算过程中的舍入误差造成的。必须深入代码才能正确地处理它。

另一个重要因素是函数本身。通过观察pi和pi/2的泰勒级数，可以看到误差传播的影响。

sin(pi+dx)=sin(pi)+cos(pi)dx+o(dx^2)=-dx+o(dx^2)
sin(pi/2+dx)=sin(pi/2)+cos(pi/2)dx+o(dx^2)=1+o(dx^2)

很明显：如果dx与eps相近，那么由于不精确输入而导致的误差将大约为eps*eps，因此与1相比几乎不可见。

- ead