Multiple Integrals and Probability : A Numerical

Exploration and an Introduction to MPI

September 8, 2014

The homework questions are due at the 17:00 on Wednesday 17 Septem-

ber

1 Riemann Integration
Recall that we can approximate integrals by Riemann sums. There are many integrals
one cannot evaluate analytically, but for which a numerical answer is required. In
this section, we shall explore a simple way of doing this on a computer. Suppose we
want to find Z 1Z 4
I2D = x2 + 2y 2 dydx.
0 0
If we do this analytically we find

I2D = 44.

Let us suppose we have forgotten how to integrate, and so we do this numerically.

We can do so by splitting integration intervals into equal sub-intervals and using sum
over integrand1 values in the grid points. The following Python code shows how to
do that:

A python program t o compute a Reimann Sum

i m p o r t numpy

nx=1000 # number o f p o i n t s i n x
xend =1.0 # l a s t d i s c r e t i z a t i o n p o i n t

1
The quantity being integrated

1
x s t a r t =0.0 # f i r s t d i s c r e t i z a t i o n point
dx=(xendx s t a r t ) / ( nx 1) # s i z e o f e a c h x subi n t e r v a l

ny=4000 # number o f p o i n t s i n y
yend =4.0 # l a s t d i s c r e t i z a t i o n p o i n t
y s t a r t =0.0 # f i r s t d i s c r e t i z a t i o n p o i n t
dy=(yendy s t a r t ) / ( ny 1) # s i z e o f e a c h y subi n t e r v a l

# c r e a t e v e c t o r s f o r p o i n t s f o r x and y
x=numpy . z e r o s ( ( nx ) , d t y p e= f l o a t ) # a r r a y for x values
y=numpy . z e r o s ( ( ny ) , d t y p e= f l o a t ) # a r r a y for y values

for i i n x r a n g e ( nx ) :
x [ i ]= x s t a r t+i dx

for j i n x r a n g e ( ny ) :
y [ j ]= y s t a r t+j dy

# Approximate t h e i n t e g r a l by a sum
I 2 d =0.00
f o r i i n x r a n g e ( nx ) :
f o r j i n x r a n g e ( ny ) :
I 2 d=I 2 d +(x [ i ] 2 + 2 . 0 y [ j ] 2 ) dydx

print I2d

Listing 1: A Python program which demonstrates how to approximate an integral

by a sum.
We can do something similar in three dimensions. Suppose we want to calculate
Z 1Z 1Z 4
I3D = x2 + 2y 2 + 3z 2 dzdydx.
0 0 0

Analytically we find that

I3D = 68
Homework question 1
a) Modify the Python code to perform the three dimensional integral.
b) Try and determine how the accuracy of either the two or three dimensional
method varies as the number of subintervals is changed.

2 Monte Carlo Integration

If we have many dimensions it may be expensive to calculate sum over all points (see
Section B). The basic idea is to take sum over randomly generated points instead of
grid points.
For example, we can use the formula in eq. (5) to approximate the volume V
under the surface z = x2 + 2y 2 over the rectangle R = (0, 1) (0, 4). Recall that the
actual volume is 44. Below is a Python code that calculates the volume using Monte
Carlo integration

2

A program t o a p p r o x i m a t e an i n t e g r a l u s i n g a Monte C a r l o method

T h i s c o u l d be made f a s t e r by u s i n g v e c t o r i z a t i o n , however i t i s
k e p t a s s i m p l e a s p o s s i b l e f o r c l a r i t y and e a s e o f t r a n s l a t i o n i n t o
other languages

i m p o r t math
i m p o r t numpy
import time
numpoints =4096 # number o f random s a m p l e p o i n t s
I 2 d =0.0 # i n i t i a l i z e value
I 2 d s q u a r e =0.0 # i n i t i a l i z e to allow f o r c a l c u l a t i o n of variance
f o r n i n x r a n g e ( numpoints ) :
x=numpy . random . u n i f o r m ( )
y =4.0numpy . random . u n i f o r m ( )
I 2 d=I 2 d+xx +2.0 yy
I 2 d s q u a r e=I 2 d s q u a r e +(xx +2.0 yy ) 2

# we s c a l e t h e i n t e g r a l by t h e t o t a l a r e a and d i v i d e by t h e number o f
# p o i n t s used
I 2 d=I 2 d 4/ numpoints
I 2 d s q u a r e=I 2 d s q u a r e 4/ numpoints
E s t i m E r r o r =4numpy . s q r t ( ( I 2 d 2 I 2 d s q u a r e ) / ( numpoints 1)) # e s t i m a t e d error
p r i n t Value : %f %I 2 d
p r i n t E r r o r e s t i m a t e : %f %E s t i m E r r o r

Listing 2: A Python program which demonstrates how to use the Monte Carlo
method to calculate the volume below z = x2 + 2y 2 , with (x, y) (0, 1) (0, 4).
The results of running this program with various numbers of random points are
shown below:
N = 16: 41.3026 +/- 30.9791
N = 256: 47.1855 +/- 9.0386
N = 4096: 43.4527 +/- 2.0280
N = 65536: 44.0026 +/- 0.5151
As you can see, the approximation is fairly good. As N , it can be shown
that the
Monte Carlo approximation converges to the actual volume (on the order
of O( N ), in computational complexity terminology).
In the above example the region R was a rectangle. To use the Monte Carlo
method for a nonrectangular (bounded) region R, only a slight modification is needed.
Pick a rectangle R that encloses R, and generate random points in that rectangle as
before. Then use those points in the calculation of f only if they are inside R. There
is no need to calculate the area of R for formula (5) in this case, since the exclusion
of points not inside R allows you to use the area of the rectangle R instead, similar
to before.
For instance, one can show that the volume under the surface z = 1 over the
nonrectangular region R = {(x, y) : 0 x2 + y 2 1} is . Since the rectangle R =
[1, 1] [1, 1] contains R, we can use a similar program to the one we used, the
largest change being a check to see if y 2 + x2 1 for a random point (x, y) in
[1, 1] [1, 1]. A Python code listing which demonstrates this is below:

3

A program t o a p p r o x i m a t e an i n t e g r a l u s i n g a Monte C a r l o method

T h i s c o u l d be made f a s t e r by u s i n g v e c t o r i z a t i o n , however i t i s
k e p t a s s i m p l e a s p o s s i b l e f o r c l a r i t y and e a s e o f t r a n s l a t i o n i n t o
other languages

i m p o r t math
i m p o r t numpy
import time
numpoints =512 # number o f random s a m p l e p o i n t s
I 2 d =0.0 # i n i t i a l i z e value
I 2 d s q u a r e =0.0 # i n i t i a l i z e to allow f o r c a l c u l a t i o n of variance
f o r n i n x r a n g e ( numpoints ) :
x =2.0numpy . random . u n i f o r m ( ) 1 . 0
y =2.0numpy . random . u n i f o r m ( ) 1 . 0
i f ( xx+yy < 1 . 0 ) :
I 2 d=I 2 d +1.0
I 2 d s q u a r e=I 2 d s q u a r e +1.0

# we s c a l e t h e i n t e g r a l by t h e t o t a l a r e a and d i v i d e by t h e number o f
# p o i n t s used
I 2 d=I 2 d 4/ numpoints
I 2 d s q u a r e=I 2 d s q u a r e 4/ numpoints
E s t i m E r r o r =4numpy . s q r t ( ( I 2 d 2 I 2 d s q u a r e ) / ( numpoints 1)) # e s t i m a t e d error
p r i n t Value : %f %I 2 d
p r i n t E r r o r e s t i m a t e : %f %E s t i m E r r o r

Listing 3: A Python program which demonstrates how to use the Monte Carlo
method to calculate the area of an irregular region and also to calculate .
The results of running the program with various numbers of random points are
shown below:

N = 16: 3.5000 +/- 2.9580

N = 256: 3.2031 +/- 0.6641
N = 4096: 3.1689 +/- 0.1639
N = 65536: 3.1493 +/- 0.0407

To use the Monte Carlo method to evaluate triple integrals, you will need to
generate random triples (x, y, z) in a parallelepiped, instead of random pairs (x, y)
in a rectangle, and use the volume of the parallelepiped instead of the area of a
rectangle in formula (5). For a more detailed discussion of numerical integration
methods, please take a further course in mathematics or scientific computing.
Homework question 2
a) Write a program
RR xy that uses the Monte Carlo method to approximate the double
integral e dA, where R = [0, 1] [0, 1]. Show the program output for
R
N = 10, 100, 1000, 10000, 100000 and 1000000 random points.
b) Write a RRR
program that uses the Monte Carlo method to approximate the triple
integral exyz dV , where S = [0, 1] [0, 1] [0, 1]. Show the program output
S
for N = 10, 100, 1000, 10000, 100000 and 1000000 random points.

4
 c) Use the Monte Carlo method to approximate the volume of a sphere of radius
1.

3 Parallel Monte Carlo Integration

As you may have noticed, the algorithms are simple, but can require very many grid
points to become accurate. It is therefore useful to run these algorithms on a parallel
computer. We will demonstrate a parallel Monte Carlo calculation of . Before we
can do this, we need to learn how to use a parallel computer2 . We shall use Rocket.
Information on this computer is available at: http://www.hpc.ut.ee/et

3.1 MPI
Before we can do Monte Carlo integration, we need to learn a little about paral-
lel programming. A copy of the current standard MPI standard can be found at
http://www.mpi-forum.org/. It allows for parallelization of Fortran, C and C++
programs. There are newer parallel programming languages such as Co-Array For-
tran (CAF) and Unified Parallel C (UPC) which allow the programmer to view mem-
ory as a single addressable space even on a distributed memory machine. However,
computer hardware limitations imply that most of the programming concepts used
when writing MPI programs will be required to write programs in CAF and UPC.
Compiler technology for these languages is also not as well developed as compiler
technology for older languages such as Fortran and C, so at the present time, For-
tran and C dominate high performance computing. An introduction to the essential
concepts required for writing and using MPI programs can be found at http://www.
shodor.org/refdesk/Resources/Tutorials/ and at http://www.citutor.org/.
More information on MPI can be found in Gropp, Lusk and Skjellum3 , Gropp, Lusk
and Thakur4 and at https://computing.llnl.gov/tutorials/mpi/. There are
many resources available online, however once the basic concepts have been mas-
tered, what is most useful is an index of MPI commands, usually a search engine will
give you sources of listings, however we have found the following sites useful:
http://publib.boulder.ibm.com/infocenter/zos/v1r13/index.jsp?topic=
%2Fcom.ibm.zos.r13.fomp200%2Fipezps00172.htm
2
Many computers and mobile telephones produced today have 2 or more cores and so can be
considered parallel, but here we mean computers with over hundreds of cores.
3
Using MPI MIT Press(1999)
4
Using MPI 2 MIT Press (1999)

5
 http://www.open-mpi.org/doc/v1.4/

Homework question 3
a) What does MPI stand for?
b) Please read the tutorials at http://www.citutor.org/ or at http://www.
shodor.org/refdesk/Resources/Tutorials/BasicMPI/ and at https://computing.
llnl.gov/tutorials/mpi/, then explain what the following commands do:
USE mpi or INCLUDE mpif.h
MPI INIT
MPI COMM SIZE
MPI COMM RANK
MPI FINALIZE
c) What is the version number of the current MPI standard?
d) Try to understand the Hello World program in listing 4. Run the program in
listing 4 on 20 and 40 MPI processes5 . Put the output of each run in your
solutions, the output will be in a file of the form
job output
An example makefile to compile this program on Rocket is in listing 5. An
example submission script is in listing 6. On Rocket, there is a maximum of 20
cores per node, so if more than 20 single threaded MPI processes are required,
one needs to change the number of nodes as well. The total number of cores
required is equal to the number of nodes multiplied by the number of processes
per node.

! --------------------------------------------------------------------
! PURPOSE
!
! This program uses MPI to print hello world from all available
! processes
!
! .. Scalars ..
! myid = process id
! numprocs = total number of MPI processes
! ierr = error code
!
! REFERENCES
! http :// en . wikipedia . org / wiki / OpenMP
!
! External libraries required
! MPI library

PROGRAM hello90
USE MPI
INTEGER ( kind =4) :: myid , numprocs , ierr

CALL MPI_INIT ( ierr )

CALL MPI_COMM_SIZE ( MPI_COMM_WORLD , numprocs , ierr )

5
One can run this program on many more than 40 processes, however, the output becomes quite
excessive.

6
 CALL MPI_COMM_RANK ( MPI_COMM_WORLD , myid , ierr )

PRINT * , Hello World from process , myid

CALL MPI_BARRIER ( MPI_COMM_WORLD , ierr )
IF ( myid == 0 ) THEN
PRINT * , There are , numprocs , MPI processes
END IF
CALL MPI_FINALIZE ( ierr )
END PROGRAM

Listing 4: A Fortran program which demonstrates parallelizm using MPI.

# define the complier

COMPILER = mpiifort
# compilation settings , optimization , precision , parallelization
FLAGS = - O0

# libraries
LIBS =
# source list for main program
SOURCES = helloworld . f90

test : $ ( SOURCES )
$ { COMPILER } -o helloworld $ ( FLAGS ) $ ( SOURCES )

clean :
rm *. o

clobber :
rm helloworld

Listing 5: An example makefile for compiling the helloworld program in listing 4.

# !/ bin / bash

# SBATCH -N 2
# SBATCH -- sockets - per - node =2
# SBATCH -- cores - per - socket =10
# SBATCH -- threads - per - core =1
# SBATCH -t 00:10:00

module load i n t e l _ c l u s t e r _ s t u d i o

export I_ MPI _PM I_L IBRA RY =/ usr / lib64 / libpmi . so

srun -- cpu_bind = cores helloworld

Listing 6: An example submission script for use on Rocket.

We now examine a Fortran program for calculating . These programs are taken
from http://chpc.wustl.edu/mpi-fortran.html, where further explanation can
be found. The original source of these programs appears to be Using MPI by Gropp,
Lusk and Skjellum.

Serial
! --------------------------------------------------------------------
! PURPOSE
!
! This program use a monte carlo method to calculate pi
!
! .. Parameters ..
! npts = total number of Monte Carlo points
! xmin = lower bound for integration region
! xmax = upper bound for integration region

7
 ! .. Scalars ..
! i = loop counter
! f = average value from summation
! sum = total sum
! randnum = random number generated from (0 ,1) uniform
! distribution
! x = current Monte Carlo location
!
! REFERENCES
! http :// chpc . wustl . edu / mpi - fortran . html
! Gropp , Lusk and Skjellum , " Using MPI " MIT press (1999)
PROGRAM monte_carlo
IMPLICIT NONE

INTEGER ( kind =8) , PARAMETER :: npts = 1 e4

REAL ( kind =8) , PARAMETER :: xmin =0.0 d0 , xmax =1.0 d0
INTEGER ( kind =8) :: i
REAL ( kind =8) :: f , sum , randnum , x

DO i =1 , npts
CALL random_number ( randnum )
x = ( xmax - xmin )* randnum + xmin
sum = sum + 4.0 d0 /(1.0 d0 + x **2)
END DO
f = sum / npts
PRINT * , PI calculated with , npts , points = ,f

STOP
END

Listing 7: A serial Fortran program which demonstrates how to calculate using a

Monte Carlo method.

# define the complier

COMPILER = mpiifort
# compilation settings , optimization , precision , parallelization
FLAGS = - O0

# libraries
LIBS =
# source list for main program
SOURCES = montecarloserial . f90

test : $ ( SOURCES )
$ { COMPILER } -o montecarloserial $ ( FLAGS ) $ ( SOURCES )

clean :
rm *. o

clobber :
rm montecarloserial

Listing 8: An example makefile for compiling the program in listing 7.

# !/ bin / bash
# number of nodes and number of processors per node requested
# in this case 1 compute nodes , 1 MPI tasks per node
# SBATCH -N 1
# SBATCH -- sockets - per - node =1
# SBATCH -- cores - per - socket =1
# SBATCH -- threads - per - core =1
# requested Wall - clock time .
# SBATCH -t 00:10:00
# name of the job , you may want to change this so it is unique to you
# SBATCH -J MPI_MCSERIAL
# Email address to send a notification to , change " youremail " appropriately
# SBATCH - mail - user = youremail@ut . ee
# send a notification for job abort , begin and end
# SBATCH - mail - type = all

# setup environment
module load i n t e l _ c l u s t e r _ s t u d i o

8
export I_ MPI _PM I_L IBRA RY =/ usr / lib64 / libpmi . so

# change to the job submission directory

cd $PBS_O_WORKDIR
# Run the job
srun -- cpu_bind = cores montecarloserial

Listing 9: An example submission script for use on Rocket.

Parallel
! --------------------------------------------------------------------
! PURPOSE
!
! This program uses MPI to do a parallel monte carlo calculation of pi
!
! .. Parameters ..
! npts = total number of Monte Carlo points
! xmin = lower bound for integration region
! xmax = upper bound for integration region
! .. Scalars ..
! mynpts = this processes number of Monte Carlo points
! myid = process id
! nprocs = total number of MPI processes
! ierr = error code
! i = loop counter
! f = average value from summation
! sum = total sum
! mysum = sum on this process
! randnum = random number generated from (0 ,1) uniform
! distribution
! x = current Monte Carlo location
! start = simulation start time
! finish = simulation end time
!
! REFERENCES
! http :// chpc . wustl . edu / mpi - fortran . html
! Gropp , Lusk and Skjellum , " Using MPI " MIT press (1999)
!
! External libraries required
! MPI library
PROGRAM monte_carlo_mpi
USE MPI
IMPLICIT NONE

INTEGER ( kind =8) , PARAMETER :: npts = 1 e8

REAL ( kind =8) , PARAMETER :: xmin =0.0 d0 , xmax =1.0 d0
INTEGER ( kind =8) :: mynpts
INTEGER ( kind =4) :: ierr , myid , nprocs
INTEGER ( kind =8) :: i
REAL ( kind =8) :: f , sum , mysum , randnum
REAL ( kind =8) :: x , start , finish

! Initialize MPI
CALL MPI_INIT ( ierr )
CALL MPI_COMM_RANK ( MPI_COMM_WORLD , myid , ierr )
CALL MPI_COMM_SIZE ( MPI_COMM_WORLD , nprocs , ierr )
start = MPI_WTIME ()

! Calculate the number of points each MPI process needs to generate

IF ( myid . eq . 0) THEN
mynpts = npts - ( nprocs -1)*( npts / nprocs )
ELSE
mynpts = npts / nprocs
ENDIF

! set initial sum to zero

mysum = 0.0 d0
! use loop on local process to generate portion of Monte Carlo integral
DO i =1 , mynpts
CALL random_number ( randnum )
x = ( xmax - xmin )* randnum + xmin
mysum = mysum + 4.0 d0 /(1.0 d0 + x **2)

9
 ENDDO

! Do a reduction and sum the results from all processes

CALL MPI_REDUCE ( mysum , sum ,1 , MPI_DOUBLE_PRECISION , MPI_SUM ,&
0 , MPI_COMM_WORLD , ierr )
finish = MPI_WTIME ()

! Get one process to output the result and running time

IF ( myid . eq . 0) THEN
f = sum / npts
PRINT * , PI calculated with , npts , points = ,f
PRINT * , Program took , finish - start , to calculate the integral
ENDIF

CALL MPI_FINALIZE ( ierr )

STOP
END PROGRAM

Listing 10: A parallel Fortran program which demonstrates how to calculate using
MPI.

# define the complier

COMPILER = mpiifort
# compilation settings , optimization , precision , parallelization
FLAGS = - O0

# libraries
LIBS =
# source list for main program
SOURCES = m on te c ar lo p ar a ll el . f90

test : $ ( SOURCES )
$ { COMPILER } -o m o nt ec a rl o pa ra l le l $ ( FLAGS ) $ ( SOURCES )

clean :
rm *. o

clobber :
rm m o nt ec a rl op a ra l le l

Listing 11: An example makefile for compiling the program in listing 10.

# !/ bin / bash
# number of nodes and number of processors per node requested
# in this case 2 compute nodes , 10 MPI tasks per node
# SBATCH -N 2
# SBATCH -- sockets - per - node =2
# SBATCH -- cores - per - socket =10
# SBATCH -- threads - per - core =1
# requested Wall - clock time .
# SBATCH -t 00:10:00
# name of the job , you may want to change this so it is unique to you
# SBATCH -J MPI_MCPARALLEL
# Email address to send a notification to , change " youremail " appropriately
# SBATCH - mail - user = youremail@ut . ee
# send a notification for job abort , begin and end
# SBATCH - mail - type = all

# setup environment
module load i n t e l _ c l u s t e r _ s t u d i o

export I_ MPI _PM I_L IBRA RY =/ usr / lib64 / libpmi . so

# change to the job submission directory

cd $PBS_O_WORKDIR
# Run the job
srun -- cpu_bind = cores m o nt e ca rl o pa r al le l

Listing 12: An example submission script for use on Rocket.

10
Homework question 4
a) Explain why using Monte Carlo to evaluate
Z 1
4
dx
0 1 + x2

allows you to find and, in your own words (hint: look up Integration table),
explain what the serial and parallel programs do.
b) Find the time it takes to run the Parallel Monte Carlo program on 32, 64, 128,
256 and 512 cores.

Homework question 5
a) Use a parallel Monte Carlo integration program to evaluate
ZZ
x2 + y 6 + exp(xy) cos(y exp(x))dA

over the unit circle.

b) Use a parallel Monte Carlo integration program to approximate the volume of
2 2 2
the ellipsoid x9 + y4 + z1 = 1.
c) Write parallel programs to find the volume of the 4 dimensional sphere
4
X
1 x2i .
i=1

Try both Monte Carlo and Riemann sum techniques.

Homework question 6
a) Why are Monte Carlo methods used in Finance?
b) Is speed of evaluation useful in such applications?
c) What is hpcwire? Consider signing up for a thrice weekly digest email from
hpcwire which will inform you about all that is happening in the HPC world..
d) Find an online job advertisement that requires high performance computing
skills. What is the expected salary? What other skills are typically desired?

11
A Probability
Definition A.1. f : U R2 R+ is a probability density function if
Z Z
f dA = 1
U

Definition A.2. If f is a probability density function which takes the set U R2 ,

then the probability of events in the set W U occurring is
Z Z
P (W ) = f dA.
W

Example A.1. The joint density for it to snow x inches tomorrow and for Immo
to win y dollars in the lottery tomorrow is given by
c
f=
(1 + x)(100 + y)

for
x, y [0, 100] [0, 100]
and f = 0 otherwise. Find c.

Definition A.3. Suppose X is a random variable with probability density function

f1 (x) and Y is a random variable with a probability density function f2 (y). Then X
and Y are independent random variables if their joint density function is

f (x, y) = f1 (x)f2 (y).

Example A.2. The probability it will snow tomorrow and the probability Immo will
win the lottery tomorrow are independent random variables.

Definition A.4. If f (x, y) is a probability density function for the random variables
X and Y , the X mean is Z Z

1 = X = xf dA

and the Y mean is Z Z

2 = Y = yf dA.

Remark A.1. The X mean and the Y mean are the expected values of X and Y.

12
Definition A.5. If f (x, y) is a probability density function for the random variables
X and Y , the X variance is
Z Z
2
1 = (X X) = 2 2 f dA
(x X)

and the Y variance is

Z Z
22 = (Y Y )2 = (y Y )2 f dA.

Definition A.6. The standard deviation is defined to be the square root of the vari-
ance.

Example A.3. Find an expression for the probability that it will snow more than
1.1 times the expected snowfall and also that Immo will win more than 1.2 times the
expected amount in the lottery.

B Monte Carlo Integration

6
It is possible to extend the Riemann integration schemes to higher and higher
dimensional integrals. This can become computationally intensive and an alternate
method of integration based on probability is often used. The method we will discuss
is called the Monte Carlo method. The idea behind it is based on the concept of the
average value of a function, which you learned in single-variable calculus. Recall that
for a continuous function f (x), the average value f of f over an interval [a, b] is
defined as
Z b
1
f = f (x) dx . (1)
ba a

The quantity b a is the length of the interval [a, b], which can be thought of
as the volume of the interval. Applying the same reasoning to functions of two or
three variables, we define the average value of f (x, y) over a region R to be
ZZ
1
f = f (x, y) dA , (2)
A(R)
R

6
This section is taken from Chapter 3 of Vector Calculus by Michael Corral which is available at
http://www.mecmath.net/ and where Java and Sage programs for doing Monte Carlo integration
can be found.

13
where A(R) is the area of the region R, and we define the average value of f (x, y, z)
over a solid S to be
ZZZ
1
f = f (x, y, z) dV , (3)
V (S)
S

where V (S) is the volume of the solid S. Thus, for example, we have
ZZ
f (x, y) dA = A(R)f . (4)
R

The average value of f (x, y) over R can be thought of as representing the sum of
all the values of f divided by the number of points in R. Unfortunately there are
an infinite number (in fact, uncountably many) points in any region, i.e. they can
not be listed in a discrete sequence. But what if we took a very large number N of
random points in the region R (which can be generated by a computer) and then
took the average of the values of f for those points, and used that average as the
value of f? This is exactly what the Monte Carlo method does. So in formula (4)
the approximation we get is
s
f 2 (f)2
ZZ
f (x, y) dA A(R)f A(R) , (5)
N 1
R

where
PN PN 2
f (xi , yi ) i=1 (f (xi , yi ))
f = i=1
and f 2 = , (6)
N N
with the sums taken over the N random points (x1 , y1 ), . . ., (xN , yN ). The error
term in formula (5) does not really provide hard bounds on the approximation.
It represents a single unbiased estimate of the standard deviation from the expected
value of the integral[9]. That is, it provides a likely bound on the error. Due to its use
of random points, the Monte Carlo method is an example of a probabilistic method (as
opposed to deterministic methods such as the Riemann sum approximation method,
which use a specific formula for generating points).

References
[1] R. Courant and F. John, Introduction to Calculus and Analysis I, II
Springer (1998,1999)

14
[2] D. Hughes-Hallett, A.M. Gleason, D.E. Flath, P.F. Lock, D.O. Lomen, D. Love-
lock, W.G. MacCallum, D. Mumford, B. G. Osgood, D. Quinney, K. Rhea, J.
Tecosky-Feldman, T.W. Tucker, and O.K. Bretscher, A. Iovita, W. Raskind, S.P.
Gordon, A. Pasquale, J.B. Thrash, Calculus, Single and Multivariable, 5th
ed. Wiley, (2008)

[3] J. Marsden and A. Tromba, Vector Calculus, W.H. Freeman, (2003).

[4] P. Lax, S. Burstein and A. Lax, Calculus with Applications and Comput-
ing, Vol. 1 Springer, (1976).

[5] W.P. Petersen and P. Arbenz, Introduction to Parallel Computing: A

Practical Guide with Examples in C Oxford University Press, (2004).

[6] W. Gropp, E. Lusk and A. Skjellum Using MPI MIT Press (1999)

[7] W. Gropp, E. Lusk and A. Skjellum Using MPI-2 MIT Press (1999)

[8] http://www.mpi-forum.org/docs/docs.html

[9] https://en.wikipedia.org/wiki/Bias_of_an_estimator

15