Blog Entry © Thursday, March 19, 2026, by James Pate Williams, Jr. Solution of Laplace’s Two-Dimensional Equation on a Square Boundary

#pragma once

class Laplace
{
private:
    static int N;

public:
    Laplace(int N);
    double** Solve(
        double L,
        double epsilon,
        int& iterations,
        double (*f)(double));
    void PrintSolution(
        double L,
        double epsilon,
        double** u,
        double** w,
        int& iterations);
};

#include "Laplace.h"
#include <float.h>
#include <math.h>
#include <iomanip>
#include <iostream>

int Laplace::N = 0;

static double fxy(double x, double y)
{
    return x * y * exp(-0.5 * (x * x + y * y));
}

static double gxy(double x, double y)
{
    return x * sin(y);
}

Laplace::Laplace(int N)
{
    this->N = N;
}

double** Laplace::Solve(
    double L, double epsilon, int& iterations,
    double (*f)(double))
{
    double h = L / (N + 1.0), uv = 0, max = 0;
    double** u = new double*[N + 2];

    if (u == nullptr)
        exit(-1);

    for (int i = 0; i < N + 2; i++)
    {
        u[i] = new double[N + 2];

        if (u[i] == nullptr)
            exit(-2);

        for (int j = 0; j < N + 2; j++)
            u[i][j] = 0.0;
    }

    double** v = new double* [N + 2];

    if (v == nullptr)
        exit(-3);

    for (int i = 0; i < N + 2; i++)
    {
        v[i] = new double[N + 2];

        if (v[i] == nullptr)
            exit(-4);

        for (int j = 0; j < N + 2; j++)
            v[i][j] = 0.0;
    }

    iterations = 0;

    for (int j = 0; j <= N + 1; j++)
        u[0][j] = 0.0;

    for (int j = 0; j <= N + 1; j++)
        u[N + 1][j] = 0.0;

    for (int i = 0; i <= N + 1; i++)
        u[i][0] = f(i * h);

    for (int i = 0; i <= N + 1; i++)
        u[i][N + 1] = 0.0;

    do
    {
        iterations++;

        for (int i = 1; i <= N; i++)
            for (int j = 1; j <= N; j++)
                u[i][j] = 0.25 * (u[i - 1][j] + u[i + 1][j] + u[i][j - 1] +
                    u[i][j + 1]);

        for (int i = 0; i <= N + 1; i++)
            for (int j = 0; j <= N + 1; j++)
                v[i][j] = u[i][j];

        for (int i = 1; i <= N; i++)
            for (int j = 1; j <= N; j++)
                v[i][j] = 0.25 * (v[i - 1][j] + v[i + 1][j] + v[i][j - 1] +
                    v[i][j + 1]);

        max = DBL_MIN;

        for (int i = 1; i <= N; i++)
        {
            for (int j = 1; j <= N; j++)
            {
                uv = fabs(v[i][j] - u[i][j]);
                max = fmax(uv, max);
            }
        }

        for (int i = 0; i <= N + 1; i++)
            for (int j = 0; j <= N + 1; j++)
                u[i][j] = v[i][j];

    } while (max > epsilon);

    for (int i = 0; i < N + 2; i++)
        delete v[i];

    return u;
}

void Laplace::PrintSolution(
    double L,
    double epsilon,
    double** u,
    double** w,
    int& iterations)
{
    double h = L / (N + 1.0);
    
    std::cout << std::endl;
    std::cout << "x\t\ty\t\tu(x, y)\t\tw(x, y)\t\tpercent error\r\n\r\n";

    for (int i = 0; i <= N + 1; i++)
    {
        double x = i * h;

        for (int j = 0; j <= N + 1; j++)
        {
            double y = j * h;
            std::cout << std::fixed << std::setprecision(5);
            std::cout << std::setw(8);
            std::cout << x << '\t';
            std::cout << std::fixed << std::setprecision(5);
            std::cout << std::setw(8);
            std::cout << y << '\t';
            std::cout << std::fixed << std::setprecision(5);
            std::cout << std::setw(8);
            std::cout << u[i][j] << '\t';
            std::cout << std::fixed << std::setprecision(5);
            std::cout << std::setw(8);
            std::cout << w[i][j] << '\t';
            double error = 0.0;
            if (w[i][j] != 0.0)
                error = 100.0 * fabs(u[i][j] - w[i][j]) /
                    fabs(w[i][j]);
            std::cout << std::fixed << std::setprecision(5);
            std::cout << std::setw(8);
            std::cout << error << std::endl;
        }
    }

    std::cout << std::endl;
    std::cout << "iterations = " << iterations;
    std::cout << std::endl;
}

// LaplaceEquationFD2d.cpp (c) Tuesday, March 17, 2026 
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// by James Pate Williams, Jr., BA, ,BS, MSwE, PhD
// Reference: Numerical Computation 2 Methods,
// Software And Analysis (c) 1995 by Chrisoph E. Ueberhuber
// Pages 392 - 393 Poisson Matrices
// https://math.libretexts.org/Bookshelves/Differential_Equations/Introduction_to_Partial_Differential_Equations_(Herman)/06%3A_Problems_in_Higher_Dimensions/6.03%3A_Laplaces_Equation_in_2D

#include <math.h>
#include <iostream>
#include "Laplace.h"

double** w;

double SimpsonsRule(double L, double a, double b,
    int Nsr, int n, double (*f)(double, double, int))
{
    double h = (b - a) / Nsr;
    double h2 = 2.0 * h;
    double s = 0.0;
    double t = 0.0;
    double x = a + h;

    for (int i = 1; i < Nsr; i += 2)
    {
        s += f(L, x, n);
        x += h2;
    }

    x = a + h2;

    for (int i = 2; i < Nsr; i += 2)
    {
        t += f(L, x, n);
        x += h2;
    }

    return
        h * (f(L, a, n) + 4 * s + 2 * t + f(L, b, n)) / 3.0;
}

double f(double x)
{
    return x * exp(-x) * cos(x);
}

double g(double L, double x, int n)
{
    double pi = 4.0 * atan(1.0);
    double fx = f(x);
    return fx * sin(n * pi * x / L);
}

double bn(double L, int Nsr, int n)
{
    double cs = 2.0 / L;
    double sr = SimpsonsRule(L, 0.0, L, Nsr, n, g);
    return cs * sr;
}

void ComputeW(double L, int N)
{
    double h = L / (N + 1.0);

    w = new double* [N + 2];

    if (w == nullptr)
        exit(-5);

    for (int i = 0; i < N + 2; i++)
    {
        w[i] = new double[N + 2];

        if (w[i] == nullptr)
            exit(-6);

        for (int j = 0; j < N + 2; j++)
            w[i][j] = 0.0;
    }

    for (int j = 0; j <= N; j++)
        w[0][j] = 0.0;

    for (int j = 0; j <= N; j++)
        w[N + 1][j] = 0.0;

    for (int i = 0; i <= N; i++)
        w[i][0] = f(i * h);

    for (int i = 0; i <= N; i++)
        w[i][N + 1] = 0.0;

    double pi = 4.0 * atan(1.0);
    double x = h, y = 0.0;

    for (int i = 1; i <= N; i++)
    {
        y = h;

        for (int j = 1; j <= N; j++)
        {
            double s = 0.0;

            for (int n = 1; n <= 64; n++)
            {
                double npi = n * pi;
                double an = bn(L, 512, n) / sinh(npi);
                s += an * sin(npi * x / L) * sinh(npi * (L - y) / L);
            }

            w[i][j] = s;
            y += h;
        }

        x += h;
    }
}

int main()
{
    while (true)
    {
        char line[128] = "";
        std::cout << "N = ";
        std::cin.getline(line, 128);
        int N = atoi(line);

        if (N < 4 || N > 128)
        {
            std::cout << "N < 4 || N > 128, try again";
            std::cout << std::endl;
        }

        std::cout << "epsilon = ";
        std::cin.getline(line, 128);

        double epsilon = atof(line);
        Laplace l(N);

        double L = 1.0;
        int iterations = 0;
        ComputeW(L, N);
        double** u = l.Solve(L, epsilon, iterations, f);
        l.PrintSolution(L, epsilon, u, w, iterations);

        for (int i = 0; i < N + 2; i++)
        {
            if (u[i] != nullptr)
                delete u[i];

            if (w[i] != nullptr)
                delete w[i];
        }

        break;
    }

    return 0;
}