Child pages
• Exercise Pi. Anyone for Pi?
Go to start of banner

Exercise Pi. Anyone for Pi?

Go to start of metadata

One way to calculate pi is using the Monte Carlo method devised by Nicholas Constantine Metropolis. It goes as follows.

Start with a circle of radius r=1 centred on the origin. The circle sits within a square of  dimension l=2r=2  also centred on the origin like so: The ratio of the areas of the circle and square is given by: LaTeX Markup problem Show/Hide latex markup \documentclass{article} \usepackage{amsmath} \usepackage{amsthm} \usepackage{amssymb} \usepackage{bm} \newcommand{\mx}{\mathbf{\bm{#1}}} % Matrix command \newcommand{\vc}{\mathbf{\bm{#1}}} % Vector command \newcommand{\T}{\text{T}} % Transpose \pagestyle{empty} \begin{document} $\rho = \frac{area of circle}{area of square} = \frac{\pi r^2}{2r}^2 = \frac{\pi}{4}$ \end{document} Show/Hide latex error Unable to find DVI conversion log file.

So to get some p all we have to do is find this ratio and multiply by 4. Simple...

But wait, how do we get the areas, and therefore ratio, without using the formulas?

This is where Monte Carlo comes in handy. If we start placing random points within the square and measure what proportion also fall within the circle then we'll get a convergent Create a Monte Carlo code to calculate p  using the intrinsic function random_number:

PROGRAM mcpi
implicit none
integer :: NTHROW,HITS,i
real    :: dist,x2,y2,x,y,ratio,pi
!  This program uses a 'hit and miss' Monte Carlo integration
!  to determine the value of pi.
!WRITE(6,*) 'Enter the number of throws desired'
!
NTHROW = 1000000
!
!WRITE(6,*) 'Enter an integer random number seed'
!
HITS=0
DO i=1,NTHROW
call random_number(x)
call random_number(y)
X2=X*X
Y2=Y*Y
DIST=SQRT(X2+Y2)
!
IF (DIST.LE.1.0) THEN
HITS=HITS+1
!write(*,"(e15.8,2x,e15.8,i4)") x,y,2
else
!write(*,"(e15.8,2x,e15.8,i4)") x,y,4
ENDIF
!
RATIO=REAL(HITS)/REAL(i)
!
PI=4.0*RATIO
!
if (mod(i,100)==0) write(10,"(i,2x,e21.12)") i,pi
!
enddo
RATIO=REAL(HITS)/REAL(NTHROW)
PI=4.0*RATIO
WRITE(*,"('The value of pi is calculated to be ',F14.8)") PI
!
END PROGRAM mcpi

Exercise: Compile and run with different numbers of NTHROW

You can print out the values of pi as a function of n and get them visualaised with gnuplot using pi-1D.plt, see