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

Maclaurin Series ExpansionsDownload
// 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

Second Virial CoefficientDownload

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

Blog Entry Wednesday October 28 2026Download
// 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

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 © Tuesday, October 20, 2026, by James Pate Williams, Jr., Three-Dimensional Heat Conduction Problem

Heat Conduction Equation in Cartesian CoordinatesDownload
// 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

Telegraph Equation SolutionDownload
// 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

Blog Entry Wednesday October 14 2026Download
// 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.

Blog Entry Sunday October 11 2026Download
#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;
}

Blog Entry © Tuesday, October 6, 2026, by James Pate Williams, Jr. Radial Basis Function Networks Results

Blog Entry Tuesday October 6 2026Download