Author: jamespatewilliamsjr

// 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;
}

// 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;
}