Category: Laplace Equation
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 © Tuesday, May 5, 2026, by James Pate Williams, Jr., Solution of the Laplace Equation on a Rectangle
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;
}
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.
Blog Entry © Tuesday, March 25, 2025, by James Pate Williams, Jr. Hydrogen Radial Wavefunctions and Related Functions
Classification of Two Mathematical Entities by James Pate Williams, Jr. Sunday, May 12, 2024
Happy Mother’s Day
Numerical Explorations 