## 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$
$I≈∑i=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)=1b–a$

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 integratedclass Function1d{    public:        virtual double operator()(double x)         {            return x;         }};// Functor for sine squared functionclass SineSquared : public Function1d{    public:        double operator()(double x)        {            return pow(sin(x), 2);        }};// Monte Carlo integrator class declarationclass MonteCarloIntegrator{    public:        double Run(int numSamples, Function1d& func, double intervalMin, double intervalMax);};// Monte Carlo integrator implementation// ::Run() method implementationdouble 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 codevoid main(){    MonteCarloIntegrator    integrator;        double output = integrator.Run(5000000, SineSquared(), 3, 5);    srand(time(NULL));    std::cout << output << endl;    getchar();}`