Numerical Explorations

// NumericalIntegrals.cpp (c) Thursday, May 14, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD

#include <iomanip>
#include <iostream>
#include <vector>
#include <stdlib.h>

static double f(double x) {
    return sin(x);
}

static double MonteCarlo(double a, double b,
    double (*f)(double), int n){ 
    double sum = 0;

    for (int i = 0; i < n; i++) {
        double x = (b - a) * (double)rand() / RAND_MAX;

        sum += f(x);
    }

    return (b - a) * sum / n;
}

static double Factorial(int n) {
    double factorial = 1.0;

    for (int i = 2; i <= n; i++)
        factorial *= i;

    return factorial;
}

static double Series(double a, double b, int n)
{
    double sumA = 0.0, sumB = 0.0;
    int sign = 1;

    for (int i = 0; i <= n; i++) {
        sumA += sign * pow(a, 2 * i + 2) /
            Factorial(2 * i + 2);

        sign *= -1;
    }

    sign = 1;

    for (int i = 0; i <= n; i++) {      
        sumB += sign * pow(b, 2 * i + 2) /
            Factorial(2 * i + 2);

        sign *= -1;
    }

    return sumB - sumA;
}

static double CompositeTrapezoidalRule(
    double a, double b, int n) {
    double pi = 4.0 * atan(1.0);
    double endPts = 0.5 * (f(a) + f(b));
    double sum = 0, xk = 0.0;
    double h = (b - a) / n;

    for (int k = 1; k <= n - 1; k++) {
        xk = a + k * h;
        sum += f(xk);
    }

    return h * (0.5 * endPts + sum);
}

static double SimpsonsRule(
    int n, double a, double b, double(*fx)(double)) {
    double h = (b - a) / n;
    double h2 = 2.0 * h;
    double s = 0.0;
    double t = 0.0;
    double x = a + h;

    for (int i = 1; i < n; i += 2) {
        s += fx(x);
        x += h2;
    }

    x = a + h2;

    for (int i = 2; i < n; i += 2) {
        t += fx(x);
        x += h2;
    }

    return h * (fx(a) + 4 * s + 2 * t + fx(b)) / 3.0;
}

static void Romberg(double a, double b,
    double (*f)(double), int mStart, int nRow,
    std::vector<std::vector<double>>& T) {
    int m = mStart;
    double h = (b - a) / m;
    double sum = 0.5 * (f(a) + f(b));

    if (m > 1) {
        for (int i = 1; i <= m - 1; i++) {
            sum += f(a + i * h);
        }
    }

    T[0][0] = sum * h;

    std::cout << "romberg t-table" << std::endl;
    std::cout << std::fixed;
    std::cout << std::setprecision(5) << T[0][0];
    std::cout << std::endl;

    if (nRow < 2)
        return;

    for (int k = 2; k <= nRow; k++) {
        h /= 2.0;
        m *= 2;
        sum = 0.0;

        for (int i = 1; i <= m; i += 2) {
            sum += f(a + i * h);
        }

        T[k][1] = 0.5 * T[k - 1LL][1] + sum * h;

        for (int j = 1; j <= k - 1; j++) {
            T[k - 1LL][j] = T[k][j] - T[k - 1LL][j];
            T[k][j + 1LL] = T[k][j] + T[k - 1LL][j] /
                (pow(4.0, j) - 1.0);
        }

        for (int j = 1; j <= k; j++) {
            std::cout << std::fixed;
            std::cout << std::setprecision(5);
            std::cout << T[k][j] << '\t';
        }

        std::cout << std::endl;
    }

    if (nRow < 3) {
        return;
    }

    std::cout << "table of ratios" << std::endl;
    double ratio = 0.0;

    for (int k = 1; k <= nRow - 2; k++) {
        for (int j = 1; j <= k; j++) {
            if (T[k + 1LL][j] == 0.0) {
                ratio = 0.0;
            }

            else {
                ratio = T[k][j] / T[k + 1LL][j];
            }

            T[k][j] = ratio;
        }

        for (int j = 1; j <= k; j++) {
            std::cout << std::fixed;
            std::cout << std::setprecision(5);
            std::cout << T[k][j] << '\t';
        }

        std::cout << std::endl;
    }
}

double MonteCarloVolume(double R, int n)
{
    double pi = 4.0 * atan(1.0), pi2 = 2.0 * pi;
    double R2 = R * R, sum = 0;

    for (int i = 0; i < n; i++)
    {
        double r = R2 * (double)rand() / RAND_MAX;
        double t = pi * (double)rand() / RAND_MAX;
        double p = pi2 * (double)rand() / RAND_MAX;
        sum += r * r * sin(t);
    }

    return R * pi * pi2 * sum / n;
}

int main()
{
    srand(1);
    std::vector<std::vector<double>> T;
    T.resize(35);
    for (int i = 0; i < 35; i++) {
        T[i].resize(35);
    }
    Romberg(0.0, 2.0, f, 2, 7, T);
    std::cout << std::setprecision(11);
    std::cout << "analytic integral of sine = " << -cos(2.0) + cos(0.0);
    std::cout << std::endl;
    std::cout << "simpson's rule integral   = " << SimpsonsRule(500, 0, 2.0, f);
    std::cout << std::endl;
    std::cout << "monte carlo integral      = " << MonteCarlo(0.0, 2.0, f, 2130);
    std::cout << std::endl;
    std::cout << "infinite series integral  = " << Series(0.0, 2.0, 16);
    std::cout << std::endl;
    double integral = CompositeTrapezoidalRule(0.0, 2.0, 175000000);
    std::cout << "romberg integral          = " << integral << std::endl;
    std::cout << "actual spherical volume   = " << 4.0 * 4.0 * atan(1.0) / 3.0;
    std::cout << std::endl;
    double volume = MonteCarloVolume(1.0, 1000000);
    std::cout << "approx spherical volume   = " << volume;
    std::cout << std::endl;
}