Numerical Solution of the
Laplace Equation
Finite Difference Method and Liebmann's Iterative Technique · Course: Numerical Methods
The Laplace Equation
A brief look at what the equation means, where it is used, and why it needs a numerical solution.
u(x, y) represents a steady-state field, such as temperature, that no longer changes with time. The equation states that its value at every interior point equals the average of the field around it.
Application
A common example is steady-state heat conduction in a thin metal plate: given fixed temperatures on all four edges, the equation gives the temperature at every interior point once the plate reaches equilibrium.
Why a numerical method is needed
Exact analytical solutions exist only for simple shapes and boundary conditions. Real plates and boundary conditions are rarely that simple, so the equation is solved using the Finite Difference Method together with Liebmann's iteration.
Finite Difference Method & Liebmann's Iteration
The plate is discretized into a grid, and the equation is replaced by a simple averaging rule applied repeatedly.
Finite difference discretization
Replacing both second derivatives with central differences on a grid of spacing h and substituting into ∇²u = 0 gives the five-point stencil below.
Every interior node equals the average of its four neighbouring nodes: top, bottom, left, and right.
Liebmann's (Gauss–Seidel) iteration
Interior nodes start at zero. The stencil formula is applied to every node in turn, always using the latest available neighbour values, and this sweep is repeated across the whole grid.
Convergence and tolerance
Each sweep records the largest change seen at any node. Once that maximum change falls below a chosen tolerance ε, the field has converged and the iteration stops.
Solution Procedure
Step-by-step algorithm executed sequentially by the solver program.
Input Collection: Read the grid size N, convergence tolerance ε, and the four boundary temperatures (top, bottom, left, and right).
Grid Initialization: Create an N × N matrix initialized with 0 for all interior nodes, and assign fixed values to boundary rows and columns.
Iterative Sweep (Liebmann's Method): Iterate sequentially through each interior node (row i from 1 to N-2, column j from 1 to N-2).
Five-Point Stencil Update: Replace each interior node value with the average of its four immediate neighbors: ui,j = ¼ (ui+1,j + ui-1,j + ui,j+1 + ui,j-1).
Convergence Check: Calculate the absolute difference between the new and old node values, tracking maximum change across the sweep. Repeat sweeps until max change < ε.
Output Results: Output the converged temperature distribution matrix, total iterations executed, and final maximum error.
Python Implementation
The core routine: build the grid, then sweep with the five-point stencil until convergence.
import numpy as np def read_inputs(): n = int(input("Grid dimension N (for an N x N mesh): ")) tol = float(input("Error tolerance, e.g. 0.001: ")) top = float(input("Top boundary value (deg C): ")) bottom = float(input("Bottom boundary value (deg C): ")) left = float(input("Left boundary value (deg C): ")) right = float(input("Right boundary value (deg C): ")) return n, tol, top, bottom, left, right def initialize_grid(n, top, bottom, left, right): u = np.zeros((n, n)) u[0, :] = top u[n - 1, :] = bottom u[:, 0] = left u[:, -1] = right return u def liebmann_sweep(u): n = u.shape[0] max_diff = 0.0 for i in range(1, n - 1): for j in range(1, n - 1): old_val = u[i, j] u[i, j] = 0.25 * (u[i + 1, j] + u[i - 1, j] + u[i, j + 1] + u[i, j - 1]) max_diff = max(max_diff, abs(u[i, j] - old_val)) return max_diff def solve_laplace(n, tol, top, bottom, left, right): u = initialize_grid(n, top, bottom, left, right) iteration = 0 while True: max_diff = liebmann_sweep(u) iteration += 1 if max_diff < tol: break return u, iteration, max_diff if __name__ == "__main__": n, tol, top, bottom, left, right = read_inputs() u, iteration, max_diff = solve_laplace(n, tol, top, bottom, left, right) print(f"Grid size : {n} x {n}") print(f"Tolerance : {tol}") print(f"Boundary values : top={top}, bottom={bottom}, left={left}, right={right}") print(f"Converged in : {iteration} iterations") print(f"Final max change : {max_diff:.6f}") print(np.round(u, 2))
Live Solver
Set the boundary conditions and watch the field relax to a steady state. Drag freely to orbit the 3D surface from any angle, including underneath.
Grid & tolerance
Boundary conditions (°C)
Conclusion
Summary of what this study has shown.
The Laplace equation was successfully solved numerically for a rectangular plate with prescribed boundary temperatures.
The Finite Difference Method converts the PDE into algebraic equations by replacing derivatives with central differences on a discrete grid.
Liebmann's iteration efficiently computes the solution by repeatedly averaging each node with its neighbours until the field converges.
The simulation validates the numerical method, visually confirming that the field settles into a smooth, physically consistent steady state.