Category: Monte Carlo Integration

Blog Entry © Friday, October 23, 2026, by James Pate Williams, Jr. Equation of State Calculations

Equation of State CalculationsDownload
// EquationOfState.cpp : This file contains the 'main' function.
// EquationOfState.cpp : This file contains the 'main' function.
// Program execution begins and ends there.
// Copyright (c) Friday, October 23, 2026
// by James Pate Williams, Jr. BA, BS, MSwE, PhD
// Reference: osti.gov/servlets/purl/4390578/			

#include <math.h>
#include <iostream>
#include <vector>

const int N = 224;
const int nRows = 14;
const int nCols = 16;
const double d = 1.0 / nRows;

typedef struct tagPoint2d
{
    double x, y;
} Point2d, *PPoint2d;

double d0 = 0.0, AA0 = 0.0, AA2 = 0.0, maxD = 0.0, pi = 0.0;
int numberSteps = 0;
std::vector<Point2d> particlePositions(N);

static double Distance(double x1, double y1, double x2, double y2)
{
    // 2-D Euclidean norm
    return sqrt(pow(x1 - x2, 2.0) + pow(y1 - y2, 2.0));
}

static double Random01()
{
    // random number in the interval [0, 1]
    double value = static_cast<double>(rand()) / RAND_MAX;
    return value;
}

static double Random11()
{
    // random number in the interval [-1, 1]
    double value = -1.0 + 2.0 * static_cast<double>(rand()) / RAND_MAX;
    return value;
}

static void InitializeData(double K, double nu)
{
    pi = 4.0 * atan(1.0);
    d0 = d * (1.0 - pow(2.0, nu - 8.0));
    AA0 = 1.0 / (sqrt(3.0) * d0 * d0 * N / 2.0);
    maxD = d - d0;
    AA2 = (K * K - 1.0) * pi * d0 * d0 / 64.0;
    AA0 = 1.0 / (sqrt(3.0) * d0 * d0 * N / 2.0);
    maxD = d - d0;
    AA2 = (K * K - 1.0) * pi * d0 * d0 / 64.0;
}

static void InitializeLattice()
{
    // initialize to a regular 2-D Cartesian lattice
    for (size_t r = 0; r < nRows; r++)
    {
        for (size_t c = 0; c < nCols; c++)
        {
            size_t address = r * nRows + c;
            particlePositions[address].x = r / 14.0;
            particlePositions[address].y = c / 16.0;
        }
    }
}

static int Move()
{
    bool result = true;
    double r = 1.0 / nRows;
    double c = 1.0 / nCols;

    InitializeLattice();
     
    for (int n = 0; n < N; n++)
    {
        double xj = particlePositions[n].x;
        double yj = particlePositions[n].y;
        double xk = 1.25 * N * r * Random11() + xj;
        double yk = 1.25 * N * c * Random11() + yj;
        double rr =
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0);

        if (rr < d)
            result = false;
    }

    return result;
}

static double A33()
{
    double a33 = 0.0;
    
    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0, f31 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 3 + 1;
        int k = rand() % 3 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a33 += 1.0;
        else if (j == 2 && k == 3 && rjk < d)
            a33 += 1.0;
        else if (j == 3 && k == 1 && rjk < d)
            a33 += 1.0;
        else if (j == 2 && k == 1 && rjk < d)
            a33 += 1.0;
        else if (j == 3 && k == 2 && rjk < d)
            a33 += 1.0;
        else if (j == 1 && k == 3 && rjk < d)
            a33 += 1.0;
    }

    a33 /= (pi * pi * pow(d, 4.0));
    a33 /= numberSteps;
    return a33;
}

static double A44()
{
    double a44 = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f41 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 4 + 1;
        int k = rand() % 4 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a44 += 1.0;
        else if (j == 2 && k == 3 && rjk < d)
            a44 += 1.0;
        else if (j == 3 && k == 4 && rjk < d)
            a44 += 1.0;
        else if (j == 4 && k == 1 && rjk < d)
            a44 += 1.0;
    }

    a44 /= (pi * pi * pow(d, 4.0));
    a44 /= numberSteps;
    return a44;
}

static double A45()
{
    double a45 = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f41 = 0.0;
        double f13 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 4 + 1;
        int k = rand() % 4 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a45 += 1.0;
        else if (j == 2 && k == 3 && rjk < d)
            a45 += 1.0;
        else if (j == 3 && k == 4 && rjk < d)
            a45 += 1.0;
        else if (j == 4 && k == 1 && rjk < d)
            a45 += 1.0;
        else if (j == 1 && k == 3 && rjk < d)
            a45 += 1.0;
    }

    a45 /= (pi * pi * pow(d, 4.0));
    a45 /= numberSteps;
    return a45;
}

static double A46()
{
    double a46 = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f41 = 0.0;
        double f13 = 0.0, f24 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 4 + 1;
        int k = rand() % 4 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a46 += 1.0;
        else if (j == 2 && k == 3 && rjk < d)
            a46 += 1.0;
        else if (j == 3 && k == 4 && rjk < d)
            a46 += 1.0;
        else if (j == 4 && k == 1 && rjk < d)
            a46 += 1.0;
        else if (j == 1 && k == 3 && rjk < d)
            a46 += 1.0;
        else if (j == 2 && k == 4 && rjk < d)
            a46 += 1.0;
    }

    a46 /= (pi * pi * pow(d, 4.0));
    a46 /= numberSteps;
    return a46;
}

static double A55()
{
    double a55 = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f45 = 0.0;
        double f51 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a55 += 1.0;
        else if (j == 2 && k == 3 && rjk < d)
            a55 += 1.0;
        else if (j == 3 && k == 4 && rjk < d)
            a55 += 1.0;
        else if (j == 4 && k == 5 && rjk < d)
            a55 += 1.0;
        else if (j == 5 && k == 1 && rjk < d)
            a55 += 1.0;
    }

    a55 /= (pi * pi * pow(d, 4.0));
    a55 /= numberSteps;
    return a55;
}

static double A56p()
{
    double a56p = 0.0;
    
    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f31 = 0.0;
        double f41 = 0.0, f35 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a56p += 1.0;
        else if (j == 2 && k == 3 && rjk < d)
            a56p += 1.0;
        else if (j == 3 && k == 4 && rjk < d)
            a56p += 1.0;
        else if (j == 3 && k == 1 && rjk < d)
            a56p += 1.0;
        else if (j == 4 && k == 1 && rjk < d)
            a56p += 1.0;
        else if (j == 3 && k == 5 && rjk < d)
            a56p += 1.0;
    }

    a56p /= (pi * pi * pow(d, 4.0));
    a56p /= numberSteps;
    return a56p;
}

static double A56pp()
{
    double a56pp = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f41 = 0.0;
        double f35 = 0.0, f51 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a56pp += 1.0;
        else if (j == 2 && k == 3 && rjk < d)
            a56pp += 1.0;
        else if (j == 3 && k == 4)
            a56pp += 1;
        else if (j == 4 && k == 1 && rjk < d)
            a56pp += 1.0;
        else if (j == 3 && k == 5 && rjk < d)
            a56pp += 1.0;
        else if (j == 5 && k == 1 && rjk < d)
            a56pp += 1.0;
    }

    a56pp /= (pi * pi * pow(d, 4.0));
    a56pp /= numberSteps;
    return a56pp;
}

static double A57p()
{
    double a57p = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f45 = 0.0;
        double f51 = 0.0, f35 = 0.0;
        double f41 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a57p += 1.0;
        else if (j == 2 && k == 3 && rjk < d)
            a57p += 1.0;
        else if (j == 3 && k == 4 && rjk < d)
            a57p += 1.0;
        else if (j == 4 && k == 5 && rjk < d)
            a57p += 1.0;
        else if (j == 5 && k == 1 && rjk < d)
            a57p += 1.0;
        else if (j == 3 && k == 5 && rjk < d)
            a57p += 1.0;
        else if (j == 4 && k == 1 && rjk < d)
            a57p += 1.0;
    }

    a57p /= (pi * pi * pow(d, 4.0));
    a57p /= numberSteps;
    return a57p;
}

static double A57pp()
{
    double a57pp = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f45 = 0.0;
        double f51 = 0.0, f31 = 0.0;
        double f41 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a57pp += 1.0;
        else if (j == 2 && k == 3 && rjk < d)
            a57pp += 1.0;
        else if (j == 3 && k == 4 && rjk < d)
            a57pp += 1.0;
        else if (j == 4 && k == 5 && rjk < d)
            a57pp += 1.0;
        else if (j == 5 && k == 1 && rjk < d)
            a57pp += 1.0;
        else if (j == 3 && k == 1 && rjk < d)
            a57pp += 1.0;
        else if (j == 4 && k == 1 && rjk < d)
            a57pp += 1.0;
    }

    a57pp /= (pi * pi * pow(d, 4.0));
    a57pp /= numberSteps;
    return a57pp;
}

static double A57ppp()
{
    double a57ppp = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f41 = 0.0;
        double f35 = 0.0, f51 = 0.0;
        double f31 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a57ppp += 1.0;
        else if (j == 2 && k == 3 && rjk < d)
            a57ppp += 1.0;
        else if (j == 3 && k == 4 && rjk < d)
            a57ppp += 1.0;
        else if (j == 4 && k == 1 && rjk < d)
            a57ppp += 1.0;
        else if (j == 3 && k == 5 && rjk < d)
            a57ppp += 1.0;
        else if (j == 5 && k == 1 && rjk < d)
            a57ppp += 1.0;
        else if (j == 3 && k == 1 && rjk < d)
            a57ppp += 1.0;
    }

    a57ppp /= (pi * pi * pow(d, 4.0));
    a57ppp /= numberSteps;
    return a57ppp;
}

static double A58p()
{
    double a58p = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f45 = 0.0;
        double f51 = 0.0, f31 = 0.0;
        double f35 = 0.0, f41 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a58p += 1.0;
        else if (j == 2 && k == 3 && rjk < d)
            a58p += 1.0;
        else if (j == 3 && k == 4 && rjk < d)
            a58p += 1.0;
        else if (j == 4 && k == 5 && rjk < d)
            a58p += 1.0;
        else if (j == 5 && k == 1 && rjk < d)
            a58p += 1.0;
        else if (j == 3 && k == 1 && rjk < d)
            a58p += 1.0;
        else if (j == 3 && k == 5 && rjk < d)
            a58p += 1.0;
        else if (j == 4 && k == 1 && rjk < d)
            a58p += 1.0;
    }

    a58p /= (pi * pi * pow(d, 4.0));
    a58p /= numberSteps;
    return a58p;
}

static double A58pp()
{
    double a58pp = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f45 = 0.0;
        double f51 = 0.0, f31 = 0.0;
        double f41 = 0.0, f25 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a58pp += 1.0;
        else if (j == 2 && k == 3 && rjk < d)
            a58pp += 1.0;
        else if (j == 3 && k == 4 && rjk < d)
            a58pp += 1.0;
        else if (j == 4 && k == 5 && rjk < d)
            a58pp += 1.0;
        else if (j == 5 && k == 1 && rjk < d)
            a58pp += 1.0;
        else if (j == 3 && k == 1 && rjk < d)
            a58pp += 1.0;
        else if (j == 4 && k == 1 && rjk < d)
            a58pp += 1.0;
        else if (j == 4 && k == 1 && rjk < d)
            a58pp += 1.0;
        else if (j == 2 && k == 5 && rjk < d)
            a58pp += 1.0;
    }

    a58pp /= (pi * pi * pow(d, 4.0));
    a58pp /= numberSteps;
    return a58pp;
}

static double A59()
{
    double a59 = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f45 = 0.0;
        double f51 = 0.0, f31 = 0.0;
        double f41 = 0.0, f25 = 0.0;
        double f35 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a59 += 1.0;
        else if (j == 2 && k == 3 && rjk < d)
            a59 += 1.0;
        else if (j == 3 && k == 4 && rjk < d)
            a59 += 1.0;
        else if (j == 4 && k == 5 && rjk < d)
            a59 += 1.0;
        else if (j == 5 && k == 1 && rjk < d)
            a59 += 1.0;
        else if (j == 3 && k == 1 && rjk < d)
            a59 += 1.0;
        else if (j == 3 && k == 5 && rjk < d)
            a59 += 1.0;
        else if (j == 4 && k == 1 && rjk < d)
            a59 += 1.0;
        else if (j == 2 && k == 5 && rjk < d)
            a59 += 1.0;
        else if (j == 3 && k == 5 && rjk < d)
            a59 += 1.0;
    }

    a59 /= (pi * pi * pow(d, 4.0));
    a59 /= numberSteps;
    return a59;
}

static double A510()
{
    double a510 = 0.0;

    for (int i = 0; i < numberSteps; i++)
    {
        double f12 = 0.0, f23 = 0.0;
        double f34 = 0.0, f45 = 0.0;
        double f51 = 0.0, f31 = 0.0;
        double f41 = 0.0, f25 = 0.0;
        double f35 = 0.0, f24 = 0.0;
        double xj = maxD * Random11();
        double yj = maxD * Random11();
        double xk = maxD * Random11();
        double yk = maxD * Random11();
        double rjk = sqrt(
            pow(xj - xk, 2.0) +
            pow(yj - yk, 2.0));
        int j = rand() % 5 + 1;
        int k = rand() % 5 + 1;

        if (j == 1 && k == 2 && rjk < d)
            a510 += 1.0;
        if (j == 2 && k == 3 && rjk < d)
            a510 += 1.0;
        if (j == 3 && k == 4 && rjk < d)
            a510 += 1.0;
        if (j == 4 && k == 5 && rjk < d)
            a510 += 1.0;
        if (j == 5 && k == 1 && rjk < d)
            a510 += 1.0;
        if (j == 3 && k == 1 && rjk < d)
            a510 += 1.0;
        if (j == 3 && k == 5 && rjk < d)
            a510 += 1.0;
        if (j == 4 && k == 1 && rjk < d)
            a510 += 1.0;
        if (j == 2 && k == 5 && rjk < d)
            a510 += 1.0;
        if (j == 2 && k == 4 && rjk < d)
            a510 += 1.0;
    }

    a510 /= (pi * pi * pow(d, 4.0));
    a510 /= numberSteps;
    return a510;
}

int main()
{
    double K[8] = { 2.5, 2.45, 2.4, 2.35, 2.3, 2.25, 2.2, 2.15 };
    double nu[8] = { 2.0, 4.0, 5.0, 5.5, 6.0, 6.25, 6.5, 7.0 };
    unsigned int seed = 1;

    srand(seed);

    for (int i = 0; i < 8; i++)
    {
        double a33 = 0.0, a44 = 0.0, a45 = 0.0, a46 = 0.0, a55 = 0.0;
        double a56p = 0.0, a56pp = 0.0, a57p = 0.0, a57pp = 0.0;
        double a57ppp = 0.0, a58p = 0.0, a58pp = 0.0;
        double N12 = 0.0, a59 = 0.0, a510 = 0.0;

        InitializeData(K[i], nu[i]);

        for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
            a33 += A33();

        for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
            a44 += A44();

        for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
            a45 += A45();

        for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
            a46 += A46();

        for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
            a55 += A55();

        for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
            a56p += A56p();

        for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
            a56pp += A56pp();

        for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
            a57p += A57p();

        for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
            a57pp += A57pp();

        for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
            a57ppp += A57ppp();

        for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
            a58p += A58p();

        for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
            a58pp += A58pp();

        for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
            a59 += A59();

        for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
            a510 += A510();

        a33 /= 4.0;
        a44 /= 4.0;
        a45 /= 4.0;
        a46 /= 4.0;
        a55 /= 4.0;
        a56p /= 4.0;
        a56pp /= 4.0;
        a57p /= 4.0;
        a57pp /= 4.0;
        a57ppp /= 4.0;
        a58p /= 4.0;
        a58pp /= 4.0;
        a59 /= 4.0;
        a510 /= 4.0;

        std::cout << "A(3, 3)    = " << a33 << std::endl;
        std::cout << "A(4, 4)    = " << a44 << std::endl;
        std::cout << "A(4, 5)    = " << a45 << std::endl;
        std::cout << "A(4, 6)    = " << a46 << std::endl;
        std::cout << "A(5, 5)    = " << a55 << std::endl;
        std::cout << "A(5, 6')   = " << a56p << std::endl;
        std::cout << "A(5, 6'')  = " << a56pp << std::endl;
        std::cout << "A(5, 7')   = " << a57p << std::endl;
        std::cout << "A(5, 7'')  = " << a57pp << std::endl;
        std::cout << "A(5, 7''') = " << a57ppp << std::endl;
        std::cout << "A(5, 8')   = " << a58p << std::endl;
        std::cout << "A(5, 8'')  = " << a58pp << std::endl;
        std::cout << "A(5, 9)    = " << a59 << std::endl;
        std::cout << "A(5, 10)   = " << a510 << std::endl;

        for (int i = 0; i < 6500; i++)
        {
            if (Move())
                N12 += 1.0;
        }

        std::cout << "N(1/2) = ";
        std::cout << N12 << std::endl;

        double C0 = 64.0 * N12 / (N * N * (K[i] * K[i] - 1.0));
        double A0 = 0.5 * sqrt(3.0) * d0 * d0 * N;

        std::cout << "C0   = " << C0 << std::endl;
        std::cout << "A0/A = " << A0 << std::endl;
        std::cout << "A/A0 = " << 1.0 / A0 << std::endl;

        double Ni = 1.0 / N12;
        double C1 = pi / sqrt(3.0);
        double C2 = (4.0 * pi * pi * A33() / 9.0) * Ni;
        double C3a = (6.0 * A45() - 3.0 * A44() - A46()) / pow(3.0, 1.5);
        double C3 = (pi * pi * pi * C3a) * Ni;
        double C40 = 8.0 * pi * pi * pi / 135.0;
        double C4a = 12.0 * a55 - 60.0 * a56p - 10.0 * a56pp;
        double C4b = 30.0 * a57p + 60.0 * a57pp + 10.0 * a57ppp;
        double C4c = -30.0 * a58p - 15.0 * a58pp;
        double C4d = 10.0 * a59 - a510;
        double C4 = (C40 * (C4a + C4b + C4c + C4d)) * Ni;

        std::cout << "K    = " << K[i] << std::endl;
        std::cout << "nu   = " << nu[i] << std::endl;
        std::cout << "C(1) = " << C1 << std::endl;
        std::cout << "C(2) = " << C2 << std::endl;
        std::cout << "C(3) = " << C3 << std::endl;
        std::cout << "C(4) = " << C4 << std::endl;

        double X3 = C1 * A0 + C2 * A0 * A0 + C3 * A0 * A0 * A0 +
            C4 * A0 * A0 * A0 * A0;

        std::cout << "X(3) = " << X3 << std::endl;
    }

    return 0;
}

Blog Entry © Thursday, May 14, 2026, by James Pate Williams, Jr. More Numerical Integration Results

Blog Entry Thursday May 14 2026Download
// 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;
}

Blog Entry © Wednesday, May 13, 2026, by James Pate Williams, Jr. Adaptive n-Quadrature Versus Monte Carlo Integration

Blog Entry Wednesday May 13 2026Download

Blog Entry © Sunday, November 9, 2025, by James Pate Williams, Jr. Hydrogenic Wavefunctions, Radial Probability Functions, Distribution Functions, and First Moment Integrals

Blog Entry Sunday November 09 2025Download

Blog Entry © Sunday, August 10, 2025, First-Order Perturbation Treatment of the Helium Atom by James Pate Williams, Jr.

Blog Entry Sunday August 10 2025Download

Blog Entry © Tuesday, July 29, 2025, Double and Triple Monte Carlo Integration by James Pate Williams, Jr.

Blog Entry Tuesday July 29 2025Download
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

static double randomRange(double lo, double hi)
{
	return (hi - lo) * (double)rand() / RAND_MAX + lo;
}

static double integrand(double r, double w)
{
	return pow(r, 4.0) * (2.0 - r) * w * w * exp(-r);
}

static double StarkEffectIntegral(double E, int N)
{
	double sum = 0.0;

	for (int i = 0; i <= N; i++)
	{
		double r = randomRange(0.0, 100.0);
		double w = randomRange(-1.0, 1.0);

		sum += integrand(r, w);
	}

	return 100.0 * 2.0 * E * sum / (16.0 * (N - 1));
}

static void firstOrderStarkEffect(double E)
{
	double exact = -3.0 * E;
	int N[9] = {
		1000000, 2000000, 3000000, 4000000,
		5000000, 6000000, 7000000, 8000000,
		9000000 };

	for (int n = 0; n < 9; n++)
	{
		int iN = N[n];
		double integ = StarkEffectIntegral(E, iN);
		double error = 100.0 * fabs(integ - exact) / fabs(exact);

		printf("N = %4ld\tintegral = %13.10lf\t%% error = %13.10lf\n",
			iN, integ, error);
	}

	printf("exact value = %13.10lf\n", exact);
}

static double ee1(int N, double R, double Z)
{
	double pi = 4.0 * atan(1.0);
	double sum = 0.0;

	for (int i = 0; i <= N; i++)
	{
		double r1 = randomRange(1.0e-25, R);
		double r2 = randomRange(0.0, r1);

		sum += R * r1 * r1 * exp(-2.0 * Z * (r1 + r2)) * r2 * r2;
	}

	return 16.0 * pi * pi * sum / (N - 1);
}

static double ee2(int N, double R, double Z)
{
	double pi = 4.0 * atan(1.0);
	double sum = 0.0;

	for (int i = 0; i <= N; i++)
	{
		double r1 = randomRange(1.0e-25, R);
		double r2 = randomRange(r1, R);
		
		sum += R * (R - r2) * r2 * exp(-2.0 * Z * (r1 + r2)) * r1 * r1;
	}

	return 16.0 * pi * pi * sum / (N - 1);
}

static void firstOrderHelium(double Z)
{
	double pi = 4.0 * atan(1.0), R = 25.0;
	double exact = 5.0 * pi * pi / (8.0 * pow(Z, 5.0));

	int N[9] = {
		1000000, 2000000, 3000000, 4000000,
		5000000, 6000000, 7000000, 8000000,
		9000000 };

	for (int n = 0; n < 9; n++)
	{
		int iN = N[n];
		double integ = ee1(iN, R, Z) + ee2(iN, R, Z);
		double error = 100.0 * fabs(integ - exact) / fabs(exact);

		printf("N = %4ld\tintegral = %13.10lf\t%% error = %13.10lf\n",
			iN, integ, error);
	}

	printf("exact value = %13.10lf\n", exact);
}

int main(void)
{
	firstOrderStarkEffect(2.0);
	firstOrderHelium(2.0);
	return 0;
}

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

static double randomRange(double lo, double hi)
{
	return (hi - lo) * (double)rand() / RAND_MAX + lo;
}

static double f(double x, double y, double z)
{
	return pow(sin(x), 2.0) + y * sin(z);
}

static double g(double x, double y, double z)
{
	return x + y * z * z;
}

static double integral(
	double x0, double x1,
	double y0, double y1,
	double z0, double z1,
	double (*f)(double, double, double),
	int N)
{
	double sum = 0.0;

	for (int n = 0; n <= N; n++)
	{
		double x = randomRange(x0, x1);
		double y = randomRange(y0, y1);
		double z = randomRange(z0, z1);

		sum += f(x, y, z);
	}

	return (x1 - x0) * (y1 - y0) * (z1 - z0) *
		sum / (N - 1);
}

int main(void)
{
	double pi = 4.0 * atan(1.0);
	double x0 = 0.0, x1 = pi;
	double y0 = 0.0, y1 = 1.0;
	double z0 = 0.0, z1 = pi;
	double exact = 0.5 * pi * (2.0 + pi);
	int N[9] = {
		1000000, 2000000, 3000000, 4000000,
		5000000, 6000000, 7000000, 8000000,
		9000000 };

	printf("integrand pow(sin(x), 2.0) + y * sin(z)\n");
	printf("x = 0 to pi, y = 0 to 1, z = 0 to pi\n");

	for (int n = 0; n < 9; n++)
	{
		int iN = N[n];
		double integ = integral(
			x0, x1, y0, y1, z0, z1, f, iN);
		double error = 100.0 * fabs(integ - exact) / fabs(exact);

		printf("N = %4ld\tintegral = %13.10lf\t%% error = %13.10lf\n",
			iN, integ, error);
	}

	printf("exact value = %13.10lf\n", exact);

	x0 = -1.0;
	x1 = 5.0;
	y0 = 2.0;
	y1 = 4.0;
	z0 = 0.0;
	z1 = 1.0;
	exact = 36.0;

	printf("integrand x + y * z * z\n");
	printf("x = -1 to 5, y = 2 to 4, z = 0 to 1\n");

	for (int n = 0; n < 9; n++)
	{
		int iN = N[n];
		double integ = integral(
			x0, x1, y0, y1, z0, z1, g, iN);
		double error = 100.0 * fabs(integ - exact) / fabs(exact);

		printf("N = %4ld\tintegral = %13.10lf\t%% error = %13.10lf\n",
			iN, integ, error);
	}

	printf("exact value = %13.10lf\n", exact);
	return 0.0;
}

Four Methods of Numerical Double Integration: Sequential Simpson’s Rule, Multitasking Simpson’s Rule, Sequential Monte Carlo and Multitasking Monte Carlo Methods © Wednesday April 16 – 18, 2025, by James Pate Williams, Jr.

Four Methods of Numerical Double IntegrationDownload

Approximation of the Ground-State Total Energy of a Beryllium Atom © Sunday, March 30 to Tuesday April 1, 2025, by James Pate Williams, Jr., BA, BS, Master of Software Engineering, PhD Computer Science

Approximation of the Ground-State Total Energy of a Beryllium AtomDownload

Blog Entry © Thursday, March 27, 2025, by James Pate Williams, Jr., BA, BS, Master of Software Engineering, PhD Lithium (Li, Z = 3) Total Ground-State Energy Numerical Experiments

Lithium Atom ExperimentsDownload

Revised Helium Ground-State Total Energy Computation (c) Sunday, March 23, 2025, by James Pate Williams, Jr.

Revised Computation of the Ground-State Energy of a Helium AtomDownload