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