Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Can't accurately calculate pi on Python

I am new member here and I'm gonna drive straight into this as I've spent my whole Sunday trying to get my head around it.

I'm new to Python, having previously learned coding on C++ to a basic-intermediate level (it was a 10-week university module).

I'm trying a couple of iterative techniques to calculate Pi but both are coming up slightly inaccurate and I'm not sure why.

The first method I was taught at university - I'm sure some of you have seen it done before.

x=0.0
y=0.0
incircle = 0.0
outcircle = 0.0
pi = 0.0
i = 0
while (i<100000):
    x = random.uniform(-1,1)
    y = random.uniform(-1,1)
    if (x*x+y*y<=1):
        incircle=incircle+1
    else:
        outcircle=outcircle+1
    i=i+1
pi = (incircle/outcircle)
print pi

It's essentially a generator for random (x,y) co-ordinates on a plane from -1 to +1 on both axes. Then if x^2+y^2 <= 1, we know the point rests inside a circle of radius 1 within the box formed by the co-ordinate axes.

Depending on the position of the point, a counter increases for incircle or outcircle.

The value for pi is then the ratio of values inside and outside the circle. The co-ordinates are randomly generated so it should be an even spread.

However, even at very high iteration values, my result for Pi is always around the 3.65 mark.

The second method is another iteration which calculates the circumference of a polygon with increasing number of sides until the polygon is almost a circle, then, Pi=Circumference/diameter. (I sort of cheated because the coding has a math.cos(Pi) term so it looks like I'm using Pi to find Pi, but this is only because you can't easily use degrees to represent angles on Python). But even for high iterations the final result seems to end around 3.20, which again is wrong. The code is here:

S = 0.0
C = 0.0
L = 1.0

n = 2.0
k = 3.0
while (n<2000):
    S = 2.0**k
    L = L/(2.0*math.cos((math.pi)/(4.0*n)))
    C = S*L
    n=n+2.0
    k=k+1.0

pi = C/math.sqrt(2.0)
print pi

I remember, when doing my C++ course, being told that the problem is a common one and it isn't due to the maths but because of something within the coding, however I can't remember exactly. It may be to do with the random number generation, or the limitations of using floating point numbers, or... anything really. It could even just be my mathematics...

Can anyone think what the issue is?

TL;DR: Trying to calculate Pi, I can get close to it but never very accurately, no matter how many iterations I do.

(Oh and another point - in the second code there's a line saying S=2.0**k. If I set 'n' to anything higher than 2000, the value for S becomes too big to handle and the code crashes. How can I fix this?)

Thanks!

like image 619
Stuart Aitken Avatar asked Jan 17 '16 16:01

Stuart Aitken


2 Answers

The algorithm for your first version should look more like this:

from __future__ import division, print_function

import sys
if sys.version_info.major < 3:
    range = xrange

import random 


incircle = 0
n = 100000
for n in range(n):
    x = random.random()
    y = random.random()
    if (x*x + y*y <= 1):
        incircle += 1
pi = (incircle / n) * 4
print(pi)

Prints:

3.14699146991

This is closer. Increase n to get even closer to pi.

The algorithm takes into account only one quarter of the unit circle, i.e. with a radius of 1.

The formula for the area of a quarter circle is:

area_c = (pi * r **2) / 4

That for the area of the square containing this circle:

area_s = r **2

where r is the radius of the circle.

Now the ratio is:

area_c / area_s

substitute the equations above, re-arange and you get:

pi = 4 * (area_c / area_s)

Going Monte Carlo, just replace both areas by a very high number that represents them. Typically, the analogy of darts thrown randomly is used here.

like image 83
Mike Müller Avatar answered Oct 29 '22 09:10

Mike Müller


For the first one, your calculation should be

pi = incircle/1000000*4  # 3.145376..

This is the number of points that landed inside of the circle over the number of total points (approximately 0.785671 on my run).

With a radius of 1 (random.uniform(-1,1)), the total area is 4, so if you multiple 4 by the ratio of points that landed inside of the circle, you get the correct answer.

like image 34
14 revs, 12 users 16% Avatar answered Oct 29 '22 11:10

14 revs, 12 users 16%