Category: Laplace Equation

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

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

Four Methods of Numerical Double Integration: Sequential Simpson’s Rule, Multitasking Simpson’s Rule, Sequential Monte Carlo and Multitasking Monte Carlo Methods © Wednesday April 16 – 18, 2025, by James Pate Williams, Jr.

Four Methods of Numerical Double IntegrationDownload

Blog Entry © Tuesday, March 25, 2025, by James Pate Williams, Jr. Hydrogen Radial Wavefunctions and Related Functions

Hydrogen Radial Wavefunctions and Related FunctionsDownload

Some Helium Coulomb Integrals over Six Dimensions by James Pate Williams, Jr. Source Code in C++ Development over December 15 – 16, 2023

some-helium-coulomb-integrals-1Download

Revised Translated Source Code from May 15, 2015, by James Pate Williams, Jr.

revisited-with-source-code-translation-and-revisionDownload

New and Corrected Ground State Energy Numerical Computation for the Helium Like Atom (Atomic Number 2) by James Pate Williams, Jr.

approximate-solution-of-the-schrodinger-equation-for-the-helium-atom-1Download
new-variational-treatment-of-the-helium-like-atom-revisedDownload

A New Calculus of Variations Solution of the Schrödinger Equation for the Lithium Like Atom’s Ground State Energy

This computation took a lot longer time to reach a much better solution than my previously published result.

new-variational-solution-for-the-lithiumDownload

A Calculus of Variations Solution to the Quantum Mechanical Schrödinger Wave Equation for the Lithium Like Atom (Atomic Number Z = 3) by James Pate Williams, Jr.

variational-solution-for-the-lithiumDownload