Category: Elementary Physics
Blog Entry (c) Monday, November 2, 2026, by James Pate Williams, Jr. Some Astrophysical Calculations
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) Friday, October 30, 2026, by James Pate Williams, Jr. One-Component Classical Gas Second Virial Coefficient Calculated Using the Lennard-Jones Potential
Blog Entry (c) Wednesday, October 28, 2016, by James Pate Williams, Jr. One-Component Classical Monatomic Gas Second Virial Coefficient Using the Lennard-Jones Potential
// LennardJonesB2.cpp : This file contains the 'main' function.
// Program execution begins and ends there.
// Copyright (c) October 28, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
#include <cmath>
#include <iomanip>
#include <iostream>
double gEpsilonKT = 0.0; // e / kT
double gRStar = 0.0; // r*
int gnSteps = 0; // number of Simpson's Rule steps
static double SimpsonsRule(
double a, double b, double(*fx)(double), int n)
{
double h = (b - a) / n;
double h2 = 2.0 * h;
double s = 0.0;
double t = 0.0;
double x = a + h;
for (int i = 1; i < n; i += 2)
{
s += fx(x);
x += h2;
}
x = a + h2;
for (int i = 2; i < n; i += 2)
{
t += fx(x);
x += h2;
}
return h * (fx(a) + 4 * s + 2 * t + fx(b)) / 3.0;
}
static double LennardJones(double epsilonKT, double y)
{
// epsilon = energy / kT
return exp(2.0 * epsilonKT * pow(y, -6.0) -
epsilonKT * pow(y, -12.0));
}
static double Integrand(double y)
{
// y = r / r*
double lj = LennardJones(1.0 / gEpsilonKT, y);
return -2.0 * 4.0 * atan(1.0) * pow(gRStar, 3.0) *
(lj - 1.0) * y * y;
}
static void PrintIntegral(double integral)
{
std::cout << std::fixed;
if (integral < 0.0)
{
std::cout << std::setprecision(10);
std::cout << integral;
}
else
{
std::cout << std::setprecision(11);
std::cout << ' ' << integral;
}
}
int main()
{
// assign global variables
gEpsilonKT = 1.0;
gRStar = 1.0;
gnSteps = 25000;
std::cout << "kT/e\tIntegral 1\tIntegral 2\t";
std::cout << "Integral 3" << std::endl;
while (gEpsilonKT <= 25.0)
{
double integral1 = SimpsonsRule(0.1, 25.0, Integrand, gnSteps);
double integral2 = SimpsonsRule(0.1, 50.0, Integrand, gnSteps);
double integral3 = SimpsonsRule(0.1, 75.0, Integrand, gnSteps);
std::cout << std::setprecision(2) << std::setw(2);
std::cout << static_cast<int>(gEpsilonKT) << '\t';
PrintIntegral(integral1);
std::cout << '\t';
PrintIntegral(integral2);
std::cout << '\t';
PrintIntegral(integral3);
std::cout << std::endl;
gEpsilonKT += 1.0;
}
return 0;
}
Blog Entry © Friday, October 23, 2026, by James Pate Williams, Jr. Equation of State Calculations
// EquationOfState.cpp : This file contains the 'main' function.
// EquationOfState.cpp : This file contains the 'main' function.
// Program execution begins and ends there.
// Copyright (c) Friday, October 23, 2026
// by James Pate Williams, Jr. BA, BS, MSwE, PhD
// Reference: osti.gov/servlets/purl/4390578/
#include <math.h>
#include <iostream>
#include <vector>
const int N = 224;
const int nRows = 14;
const int nCols = 16;
const double d = 1.0 / nRows;
typedef struct tagPoint2d
{
double x, y;
} Point2d, *PPoint2d;
double d0 = 0.0, AA0 = 0.0, AA2 = 0.0, maxD = 0.0, pi = 0.0;
int numberSteps = 0;
std::vector<Point2d> particlePositions(N);
static double Distance(double x1, double y1, double x2, double y2)
{
// 2-D Euclidean norm
return sqrt(pow(x1 - x2, 2.0) + pow(y1 - y2, 2.0));
}
static double Random01()
{
// random number in the interval [0, 1]
double value = static_cast<double>(rand()) / RAND_MAX;
return value;
}
static double Random11()
{
// random number in the interval [-1, 1]
double value = -1.0 + 2.0 * static_cast<double>(rand()) / RAND_MAX;
return value;
}
static void InitializeData(double K, double nu)
{
pi = 4.0 * atan(1.0);
d0 = d * (1.0 - pow(2.0, nu - 8.0));
AA0 = 1.0 / (sqrt(3.0) * d0 * d0 * N / 2.0);
maxD = d - d0;
AA2 = (K * K - 1.0) * pi * d0 * d0 / 64.0;
AA0 = 1.0 / (sqrt(3.0) * d0 * d0 * N / 2.0);
maxD = d - d0;
AA2 = (K * K - 1.0) * pi * d0 * d0 / 64.0;
}
static void InitializeLattice()
{
// initialize to a regular 2-D Cartesian lattice
for (size_t r = 0; r < nRows; r++)
{
for (size_t c = 0; c < nCols; c++)
{
size_t address = r * nRows + c;
particlePositions[address].x = r / 14.0;
particlePositions[address].y = c / 16.0;
}
}
}
static int Move()
{
bool result = true;
double r = 1.0 / nRows;
double c = 1.0 / nCols;
InitializeLattice();
for (int n = 0; n < N; n++)
{
double xj = particlePositions[n].x;
double yj = particlePositions[n].y;
double xk = 1.25 * N * r * Random11() + xj;
double yk = 1.25 * N * c * Random11() + yj;
double rr =
pow(xj - xk, 2.0) +
pow(yj - yk, 2.0);
if (rr < d)
result = false;
}
return result;
}
static double A33()
{
double a33 = 0.0;
for (int i = 0; i < numberSteps; i++)
{
double f12 = 0.0, f23 = 0.0, f31 = 0.0;
double xj = maxD * Random11();
double yj = maxD * Random11();
double xk = maxD * Random11();
double yk = maxD * Random11();
double rjk = sqrt(
pow(xj - xk, 2.0) +
pow(yj - yk, 2.0));
int j = rand() % 3 + 1;
int k = rand() % 3 + 1;
if (j == 1 && k == 2 && rjk < d)
a33 += 1.0;
else if (j == 2 && k == 3 && rjk < d)
a33 += 1.0;
else if (j == 3 && k == 1 && rjk < d)
a33 += 1.0;
else if (j == 2 && k == 1 && rjk < d)
a33 += 1.0;
else if (j == 3 && k == 2 && rjk < d)
a33 += 1.0;
else if (j == 1 && k == 3 && rjk < d)
a33 += 1.0;
}
a33 /= (pi * pi * pow(d, 4.0));
a33 /= numberSteps;
return a33;
}
static double A44()
{
double a44 = 0.0;
for (int i = 0; i < numberSteps; i++)
{
double f12 = 0.0, f23 = 0.0;
double f34 = 0.0, f41 = 0.0;
double xj = maxD * Random11();
double yj = maxD * Random11();
double xk = maxD * Random11();
double yk = maxD * Random11();
double rjk = sqrt(
pow(xj - xk, 2.0) +
pow(yj - yk, 2.0));
int j = rand() % 4 + 1;
int k = rand() % 4 + 1;
if (j == 1 && k == 2 && rjk < d)
a44 += 1.0;
else if (j == 2 && k == 3 && rjk < d)
a44 += 1.0;
else if (j == 3 && k == 4 && rjk < d)
a44 += 1.0;
else if (j == 4 && k == 1 && rjk < d)
a44 += 1.0;
}
a44 /= (pi * pi * pow(d, 4.0));
a44 /= numberSteps;
return a44;
}
static double A45()
{
double a45 = 0.0;
for (int i = 0; i < numberSteps; i++)
{
double f12 = 0.0, f23 = 0.0;
double f34 = 0.0, f41 = 0.0;
double f13 = 0.0;
double xj = maxD * Random11();
double yj = maxD * Random11();
double xk = maxD * Random11();
double yk = maxD * Random11();
double rjk = sqrt(
pow(xj - xk, 2.0) +
pow(yj - yk, 2.0));
int j = rand() % 4 + 1;
int k = rand() % 4 + 1;
if (j == 1 && k == 2 && rjk < d)
a45 += 1.0;
else if (j == 2 && k == 3 && rjk < d)
a45 += 1.0;
else if (j == 3 && k == 4 && rjk < d)
a45 += 1.0;
else if (j == 4 && k == 1 && rjk < d)
a45 += 1.0;
else if (j == 1 && k == 3 && rjk < d)
a45 += 1.0;
}
a45 /= (pi * pi * pow(d, 4.0));
a45 /= numberSteps;
return a45;
}
static double A46()
{
double a46 = 0.0;
for (int i = 0; i < numberSteps; i++)
{
double f12 = 0.0, f23 = 0.0;
double f34 = 0.0, f41 = 0.0;
double f13 = 0.0, f24 = 0.0;
double xj = maxD * Random11();
double yj = maxD * Random11();
double xk = maxD * Random11();
double yk = maxD * Random11();
double rjk = sqrt(
pow(xj - xk, 2.0) +
pow(yj - yk, 2.0));
int j = rand() % 4 + 1;
int k = rand() % 4 + 1;
if (j == 1 && k == 2 && rjk < d)
a46 += 1.0;
else if (j == 2 && k == 3 && rjk < d)
a46 += 1.0;
else if (j == 3 && k == 4 && rjk < d)
a46 += 1.0;
else if (j == 4 && k == 1 && rjk < d)
a46 += 1.0;
else if (j == 1 && k == 3 && rjk < d)
a46 += 1.0;
else if (j == 2 && k == 4 && rjk < d)
a46 += 1.0;
}
a46 /= (pi * pi * pow(d, 4.0));
a46 /= numberSteps;
return a46;
}
static double A55()
{
double a55 = 0.0;
for (int i = 0; i < numberSteps; i++)
{
double f12 = 0.0, f23 = 0.0;
double f34 = 0.0, f45 = 0.0;
double f51 = 0.0;
double xj = maxD * Random11();
double yj = maxD * Random11();
double xk = maxD * Random11();
double yk = maxD * Random11();
double rjk = sqrt(
pow(xj - xk, 2.0) +
pow(yj - yk, 2.0));
int j = rand() % 5 + 1;
int k = rand() % 5 + 1;
if (j == 1 && k == 2 && rjk < d)
a55 += 1.0;
else if (j == 2 && k == 3 && rjk < d)
a55 += 1.0;
else if (j == 3 && k == 4 && rjk < d)
a55 += 1.0;
else if (j == 4 && k == 5 && rjk < d)
a55 += 1.0;
else if (j == 5 && k == 1 && rjk < d)
a55 += 1.0;
}
a55 /= (pi * pi * pow(d, 4.0));
a55 /= numberSteps;
return a55;
}
static double A56p()
{
double a56p = 0.0;
for (int i = 0; i < numberSteps; i++)
{
double f12 = 0.0, f23 = 0.0;
double f34 = 0.0, f31 = 0.0;
double f41 = 0.0, f35 = 0.0;
double xj = maxD * Random11();
double yj = maxD * Random11();
double xk = maxD * Random11();
double yk = maxD * Random11();
double rjk = sqrt(
pow(xj - xk, 2.0) +
pow(yj - yk, 2.0));
int j = rand() % 5 + 1;
int k = rand() % 5 + 1;
if (j == 1 && k == 2 && rjk < d)
a56p += 1.0;
else if (j == 2 && k == 3 && rjk < d)
a56p += 1.0;
else if (j == 3 && k == 4 && rjk < d)
a56p += 1.0;
else if (j == 3 && k == 1 && rjk < d)
a56p += 1.0;
else if (j == 4 && k == 1 && rjk < d)
a56p += 1.0;
else if (j == 3 && k == 5 && rjk < d)
a56p += 1.0;
}
a56p /= (pi * pi * pow(d, 4.0));
a56p /= numberSteps;
return a56p;
}
static double A56pp()
{
double a56pp = 0.0;
for (int i = 0; i < numberSteps; i++)
{
double f12 = 0.0, f23 = 0.0;
double f34 = 0.0, f41 = 0.0;
double f35 = 0.0, f51 = 0.0;
double xj = maxD * Random11();
double yj = maxD * Random11();
double xk = maxD * Random11();
double yk = maxD * Random11();
double rjk = sqrt(
pow(xj - xk, 2.0) +
pow(yj - yk, 2.0));
int j = rand() % 5 + 1;
int k = rand() % 5 + 1;
if (j == 1 && k == 2 && rjk < d)
a56pp += 1.0;
else if (j == 2 && k == 3 && rjk < d)
a56pp += 1.0;
else if (j == 3 && k == 4)
a56pp += 1;
else if (j == 4 && k == 1 && rjk < d)
a56pp += 1.0;
else if (j == 3 && k == 5 && rjk < d)
a56pp += 1.0;
else if (j == 5 && k == 1 && rjk < d)
a56pp += 1.0;
}
a56pp /= (pi * pi * pow(d, 4.0));
a56pp /= numberSteps;
return a56pp;
}
static double A57p()
{
double a57p = 0.0;
for (int i = 0; i < numberSteps; i++)
{
double f12 = 0.0, f23 = 0.0;
double f34 = 0.0, f45 = 0.0;
double f51 = 0.0, f35 = 0.0;
double f41 = 0.0;
double xj = maxD * Random11();
double yj = maxD * Random11();
double xk = maxD * Random11();
double yk = maxD * Random11();
double rjk = sqrt(
pow(xj - xk, 2.0) +
pow(yj - yk, 2.0));
int j = rand() % 5 + 1;
int k = rand() % 5 + 1;
if (j == 1 && k == 2 && rjk < d)
a57p += 1.0;
else if (j == 2 && k == 3 && rjk < d)
a57p += 1.0;
else if (j == 3 && k == 4 && rjk < d)
a57p += 1.0;
else if (j == 4 && k == 5 && rjk < d)
a57p += 1.0;
else if (j == 5 && k == 1 && rjk < d)
a57p += 1.0;
else if (j == 3 && k == 5 && rjk < d)
a57p += 1.0;
else if (j == 4 && k == 1 && rjk < d)
a57p += 1.0;
}
a57p /= (pi * pi * pow(d, 4.0));
a57p /= numberSteps;
return a57p;
}
static double A57pp()
{
double a57pp = 0.0;
for (int i = 0; i < numberSteps; i++)
{
double f12 = 0.0, f23 = 0.0;
double f34 = 0.0, f45 = 0.0;
double f51 = 0.0, f31 = 0.0;
double f41 = 0.0;
double xj = maxD * Random11();
double yj = maxD * Random11();
double xk = maxD * Random11();
double yk = maxD * Random11();
double rjk = sqrt(
pow(xj - xk, 2.0) +
pow(yj - yk, 2.0));
int j = rand() % 5 + 1;
int k = rand() % 5 + 1;
if (j == 1 && k == 2 && rjk < d)
a57pp += 1.0;
else if (j == 2 && k == 3 && rjk < d)
a57pp += 1.0;
else if (j == 3 && k == 4 && rjk < d)
a57pp += 1.0;
else if (j == 4 && k == 5 && rjk < d)
a57pp += 1.0;
else if (j == 5 && k == 1 && rjk < d)
a57pp += 1.0;
else if (j == 3 && k == 1 && rjk < d)
a57pp += 1.0;
else if (j == 4 && k == 1 && rjk < d)
a57pp += 1.0;
}
a57pp /= (pi * pi * pow(d, 4.0));
a57pp /= numberSteps;
return a57pp;
}
static double A57ppp()
{
double a57ppp = 0.0;
for (int i = 0; i < numberSteps; i++)
{
double f12 = 0.0, f23 = 0.0;
double f34 = 0.0, f41 = 0.0;
double f35 = 0.0, f51 = 0.0;
double f31 = 0.0;
double xj = maxD * Random11();
double yj = maxD * Random11();
double xk = maxD * Random11();
double yk = maxD * Random11();
double rjk = sqrt(
pow(xj - xk, 2.0) +
pow(yj - yk, 2.0));
int j = rand() % 5 + 1;
int k = rand() % 5 + 1;
if (j == 1 && k == 2 && rjk < d)
a57ppp += 1.0;
else if (j == 2 && k == 3 && rjk < d)
a57ppp += 1.0;
else if (j == 3 && k == 4 && rjk < d)
a57ppp += 1.0;
else if (j == 4 && k == 1 && rjk < d)
a57ppp += 1.0;
else if (j == 3 && k == 5 && rjk < d)
a57ppp += 1.0;
else if (j == 5 && k == 1 && rjk < d)
a57ppp += 1.0;
else if (j == 3 && k == 1 && rjk < d)
a57ppp += 1.0;
}
a57ppp /= (pi * pi * pow(d, 4.0));
a57ppp /= numberSteps;
return a57ppp;
}
static double A58p()
{
double a58p = 0.0;
for (int i = 0; i < numberSteps; i++)
{
double f12 = 0.0, f23 = 0.0;
double f34 = 0.0, f45 = 0.0;
double f51 = 0.0, f31 = 0.0;
double f35 = 0.0, f41 = 0.0;
double xj = maxD * Random11();
double yj = maxD * Random11();
double xk = maxD * Random11();
double yk = maxD * Random11();
double rjk = sqrt(
pow(xj - xk, 2.0) +
pow(yj - yk, 2.0));
int j = rand() % 5 + 1;
int k = rand() % 5 + 1;
if (j == 1 && k == 2 && rjk < d)
a58p += 1.0;
else if (j == 2 && k == 3 && rjk < d)
a58p += 1.0;
else if (j == 3 && k == 4 && rjk < d)
a58p += 1.0;
else if (j == 4 && k == 5 && rjk < d)
a58p += 1.0;
else if (j == 5 && k == 1 && rjk < d)
a58p += 1.0;
else if (j == 3 && k == 1 && rjk < d)
a58p += 1.0;
else if (j == 3 && k == 5 && rjk < d)
a58p += 1.0;
else if (j == 4 && k == 1 && rjk < d)
a58p += 1.0;
}
a58p /= (pi * pi * pow(d, 4.0));
a58p /= numberSteps;
return a58p;
}
static double A58pp()
{
double a58pp = 0.0;
for (int i = 0; i < numberSteps; i++)
{
double f12 = 0.0, f23 = 0.0;
double f34 = 0.0, f45 = 0.0;
double f51 = 0.0, f31 = 0.0;
double f41 = 0.0, f25 = 0.0;
double xj = maxD * Random11();
double yj = maxD * Random11();
double xk = maxD * Random11();
double yk = maxD * Random11();
double rjk = sqrt(
pow(xj - xk, 2.0) +
pow(yj - yk, 2.0));
int j = rand() % 5 + 1;
int k = rand() % 5 + 1;
if (j == 1 && k == 2 && rjk < d)
a58pp += 1.0;
else if (j == 2 && k == 3 && rjk < d)
a58pp += 1.0;
else if (j == 3 && k == 4 && rjk < d)
a58pp += 1.0;
else if (j == 4 && k == 5 && rjk < d)
a58pp += 1.0;
else if (j == 5 && k == 1 && rjk < d)
a58pp += 1.0;
else if (j == 3 && k == 1 && rjk < d)
a58pp += 1.0;
else if (j == 4 && k == 1 && rjk < d)
a58pp += 1.0;
else if (j == 4 && k == 1 && rjk < d)
a58pp += 1.0;
else if (j == 2 && k == 5 && rjk < d)
a58pp += 1.0;
}
a58pp /= (pi * pi * pow(d, 4.0));
a58pp /= numberSteps;
return a58pp;
}
static double A59()
{
double a59 = 0.0;
for (int i = 0; i < numberSteps; i++)
{
double f12 = 0.0, f23 = 0.0;
double f34 = 0.0, f45 = 0.0;
double f51 = 0.0, f31 = 0.0;
double f41 = 0.0, f25 = 0.0;
double f35 = 0.0;
double xj = maxD * Random11();
double yj = maxD * Random11();
double xk = maxD * Random11();
double yk = maxD * Random11();
double rjk = sqrt(
pow(xj - xk, 2.0) +
pow(yj - yk, 2.0));
int j = rand() % 5 + 1;
int k = rand() % 5 + 1;
if (j == 1 && k == 2 && rjk < d)
a59 += 1.0;
else if (j == 2 && k == 3 && rjk < d)
a59 += 1.0;
else if (j == 3 && k == 4 && rjk < d)
a59 += 1.0;
else if (j == 4 && k == 5 && rjk < d)
a59 += 1.0;
else if (j == 5 && k == 1 && rjk < d)
a59 += 1.0;
else if (j == 3 && k == 1 && rjk < d)
a59 += 1.0;
else if (j == 3 && k == 5 && rjk < d)
a59 += 1.0;
else if (j == 4 && k == 1 && rjk < d)
a59 += 1.0;
else if (j == 2 && k == 5 && rjk < d)
a59 += 1.0;
else if (j == 3 && k == 5 && rjk < d)
a59 += 1.0;
}
a59 /= (pi * pi * pow(d, 4.0));
a59 /= numberSteps;
return a59;
}
static double A510()
{
double a510 = 0.0;
for (int i = 0; i < numberSteps; i++)
{
double f12 = 0.0, f23 = 0.0;
double f34 = 0.0, f45 = 0.0;
double f51 = 0.0, f31 = 0.0;
double f41 = 0.0, f25 = 0.0;
double f35 = 0.0, f24 = 0.0;
double xj = maxD * Random11();
double yj = maxD * Random11();
double xk = maxD * Random11();
double yk = maxD * Random11();
double rjk = sqrt(
pow(xj - xk, 2.0) +
pow(yj - yk, 2.0));
int j = rand() % 5 + 1;
int k = rand() % 5 + 1;
if (j == 1 && k == 2 && rjk < d)
a510 += 1.0;
if (j == 2 && k == 3 && rjk < d)
a510 += 1.0;
if (j == 3 && k == 4 && rjk < d)
a510 += 1.0;
if (j == 4 && k == 5 && rjk < d)
a510 += 1.0;
if (j == 5 && k == 1 && rjk < d)
a510 += 1.0;
if (j == 3 && k == 1 && rjk < d)
a510 += 1.0;
if (j == 3 && k == 5 && rjk < d)
a510 += 1.0;
if (j == 4 && k == 1 && rjk < d)
a510 += 1.0;
if (j == 2 && k == 5 && rjk < d)
a510 += 1.0;
if (j == 2 && k == 4 && rjk < d)
a510 += 1.0;
}
a510 /= (pi * pi * pow(d, 4.0));
a510 /= numberSteps;
return a510;
}
int main()
{
double K[8] = { 2.5, 2.45, 2.4, 2.35, 2.3, 2.25, 2.2, 2.15 };
double nu[8] = { 2.0, 4.0, 5.0, 5.5, 6.0, 6.25, 6.5, 7.0 };
unsigned int seed = 1;
srand(seed);
for (int i = 0; i < 8; i++)
{
double a33 = 0.0, a44 = 0.0, a45 = 0.0, a46 = 0.0, a55 = 0.0;
double a56p = 0.0, a56pp = 0.0, a57p = 0.0, a57pp = 0.0;
double a57ppp = 0.0, a58p = 0.0, a58pp = 0.0;
double N12 = 0.0, a59 = 0.0, a510 = 0.0;
InitializeData(K[i], nu[i]);
for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
a33 += A33();
for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
a44 += A44();
for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
a45 += A45();
for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
a46 += A46();
for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
a55 += A55();
for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
a56p += A56p();
for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
a56pp += A56pp();
for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
a57p += A57p();
for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
a57pp += A57pp();
for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
a57ppp += A57ppp();
for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
a58p += A58p();
for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
a58pp += A58pp();
for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
a59 += A59();
for (numberSteps = 200000; numberSteps <= 800000; numberSteps += 200000)
a510 += A510();
a33 /= 4.0;
a44 /= 4.0;
a45 /= 4.0;
a46 /= 4.0;
a55 /= 4.0;
a56p /= 4.0;
a56pp /= 4.0;
a57p /= 4.0;
a57pp /= 4.0;
a57ppp /= 4.0;
a58p /= 4.0;
a58pp /= 4.0;
a59 /= 4.0;
a510 /= 4.0;
std::cout << "A(3, 3) = " << a33 << std::endl;
std::cout << "A(4, 4) = " << a44 << std::endl;
std::cout << "A(4, 5) = " << a45 << std::endl;
std::cout << "A(4, 6) = " << a46 << std::endl;
std::cout << "A(5, 5) = " << a55 << std::endl;
std::cout << "A(5, 6') = " << a56p << std::endl;
std::cout << "A(5, 6'') = " << a56pp << std::endl;
std::cout << "A(5, 7') = " << a57p << std::endl;
std::cout << "A(5, 7'') = " << a57pp << std::endl;
std::cout << "A(5, 7''') = " << a57ppp << std::endl;
std::cout << "A(5, 8') = " << a58p << std::endl;
std::cout << "A(5, 8'') = " << a58pp << std::endl;
std::cout << "A(5, 9) = " << a59 << std::endl;
std::cout << "A(5, 10) = " << a510 << std::endl;
for (int i = 0; i < 6500; i++)
{
if (Move())
N12 += 1.0;
}
std::cout << "N(1/2) = ";
std::cout << N12 << std::endl;
double C0 = 64.0 * N12 / (N * N * (K[i] * K[i] - 1.0));
double A0 = 0.5 * sqrt(3.0) * d0 * d0 * N;
std::cout << "C0 = " << C0 << std::endl;
std::cout << "A0/A = " << A0 << std::endl;
std::cout << "A/A0 = " << 1.0 / A0 << std::endl;
double Ni = 1.0 / N12;
double C1 = pi / sqrt(3.0);
double C2 = (4.0 * pi * pi * A33() / 9.0) * Ni;
double C3a = (6.0 * A45() - 3.0 * A44() - A46()) / pow(3.0, 1.5);
double C3 = (pi * pi * pi * C3a) * Ni;
double C40 = 8.0 * pi * pi * pi / 135.0;
double C4a = 12.0 * a55 - 60.0 * a56p - 10.0 * a56pp;
double C4b = 30.0 * a57p + 60.0 * a57pp + 10.0 * a57ppp;
double C4c = -30.0 * a58p - 15.0 * a58pp;
double C4d = 10.0 * a59 - a510;
double C4 = (C40 * (C4a + C4b + C4c + C4d)) * Ni;
std::cout << "K = " << K[i] << std::endl;
std::cout << "nu = " << nu[i] << std::endl;
std::cout << "C(1) = " << C1 << std::endl;
std::cout << "C(2) = " << C2 << std::endl;
std::cout << "C(3) = " << C3 << std::endl;
std::cout << "C(4) = " << C4 << std::endl;
double X3 = C1 * A0 + C2 * A0 * A0 + C3 * A0 * A0 * A0 +
C4 * A0 * A0 * A0 * A0;
std::cout << "X(3) = " << X3 << std::endl;
}
return 0;
}
Blog Entry © Tuesday, October 20, 2026, by James Pate Williams, Jr., Three-Dimensional Heat Conduction Problem
// HeatConductionCube.cpp : This file contains the 'main' function.
// Program execution begins and ends there.
// (c) Monday, October 19, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// References:
// https://mathworld.wolfram.com/HeatConductionEquation.html
// https://www.wolframalpha.com/calculators/mathematics-calculus-integrals-definite-integral-calculator
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>
using namespace std;
double a = 2, b = 3, c = 4, kappa = 1.0;
int l = 5, m = 5, n = 5, N = 5;
vector<double> C, E, G, U;
static double U0(
double a, double b, double c,
double x, double y, double z,
vector<double> C, vector<double> E,
vector<double> G, int nn)
{
double pi = 4.0 * atan(1.0), sum = 0.0;
double fl = 0.0, fm = 0.0, fn = 0.0;
int l = 0, m = 0, n = 0;
for (l = 1; l <= nn; l++)
{
fl = C[l] * sin(l * pi * x / a);
for (m = 1; m <= nn; m++)
{
fm = E[m] * sin(m * pi * y / b);
for (n = 1; n <= nn; n++)
{
fn = G[n] * sin(n * pi * z / c);
sum += fl * fm * fn;
}
}
}
return sum;
}
static double DefiniteIntegral(double a, int p)
{
double pi = 4.0 * atan(1.0), pi2 = pi * pi;
double pi3 = pi * pi2, a3 = a * a * a;
double integral = 0.0;
double p2 = p * p, p3 = p * p2;
integral = (a3 * (2.0 - pi2 * p2) * cos(pi * p) +
2.0 * pi * p * sin(pi * p) - 2.0) / (pi3 * p3);
return integral;
}
static void ComputeVectors()
{
double x = 0.0, y = 0.0, z = 0.0;
double xs = a / N, ys = b / N, zs = c / N;
C.resize(N + 1LL);
E.resize(N + 1LL);
G.resize(N + 1LL);
for (int i = 1; i <= N; i++)
{
for (int j = 1; j <= N; j++)
{
for (int k = 1; k <= N; k++)
{
for (int p = 1; p <= N; p++)
C[p] = DefiniteIntegral(a, p);
for (int q = 1; q <= N; q++)
E[q] = DefiniteIntegral(b, q);
for (int r = 1; r <= N; r++)
G[r] = DefiniteIntegral(c, r);
z += zs;
}
y += ys;
}
x += xs;
}
}
static void ComputeU()
{
double pi = 4.0 * atan(1.0);
double ts = a / N, xs = a / N, ys = b / N, zs = c / N;
double t = ts, x = xs, y = ys, z = zs;
U.clear();
for (int h = 1; h <= N; h++)
{
for (int i = 1; i <= N; i++)
{
for (int j = 1; j <= N; j++)
{
for (int k = 1; k <= N; k++)
{
for (int l = 1; l <= N; l++)
{
for (int m = 1; m <= N; m++)
{
for (int n = 1; n <= N; n++)
{
double CEG = C[l] * E[m] * G[n];
double sl = sin(l * pi * x / a);
double sm = sin(m * pi * y / b);
double sn = sin(n * pi * z / c);
double el = exp(-kappa * pow(l * pi / a, 2.0) * t);
double em = exp(-kappa * pow(m * pi / b, 2.0) * t);
double en = exp(-kappa * pow(n * pi / c, 2.0) * t);
double uu = CEG * sl * sm * sn * el * em * en;
U.push_back(uu);
if (fabs(uu) < 1.0e-15)
return;
}
}
}
x += ys;
}
y += ys;
}
z += zs;
}
t += ts;
}
}
static void ReadData()
{
cout << "a = ";
cin >> a;
cout << "b = ";
cin >> b;
cout << "c = ";
cin >> c;
cout << "k = ";
cin >> kappa;
cout << "N = ";
cin >> N;
}
static void WriteResults()
{
ComputeVectors();
cout << "C\r\n";
for (int i = 1; i <= N; i++)
{
cout << fixed << right << setw(16);
cout << setprecision(13);
cout << C[i] << endl;
}
cout << "E\r\n";
for (int i = 1; i <= N; i++)
{
cout << fixed << right << setw(16);
cout << setprecision(13);
cout << E[i] << endl;
}
cout << "G\r\n";
for (int i = 1; i <= N; i++)
{
cout << fixed << right << setw(16);
cout << setprecision(13);
cout << G[i] << endl;
}
ComputeU();
cout << "U\r\n";
for (size_t i = 0; i < U.size(); i++)
{
cout << fixed << right << setw(16);
cout << setprecision(13);
cout << U[i] << endl;
}
}
int main()
{
ReadData();
WriteResults();
return 0;
}
Blog Entry @ Thursday, October 15, 2026, by James Pate Williams, Jr. Solution of the Telegraph Equation
// TelegraphEquation.cpp
// Copyright (c) Wednesday, October 14, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
#include <iomanip>
#include <iostream>
#include <vector>
static double omega(
double alpha, double beta, double c, double l, int n)
{
double pi = 4.0 * atan(1.0);
double term1 = pow(n * pi * c / l, 2.0);
double term2 = 0.25 * pow(alpha - beta, 2.0);
return sqrt(term1 - term2);
}
static double phi(
double alpha, double beta, double c, double l, int n)
{
double d = 0.5 * (alpha + beta);
return atan(d / omega(alpha, beta, c, l, n));
}
static double An(
double alpha, double beta, double a, double c, double l, int n)
{
double factor1 = 2.0 / (l, cos(phi(alpha, beta, c, l, n)));
double factor2 = sin(n * 4.0 * atan(1.0) * a / l);
return factor1 * factor2;
}
static double analytic_u(
double alpha, double beta, double a,
double c, double l, double x, double t,
int N)
{
double sum = 0.0;
for (int n = 1; n <= N; n++)
{
double d = 0.5 * (alpha + beta);
double an = An(alpha, beta, a, c, l, n);
double factor1 = an * exp(-d * t);
double omegaN = omega(alpha, beta, c, l, n);
double phiN = phi(alpha, beta, c, l, n);
double factor2 = cos(omegaN * t - phiN);
double factor3 = sin(n * 4.0 * atan(1.0) * x / l);
sum += factor1 * factor2 * factor3;
}
return sum;
}
int main()
{
char line[128] = "";
double alpha = 0.0, beta = 0.0, a = 0.0;
double c = 0.0, l = 0.0;
int N = 0, nPoints = 0;
std::cout << "alpha = ";
std::cin.getline(line, 128);
alpha = atof(line);
std::cout << "beta = ";
std::cin.getline(line, 128);
beta = atof(line);
std::cout << "a = ";
std::cin.getline(line, 128);
a = atof(line);
std::cout << "c = ";
std::cin.getline(line, 128);
c = atof(line);
std::cout << "l = ";
std::cin.getline(line, 128);
l = atof(line);
std::cout << "N = ";
std::cin.getline(line, 128);
N = atoi(line);
std::cout << "nPoints = ";
std::cin.getline(line, 128);
nPoints = atoi(line);
double deltaX = 1.0 / nPoints;
double deltaT = 1.0 / nPoints;
std::vector<double> u(nPoints), x(nPoints), t(nPoints);
for (int i = 0; i < nPoints; i++)
{
x[i] = i * deltaX;
t[i] = i * deltaT;
u[i] = analytic_u(alpha, beta, a, c, l, x[i], t[i], N);
std::cout << std::setprecision(10);
std::cout << std::fixed;
std::cout << x[i] << '\t' << t[i] << '\t';
if (u[i] >= 0.0)
std::cout << ' ';
std::cout << u[i] << std::endl;
}
return 0;
}
Blog Entry © 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;
}
