Calculate pi with MATLAB, bad results

Question

I'm using this .m file to calculate pi with 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)

When n <= 8, it's alright. But when n >= 9, the result turns exactly into the actual pi. All I want is to get the actual result to analyse this method.

>> calpi(9)
S =
3.1415926535897932384626433832795028841971693993751

>> vpa(pi,50)
ans =
3.1415926535897932384626433832795028841971693993751

What's wrong with MATLAB?

Answer 1

First of all: You need to specify to MATLAB that you want variable precision before you actually do the calculation.

What your function calpi does is: Compute an approximation using the built-in double-precision datatype and then convert it to a symbolic vpa afterwards. You can't get an approximation that is more accurate than a double if you do it this way.

So a first step would be to use vpa before the calculation:

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;

If you have a look at the output, you will see that still only double-precision results are returned. This is because when the calculation 1/(2*i-1)*... is done, it is still evaluated using only double-precision. You can fix this by using the symbolic evaluation

vpa(subs('1/(2*i-1)*(4*((1/5)^(2*i-1))-(1/239)^(2*i-1))','i',i),50)

instead:

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;

---

As for why your calculation is actually more exact than you would expect: When you call vpa(x,50) with a normal double variable x, the default conversion sym(x,'r') happens. (See 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.

So when you call vpa(S,50), MATLAB checks if your double value is near some fraction p/q or a fraction p*pi/q, etc., and if it is, it rounds the value towards that. In our case if S is a sufficiently good approximation of pi, vpa(S,50) is rounded towards sym('pi'). If we didn't want to round, we could use something like vpa(sym(S,'f'), 50) instead.