Poisson equation finite-difference convergence analysis

What this research found

Verification asks a narrow but essential question of any numerical solver: does the implementation actually achieve the accuracy the theory promises? K-Dense built a five-point finite-difference solver for the two-dimensional Poisson equation on the unit square and tested it with the method of manufactured solutions, choosing a smooth target field whose exact source term can be derived analytically. Across five grids up to 65,536 unknowns, the measured convergence rates climbed to 2.0000 in the averaged norm and 1.9999 in the worst-case norm, confirming clean second-order accuracy with no order-degrading defects.

• Errors fell by very nearly a factor of four with each halving of the grid spacing, the signature of a second-order method. The averaged error dropped from 1.4254 x 10^-3 on the coarsest grid to 6.2262 x 10^-6 on the finest, and the worst-case pointwise error from 2.8265 x 10^-3 to 1.2452 x 10^-5.
• Observed convergence rates approached the theoretical value monotonically and from a consistent direction: the averaged-norm rate went 2.0019, 2.0005, 2.0001, 2.0000, and the maximum-norm rate 1.9924, 1.9980, 1.9995, 1.9999. A first-order boundary treatment, a mis-scaled grid factor, or an off-by-one indexing error at the domain edge would have stalled these rates near one instead.
• The ratio between the two error norms stayed essentially constant at roughly 1.98 to 2.00 across all five grids, showing the error keeps a single self-similar shape under refinement rather than developing a boundary layer.
• The spatial error is a smooth, nearly radially symmetric mound peaking at 1.245 x 10^-5 at the centre of the domain and vanishing on the boundary, with no oscillations or grid-aligned artefacts. This matches the truncation analysis, which for this manufactured field makes the leading error proportional to the solution itself and therefore largest at the central antinode.
• Magnitudes agree as well as shapes. The predicted leading truncation error at the centre on the finest grid is about 2.460 x 10^-5, while the measured maximum global error of 1.245 x 10^-5 is the same order and slightly smaller, reflecting the smoothing action of the discrete Green's function that maps local truncation error into global solution error.

How it was done

The manufactured solution was taken as the product of two sine half-waves, which is infinitely differentiable, vanishes exactly on all four edges, and yields a matching source term in closed form. The discrete negative Laplacian was assembled not entry by entry but as a Kronecker sum of the one-dimensional second-difference operator with the identity, which reproduces the five-point coupling exactly and guarantees a symmetric positive-definite, block-tridiagonal matrix with at most five nonzeros per row. Each system was solved by direct sparse LU factorization rather than an iterative method, so the algebraic error sits at machine round-off — many orders of magnitude below the discretization error being measured — and no iteration tolerance can bias the observed rate. Errors were reported in both a grid-normalized averaged norm and a maximum norm on grids of 16, 32, 64, 128, and 256 interior nodes per dimension, with the order of convergence estimated as the log-log slope between consecutive grids.

Data sources

• Roache, Journal of Fluids Engineering 124:4 (2002) — code verification by the method of manufactured solutions
• LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations, SIAM (2007)
• Oberkampf & Roy, Verification and Validation in Scientific Computing, Cambridge (2010)
• Davis, ACM Transactions on Mathematical Software 30:196 (2004) — UMFPACK multifrontal sparse LU