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 (c) Saturday, October 31, 2026, by James Pate Williams, Jr. Simple Maclaurin and Taylor Series Expansions
// MaclaurinSeries.cpp : This file contains the 'main' function.
// Program execution begins and ends there.
// Copyright (c) Saturday, October 31, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
#include <cmath>
#include <iomanip>
#include <iostream>
// iterative n-factorial function
static double iFactorial(int n)
{
double factorial = 1.0;
for (int i = 2; i <= n; i++)
factorial *= i;
return factorial;
}
static double Series1(double x, int n)
{
double x2 = x * x, x4 = x2 * x2, x6 = x2 * x4;
switch (n)
{
case 0:
return 1.0;
case 1:
return 1.0 - x2;
case 2:
return 1.0 - x2 + 9 * x4 / iFactorial(4);
case 3:
return 1.0 - x2 + 9 * x4 / iFactorial(4) - 225 * x6 / iFactorial(6);
case 4:
return 1.0 - x2 + 9 * x4 / iFactorial(4) - 225 * x6 / iFactorial(6) +
10825 * x6 * x2 / iFactorial(8);
default:
return 0.0;
}
}
static double Series2(double x, int n)
{
double y = x - 1;
double y2 = y * y, y3 = y * y2, y4 = y * y3, y5 = y * y4;
double y6 = y * y5, y7 = y * y6, y8 = y * y7, y9 = y * y8;
double y10 = y * y9;
switch (n)
{
case 0:
return log(1.0);
case 1:
return log(1.0) + y;
case 2:
return log(1.0) + y - y2 / iFactorial(2);
case 3:
return log(1.0) + y - y2 / iFactorial(2) + 2.0 * y3 / iFactorial(3);
case 4:
return log(1.0) + y - y2 / iFactorial(2) + 2.0 * y3 / iFactorial(3)
- 6.0 * y4 / iFactorial(4);
case 5:
return log(1.0) + y - y2 / iFactorial(2) + 2.0 * y3 / iFactorial(3)
- 6.0 * y4 / iFactorial(4) + 24.0 * y5 / iFactorial(5);
case 6:
return log(1.0) + y - y2 / iFactorial(2) + 2.0 * y3 / iFactorial(3)
- 6.0 * y4 / iFactorial(4) + 24.0 * y5 / iFactorial(5)
- 120.0 * y6 / iFactorial(6);
case 7:
return log(1.0) + y - y2 / iFactorial(2) + 2.0 * y3 / iFactorial(3)
- 6.0 * y4 / iFactorial(4) + 24.0 * y5 / iFactorial(5)
- 120.0 * y6 / iFactorial(6) + 720.0 * y7 / iFactorial(7);
case 8:
return log(1.0) + y - y2 / iFactorial(2) + 2.0 * y3 / iFactorial(3)
- 6.0 * y4 / iFactorial(4) + 24.0 * y5 / iFactorial(5)
- 120.0 * y6 / iFactorial(6) + 720.0 * y7 / iFactorial(7)
- 5040.0 * y8 / iFactorial(8);
case 9:
return log(1.0) + y - y2 / iFactorial(2) + 2.0 * y3 / iFactorial(3)
- 6.0 * y4 / iFactorial(4) + 24.0 * y5 / iFactorial(5)
- 120.0 * y6 / iFactorial(6) + 720.0 * y7 / iFactorial(7)
- 5040.0 * y8 / iFactorial(8) + 40320 * y9 / iFactorial(9);
case 10:
return log(1.0) + y - y2 / iFactorial(2) + 2.0 * y3 / iFactorial(3)
- 6.0 * y4 / iFactorial(4) + 24.0 * y5 / iFactorial(5)
- 120.0 * y6 / iFactorial(6) + 720.0 * y7 / iFactorial(7)
- 5040.0 * y8 / iFactorial(8) + 40320 * y9 / iFactorial(9)
- 362880 * y10 / iFactorial(10);
default:
{
double sign = 1.0, sum = log(1.0);
for (int i = 1; i <= n; i++)
{
sum += sign * pow(y, i) / i;
sign *= -1.0;
}
return sum;
}
}
}
int main()
{
std::cout << "n" << '\t' << "n!" << std::endl;
for (int n = 0; n <= 10; n++)
std::cout << n << '\t' << iFactorial(n) << std::endl;
double x = 0.015;
std::cout << std::fixed << std::setprecision(10);
std::cout << "cos(arctan(x)) = (1 + x^2)^-1/2" << std::endl;
std::cout << "Maclaurin Series Expansion of (1 + x^2)^-1/2" << std::endl;
std::cout << "# is the number + 1 of nonzero terms in the summation\r\n";
std::cout << "x = " << x << std::endl;
std::cout << "#" << '\t' << "Series" << "\t\t" << "(1 + x^2)^-1/2";
std::cout << std::endl;
for (int n = 0; n <= 4; n++)
{
std::cout << n << '\t';
std::cout << std::fixed << std::setprecision(10);
std::cout << Series1(x, n) << '\t';
std::cout << pow(1.0 + x * x, -0.5) << std::endl;
}
x = 0.125;
std::cout << "Maclaurin Series Expansion of (1 + x^2)^-1/2" << std::endl;
std::cout << "x = " << x << std::endl;
std::cout << "#" << '\t' << "Series" << "\t\t" << "(1 + x^2)^-1/2";
std::cout << std::endl;
for (int n = 0; n <= 4; n++)
{
std::cout << n << '\t';
std::cout << std::fixed << std::setprecision(10);
std::cout << Series1(x, n) << '\t';
std::cout << pow(1.0 + x * x, -0.5) << std::endl;
}
x = 0.5;
std::cout << "Maclaurin Series Expansion of (1 + x^2)^-1/2" << std::endl;
std::cout << "x = " << x << std::endl;
std::cout << "#" << '\t' << "Series" << "\t\t" << "(1 + x^2)^-1/2";
std::cout << std::endl;
for (int n = 0; n <= 4; n++)
{
std::cout << n << '\t';
std::cout << std::fixed << std::setprecision(10);
std::cout << Series1(x, n) << '\t';
std::cout << pow(1.0 + x * x, -0.5) << std::endl;
}
x = 1.25;
std::cout << "Taylor Series Expansion of log(x)" << std::endl;
std::cout << "About x = 1" << std::endl;
std::cout << "x = " << x << std::endl;
std::cout << "#" << '\t' << "Series" << "\t\t" << "ln(x)";
std::cout << std::endl;
for (int n = 0; n <= 15; n++)
{
std::cout << n << '\t';
std::cout << std::fixed << std::setprecision(10);
std::cout << Series2(x, n) << '\t';
std::cout << std::log(x) << std::endl;
}
x = 1.5;
std::cout << "Taylor Series Expansion of log(x)" << std::endl;
std::cout << "About x = 1" << std::endl;
std::cout << "x = " << x << std::endl;
std::cout << "#" << '\t' << "Series" << "\t\t" << "ln(x)";
std::cout << std::endl;
for (int n = 0; n <= 30; n++)
{
std::cout << n << '\t';
std::cout << std::fixed << std::setprecision(10);
std::cout << Series2(x, n) << '\t';
std::cout << std::log(x) << std::endl;
}
return 0;
}
Blog Entry (c) Friday, October 30, 2026, by James Pate Williams, Jr. One-Component Classical Gas Second Virial Coefficient Calculated Using the Lennard-Jones Potential
Blog Entry (c) Wednesday, October 28, 2016, by James Pate Williams, Jr. One-Component Classical Monatomic Gas Second Virial Coefficient Using the Lennard-Jones Potential
// LennardJonesB2.cpp : This file contains the 'main' function.
// Program execution begins and ends there.
// Copyright (c) October 28, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
#include <cmath>
#include <iomanip>
#include <iostream>
double gEpsilonKT = 0.0; // e / kT
double gRStar = 0.0; // r*
int gnSteps = 0; // number of Simpson's Rule steps
static double SimpsonsRule(
double a, double b, double(*fx)(double), int n)
{
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 double LennardJones(double epsilonKT, double y)
{
// epsilon = energy / kT
return exp(2.0 * epsilonKT * pow(y, -6.0) -
epsilonKT * pow(y, -12.0));
}
static double Integrand(double y)
{
// y = r / r*
double lj = LennardJones(1.0 / gEpsilonKT, y);
return -2.0 * 4.0 * atan(1.0) * pow(gRStar, 3.0) *
(lj - 1.0) * y * y;
}
static void PrintIntegral(double integral)
{
std::cout << std::fixed;
if (integral < 0.0)
{
std::cout << std::setprecision(10);
std::cout << integral;
}
else
{
std::cout << std::setprecision(11);
std::cout << ' ' << integral;
}
}
int main()
{
// assign global variables
gEpsilonKT = 1.0;
gRStar = 1.0;
gnSteps = 25000;
std::cout << "kT/e\tIntegral 1\tIntegral 2\t";
std::cout << "Integral 3" << std::endl;
while (gEpsilonKT <= 25.0)
{
double integral1 = SimpsonsRule(0.1, 25.0, Integrand, gnSteps);
double integral2 = SimpsonsRule(0.1, 50.0, Integrand, gnSteps);
double integral3 = SimpsonsRule(0.1, 75.0, Integrand, gnSteps);
std::cout << std::setprecision(2) << std::setw(2);
std::cout << static_cast<int>(gEpsilonKT) << '\t';
PrintIntegral(integral1);
std::cout << '\t';
PrintIntegral(integral2);
std::cout << '\t';
PrintIntegral(integral3);
std::cout << std::endl;
gEpsilonKT += 1.0;
}
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 © Tuesday, October 20, 2026, by James Pate Williams, Jr., Three-Dimensional Heat Conduction Problem
// HeatConductionCube.cpp : This file contains the 'main' function.
// Program execution begins and ends there.
// (c) Monday, October 19, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// References:
// https://mathworld.wolfram.com/HeatConductionEquation.html
// https://www.wolframalpha.com/calculators/mathematics-calculus-integrals-definite-integral-calculator
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>
using namespace std;
double a = 2, b = 3, c = 4, kappa = 1.0;
int l = 5, m = 5, n = 5, N = 5;
vector<double> C, E, G, U;
static double U0(
double a, double b, double c,
double x, double y, double z,
vector<double> C, vector<double> E,
vector<double> G, int nn)
{
double pi = 4.0 * atan(1.0), sum = 0.0;
double fl = 0.0, fm = 0.0, fn = 0.0;
int l = 0, m = 0, n = 0;
for (l = 1; l <= nn; l++)
{
fl = C[l] * sin(l * pi * x / a);
for (m = 1; m <= nn; m++)
{
fm = E[m] * sin(m * pi * y / b);
for (n = 1; n <= nn; n++)
{
fn = G[n] * sin(n * pi * z / c);
sum += fl * fm * fn;
}
}
}
return sum;
}
static double DefiniteIntegral(double a, int p)
{
double pi = 4.0 * atan(1.0), pi2 = pi * pi;
double pi3 = pi * pi2, a3 = a * a * a;
double integral = 0.0;
double p2 = p * p, p3 = p * p2;
integral = (a3 * (2.0 - pi2 * p2) * cos(pi * p) +
2.0 * pi * p * sin(pi * p) - 2.0) / (pi3 * p3);
return integral;
}
static void ComputeVectors()
{
double x = 0.0, y = 0.0, z = 0.0;
double xs = a / N, ys = b / N, zs = c / N;
C.resize(N + 1LL);
E.resize(N + 1LL);
G.resize(N + 1LL);
for (int i = 1; i <= N; i++)
{
for (int j = 1; j <= N; j++)
{
for (int k = 1; k <= N; k++)
{
for (int p = 1; p <= N; p++)
C[p] = DefiniteIntegral(a, p);
for (int q = 1; q <= N; q++)
E[q] = DefiniteIntegral(b, q);
for (int r = 1; r <= N; r++)
G[r] = DefiniteIntegral(c, r);
z += zs;
}
y += ys;
}
x += xs;
}
}
static void ComputeU()
{
double pi = 4.0 * atan(1.0);
double ts = a / N, xs = a / N, ys = b / N, zs = c / N;
double t = ts, x = xs, y = ys, z = zs;
U.clear();
for (int h = 1; h <= N; h++)
{
for (int i = 1; i <= N; i++)
{
for (int j = 1; j <= N; j++)
{
for (int k = 1; k <= N; k++)
{
for (int l = 1; l <= N; l++)
{
for (int m = 1; m <= N; m++)
{
for (int n = 1; n <= N; n++)
{
double CEG = C[l] * E[m] * G[n];
double sl = sin(l * pi * x / a);
double sm = sin(m * pi * y / b);
double sn = sin(n * pi * z / c);
double el = exp(-kappa * pow(l * pi / a, 2.0) * t);
double em = exp(-kappa * pow(m * pi / b, 2.0) * t);
double en = exp(-kappa * pow(n * pi / c, 2.0) * t);
double uu = CEG * sl * sm * sn * el * em * en;
U.push_back(uu);
if (fabs(uu) < 1.0e-15)
return;
}
}
}
x += ys;
}
y += ys;
}
z += zs;
}
t += ts;
}
}
static void ReadData()
{
cout << "a = ";
cin >> a;
cout << "b = ";
cin >> b;
cout << "c = ";
cin >> c;
cout << "k = ";
cin >> kappa;
cout << "N = ";
cin >> N;
}
static void WriteResults()
{
ComputeVectors();
cout << "C\r\n";
for (int i = 1; i <= N; i++)
{
cout << fixed << right << setw(16);
cout << setprecision(13);
cout << C[i] << endl;
}
cout << "E\r\n";
for (int i = 1; i <= N; i++)
{
cout << fixed << right << setw(16);
cout << setprecision(13);
cout << E[i] << endl;
}
cout << "G\r\n";
for (int i = 1; i <= N; i++)
{
cout << fixed << right << setw(16);
cout << setprecision(13);
cout << G[i] << endl;
}
ComputeU();
cout << "U\r\n";
for (size_t i = 0; i < U.size(); i++)
{
cout << fixed << right << setw(16);
cout << setprecision(13);
cout << U[i] << endl;
}
}
int main()
{
ReadData();
WriteResults();
return 0;
}
Blog Entry @ Thursday, October 15, 2026, by James Pate Williams, Jr. Solution of the Telegraph Equation
// TelegraphEquation.cpp
// Copyright (c) Wednesday, October 14, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
#include <iomanip>
#include <iostream>
#include <vector>
static double omega(
double alpha, double beta, double c, double l, int n)
{
double pi = 4.0 * atan(1.0);
double term1 = pow(n * pi * c / l, 2.0);
double term2 = 0.25 * pow(alpha - beta, 2.0);
return sqrt(term1 - term2);
}
static double phi(
double alpha, double beta, double c, double l, int n)
{
double d = 0.5 * (alpha + beta);
return atan(d / omega(alpha, beta, c, l, n));
}
static double An(
double alpha, double beta, double a, double c, double l, int n)
{
double factor1 = 2.0 / (l, cos(phi(alpha, beta, c, l, n)));
double factor2 = sin(n * 4.0 * atan(1.0) * a / l);
return factor1 * factor2;
}
static double analytic_u(
double alpha, double beta, double a,
double c, double l, double x, double t,
int N)
{
double sum = 0.0;
for (int n = 1; n <= N; n++)
{
double d = 0.5 * (alpha + beta);
double an = An(alpha, beta, a, c, l, n);
double factor1 = an * exp(-d * t);
double omegaN = omega(alpha, beta, c, l, n);
double phiN = phi(alpha, beta, c, l, n);
double factor2 = cos(omegaN * t - phiN);
double factor3 = sin(n * 4.0 * atan(1.0) * x / l);
sum += factor1 * factor2 * factor3;
}
return sum;
}
int main()
{
char line[128] = "";
double alpha = 0.0, beta = 0.0, a = 0.0;
double c = 0.0, l = 0.0;
int N = 0, nPoints = 0;
std::cout << "alpha = ";
std::cin.getline(line, 128);
alpha = atof(line);
std::cout << "beta = ";
std::cin.getline(line, 128);
beta = atof(line);
std::cout << "a = ";
std::cin.getline(line, 128);
a = atof(line);
std::cout << "c = ";
std::cin.getline(line, 128);
c = atof(line);
std::cout << "l = ";
std::cin.getline(line, 128);
l = atof(line);
std::cout << "N = ";
std::cin.getline(line, 128);
N = atoi(line);
std::cout << "nPoints = ";
std::cin.getline(line, 128);
nPoints = atoi(line);
double deltaX = 1.0 / nPoints;
double deltaT = 1.0 / nPoints;
std::vector<double> u(nPoints), x(nPoints), t(nPoints);
for (int i = 0; i < nPoints; i++)
{
x[i] = i * deltaX;
t[i] = i * deltaT;
u[i] = analytic_u(alpha, beta, a, c, l, x[i], t[i], N);
std::cout << std::setprecision(10);
std::cout << std::fixed;
std::cout << x[i] << '\t' << t[i] << '\t';
if (u[i] >= 0.0)
std::cout << ' ';
std::cout << u[i] << std::endl;
}
return 0;
}
Blog Entry © Wednesday, October 14, 2026, by James Pate Williams Jr. Three Pi Related Integrals
// ThreeIntegrals.cpp
// Copyright (c) Sunday, October 11, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Integrals due to Gauss, Cauchy, and Dirichlet
#include <chrono>
#include <iomanip>
#include <iostream>
#include <vector>
static double Gauss(double x)
{
return exp(-x * x);
}
static double Cauchy(double x)
{
return 1.0 / (x * x + 1.0);
}
static double Dirichlet(double x)
{
return sin(x) / x;
}
static double SimpsonsRule(
double a, double b, double(*fx)(double), int nSteps)
{
double h = (b - a) / nSteps;
double h2 = 2.0 * h;
double s = 0.0;
double t = 0.0;
double x = a + h;
for (int i = 1; i < nSteps; i += 2)
{
s += fx(x);
x += h2;
}
x = a + h2;
for (int i = 2; i < nSteps; i += 2)
{
t += fx(x);
x += h2;
}
return h * (fx(a) + 4 * s + 2 * t + fx(b)) / 3.0;
}
static double ComputePi1(double bound)
{
return SimpsonsRule(-bound, bound, Gauss, 16384);
}
static double ComputePi2(double bound, int nSteps)
{
return 2.0 * SimpsonsRule(0.0, bound, Cauchy, nSteps);
}
static double ComputePi3(double bound, int nSteps)
{
double zero = 1.0e-32;
double integral1 = SimpsonsRule(-bound, -zero, Dirichlet, nSteps);
double integral2 = SimpsonsRule(+zero, +bound, Dirichlet, nSteps);
return integral1 + integral2;
}
static double nFactorial(int n)
{
double factorial = 1.0;
for (int i = 1; i <= n; i++)
factorial *= i;
return factorial;
}
typedef struct tagGene {
int R; int nSteps;
} Gene, * PGene;
static double HillClimber(
int fNumber, int generations, int population, Gene& gene)
{
double pi = 4.0 * atan(1.0), pi2 = 0.0, pi3 = 0.0;
std::vector<double> fitness(population);
std::vector<std::vector<Gene>> genome;
genome.resize(population);
for (int i = 0; i < population; i++)
{
genome[i].resize(1);
Gene g = { 0 };
if (fNumber == 2)
{
g.R = rand() % 100000;
g.nSteps = rand() % 100000;
double pi2 = ComputePi2(g.R, g.nSteps);
fitness[i] = fabs(pi - pi2);
if (fitness[i] <= 1.0e-15)
{
gene = g;
return pi2;
}
genome[i][0] = g;
}
else if (fNumber == 3)
{
g.R = rand() % 10000000;
g.nSteps = rand() % 10000000;
double pi3 = ComputePi3(g.R, g.nSteps);
fitness[i] = fabs(pi - pi3);
if (fitness[i] <= 1.0e-15)
{
gene = g;
return pi3;
}
genome[i][0] = g;
}
}
for (int g = 0; g < generations; g++)
{
double childFitness = DBL_MAX;
int parent1 = rand() % population;
int parent2 = rand() % population;
int child = parent1 > parent2 ? parent1 : parent2;
Gene childGene = genome[child][0];
if (fNumber == 2)
{
pi2 = ComputePi2(childGene.R, childGene.nSteps);
childFitness = fabs(pi - pi2);
if (childFitness <= 1.0e-15)
{
gene = childGene;
return childFitness;
}
}
else if (fNumber == 3)
{
pi3 = ComputePi3(childGene.R, childGene.nSteps);
childFitness = fabs(pi - pi3);
if (childFitness <= 1.0e-15)
{
gene = childGene;
return childFitness;
}
}
double maxFitness = fitness[0];
for (int i = 1; i < population; i++)
if (fitness[i] > maxFitness)
maxFitness = fitness[i];
std::vector<int> maxIndices;
for (int i = 0; i < population; i++)
if (fitness[i] == maxFitness)
maxIndices.push_back(i);
int selection = maxIndices[rand() % maxIndices.size()];
fitness[selection] = childFitness;
genome[selection][0] = childGene;
}
double minFitness = fitness[0];
for (int i = 1; i < population; i++)
if (fitness[i] < minFitness)
minFitness = fitness[i];
std::vector<int> minIndices;
for (int i = 0; i < population; i++)
if (fitness[i] == minFitness)
minIndices.push_back(i);
int selection = minIndices[rand() % minIndices.size()];
gene = genome[selection][0];
if (fNumber == 2)
return ComputePi2(gene.R, gene.nSteps);
else if (fNumber == 3)
return ComputePi3(gene.R, gene.nSteps);
else
return 0.0;
}
int main()
{
Gene gene = { 0 };
std::cout << std::fixed << std::setprecision(12);
std::cout << "Function\t\tValue of Pi\tRuntime (mu-s)" << std::endl;
auto start0 = std::chrono::high_resolution_clock::now();
double pi0 = 4.0 * atan(1.0);
auto stop0 = std::chrono::high_resolution_clock::now();
auto duration0 = std::chrono::duration_cast<std::chrono::microseconds>(stop0 - start0);
std::cout << "Pi [4 * atan(1)] = " << pi0 << '\t' << duration0.count() << std::endl;
// Gauss Integral
auto start1 = std::chrono::high_resolution_clock::now();
double pi1 = pow(ComputePi1(6.0e2), 2.0);
auto stop1 = std::chrono::high_resolution_clock::now();
auto duration1 = std::chrono::duration_cast<std::chrono::microseconds>(stop1 - start1);
std::cout << "Gauss Integral = " << pi1 << '\t' << duration1.count() << std::endl;
// Cauchy Integral
auto start2 = std::chrono::high_resolution_clock::now();
double pi2 = HillClimber(2, 200000, 3000, gene);
auto stop2 = std::chrono::high_resolution_clock::now();
auto duration2 = std::chrono::duration_cast<std::chrono::microseconds>(stop2 - start2);
std::cout << "Cauchy Integral = " << pi2 << '\t' << duration2.count() << std::endl;
// Dirichlet integral
auto start3 = std::chrono::high_resolution_clock::now();
double pi3 = HillClimber(3, 50000, 1000, gene);
auto stop3 = std::chrono::high_resolution_clock::now();
auto duration3 = std::chrono::duration_cast<std::chrono::microseconds>(stop3 - start3);
std::cout << "Dirichlet Integral = " << pi3 << '\t' << duration3.count() << std::endl;
return 0;
}
Blog Entry © Sunday, October 11, 2026, by James Pate Williams, Jr.
#pragma once
class Potential
{
public:
static double Slot1SimpsonsRule(
double a, double A, double B, int n, int nSteps,
double(*fx)(double, double, int));
static double SquareExact(
double a, double b,
double x, double y,
int N);
static void ComputeIntegrals(
double a, int N, int nSteps, std::vector<double>& c);
static double Slot1Exact(
double a, double x, double y,
int N, std::vector<double>& c);
static void SquareFixedPoint(
double a, double b,
int Nx, int Ny, int iterations,
std::vector<std::vector<double>>& U);
static void Slot1FixedPoint(
double a, double b,
int Nx, int Ny, int iterations,
std::vector<std::vector<double>>& U);
static double Disk_f(double theta);
static double Disk_a0(double theta);
static double Disk_an(double c, double theta, int n);
static double Disk_bn(double c, double theta, int n);
static double SimpsonsRule_Disk_a0(int nSteps);
static double SimpsonsRule_Disk_an(
double c, int n, int nSteps);
static double SimpsonsRule_Disk_bn(
double c, int n, int nSteps);
static double Disk_u(
double c, double r, double theta,
double a0, std::vector<double>& an,
std::vector<double>& bn, int Nt,
int nSteps);
};
#include "pch.h"
#include "Potential.h"
double Potential::Slot1SimpsonsRule(
double a, double A, double B, int n, int nSteps,
double(*fx)(double, double, int))
{
double h = (B - A) / nSteps;
double h2 = 2.0 * h;
double s = 0.0;
double t = 0.0;
double x = A + h;
for (int i = 1; i < nSteps; i += 2)
{
s += fx(a, x, n);
x += h2;
}
x = A + h2;
for (int i = 2; i < nSteps; i += 2)
{
t += fx(a, x, n);
x += h2;
}
return h * (fx(a, A, n) + 4 * s + 2 * t + fx(a, B, n)) / 3.0;
}
double Potential::SquareExact(
double a, double b,
double x, double y,
int N)
{
double pi = 4.0 * atan(1.0);
double sum = 0.0;
for (int n = 1; n <= N; n++)
{
double l2 = pi * n / a;
double s1 = sin(0.5 * n * pi) / (n * n);
double s2 = sinh(l2 * y);
double s3 = sinh(l2 * (b - y));
double s4 = sinh(l2 * b);
double s5 = sin(l2 * x);
sum += s1 * ((s2 + s3) / s4) * s5;
}
return 8.0 * sum / (pi * pi);
}
double Slot1Integrand(double a, double x, int n)
{
double pi = 4.0 * atan(1.0);
return exp(-x) * sin(n * pi * x / a);
}
void Potential::ComputeIntegrals(
double a, int N, int nSteps, std::vector<double>& c)
{
c[0] = 0.0;
for (int n = 1; n <= N; n++)
{
c[n] = 2.0 * Slot1SimpsonsRule(a, 0.0, a,
n, nSteps, Slot1Integrand) / a;
}
}
double Potential::Slot1Exact(
double a, double x, double y,
int N, std::vector<double>& c)
{
double pi = 4.0 * atan(1.0), uxy = 0.0;
for (int n = 1; n <= N; n++)
uxy += c[n] * sin(n * pi * x / a) * exp(-n * pi * y / a);
return uxy;
}
static double SquareF(double x)
{
if (x > 0.0 && x < 0.5)
return 2.0 * x;
else if (x >= 0.5 && x < 1.0)
return 2.0 - 2.0 * x;
else
return 0.0;
}
void Potential::SquareFixedPoint(
double a, double b,
int Nx, int Ny, int iterations,
std::vector<std::vector<double>>& U)
{
double deltaX = a / Nx, deltaY = b / Ny;
U.resize(Nx + 1LL);
for (int i = 0; i <= Nx; i++)
U[i].resize(Ny + 1LL);
for (int i = 0; i <= Ny; i++)
U[0][i] = 0.0;
for (int i = 1; i <= Ny; i++)
U[Nx][i] = 0.0;
for (int i = 0; i <= Nx; i++)
U[i][0] = SquareF(i * deltaX);
for (int i = 0; i <= Nx; i++)
U[i][Ny] = SquareF(i * deltaX);
for (int m = 0; m < iterations; m++)
{
for (int j = 1; j < Nx; j++)
{
for (int k = 1; k < Ny; k++)
U[j][k] = 0.25 *
(U[j + 1][k] + U[j - 1][k] +
U[j][k + 1] + U[j][k - 1]);
}
}
}
void Potential::Slot1FixedPoint(
double a, double b,
int Nx, int Ny, int iterations,
std::vector<std::vector<double>>& U)
{
double deltaX = a / Nx, deltaY = b / Ny;
U.resize(Nx + 1);
for (int i = 0; i <= Nx; i++)
U[i].resize(Ny + 1, 0.0);
for (int i = 0; i <= Nx; i++)
U[i][0] = exp(-i * deltaX);
for (int i = 0; i < iterations; i++)
{
for (int j = 1; j < Nx; j++)
{
for (int k = 1; k < Ny; k++)
U[j][k] = 0.25 *
(U[j + 1][k] + U[j - 1][k] +
U[j][k + 1] + U[j][k - 1]);
}
}
}
double Potential::Disk_f(double theta)
{
return theta * theta;
}
double Potential::Disk_a0(double theta)
{
double pi = 4.0 * atan(1.0);
return Disk_f(theta);
}
double Potential::Disk_an(double c, double theta, int n)
{
double pi = 4.0 * atan(1.0);
return Disk_f(theta) * cos(n * theta) / (pi * pow(c, n));
}
double Potential::Disk_bn(double c, double theta, int n)
{
double pi = 4.0 * atan(1.0);
return Disk_f(theta) * cos(n * theta) / (pi * pow(c, n));
}
double Potential::SimpsonsRule_Disk_a0(int nSteps)
{
double pi = 4.0 * atan(1.0);
double h = (pi + pi) / nSteps;
double h2 = 2.0 * h;
double s = 0.0;
double t = 0.0;
double theta = -pi + h;
for (int i = 1; i < nSteps; i += 2)
{
s += Disk_a0(theta);
theta += h2;
}
theta = -pi + h2;
for (int i = 2; i < nSteps; i += 2)
{
t += Disk_a0(theta);
theta += h2;
}
return h * (Disk_a0(0.0) + 4 * s + 2 * t + Disk_a0(pi)) / 3.0;
}
double Potential::SimpsonsRule_Disk_an(
double c, int n, int nSteps)
{
double pi = 4.0 * atan(1.0);
double h = (pi + pi) / nSteps;
double h2 = 2.0 * h;
double s = 0.0;
double t = 0.0;
double theta = -pi + h;
for (int i = 1; i < nSteps; i += 2)
{
s += Disk_an(c, theta, n);
theta += h2;
}
theta = -pi + h2;
for (int i = 2; i < nSteps; i += 2)
{
t += Disk_an(c, theta, n);
theta += h2;
}
double integral =
h * (Disk_an(c, 0.0, n) + 4 * s + 2 * t +
Disk_an(c, pi, n)) / 3.0;
return integral / (pi * pow(c, n));
}
double Potential::SimpsonsRule_Disk_bn(
double c, int n, int nSteps)
{
double pi = 4.0 * atan(1.0);
double h = (pi + pi) / nSteps;
double h2 = 2.0 * h;
double s = 0.0;
double t = 0.0;
double theta = -pi + h;
for (int i = 1; i < nSteps; i += 2)
{
s += Disk_bn(c, theta, n);
theta += h2;
}
theta = -pi + h2;
for (int i = 2; i < nSteps; i += 2)
{
t += Disk_bn(c, theta, n);
theta += h2;
}
double integral =
h * (Disk_bn(c, 0.0, n) + 4 * s + 2 * t +
Disk_bn(c, pi, n)) / 3.0;
return integral / (pi * pow(c, n));
}
double Potential::Disk_u(
double c, double r, double theta,
double a0, std::vector<double>& an,
std::vector<double>& bn, int Nt,
int nSteps)
{
double urt = a0;
for (int n = 1; n <= Nt; n++)
{
urt += pow(r, n) * (
an[n] * cos(n * theta) +
bn[n] * sin(n * theta));
}
return urt;
}
