MATLAB: Using a For Loop to calculate the pi for a taylor series

forloops

Write a program (using a loop) that determines for a given n. Run the program with n = 10, n = 100, and n = 1,000. Compare the result with pi. (Use format long.)

This is my code thus far,

clear;clc;

format Long

n=input('Enter the Value of n ');

Sum_n=0;

for ii=1:length(n)

Sum_n=Sum_n+((-1).^n)/((2.*n+1).^3);

end

value= nthroot(32,3)*Sum_n;

I know what I want to do with the code, I just dont know how to input it. I want the the FOR loop to add each increment of n+1 until it reaches the n the user inputs. I want the the sum of all those n's and in the end it should give me something that outpts a number that gives me pi.

Best Answer

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.

https://en.wikipedia.org/wiki/List_of_formulae_involving_π

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

Related Solutions

MATLAB: MATLAB program in a script file that determines

Ok, you made a decent start, so lets see what went wrong.

First, try not to use variables named error, or any other existing, useful function names. I'd suggest a different name, perhaps errortol, that indicates better what the variable means, here an error tolerance.

Next, you have a test in the while loop definition. But your test is p>=pivalue. Since this is an alternating series, I expect the result to flip-flop over and under pi. So >= is just abad idea. Anyway, you don't need it, as you then test inside the loop, and use a break.

Next, you are trying to stuff the terms into a vector sigma that will grow with each iteration. growing a vector dynamically is usually a bad idea, as that forces MATLAB to reallocate memory for the vector at every iteration. Luckily, we know this cannot take too many iterations. But, better is to preallocate sigma if you want to keep the individual terms in the sum. I'm not sure why you would bother though. Anyway, that way you are forced to re-sum all of the terms at EVERY iteration. So that is quite inefficient.

A big problem was you were computing pivalue wrong. Here is the offending line:

pivalue= abs(sum(sigma(n)).*32).^(1/3);

See that you formed this subexpression:

sum(sigma(n))

n is a SCALAR! so sigma(n) is just the last term in that vector. If you wanted to form the sum, sum(sigma) would have been better.

Finally, you had one other problem, and this is why I wanted to change the name of error to errortol. It makes it clearer to you what the test inside the loop should be.

    if errortol >= abs(pi - pivalue)

You used <=, but really, you need to stop only when the tolerance is greater than the error from the true value of pi.

I might use nthroot, instead of taking a cuberoot using .^(1/3). nthroot is a bit safer. The point is, if the number is negative, then the form you used will compute ONE of the nth roots of a negative number, which will just happen to be a complex number.

nthroot(-1,3)
ans =
    -1
(-1).^(1/3)
ans =
                        0.5 +     0.866025403784439i

As I said, nthroot will not surprise you.

oops. One other thing. In that series, the first term starts at n=0. You increment n before computing the term,so you are effectively starting at n=1.

There are many things to think of when you write good code. You will learn to be careful, but that takes some practice.

So here is how I would have changed your code:

errortol=1.0e-6;
pivalue=0;
sigma=0;
n=0;
% the while loop could run forever like this, but we
% know that it won't, because of that test inside.
while true
  sigma = sigma + (-1).^(n)/(2.*n+1).^3;
  pivalue = nthroot(sigma.*32,3);
  if errortol >= abs(pi - pivalue)
    break
  end
  n=n+1;
end

I will point out that my version is very little different from yours. You had the right idea all along. You just needed some minor fixes in the implementation.

When I run this, I see it has stopped at n=40;

n
n =
    40
pivalue
pivalue =
          3.14159363278731

Yes, there were a bunch of things I had to fix, but they were all subtle changes, worth thinking about why I did it that way.

Having said all of that, how might I have written this in a more classical MATLAB style? Yes, I know that you need to use a loop here. That teaches you much. But you can also learn by comparison.

n = 0:100;
piest = nthroot(cumsum((-1).^n./(2*n+1).^3)*32,3);
semilogy(n,abs(pi-piest))
grid on

Here you can see clearly at what term the approximation crosses below 1e-6.

It is worth looking at what I wrote, and understanding that too.

MATLAB: Having trouble with for loops

out = 1;
for ii = 2:1000
   out = out + ii^2;
end

Best Answer

Related Solutions

MATLAB: MATLAB program in a script file that determines

MATLAB: Having trouble with for loops

Related Question