Category: Monte Carlo Integration
Blog Entry (c) Sunday, November 1, 2026, by James Pate Williams, Jr. Some Hard Sphere Calculations
References: “Introduction to Statistical Thermodynamics” (c) 1960 by Terrell L. Hill 8.14.2: Properties of Nobel Gases – Chemistry LibreTexts
First, we approximate the hard sphere second virial coefficient for the Noble Gases using Monte Carlo integration. Analytic integration yields:
B2(T) = 2 * pi * a ^ 3 / 3
Where a is the diameter of the atom. The volume utilized was V = X ^ 3 where X = 25.0 * 100.0 * a. Then we compute four cluster integrals.
// HardSphereZ3.cpp : This file contains the 'main' function.
// Program execution begins and ends there.
// Copyright (c) Sunday, November 1, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// all potential energy is for a hard sphere
#include <cmath>
#include <iomanip>
#include <iostream>
#include <string>
// second virial coefficient
static double B2(double X, double a, int N)
{
double I = 0.0, V = X * X * X;
for (int n = 0; n < N; n++)
{
double x1 = (X * rand()) / RAND_MAX;
double x2 = (X * rand()) / RAND_MAX;
double y1 = (X * rand()) / RAND_MAX;
double y2 = (X * rand()) / RAND_MAX;
double z1 = (X * rand()) / RAND_MAX;
double z2 = (X * rand()) / RAND_MAX;
double r12 = std::sqrt(
std::pow(x1 - x2, 2.0) +
std::pow(y1 - y2, 2.0) +
std::pow(z1 - z2, 2.0));
if (r12 >= a)
I += 1.0;
}
return V * I / N;
}
// three particle cluster integral
static double Z30(double X, double a, int N)
{
double I = 0.0, V = X * X * X;
for (int n = 0; n < N; n++)
{
double x1 = (X * rand()) / RAND_MAX;
double x2 = (X * rand()) / RAND_MAX;
double x3 = (X * rand()) / RAND_MAX;
double y1 = (X * rand()) / RAND_MAX;
double y2 = (X * rand()) / RAND_MAX;
double y3 = (X * rand()) / RAND_MAX;
double z1 = (X * rand()) / RAND_MAX;
double z2 = (X * rand()) / RAND_MAX;
double z3 = (X * rand()) / RAND_MAX;
double r12 = std::sqrt(
std::pow(x1 - x2, 2.0) +
std::pow(y1 - y2, 2.0) +
std::pow(z1 - z2, 2.0));
double r13 = std::sqrt(
std::pow(x1 - x3, 2.0) +
std::pow(y1 - y3, 2.0) +
std::pow(z1 - z3, 2.0));
double r23 = std::sqrt(
std::pow(x2 - x3, 2.0) +
std::pow(y2 - y3, 2.0) +
std::pow(z2 - z3, 2.0));
double u = 0.0;
if (r12 >= a)
u += 1.0;
if (r13 >= a)
u += 1.0;
if (r23 >= a)
u += 1.0;
if (u > 0.0)
I += u;
}
return V * I / N;
}
// four particle cluster integral
static double Z40(double X, double a, int N)
{
double I = 0.0, V = X * X * X;
for (int n = 0; n < N; n++)
{
double x1 = (X * rand()) / RAND_MAX;
double x2 = (X * rand()) / RAND_MAX;
double x3 = (X * rand()) / RAND_MAX;
double x4 = (X * rand()) / RAND_MAX;
double y1 = (X * rand()) / RAND_MAX;
double y2 = (X * rand()) / RAND_MAX;
double y3 = (X * rand()) / RAND_MAX;
double y4 = (X * rand()) / RAND_MAX;
double z1 = (X * rand()) / RAND_MAX;
double z2 = (X * rand()) / RAND_MAX;
double z3 = (X * rand()) / RAND_MAX;
double z4 = (X * rand()) / RAND_MAX;
double r12 = std::sqrt(
std::pow(x1 - x2, 2.0) +
std::pow(y1 - y2, 2.0) +
std::pow(z1 - z2, 2.0));
double r23 = std::sqrt(
std::pow(x2 - x3, 2.0) +
std::pow(y2 - y3, 2.0) +
std::pow(z2 - z3, 2.0));
double r34 = std::sqrt(
std::pow(x3 - x4, 2.0) +
std::pow(y3 - y4, 2.0) +
std::pow(z3 - z4, 2.0));
double r41 = std::sqrt(
std::pow(x4 - x1, 2.0) +
std::pow(y4 - y1, 2.0) +
std::pow(z4 - z1, 2.0));
double u = 0.0;
if (r12 >= a)
u += 1.0;
if (r23 >= a)
u += 1.0;
if (r34 >= a)
u += 1.0;
if (r41 >= a)
u += 1.0;
if (u > 0.0)
I += u;
}
return V * I / N;
}
static double Z41(double X, double a, int N)
{
double I = 0.0, V = X * X * X;
for (int n = 0; n < N; n++)
{
double x1 = (X * rand()) / RAND_MAX;
double x2 = (X * rand()) / RAND_MAX;
double x3 = (X * rand()) / RAND_MAX;
double x4 = (X * rand()) / RAND_MAX;
double y1 = (X * rand()) / RAND_MAX;
double y2 = (X * rand()) / RAND_MAX;
double y3 = (X * rand()) / RAND_MAX;
double y4 = (X * rand()) / RAND_MAX;
double z1 = (X * rand()) / RAND_MAX;
double z2 = (X * rand()) / RAND_MAX;
double z3 = (X * rand()) / RAND_MAX;
double z4 = (X * rand()) / RAND_MAX;
double r12 = std::sqrt(
std::pow(x1 - x2, 2.0) +
std::pow(y1 - y2, 2.0) +
std::pow(z1 - z2, 2.0));
double r23 = std::sqrt(
std::pow(x2 - x3, 2.0) +
std::pow(y2 - y3, 2.0) +
std::pow(z2 - z3, 2.0));
double r34 = std::sqrt(
std::pow(x3 - x4, 2.0) +
std::pow(y3 - y4, 2.0) +
std::pow(z3 - z4, 2.0));
double r41 = std::sqrt(
std::pow(x4 - x1, 2.0) +
std::pow(y4 - y1, 2.0) +
std::pow(z4 - z1, 2.0));
double r24 = std::sqrt(
std::pow(x2 - x4, 2.0) +
std::pow(y2 - y4, 2.0) +
std::pow(z2 - z4, 2.0));
double u = 0.0;
if (r12 >= a)
u += 1.0;
if (r23 >= a)
u += 1.0;
if (r34 >= a)
u += 1.0;
if (r41 >= a)
u += 1.0;
if (r24 >= a)
u += 1.0;
if (u > 0.0)
I += u;
}
return V * I / N;
}
static double Z42(double X, double a, int N)
{
double I = 0.0, V = X * X * X;
for (int n = 0; n < N; n++)
{
double x1 = (X * rand()) / RAND_MAX;
double x2 = (X * rand()) / RAND_MAX;
double x3 = (X * rand()) / RAND_MAX;
double x4 = (X * rand()) / RAND_MAX;
double y1 = (X * rand()) / RAND_MAX;
double y2 = (X * rand()) / RAND_MAX;
double y3 = (X * rand()) / RAND_MAX;
double y4 = (X * rand()) / RAND_MAX;
double z1 = (X * rand()) / RAND_MAX;
double z2 = (X * rand()) / RAND_MAX;
double z3 = (X * rand()) / RAND_MAX;
double z4 = (X * rand()) / RAND_MAX;
double r12 = std::sqrt(
std::pow(x1 - x2, 2.0) +
std::pow(y1 - y2, 2.0) +
std::pow(z1 - z2, 2.0));
double r23 = std::sqrt(
std::pow(x2 - x3, 2.0) +
std::pow(y2 - y3, 2.0) +
std::pow(z2 - z3, 2.0));
double r34 = std::sqrt(
std::pow(x3 - x4, 2.0) +
std::pow(y3 - y4, 2.0) +
std::pow(z3 - z4, 2.0));
double r41 = std::sqrt(
std::pow(x4 - x1, 2.0) +
std::pow(y4 - y1, 2.0) +
std::pow(z4 - z1, 2.0));
double r24 = std::sqrt(
std::pow(x2 - x4, 2.0) +
std::pow(y2 - y4, 2.0) +
std::pow(z2 - z4, 2.0));
double r31 = std::sqrt(
std::pow(x3 - x1, 2.0) +
std::pow(y3 - y1, 2.0) +
std::pow(z3 - z1, 2.0));
double u = 0.0;
if (r12 >= a)
u += 1.0;
if (r23 >= a)
u += 1.0;
if (r34 >= a)
u += 1.0;
if (r41 >= a)
u += 1.0;
if (r24 >= a)
u += 1.0;
if (r31 >= a)
u += 1.0;
if (u > 0.0)
I += u;
}
return V * I / N;
}
int main()
{
double pi = 4.0 * atan(1.0);
// radii of Noble Gases in picometers (10^-12 m)
const double Radii[6] = { 31, 38, 71, 88, 108, 120 };
// spherically symmetric atom's name
const std::string Names[6] = {
"Helium", "Neon", "Argon", "Krypton", "Xenon", "Radon" };
std::cout << "Hard Sphere Second Virial Coefficients for Noble Gases";
std::cout << std::endl;
std::cout << "Monte Carlo Integration is Used";
std::cout << std::endl;
std::cout << "B2(T)" << std::endl;
std::cout << "Atom\tExperiment\tTheory" << std::endl;
srand(1);
for (int i = 0; i < 6; i++)
{
double a = 2.0 * Radii[i] * 1.0e-12;
double X = 25.0 * 100.0 * a;
double b2 = B2(X, a, 1000000);
double bt = 2.0 * pi * std::pow(a, 3.0) / 3.0;
std::cout << Names[i] << "\t";
std::cout << std::scientific;
std::cout << std::setprecision(6);
std::cout << std::right;
std::cout << b2 << '\t';
std::cout << bt << std::endl;
}
std::cout << "Hard Sphere Configuration Integrals";
std::cout << std::endl;
std::cout << "Atom\tZ30\t\tZ40\t\tZ41\t\tZ42";
std::cout << std::endl;
for (int i = 0; i < 6; i++)
{
double a = 2.0 * Radii[i] * 1.0e-12;
double X = 25.0 * 100.0 * a;
double z30 = Z30(X, a, 1000000);
double z40 = Z40(X, a, 1000000);
double z41 = Z41(X, a, 1000000);
double z42 = Z42(X, a, 1000000);
std::cout << Names[i] << '\t';
std::cout << std::scientific;
std::cout << std::setprecision(6);
std::cout << std::right;
std::cout << z30 << '\t';
std::cout << z40 << '\t';
std::cout << z41 << '\t';
std::cout << z42 << '\t';
std::cout << std::endl;
}
return 0;
}
Blog Entry © Friday, October 23, 2026, by James Pate Williams, Jr. Equation of State Calculations
// 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
// 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 © Sunday, November 9, 2025, by James Pate Williams, Jr. Hydrogenic Wavefunctions, Radial Probability Functions, Distribution Functions, and First Moment Integrals
Blog Entry © Sunday, August 10, 2025, First-Order Perturbation Treatment of the Helium Atom by James Pate Williams, Jr.
Blog Entry © Tuesday, July 29, 2025, Double and Triple Monte Carlo Integration by James Pate Williams, Jr.
#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;
}
