UCL WIKI

UCL Logo
Child pages
  • Exercise Pi. Anyone for Pi?
Skip to end of metadata
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

Show/Hide latex error

 

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'
      !READ(5,*) NTHROW
      !
      NTHROW = 1000000
      !
      !WRITE(6,*) 'Enter an integer random number seed'
      !READ(5,*) idum
      !
      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

https://wiki.ucl.ac.uk/display/GradPrgCrs/gnuplot+scripts

gnuplot -persist pi-1D.plt

or with xmgrace 

xmgrace fort.10

Exercise: Convert this code into a 3D code by using the ratio between the box and sphere volumes. 

 

  • No labels