Child pages
• Question 13 - Parallel Monte Carlo Method for calculating an integral
Go to start of banner

Question 13 - Parallel Monte Carlo Method for calculating an integral

Objective: to calculate the area under the curve f(x) > 0 in the interval [a,b], which corresponds to the integral from a to b of f(x)
Approach: Consider a rectangle of area A =(b -a)*sup(f(x)), if we throw two random numbers (x,y) on this rectangle , the probability that this pair of numbers fall under the curve f(x) is proportional to the area of the curve under f(x) which is the integral of f(x) between a and b. And the probability is the number of counts that fall under the curve divided by the total number of throws. And the the area under the curve would be the total area of the rectangle multiplied by the probability.

Output: Estimate of the integral of f(x) > 0 on [a,b]

Algorithm:

1. Each process generates two random numbers (x,y) between [a,b] for x, and between [0,sup(f(x))] for y. This is done for a given specified number of times.

2. Each process checks if the generated random pairs satisfy

y < f(x) , i.e. if the pairs are below the curve of f(x)

For each successful throw increase the counter by one, this loop continues until the given specified number of times.

3. Each process, except process 0 , sends its count to process 0.

4. Process 0 sums all the counts received from other process and calculates the area and prints out the result

Code:

```//This program calcualtes the area above the
//function y = x^3 between x = 0 and x = 2

#include "mpi.h"
#include <math.h>
#include <cstdlib>
#include <iostream>
#include <time.h>

int main( int argc, char *argv[] )
{
int n, myid, numprocs,source, tag, dest;
int count = 0;
int mycount = 0;
MPI_Status status;
MPI_Init(&argc,&argv);
MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
MPI_Comm_rank(MPI_COMM_WORLD,&myid);
int m=100000;
dest =0;
tag =100;
//to make sure we do not have same consecutive random numbers
srand ( time(NULL) );

for(int i = 0;i<m;i++)
{
//to make sure x is between  between 0 and 2
double x = (double)2*rand()/(RAND_MAX + 1.0);
//to make sure y is between   0 and 8
double y = (double)8*rand()/(RAND_MAX + 1.0);
if (y < x*x*x)
{
mycount++;
}
}
//Collecting the counts from each process by the master node
if(myid ==0)
{
count = mycount;
for(int i =1;i<numprocs;i++)
{
source = i;
MPI_Recv(&mycount,1,MPI_DOUBLE,source,tag, MPI_COMM_WORLD,&status);
count+= mycount;
}

}

//process other than the master sent their count to process 0
else
MPI_Send(&mycount, 1, MPI_DOUBLE, dest, tag,MPI_COMM_WORLD);

//calculating the area from process 0
if(myid==0)
{
//calculating the area above y = x^3
//the probability multipled by the total area of the rectangle

double area = (double)16*count/((double)(numprocs*m));
std::cout<<"area = "<<area<<"\n";
}

MPI_Finalize();
return 0;
}
```
• No labels