// 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>
using namespace std;
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;
vector<double> Nm(65);
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;
}
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()
{
double r = 1.0 / nRows;
double c = 1.0 / nCols;
int collisions = 0;
for (int n = 0; !collisions && n < N; n++)
{
double xj = particlePositions[n].x;
double yj = particlePositions[n].y;
double xk = 100.0 * r * Random11() + xj;
double yk = 100.0 * c * Random11() + yj;
double rr =
pow(xj - xk, 2.0) +
pow(yj - yk, 2.0);
if (rr < d)
collisions++;
if (!collisions)
{
particlePositions[n].x = xk;
particlePositions[n].y = yk;
}
}
return collisions;
}
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;
}
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;
if (j == 2 && k == 3 && rjk < d)
a45 += 1.0;
if (j == 3 && k == 4 && rjk < d)
a45 += 1.0;
if (j == 4 && k == 1 && rjk < d)
a45 += 1.0;
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 f24 = 0.0, f31 = 0.0;
double f14 = 0.0, f34 = 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;
if (j == 2 && k == 3 && rjk < d)
a55 += 1.0;
if (j == 3 && k == 4 && rjk < d)
a55 += 1.0;
if (j == 4 && 1 == 5 && rjk < d)
a55 += 1.0;
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 a55p = 0.0;
for (int i = 0; i < numberSteps; i++)
{
double f12 = 0.0, f23 = 0.0;
double f24 = 0.0, f31 = 0.0;
double f14 = 0.0, f34 = 0.0;
double f53 = 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)
a55p += 1.0;
if (j == 2 && k == 3 && rjk < d)
a55p += 1.0;
if (j == 3 && k == 4 && rjk < d)
a55p += 1.0;
if (j == 4 && 1 == 5 && rjk < d)
a55p += 1.0;
if (j == 5 && k == 1 && rjk < d)
a55p += 1.0;
}
a55p /= (pi * pi * pow(d, 4.0));
a55p /= numberSteps;
return a55p;
}
static double A56pp()
{
double a55pp = 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)
a55pp += 1.0;
if (j == 2 && k == 3 && rjk < d)
a55pp += 1.0;
if (j == 4 && k == 1 && rjk < d)
a55pp += 1.0;
if (j == 3 && 1 == 5 && rjk < d)
a55pp += 1.0;
if (j == 5 && k == 1 && rjk < d)
a55pp += 1.0;
}
a55pp /= (pi * pi * pow(d, 4.0));
a55pp /= numberSteps;
return a55pp;
}
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, f35 = 0.0;
double f45 = 0.0, f51 = 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;
if (j == 2 && k == 3 && rjk < d)
a57p += 1.0;
if (j == 3 && k == 4 && rjk < d)
a57p += 1.0;
if (j == 3 && k == 5 && rjk < d)
a57p += 1.0;
if (j == 5 && k == 1 && rjk < d)
a57p += 1.0;
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;
if (j == 2 && k == 3 && rjk < d)
a57pp += 1.0;
if (j == 3 && k == 4 && rjk < d)
a57pp += 1.0;
if (j == 4 && k == 5 && rjk < d)
a57pp += 1.0;
if (j == 5 && k == 1 && rjk < d)
a57pp += 1.0;
if (j == 3 && k == 1 && rjk < d)
a57pp += 1.0;
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;
if (j == 2 && k == 3 && rjk < d)
a57ppp += 1.0;
if (j == 3 && k == 4 && rjk < d)
a57ppp += 1.0;
if (j == 4 && k == 1 && rjk < d)
a57ppp += 1.0;
if (j == 3 && k == 5 && rjk < d)
a57ppp += 1.0;
if (j == 5 && k == 1 && rjk < d)
a57ppp += 1.0;
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;
if (j == 2 && k == 3 && rjk < d)
a58p += 1.0;
if (j == 3 && k == 4 && rjk < d)
a58p += 1.0;
if (j == 4 && k == 5 && rjk < d)
a58p += 1.0;
if (j == 5 && k == 1 && rjk < d)
a58p += 1.0;
if (j == 3 && k == 1 && rjk < d)
a58p += 1.0;
if (j == 3 && k == 5 && rjk < d)
a58p += 1.0;
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;
if (j == 2 && k == 3 && rjk < d)
a58pp += 1.0;
if (j == 3 && k == 4 && rjk < d)
a58pp += 1.0;
if (j == 4 && k == 5 && rjk < d)
a58pp += 1.0;
if (j == 5 && k == 1 && rjk < d)
a58pp += 1.0;
if (j == 3 && k == 1 && rjk < d)
a58pp += 1.0;
if (j == 3 && k == 5 && rjk < d)
a58pp += 1.0;
if (j == 4 && k == 1 && rjk < d)
a58pp += 1.0;
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;
if (j == 2 && k == 3 && rjk < d)
a59 += 1.0;
if (j == 3 && k == 4 && rjk < d)
a59 += 1.0;
if (j == 4 && k == 5 && rjk < d)
a59 += 1.0;
if (j == 5 && k == 1 && rjk < d)
a59 += 1.0;
if (j == 3 && k == 1 && rjk < d)
a59 += 1.0;
if (j == 3 && k == 5 && rjk < d)
a59 += 1.0;
if (j == 4 && k == 1 && rjk < d)
a59 += 1.0;
if (j == 2 && k == 5 && rjk < d)
a59 += 1.0;
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 = 0.0, nu = 0.0;
int N12 = 0;
cout << "K = ";
cin >> K;
cout << "nu = ";
cin >> nu;
InitializeData(K, nu);
for (int i = 0; i < 10000; i++)
{
if (Move() == 0)
N12++;
}
cout << "N(1/2) = ";
cout << N12 << endl;
double C0 = 64.0 * N12 / (N * N * (K * K - 1.0));
double A0 = 0.5 * sqrt(3.0) * d0 * d0 * N;
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 a59 = 0.0, a510 = 0.0;
cout << "C0 = " << C0 << endl;
cout << "A0 = " << A0 << '\t';
cout << "1 / A0 = " << 1.0 / A0 << endl;
for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
a33 += A33();
for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
a44 += A44();
for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
a45 += A45();
for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
a46 += A46();
for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
a55 += A55();
for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
a56p += A56p();
for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
a56pp += A56pp();
for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
a57p += A57p();
for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
a57pp += A57pp();
for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
a57ppp += A57ppp();
for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
a58p += A58p();
for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
a58pp += A58pp();
for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
a59 += A59();
for (numberSteps = 20000; numberSteps <= 80000; numberSteps += 20000)
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;
cout << "A(3, 3) = " << a33 << endl;
cout << "A(4, 4) = " << a44 << endl;
cout << "A(4, 5) = " << a45 << endl;
cout << "A(4, 6) = " << a46 << endl;
cout << "A(5, 5) = " << a55 << endl;
cout << "A(5, 6') = " << a56p << endl;
cout << "A(5, 6'') = " << a56pp << endl;
cout << "A(5, 7') = " << a57p << endl;
cout << "A(5, 7'') = " << a57pp << endl;
cout << "A(5, 7''') = " << a57ppp << endl;
cout << "A(5, 8') = " << a58p << endl;
cout << "A(5, 8'') = " << a58pp << endl;
cout << "A(5, 9) = " << a59 << endl;
cout << "A(5, 10) = " << a510 << endl;
double C1 = pi / sqrt(3.0);
double C2 = N12 / ((4.0 * pi * pi * A33()) / 9.0);
double C3a = (6.0 * A45() - 3.0 * A44() - A46()) / pow(3.0, 1.5);
double C3b = N12 / (pi * pi * pi * C3a);
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 = C4a + C4b + C4c + C4d;
C4 = C40 * N12 / C4;
cout << "C(1) = " << C1 << endl;
cout << "C(2) = " << 2.0 * C2 << endl;
cout << "C(3) = " << 9.0 * C3b << endl;
cout << "C(4) = " << C4 << endl;
double X3 = C1 * A0 + 2.0 * C2 * A0 * A0 + 9.0 * C3b * A0 * A0 * A0;
X3 += C4 * pow(A0, 4.0);
X3 += 1.0;
cout << "X(3) = " << X3 << endl;
return 0;
}
Tag: C/C++
The programming language I utilized.
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;
}
Blog Entry © Tuesday, October 6, 2026, by James Pate Williams, Jr. Radial Basis Function Networks Results
Blog Entry © Sunday, October 4, 2026, by James Pate Williams, Jr. Riemannian Integration Tests
// Integration1dTests.cpp
// Copyright Saturday, October 3, 2026
// James Pate Williams, Jr. BA, BS, MSwE, PhD
// One-dimensional integration tests
// Our canonical integrands are x^3 exp(x) for
// x = 0 to x = 1 and x^3 sin(x) for x = 0 to x = pi
#include <iomanip>
#include <iostream>
using namespace std;
const int steps = 1024;
static double f1(double x)
{
return x * x * x * exp(x);
}
static double f2(double x)
{
return x * x * x * sin(x);
}
static double g1()
{
return -2.0 * exp(1.0) + 6;
}
static double g2()
{
double pi = 4.0 * atan(1.0);
return -pi * pi * pi * cos(pi) + 6.0 * pi * cos(pi);
}
double TrapezoidalRule(
int n, double a, double b, double(*fx)(double))
{
double h = (b - a) / n;
double s = 0.5 * (fx(a) + fx(b));
double x = a + h;
for (int i = 1; i < n; i++)
{
s += fx(x);
x += h;
}
return h * s;
}
double SimpsonsRule(
int n, double a, double b, double(*fx)(double))
{
double h = (b - a) / n;
double h2 = 2.0 * h;
double s = 0.0;
double t = 0.0;
double x = a + h;
for (int i = 1; i < n; i += 2)
{
s += fx(x);
x += h2;
}
x = a + h2;
for (int i = 2; i < n; i += 2)
{
t += fx(x);
x += h2;
}
return h * (fx(a) + 4 * s + 2 * t + fx(b)) / 3.0;
}
static void PrintResult(double x)
{
cout << setw(15) << fixed;
cout << right << setprecision(12);
cout << setfill(' ');
cout << x;
}
static void PrintError(double x, double exact)
{
cout << '\t';
cout << setw(15) << fixed;
cout << right << setprecision(12);
cout << setfill(' ');
cout << 100.0 * fabs(x - exact) /
fabs(exact) << endl;
}
int main()
{
double gAct1 = g1();
double gAct2 = g2();
double pi = 4.0 * atan(1.0);
double trap1 = TrapezoidalRule(steps, 0.0, 1.0, f1);
double simp1 = SimpsonsRule(steps, 0.0, 1.0, f1);
double trap2 = TrapezoidalRule(steps, 0.0,pi, f2);
double simp2 = SimpsonsRule(steps, 0.0, pi, f2);
cout << "Functions are: " << endl;
cout << "f1(x) = x ^ 3 * exp(x)" << endl;
cout << "f2(x) = x ^ 3 * sin(x)" << endl;
cout << "Trapezoidal and Simpson Use 1024 Steps" << endl;
cout << " Value % Error";
cout << endl;
cout << "trap1 = ";
PrintResult(trap1);
PrintError(trap1, gAct1);
cout << "simp1 = ";
PrintResult(simp1);
PrintError(simp1, gAct1);
cout << "true1 = ";
PrintResult(gAct1);
cout << endl;
cout << "trap2 = ";
PrintResult(trap2);
PrintError(trap2, gAct2);
cout << "simp2 = ";
PrintResult(simp2);
PrintError(simp2, gAct2);
cout << "true2 = ";
PrintResult(gAct2);
cout << endl;
return 0;
}
Blog Entry © Saturday, October 3, 2026, by James Pate Williams, Jr. Planck’s Law Black Body Computations (C/C++ Code Translated from a 2021 C# App)
Blog Entry © Thursday, October 1, 2026, Classical Ballistics by James Pate Williams, Jr.
Blog Entry © Monday, September 21, 2026, by James Pate Williams, Jr. Using an Evolutionary Hill Climber for Function Optimization
// Minima.cpp : This file contains the 'main' function.
// Program execution begins and ends there.
// Copyright (c) Sunday, September 20, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Evolutionary hill climber (c) 1999 by Pate
#include <iostream>
#include <time.h>
#include "Algorithms.h"
int main()
{
while (true)
{
double minimum = 0.0, x0 = 0.0, y0 = 0.0;
double f1 = 0.0, x = 0.0, y = 0.0;
int evaluations = 0, maxIterations = 10000;
int generation = 0, population = 0, option = 0;
unsigned int seed = 1;
std::cout << "== Menu == \r\n";
std::cout << "1 Rastrigin Function\r\n";
std::cout << "2 Ackley Function\r\n";
std::cout << "3 Sphere Function\r\n";
std::cout << "4 Rosenbrock Function\r\n";
std::cout << "5 Beale Function\r\n";
std::cout << "6 Exit\r\n";
std::cout << "Option (1 - 6): ";
std::cin >> option;
if (option == 6)
break;
std::cout << "population = ";
std::cin >> population;
std::cout << "generation = ";
std::cin >> generation;
std::cout << "seed = ";
std::cin >> seed;
srand(seed);
x0 = y0 = 10.0;
clock_t clock0 = clock();
if (option == 1)
{
Algorithms::HillClimber(
population,
generation,
Algorithms::Rastrigin,
x, y, minimum,
evaluations,
seed);
}
else if (option == 2)
{
Algorithms::HillClimber(
population,
generation,
Algorithms::Ackley,
x, y, minimum,
evaluations,
seed);
}
else if (option == 3)
{
Algorithms::HillClimber(
population,
generation,
Algorithms::Sphere,
x, y, minimum,
evaluations,
seed);
}
else if (option == 4)
{
Algorithms::HillClimber(
population,
generation,
Algorithms::Rosenbrock,
x, y, minimum,
evaluations,
seed);
}
else if (option == 5)
{
Algorithms::HillClimber(
population,
generation,
Algorithms::Beale,
x, y, minimum,
evaluations,
seed);
}
clock_t clock1 = clock();
double runtime = (double)(clock1 - clock0) /
CLOCKS_PER_SEC;
std::cout << "x = " << x << std::endl;
std::cout << "y = " << y << std::endl;
std::cout << "f = " << minimum << std::endl;
std::cout << "t = " << runtime << " seconds";
std::cout << std::endl;
}
return 0;
}
#pragma once
class Algorithms
{
public:
static double Rastrigin(double x, double y);
static double Ackley(double x, double y);
static double Sphere(double x, double y);
static double Rosenbrock(double x, double y);
static double Beale(double x, double y);
static void HillClimber(
int population,
int generation,
double(*f)(double, double),
double& x, double& y,
double& minimum,
int& evaluations,
unsigned int seed);
};
#include "Algorithms.h"
#include <vector>
double Algorithms::Rastrigin(double x, double y)
{
double A = 10.0, pi = 4.0 * atan(1.0);
return
2.0 * A +
x * x - A * cos(2.0 * pi * x) +
y * y - A * cos(2.0 * pi * y);
}
double Algorithms::Ackley(double x, double y)
{
double pi = 4.0 * atan(1.0);
double xy = sqrt(0.5 * (x * x + y * y));
double cx = cos(2.0 * pi * x);
double cy = cos(2.0 * pi * y);
double expf1 = -20.0 * exp(-0.2 * xy);
double expf2 = -exp(0.5 * (cx + cy));
return expf1 + expf2 + exp(1.0) + 20.0;
}
double Algorithms::Sphere(double x, double y)
{
return x * x + y * y;
}
double Algorithms::Rosenbrock(double x, double y)
{
return 100.0 * pow(y - x * x, 2.0) + pow(1.0 - x, 2.0);
}
double Algorithms::Beale(double x, double y)
{
return
pow(1.500 - x + x * y, 2.0) +
pow(2.250 - x + x * y * y, 2.0) +
pow(2.625 - x + x * y * y * y, 2.0);
}
void Algorithms::HillClimber(
int population,
int generation,
double(*f)(double, double),
double& x, double& y,
double& minimum,
int& evaluations,
unsigned int seed)
{
srand(seed);
std::vector<double> fitness(population);
std::vector<std::vector<double>> genes(
population, std::vector<double>(2));
for (int i = 0; i < population; i++)
{
genes[i][0] = (10.0 * rand()) / RAND_MAX;
genes[i][1] = (10.0 * rand()) / RAND_MAX;
fitness[i] = f(genes[i][0], genes[i][1]);
}
for (int i = 0; i < generation; i++)
{
int parent1 = rand() % population;
int parent2 = rand() % population;
int parent0 = parent1 > parent2 ? parent1 : parent2;
int sign = rand() % 2 == 0 ? -1 : +1;
double childGeneX = (10.0 * rand()) / RAND_MAX;
double childGeneY = (10.0 * rand()) / RAND_MAX;
double childFitness = f(childGeneX, childGeneY);
double maxFitness = fitness[0];
for (int i = 1; i < population; i++)
if (fitness[i] > maxFitness)
maxFitness = fitness[i];
std::vector<int> maxIndex;
for (int i = 0; i < population; i++)
if (fitness[i] == maxFitness)
maxIndex.push_back(i);
int replace = maxIndex[rand() % maxIndex.size()];
fitness[replace] = childFitness;
genes[replace][0] = childGeneX;
genes[replace][1] = childGeneY;
}
evaluations = generation + population;
double minFitness = fitness[0];
for (int i = 1; i < population; i++)
if (fitness[i] < minFitness)
minFitness = fitness[i];
std::vector<int> minIndex;
for (int i = 0; i < population; i++)
if (fitness[i] == minFitness)
minIndex.push_back(i);
int replace = minIndex[rand() % minIndex.size()];
minimum = fitness[replace];
x = genes[replace][0];
y = genes[replace][1];
}
