Posts Tagged ‘numerical integration’

Monte Carlo integration

I was reading a bit about random numbers and remembered that I wrote a simple Monte Carlo integrator in C++ a few years ago. I took a few minutes to cleanup and comment the code, which is presented below.

Monte Carlo integration is simple, but surprisingly powerful:

I=abf(x)dx
Ii=1nf(xi)p(xi)

xi is a random value within the range [a,b]

p(x) represents the distribution of random values, for a uniform distribution:

p(x)=1ba

This presentation by Fabio Pellacini provides a lot more details.

The test code in the main() method computes the integral of sin2(x) in the interval [3,5].

#include <iostream>
#include <cmath>
#include <ctime>
using namespace std;

// Functor base class for encapsulating 1-dimensional function to be integrated
class Function1d
{
    
public:
        
virtual double operator()(double x)
        {
            
return x;
        }
};

// Functor for sine squared function
class SineSquared : public Function1d
{
    
public:
        
double operator()(double x)
        {
            
return pow(sin(x), 2);
        }
};

// Monte Carlo integrator class declaration
class MonteCarloIntegrator
{
    
public:
        
double Run(int numSamples, Function1d& func, double intervalMin, double intervalMax);
};

// Monte Carlo integrator implementation
// ::Run() method implementation
double MonteCarloIntegrator::Run(int numSamples, Function1d& func, double intervalMin, double intervalMax)
{
    
double sum = 0;
    
double div = 1.0 / (double)RAND_MAX;
    
double intervalScale = 1.0 / (intervalMax-intervalMin);

    
for(int i=0; i<numSamples; i++)
    {
        
double rnd1 = intervalMin + ( ((double)rand() * div) * (intervalMax-intervalMin) );
        sum += (func)(rnd1) / intervalScale;
    }

    
return (1.0/(double)numSamples) * sum;
}

// main() function with test code
void main()
{
    MonteCarloIntegrator    integrator;    
    
double output = integrator.Run(5000000, SineSquared(), 3, 5);

    srand(time(NULL));
    std::cout << output << endl;

    getchar();
}