用MATLAB计算圆周率，结果不好

Question

我正在使用这个 .m 文件通过 MATLAB 计算圆周率

function calpi(n)
    S = 0;
    for i=1:n
        if mod(i,2) == 0
            S=S-1/(2*i-1)*(4*((1/5)^(2*i-1))-(1/239)^(2*i-1));
        else
            S=S+1/(2*i-1)*(4*((1/5)^(2*i-1))-(1/239)^(2*i-1));
        end    
    end
    S = 4*S;
    S=vpa(S,50)

n <= 8的时候就可以了。但是当n >= 9时，结果正好变成实际的圆周率。我只想得到实际结果来分析这个方法。

>> calpi(9)
S =
3.1415926535897932384626433832795028841971693993751

>> vpa(pi,50)
ans =
3.1415926535897932384626433832795028841971693993751

MATLAB 有什么问题？

Answer 1

首先：您需要在实际进行计算之前向 MATLAB 指定您需要可变精度。

您的函数 calpi 所做的是：使用内置 double 精度数据类型计算近似值，然后将其转换为符号 vpa 。如果你这样做，你无法获得比 double 更准确的近似值。

所以第一步是在计算之前使用vpa：

function S = calpi(n) S = vpa(0,50); for i = 1:n if mod(i,2) == 0 S=S-1/(2*i-1)*(4*((1/5)^(2*i-1))-(1/239)^(2*i-1)); else S=S+1/(2*i-1)*(4*((1/5)^(2*i-1))-(1/239)^(2*i-1)); end end S = 4*S;

如果您查看输出，您会发现仍然只返回双精度结果。这是因为当计算 1/(2*i-1)*... 完成时，它仍然仅使用双精度进行计算。您可以使用符号求值
来解决这个问题
vpa(subs('1/(2*i-1)*(4*((1/5)^(2*i-1))-(1/239)^(2*i-1))','i',i),50)

改为：

function S = calpi(n) S = vpa(0,50); for i = 1:n tmp = vpa(subs('1/(2*i-1)*(4*((1/5)^(2*i-1))-(1/239)^(2*i-1))','i',i),50); if mod(i,2) == 0 S = S-tmp; else S = S+tmp; end end S = 4*S;

至于为什么您的计算实际上 比您预期的更精确：当您使用普通 double 变量 x 调用 vpa(x,50) 时，会发生默认转换 sym(x,'r')。（参见 help sym）。

'r' stands for 'rational'. Floating point numbers obtained by evaluating expressions of the form p/q, p*pi/q, sqrt(p), 2^q and 10^q for modest sized integers p and q are converted to the corresponding symbolic form.

因此，当您调用 vpa(S,50) 时，MATLAB 会检查您的 double 值是否接近某个小数 p/q 或小数 p*pi/q 等，如果是，它会舍入对此的价值。在我们的例子中，如果 S 是 pi 的足够好的近似值，则 vpa(S,50) 向 sym('pi') 舍入。如果我们不想四舍五入，我们可以使用 vpa(sym(S,'f'), 50) 之类的东西来代替。

用MATLAB计算圆周率，结果不好

Calculate pi with MATLAB, bad results

matlab

symbolic-math

arbitrary-precision