compute the approximation of pi

Write a program that approximates PI by computing the sum

.

The more terms you keep in the summation, the more accurate your answer will be. (In fact, the series converges to PI as m goes to infinity.) See how many terms you need to approximate PI with 5 decimals. (Note: This is by no means the most efficient way to approximate PI, but the formula is quite beautiful...)

Where have I gone wrong in this code, I do not get anywhere near the approcimation of pi?

for k = 0:25
    
    (-1)^(-k)/(2*k +1);
    sum = 0;
    
    for m = 0:10
        
        sum = 4*(sum + (-1)^(-k)/(2*k +1));
    end
    
    disp(sum);
end

Answer 1

[Steven Lord]

Multiply by 4 after summing the terms not while you're summing the terms. Otherwise the first term will be multiplied by 4 many times. Also you don't need a nested pair of loops; you're summing over k (with an upper limit of m = 25) so you don't need to loop over m. If you wanted to show the effect of increasing the upper limit of summation, you could make a vector in which you store the current result for each value of k.

Finally, I'd avoid calling your variable sum. While that variable exists you will not be able to call this somewhat useful function.

[Bradley Johnson]

I changed it to this, is this what you were saying to do.My problem is that i can not even come close to approximating pi which is the end goal.

for k = 0:20
    
    (-1)^(k)/(2*k +1);
    A = 0;
    A = (A + (-1)^(k)/(2*k +1));
    B = 4 * A
   
    
    disp(B);
end

[Stephen23]

"...is this what you were saying to do."

No.

B = 4 * A
%   ^ This must not be INSIDE the loop. It must be AFTER the loop.

Compare what your code does with the formula you were given: you multiply by four on every loop iteration... whereas the formula only multiplies by four once after the summation. You should do the same.

[Steven Lord]

You're getting closer, but you don't want to reset A to 0 each time through the loop. Initialize it to 0 before the loop, add each term to it inside the loop, and multiply by 4 after the loop as Stephen Cobeldick said.

If I gave you a list of numbers to add together, one way to do that is to (implicitly) start a running total at 0, add each number in the list to the running total, and tell me the total once you've finished the list. This is the same idea.

[Bradley Johnson]

I got it, Thank you

Answer 2

[John D'Errico]

Now that Steven has had his answer accepted, and you should have finished your work, I'll add an answer to expand on the ideas in such an approximation for pi. First, using a while loop, to see how far we need to go.

One trick in such a series is to recognize that you will need to sum terms for at roughly as long as how long it takes for any term in the series to be smaller than the tolerance. Given a tolerance of 1e-5, I would expect to need at least 4/1e-5 terms, so perhaps 25000 terms or more.

% I'll multiply by 4 at the end, so the tolerance should be decreased
% by a factor of 4 in the while loop.
pitol = 1e-5/4;
pisum = 0;
pierr = inf;
S = 1;
K = -1;
while pierr > pitol
  K = K + 1;
  pisum = pisum + S/(2*K+1);
  S = -S;
  pierr = abs(pi/4 - pisum);
end
pisum = pisum*4;

disp([K,pisum,pi - pisum])
        57648          3.14160999994444     -9.99994444184082e-06

This is an alternating series, so the approximants jump back and forth across the final value. That also explains why it tooks roughly twice as many terms as I might have guessed.

Can we write the series in a vectorized form? This next generates the entire series, so we can plot the convergence.

N = 100000;
K = 0:N;
piapprox = 4*cumsum((1 - mod(K,2)*2)./(2*K+1));

An important thing to learn is to plot everything. Look at what you see. Think about what you see. Does it make sense?

plot(K,pi - piapprox,'.')

Honestly that plot seems pretty boring. But it should give you a hint as to how to view it better

semilogy(K,abs(pi - piapprox),'b.')
xlabel K
ylabel 'Absolute error'

This was more interesting. But it begs the question as to why there appear to be TWO curves there.

semilogy(K(1:2:end),abs(pi - piapprox(1:2:end)),'b.')
hold on
semilogy(K(2:2:end),abs(pi - piapprox(2:2:end)),'r.')
xlabel K
ylabel 'Absolute error'
legend('Even numbered cumulants','Odd numbered cumlants')

Now does it make a little more sense? Looking at the series, the odd numbered cumulants in this series are closer to pi than are the even numbered ones. In fact, that sub-sequence appears to be converging more quickly to pi. What is happening?

To understand that behavior, we need to recall this is an alternating series. Every term switches sign from the previous term. At the same time, every term is quite close in absolute value to the previous term.

So we might break this series up into a new series. I'll do that by grouping pairs of terms.

4*(1 + (-1/3 + 1/5) + (-1/7 + 1/9) + (-1/11 + 1/13) + (-1/15 + 1/17) + (-1/19 + 1/21) + ...)

Combining those pairs of terms, we would get a series as

4*(1 + sum(-1/(4*m+3) + 1/(4*m+5))

where the sum index m goes from 0 to infinity.

N = 10000;
m = 0:N
piapprox2 = 4*(1 - cumsum(2./((4*m+3).*(4*m+5))));
loglog(m,abs(pi - piapprox2),'b.')

This time plotted on log-log axes.