Monte Carlo Integration

Numerical Methods

Monte Carlo integration is a numerical method that uses random sampling to estimate integrals. While not as accurate as deterministic methods for low-dimensional integrals, it becomes essential for high-dimensional problems.

Basic Monte Carlo Integration

For an integral:

I = \int_{a}^{b} f (x) d x

The Monte Carlo estimate is:

I \approx \frac{b - a}{N} i = 1 \sum N f (x_{i})

where $x_{i}$ are uniformly distributed random numbers in $[a, b]$ .

Convergence

The error in Monte Carlo integration scales as:

Error \sim \frac{σ}{N}

where $σ$ is the standard deviation of $f (x)$ . This $1/ N$ convergence is independent of dimension, making Monte Carlo superior for high-dimensional integrals.

C++ Implementation

The following code implements Monte Carlo integration:

#include <functional>
#include <random>
#include <cmath>

double monteCarloIntegral(double a, double b, int N, 
                          std::function<double(double)> func) {
    // Random number generator
    std::random_device dev;
    std::mt19937 rng(dev());
    std::uniform_real_distribution<double> dis{a, b};

    double result = 0.0;
    for(int i = 0; i < N; i++) {
        double random_num = dis(rng);
        result += func(random_num);
    }
    
    result = (1.0 / static_cast<double>(N)) * result * (b - a);
    return result;
}

// Example: Integrate exp(x) from -10 to 1
int main() {
    auto exp_func = static_cast<double(*)(double)>(std::exp);
    
    double integral = monteCarloIntegral(-10.0, 1.0, 1000000, exp_func);
    double exact = std::exp(1.0);
    
    std::cout << "Monte Carlo result: " << integral << std::endl;
    std::cout << "Exact result: " << exact << std::endl;
    std::cout << "Error: " << std::abs(integral - exact) << std::endl;
    
    return 0;
}

// Example: Integrate Gaussian (tests numerical precision)
void testGaussianIntegration() {
    double alpha = 1.0;
    auto gauss = [alpha](double x) { 
        return std::exp(-alpha * x * x); 
    };
    
    double analytical = std::sqrt(std::numbers::pi_v<double> / alpha);
    double integral = monteCarloIntegral(-5.0, 5.0, 1000000, gauss);
    
    std::cout << "Gaussian integral (MC): " << integral << std::endl;
    std::cout << "Analytical: " << analytical << std::endl;
    std::cout << "Error: " << std::abs(integral - analytical) << std::endl;
}

Key Features

Simple implementation using C++11 random number generators
Generic function interface
Convergence independent of dimension
Useful for testing numerical precision

Advantages and Limitations

Advantages:

Dimension-independent convergence rate
Easy to implement
Works for complex integration domains
Parallelizable

Limitations:

Slow convergence ( $O (1/ N)$ )
Not competitive for 1D or 2D integrals
Requires good random number generator
Variance reduction techniques often needed