compute the approximation of pi
58 次查看(过去 30 天)
显示 更早的评论
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...)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/352445/image.gif)
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
0 个评论
采纳的回答
Steven Lord
2020-8-27
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.
4 个评论
Steven Lord
2020-8-27
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.
更多回答(1 个)
John D'Errico
2020-8-27
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,'.')
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/352544/image.jpeg)
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'
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/352547/image.jpeg)
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')
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/352550/image.jpeg)
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.')
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/352553/image.jpeg)
This time plotted on log-log axes.
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!