Category: Laplace Equation

Blog Entry © Sunday, October 11, 2026, by James Pate Williams, Jr.

Blog Entry Sunday October 11 2026Download
#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;
}

Blog Entry © Saturday July 4, 2026, by James Pate Williams, Jr. Permutations Etc.

Blog Entry Saturday July 4 2026Download

Blog Entry © Thursday, May 28, 2026, by James Pate Williams, Jr. and Microsoft’s Copilot Solution of the Potential Equation in Rectangle using Fixed Point Iteration in Python

# NOTE:
# This implementation prioritizes clarity and correctness over optimization.
# Further performance improvements can be made if needed.
# (c) May 26, 2026 by James Pate Willims, Jr.
# I had some help from the Microsoft Copilot
# to calculate runtimes and define matrices
# Computes the potential in a rectangle
# Reference: "Boundary Value Problems
# Second Edition" by David L. Powers
# See pages 179 to 182 for the analytic
# solution of this Laplace Equation
# Stand alone application using
# Microsoft Visual Studio 2022
# Community Version

import math
import time

xi = yi = 10
u = [[0.0 for _ in range(xi + 2)] for _ in range(yi + 2)]
v = [[0.0 for _ in range(xi + 2)] for _ in range(yi + 2)]

def ComputeBoundaryValues(x, y):
    if x == 0:
        return 0
    if x == 1:
        return 0
    if y == 0 or y == 1:
        if x > 0.0 and x < 0.5:
            return 2.0 * x
        elif x >= 0.5 and x < 1.0:
            return 2.0 - 2.0 * x
                 
    return 0.0

def ComputeParams(its, norm, params):
    params['iterations'] = its
    params['norm'] = norm

def Compute(h, k, xi, yi, maxIts, params):
    # Use a simple fixed-point iteration to
    # compute an approximate solution
    for i in range(0, xi + 1):
        for j in range(0, yi + 1):
            u[i][j] = ComputeBoundaryValues(i * h, j * k)

    for its in range(1, maxIts + 1):
        for i in range(1, xi):
            for j in range(1, yi): 
                u[i][j] = 0.25 * (u[i + 1][j] + u[i - 1][j] + u[i][j + 1] + u[i][j - 1]);

    norm = 0

    for i in range(0, xi + 1):
        for j in range(0, yi + 1):
            norm += math.fabs(u[i][j] * u[i][j])

    norm = math.sqrt(norm)
    params['iterations'] = its
    params['norm'] = norm

def f(x, y):
    # Analytic solution series expansion n = 1 to 100 
    sum = 0.0

    for n in range(1, 101):
        factor1 = math.sin(n * math.pi / 2.0) / (n * n)
        factor2 = math.sinh(n * math.pi * y)
        factor3 = math.sinh(n * math.pi * (1 - y))
        factor4 = math.sin(n * math.pi * x)
        term = (factor2 + factor3) / math.sinh(n * math.pi)
        sum += factor1 * term * factor4
    return 8.0 * sum / (math.pi * math.pi)

avgPE = 0
deltaX = 1.0 / xi
deltaY = 1.0 / yi
maxIts = 50
start_time = time.perf_counter()

for i in range(0, xi + 1):
    for j in range(0, yi + 1):
        v[i][j] = f(i * deltaX, j * deltaY)

minPE = +1000000000
maxPE = -1000000000
params = {}

Compute(deltaX, deltaY, xi, yi, maxIts, params)
print("Approximate\tAnalytic\tPercent Error")

for i in range(0, xi + 1):
    for j in range(0, yi + 1):
        if (math.fabs(u[i][j]) > 1.0e-12 and
            math.fabs(v[i][j]) > 1.0e-12):
            pe = 100.0 * math.fabs((v[i][j] - u[i][j]) / v[i][j])
        else:
            pe = 0.0

        avgPE += pe

        if (pe < minPE):
            minPE = pe

        if (pe > maxPE):
            maxPE = pe

        if math.fabs(pe) != 0.0:
            print("{:10.8f}".format(u[i][j]), "\t", "{:10.8f}".format(v[i][j]), "\t", "{:10.8f}".format(pe))

avgPE /= (xi * yi)
end_time = time.perf_counter()
# Calculate elapsed time in milliseconds
elapsed_ms = (end_time - start_time) * 1000

print("Iterations = ", params['iterations'])
print("Frobenius Norm = ", params['norm'])
print("Minimum Percent Error = ", "{:10.8f}".format(minPE))
print("Average Percent Error = ", "{:10.8f}".format(avgPE))
print("Maximum Percent Error = ", "{:10.8f}".format(maxPE))
print("Elapsed Milliseconds  = ", "{:10.8f}".format(elapsed_ms))


Approximate     Analytic        Percent Error
0.20000000       0.19999972      0.00013831
0.16633455       0.16663592      0.18085704
0.13739427       0.13768928      0.21425509
0.11591159       0.11605132      0.12040154
0.10292732       0.10291871      0.00836375
0.09864668       0.09854114      0.10710198
0.10305612       0.10291871      0.13351975
0.11613103       0.11605132      0.06868445
0.13763395       0.13768928      0.04018144
0.16650309       0.16663592      0.07971563
0.20000000       0.19999972      0.00013831
0.40000000       0.39999927      0.00018169
0.32813882       0.32854798      0.12453472
0.26763492       0.26776105      0.04710489
0.22369849       0.22344522      0.11334607
0.19755216       0.19705465      0.25246951
0.18899096       0.18834090      0.34515051
0.19778524       0.19705465      0.37075402
0.22409557       0.22344522      0.29105301
0.26806863       0.26776105      0.11486989
0.32844379       0.32854798      0.03171204
0.40000000       0.39999927      0.00018169
0.60000000       0.59999811      0.00031532
0.47888974       0.47875999      0.02710181
0.38176991       0.38059768      0.30799710
0.31425594       0.31267225      0.50650231
0.27518847       0.27354681      0.60013876
0.26255192       0.26082096      0.66365837
0.27549367       0.27354681      0.71170992
0.31477587       0.31267225      0.67278721
0.38233779       0.38059768      0.45720485
0.47928905       0.47875999      0.11050679
0.60000000       0.59999811      0.00031532
0.80000000       0.79999202      0.00099701
0.60602379       0.60222488      0.63081180
0.46685956       0.46199176      1.05365615
0.37704317       0.37308607      1.06063957
0.32711029       0.32392135      0.98448025
0.31121899       0.30818168      0.98555933
0.32745161       0.32392135      1.08985035
0.37762462       0.37308607      1.21648870
0.46749464       0.46199176      1.19112076
0.60647034       0.60222488      0.70496220
0.80000000       0.79999202      0.00099701
1.00000000       0.99594729      0.40692036
0.67874673       0.65811281      3.13531720
0.50319768       0.49282441      2.10486107
0.40066316       0.39470092      1.51057096
0.34574699       0.34157202      1.22228048
0.32848272       0.32468552      1.16950228
0.34608839       0.34157202      1.32223096
0.40124475       0.39470092      1.65792212
0.50383291       0.49282441      2.23375660
0.67919339       0.65811281      3.20318626
1.00000000       0.99594729      0.40692036
0.80000000       0.79999202      0.00099701
0.60615279       0.60222488      0.65223206
0.46709299       0.46199176      1.10418350
0.37734885       0.37308607      1.14257161
0.32745218       0.32392135      1.09002599
0.31156100       0.30818168      1.09653477
0.32776105       0.32392135      1.18538026
0.37787503       0.37308607      1.28360546
0.46766769       0.46199176      1.22857864
0.60655688       0.60222488      0.71933113
0.80000000       0.79999202      0.00099701
0.60000000       0.59999811      0.00031532
0.47910949       0.47875999      0.07300025
0.38216756       0.38059768      0.41247566
0.31477665       0.31267225      0.67303766
0.27577086       0.27354681      0.81304189
0.26313452       0.26082096      0.88702848
0.27602079       0.27354681      0.90440981
0.31520243       0.31267225      0.80920943
0.38263258       0.38059768      0.53465917
0.47943646       0.47875999      0.14129609
0.60000000       0.59999811      0.00031532
0.40000000       0.39999927      0.00018169
0.32837881       0.32854798      0.05149022
0.26806920       0.26776105      0.11508237
0.22426717       0.22344522      0.36785094
0.19818820       0.19705465      0.57524397
0.18962723       0.18834090      0.68297865
0.19836093       0.19705465      0.66290001
0.22456142       0.22344522      0.49953914
0.26839057       0.26776105      0.23510691
0.32860478       0.32854798      0.01728752
0.40000000       0.39999927      0.00018169
0.20000000       0.19999972      0.00013831
0.16650327       0.16663592      0.07960469
0.13769959       0.13768928      0.00748883
0.11631140       0.11605132      0.22411146
0.10337449       0.10291871      0.44285504
0.09909401       0.09854114      0.56105764
0.10346087       0.10291871      0.52678322
0.11645855       0.11605132      0.35090580
0.13786030       0.13768928      0.12420922
0.16661627       0.16663592      0.01179315
0.20000000       0.19999972      0.00013831
Iterations =  50
Frobenius Norm =  4.028216200275417
Minimum Percent Error =  0.00000000
Average Percent Error =  0.54286140
Maximum Percent Error =  3.20318626
Elapsed Milliseconds  =  36.23520000
Press any key to continue . . .

Blog Entry © Monday, May 11, 2026, by James Pate Williams, Jr., Laplace Equation in a Solid Cylinder

Blog Entry Monday May 11 2026Download
Laplace Equation in a Solid CylinderDownload

Blog Entry © Tuesday, May 5, 2026, by James Pate Williams, Jr., Solution of the Laplace Equation on a Rectangle

Blog Entry Tuesday May 5 2026Download

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