I had to look online to find a formula for pi that used an alternating sum of reciprocals of cubes of odd numbers. There had to be one, and Wikipedia has it, though I had to check a few different sites before I saw that series.
In there, I see that we have
pi^3/32 = 1 - 1/3^3 + 1/5^3 -1/7^3 + ...
That means, you solve for the series, then to compute the current approximation to pi, you must multiply by 32, and then take the cube root. I think you got confused there, because you were multiplying by the cube root of 32 in the code you wrote.
A simple code to implement this for 21 terms is:
n = (0:20)';
[n,nthroot(32*cumsum((-1).^n.*1./(2*n+1).^3),3)]
ans =
0 3.1748021039364
1 3.13511291696101
2 3.14377083641878
3 3.14062114485715
4 3.14210388509366
5 3.14129194905678
6 3.14178389122763
7 3.14146367259822
8 3.14168365475933
9 3.14152608792951
10 3.14164278860378
11 3.1415539618301
12 3.14162313060566
13 3.14156822245079
14 3.14161253590592
15 3.14157625789938
16 3.14160633164585
17 3.14158112444928
18 3.14160246099173
19 3.14158424155425
20 3.14159992269186
Of course, this is not what Jose wants to use, but that is how I would write the series sums, using MATLAB as it might be used. As you can see, it does reasonably well in converging to about 6 digits after only 21 terms.
A simple looped code might be:
N = [10;100;1000];
sum_n = zeros(size(N));
for iN = 1:numel(N)
for k = 0:N(iN)
sum_n(iN) = sum_n(iN) + (-1)^k/(2*k+1)^3;
end
end
piapprox = nthroot(sum_n*32,3);
[N,piapprox,piapprox - pi]
ans =
10 3.14164278860378 5.01350139914258e-05
100 3.14159271914105 6.55512568670247e-08
1000 3.14159265365714 6.7346128673762e-11
Best Answer