#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;
}
Numerical Explorations 