Monte Carlo Method

Integration

Integration using Monte Carlo method:

src/cpp/coding-practice/01-monte-carlo-integration.cpp

// Monte carlo method of integration
// calculate area under the curve y = x^2 within x limits [1, 3]
#include <iostream>
#include <cmath>   // pow
#include <cstdlib> // rand
using namespace std;

int main()
{
    double x_lower = 1;
    double x_upper = 3;
    int iterations = 1'000'000;

    double y_lower = 0;
    double y_upper = pow(x_upper, 2);
    double area_rect = y_upper * (x_upper - x_lower);

    int total_area = 0;
    int monte_area = 0;
    double result = 0;

    for (int i = 0; i < iterations; i++)
    {
        double x, y;
        x = rand() % 1000; // generate random integer upto 1000
        y = rand() % 1000;

        x = x / 1000; // scale the numbers to [0, 1]
        y = y / 1000;

        x = x_lower + x * (x_upper - x_lower);
        y = y_lower + y * (y_upper - y_lower);

        // this check is unnecessary since all points are bound be in this range
        if (x >= x_lower && x <= x_upper)
        {
            total_area += 1;

            if (y <= pow(x, 2))
            {
                monte_area += 1;
            }
        }
    }

    if (total_area != 0)
    {
        result = area_rect * monte_area / total_area;
    }

    cout << "Area = " << result << endl;
    cout << "Analytical result = " << (27.0 - 1) / 3 << endl;

    return 0;
}

Value of pi

Approximate value of pi using Monte Carlo method:

src/cpp/coding-practice/02-monte-carlo-pi.cpp

// find out value of pi by calculating area of circle by monte carlo method
// area of circle = pi * r^2
// area of enclosing rectangle = (2 * r)^2
#include <iostream>
#include <cstdlib>
#include <cmath>
using namespace std;

int main()
{
    int iterations = 1'000'000;
    double x_lower = -1.0;
    double y_lower = -1.0;
    double x_upper = 1.0;
    double y_upper = 1.0;
    double monte_area = 0;

    for (int i = 0; i < iterations; i++)
    {
        double x, y;
        x = rand() % 1000; // generate random integer upto 1000
        y = rand() % 1000;

        x = x / 1000; // scale the numbers to [0, 1]
        y = y / 1000;

        x = x_lower + x * (x_upper - x_lower);
        y = y_lower + y * (y_upper - y_lower);

        if (x * x + y * y < 1)
        {
            monte_area += 1;
        }
    }

    // pi = 4 * area_circ / area_rect = 4 * monte_area / iterations
    cout << "pi = " << 4 * monte_area / iterations << endl;
    cout << "actual value = " << 4 * atan(1) << endl;

    return 0;
}