Laplace Boundary Problems

Charge-free regions and the boundary-value problem

Laplace’s equation determines electrostatic potential inside a region that contains no volume charge. Begin with Gauss’s law and the potential-gradient relation:

\nabla \cdot E = \frac{ρ}{ε _{0}}, E = - \nabla V .

Substitution gives Poisson’s equation,

\nabla^{2} V = - \frac{ρ}{ε _{0}} .

In a homogeneous charge-free region, $ρ = 0$ , so

\nabla^{2} V = 0.

The charge-free condition applies only inside the chosen domain. Charges can lie outside it. Conductors, electrodes, and remote sources affect the solution through boundary conditions. A narrow vacuum gap between charged metal objects therefore often satisfies Laplace’s equation even though metal surface charge created the surrounding electric field.

A charge-free region enclosed by four electrodes. Surface charge sits on the metal boundaries; every interior vacuum point satisfies Laplace's equation

\nabla^{2} V = 0

, and the four boundary potentials fix the interior solution.

A boundary-value problem requires three pieces of information.

Domain: the spatial region $Ω$ where the differential equation is required.
Boundary data: potential values, normal derivatives, total conductor charge, or a stated combination on the enclosing boundary $\partial Ω$ .
Admissible solution: a potential satisfying the equation in $Ω$ and every stated condition on $\partial Ω$ .

The boundary is part of the physics. Replacing a grounded electrode with an isolated electrode changes the problem even when the metal has the same shape. Grounding specifies a potential. Isolation commonly specifies total charge or an electrical connection constraint, leaving the conductor potential to be determined with the rest of the solution.

Laplace’s equation is local and second order. In Cartesian coordinates it is

\frac{\partial ^{2} V}{\partial x ^{2}} + \frac{\partial ^{2} V}{\partial y ^{2}} + \frac{\partial ^{2} V}{\partial z ^{2}} = 0.

Each second derivative measures local curvature along one coordinate. Their sum must vanish. Positive curvature along one direction must be balanced by negative curvature along another, except when every term is zero. A linear potential satisfies Laplace’s equation because every second derivative vanishes. A potential that curves upward in both $x$ and $y$ fails the charge-free condition in a homogeneous medium.

Harmonic potentials have a mean-value property. In two dimensions, the value at a point equals the average around any sufficiently small circle centered at that point:

V (r_{0}) = \frac{1}{2 π} 0 \int 2 π V (r_{0} + a \overset{e}{^}_{θ}) d θ .

The corresponding three-dimensional statement uses the average over a sphere. A nonconstant potential has extrema on the boundary of the connected region. The maximum principle is both a physical statement about charge-free space and a numerical diagnostic: a relaxation calculation that creates a new interior maximum from bounded Dirichlet data contains an error.

Potential values require locations and boundary interpretation before they identify a physical solution. A map sampled on a square may represent a conducting cavity, a cross-section through a long apparatus, or a mathematical model with translational symmetry. The governing equation and boundary data must match that interpretation. A two-dimensional model assumes that omitted variation is absent or has been accounted for by an effective boundary condition.

Boundary data, conductors, and uniqueness

Dirichlet data specify potential on a boundary.

V_{\partial Ω} = g .

Grounded metal has $g = 0$ . An electrode held at a known voltage has constant $g = V_{c}$ on its connected surface. Dirichlet data are natural for apparatuses driven by calibrated supplies because the voltage is measured directly at the electrode connection.

Neumann data specify the outward normal derivative.

\frac{\partial V}{\partial n}_{\partial Ω} = h .

Because $E = - \nabla V$ , the normal electric-field component is $E_{n} = - \partial V / \partial n$ . A given normal derivative is therefore a given normal flux density after multiplication by $ε_{0}$ . Neumann data appear when surface charge is known, when symmetry forces zero normal derivative, or when a boundary represents an insulating symmetry plane in a reduced model.

A charge-free domain with pure Neumann data requires a compatibility condition.

\oint_{\partial Ω} \frac{\partial V}{\partial n} d a = 0.

More generally, integrating Poisson’s equation gives

\oint_{\partial Ω} \frac{\partial V}{\partial n} d a = - \frac{1}{ε _{0}} Ω \int ρ d τ .

The net outward derivative must match enclosed charge. In addition, Neumann data determine potential only up to one additive constant. Voltage differences and electric field remain determined. One reference potential removes that constant; a grounded point or an imposed mean potential is enough.

Mixed problems use Dirichlet conditions on one part of the boundary and Neumann conditions on another. An isolated conductor introduces another common condition: potential is constant over its surface, but the constant is unknown; the integral of exterior normal field fixes the conductor’s stated total charge. Solving for that unknown conductor potential together with the interior potential preserves the physical charge constraint.

At a boundary between two ordinary dielectric materials, potential remains continuous in electrostatics:

V_{1} = V_{2} .

The normal displacement condition is

ε_{2} \frac{\partial V _{2}}{\partial n} - ε_{1} \frac{\partial V _{1}}{\partial n} = - σ_{f},

where $\overset{n}{^}$ points from material 1 to material 2 and $σ_{f}$ is free surface charge at the interface. Homogeneous-vacuum Laplace problems avoid this extra material condition. It becomes essential when a numerical model includes dielectric pieces or an interface that the grid crosses.

The uniqueness theorem makes boundary-value methods reliable. Suppose $V_{1}$ and $V_{2}$ satisfy the same Laplace equation and the same Dirichlet boundary values. Their difference $u = V_{1} - V_{2}$ satisfies

\nabla^{2} u = 0, u_{\partial Ω} = 0.

Green’s first identity gives

Ω \int ∣\nabla u ∣^{2} d τ = \oint_{\partial Ω} u \frac{\partial u}{\partial n} d a - Ω \int u \nabla^{2} u d τ = 0.

The integrand $∣\nabla u ∣^{2}$ is nonnegative, so it vanishes throughout the connected domain. Thus $u$ is constant; the zero boundary value makes that constant zero. The two solutions are identical. The proof for Neumann data leaves one arbitrary additive constant, matching the earlier physical interpretation.

Uniqueness changes the role of an approximate construction. A separated solution, a converged finite-difference solution, and an experimental model all describe the same physical potential when they satisfy the same equation and boundary data within their stated uncertainty. Agreement between methods detects boundary transcription mistakes, insufficient grid resolution, and measurements that perturb the apparatus.

One-dimensional solutions and their discrete form

Two broad parallel electrodes with negligible variation along their faces provide a one-dimensional region. Let $0 < x < d$ , with $V (0) = V_{A}$ and $V (d) = V_{B}$ . Laplace’s equation reduces to

\frac{d ^{2} V}{d x ^{2}} = 0.

Integrating twice gives

V (x) = V_{A} + \frac{V _{B} - V _{A}}{d} x, E_{x} = - \frac{V _{B} - V _{A}}{d} .

The potential changes linearly and the normal electric field is constant away from electrode edges. The result applies to the stated one-dimensional boundary geometry. Fringing restores two- and three-dimensional variation close to the ends of finite plates.

One-dimensional Laplace solution between broad parallel electrodes. The potential falls linearly from one conductor value to the other, and its constant slope sets the uniform interior field away from the plate edges.

A finite-difference grid reproduces this result without symbolic integration. For equally spaced interior nodes separated by $h$ , approximate the second derivative by

\frac{d ^{2} V}{d x ^{2}}_{i} \approx \frac{V _{i + 1} - 2 V _{i} + V _{i - 1}}{h ^{2}} .

The discrete Laplace equation is therefore

V_{i} = \frac{V _{i - 1} + V _{i + 1}}{2} .

Each interior value is the average of its two neighbors. Starting from the boundary values and applying this relation produces an arithmetic sequence. The exact linear continuum solution is reproduced at grid nodes for any uniform spacing, a direct test for a one-dimensional numerical code.

One-dimensional models also provide a calibration standard for field reconstruction. A central difference of a linear numerical profile gives the same field at every interior node. A reconstructed field that varies across a uniform one-dimensional reference identifies a coordinate error, an unequal grid spacing error, or a voltage conversion error before a more complicated two-dimensional solve is trusted.

Separation of variables in a rectangle

Separation of variables constructs analytic solutions when the domain and boundary data align with a coordinate system. Consider the rectangle

0 < x < a, 0 < y < b,

with three grounded edges and a prescribed top-edge potential:

V (0, y) = V (a, y) = V (x, 0) = 0, V (x, b) = f (x) .

The endpoint values of $f$ must agree with the side boundaries at the upper corners. A smooth physical electrode geometry rounds those corners, but an ideal mathematical rectangle can have abrupt boundary changes. The solution remains accurate away from the corner; numerical resolution and Fourier convergence require extra care near a discontinuity.

Rectangular Dirichlet problem for separation of variables. Three edges are grounded and the top edge carries a prescribed profile; the sine basis is chosen because it vanishes on the two vertical edges.

Seek a product form

V (x, y) = X (x) Y (y) .

Substitution into $\partial^{2} V / \partial x^{2} + \partial^{2} V / \partial y^{2} = 0$ and division by $X Y$ gives

\frac{X ^{''}}{X} = - \frac{Y ^{''}}{Y} = - k^{2} .

The choice $- k^{2}$ produces oscillatory $x$ dependence and hyperbolic $y$ dependence:

X^{''} + k^{2} X = 0, Y^{''} - k^{2} Y = 0.

The side conditions $X (0) = X (a) = 0$ allow nonzero solutions only at

k_{n} = \frac{nπ}{a}, X_{n} (x) = sin \frac{nπ x}{a}, n = 1, 2, 3, \dots

The lower grounded edge selects $Y_{n} (0) = 0$ , hence

Y_{n} (y) = sinh \frac{nπ y}{a} .

Each admissible product solution is zero on the three grounded edges. Linear superposition constructs the general top-boundary solution:

V (x, y) = n = 1 \sum \infty b_{n} sin \frac{nπ x}{a} \frac{sinh ( nπ y / a )}{sinh ( nπ b / a )} .

The coefficients are the sine-series coefficients of the imposed top-edge function:

b_{n} = \frac{2}{a} 0 \int a f (x) sin \frac{nπ x}{a} d x .

The denominator normalizes every mode to unit amplitude at $y = b$ . It prevents large hyperbolic factors from obscuring the boundary amplitude, and it makes the top-edge match transparent:

V (x, b) = n = 1 \sum \infty b_{n} sin \frac{nπ x}{a} = f (x) .

The first three lateral sine modes on the driven edge. Each mode vanishes at both side walls; higher mode number packs more half-waves across the width and decays faster into the interior.

The mode number controls penetration into the rectangle. At a distance $b - y$ below the driven edge, the factor

\frac{sinh ( nπ y / a )}{sinh ( nπ b / a )}

decreases approximately as $exp [- nπ (b - y) / a]$ when the region is tall enough for the exponential approximation. Fine lateral detail in the boundary pattern therefore remains close to its source edge. Broad features penetrate farther into the interior. This scale dependence explains why a coarse numerical grid can reproduce a smooth central potential while missing sharp boundary structure near an electrode gap.

Decay of separated modes with depth below the driven edge. Higher lateral mode number carries a shorter transverse scale and falls off faster, so only the lowest modes reach the deep interior.

A single-mode boundary illustrates the calculation without truncation.

f (x) = V_{0} sin \frac{π x}{a} .

Only $b_{1} = V_{0}$ is nonzero, so

V (x, y) = V_{0} sin \frac{π x}{a} \frac{sinh ( π y / a )}{sinh ( π b / a )} .

Differentiate the potential to obtain the electric field.

E_{x} = - \frac{π V _{0}}{a} cos \frac{π x}{a} \frac{sinh ( π y / a )}{sinh ( π b / a )},

E_{y} = - \frac{π V _{0}}{a} sin \frac{π x}{a} \frac{cosh ( π y / a )}{sinh ( π b / a )} .

The signs and component directions change across the rectangle. Along the central line $x = a /2$ , $E_{x} = 0$ by symmetry for this single mode, while $E_{y}$ is negative for positive $V_{0}$ . The top electrode is the highest-potential part of the boundary, so the field points away from it into lower potential.

Superposition handles a top-edge function with several Fourier components. For example, a boundary pattern containing mode 1 and mode 4 has one broad contribution and one narrow contribution. The mode-4 structure can dominate near the driven electrode while contributing little in the middle of a tall region. A measured interior map that retains large short-wavelength oscillations far from the electrode signals either a short domain, a different boundary condition, or measurement noise.

Separation of variables also clarifies what an analytic solution assumes. The boundary functions must be specified over the full boundary, coordinates must align with the domain shape, and the series must be evaluated to enough terms near sharp changes. A finite-difference solve is usually more efficient when conductor shapes are irregular, while the separated solution remains valuable as a reference case and as a check on numerical convergence.

Two-dimensional finite differences and conductor masks

A finite-difference method replaces the continuous domain by nodes $(x_{i}, y_{j})$ separated by $h_{x}$ and $h_{y}$ . The potential at each unknown interior node becomes one number in a coupled linear system. For equal spacing $h_{x} = h_{y} = h$ , centered second differences give

\frac{\partial ^{2} V}{\partial x ^{2}}_{i, j} \approx \frac{V _{i + 1, j} - 2 V _{i, j} + V _{i - 1, j}}{h ^{2}},

\frac{\partial ^{2} V}{\partial y ^{2}}_{i, j} \approx \frac{V _{i, j + 1} - 2 V _{i, j} + V _{i, j - 1}}{h ^{2}} .

Applying Laplace’s equation produces the five-point stencil.

V_{i, j} = \frac{V _{i + 1, j} + V _{i - 1, j} + V _{i, j + 1} + V _{i, j - 1}}{4} .

Every unknown interior node equals the average of its four cardinal neighbors. The local relations form one coupled system. Boundary nodes supply known values, and an iterative or direct linear-system method resolves the coupling.

The five-point Laplace stencil on an equal-spacing grid. The central unknown equals the average of its north, south, east, and west neighbors.

The average form is a discrete version of the mean-value property. If all boundary values lie between $V_{min}$ and $V_{max}$ , a correctly solved Dirichlet grid keeps every interior value inside that interval. Checking the extrema after every numerical solve catches swapped boundary signs, unintended node updates on a conductor, and array-index errors. The check establishes boundedness; an incorrect boundary transcription can still yield a bounded harmonic grid.

Grid spacing changes the weighting. With unequal rectangular spacings, the discrete equation becomes

V_{i, j} = \frac{( V _{i + 1, j} + V _{i - 1, j} ) / h _{x}^{2} + ( V _{i, j + 1} + V _{i, j - 1} ) / h _{y}^{2}}{2/ h _{x}^{2} + 2/ h _{y}^{2}} .

Neighbors in the direction with smaller spacing receive greater weight. A formula that averages all four neighbors equally on a stretched grid solves a different discrete problem and gives the wrong spatial scale. Store the coordinates or spacings with the potential array; raw node indices lack the physical geometry.

Second-order centered differences have truncation error proportional to $h_{x}^{2}$ and $h_{y}^{2}$ when the potential is sufficiently smooth. The actual solution error also contains geometry error, boundary-data error, iteration error, and measurement error. Refining a grid while retaining a blocky stair-step representation of a curved conductor improves one part of the model while leaving another part limited by boundary geometry. A convergence claim must identify which error is being reduced.

Dirichlet boundaries are imposed by writing their known potential into every boundary node and excluding those nodes from the update. An interior conductor is treated the same way: mark every node inside the conductor and every node intersecting its surface as prescribed-potential nodes. The neighboring vacuum nodes remain unknown. The mask must follow the physical conductor shape at the chosen resolution, including narrow gaps that can disappear when $h$ is too large.

A stair-step mask approximates a smooth boundary by grid-aligned segments. The resulting error is often first order near the physical surface even though the interior stencil is second order. Finer spacing reduces the stair-step size. Embedded-boundary and finite-element methods improve curved geometry when a Cartesian mask becomes inadequate, but their boundary data and convergence checks follow the same physical logic.

Neumann boundaries need a derivative equation in addition to the interior stencil. At a left edge with zero outward normal derivative, a centered ghost node construction gives

\frac{V _{1, j} - V _{- 1, j}}{2 h _{x}} = 0, V_{- 1, j} = V_{1, j} .

Substituting the ghost value into the boundary stencil produces an equation for the edge node. A one-sided derivative can also be used:

\frac{- 3 V _{0, j} + 4 V _{1, j} - V _{2, j}}{2 h _{x}} = 0.

The one-sided formula is second order for a smooth boundary. Sharp corners and incompatible meeting conditions reduce its accuracy. State the discrete form used at each boundary; a label such as “insulated edge” must be accompanied by the numerical equation.

Grid refinement must preserve physical dimensions. A coarse grid with $h = 2 mm$ and a refined grid with $h = 1 mm$ should use the same electrode coordinates, electrode potentials, and outer enclosure. Moving an electrode by half a cell during refinement changes the physical problem and can be mistaken for numerical error. Compare potentials at common physical locations after interpolation.

Relaxation, residuals, and convergence

The discrete Laplace equations form a sparse linear system. Direct matrix factorization is practical for modest grids, but relaxation methods use the local stencil repeatedly and scale well for large structured domains. Their output remains a numerical approximation. Convergence must be assessed with a residual and grid refinement. A visually smooth contour plot is insufficient evidence.

Jacobi iteration updates every unknown node from values stored during the previous sweep:

V_{i, j}^{(k + 1)} = \frac{V _{i + 1, j}^{(k)} + V _{i - 1, j}^{(k)} + V _{i, j + 1}^{(k)} + V _{i, j - 1}^{(k)}}{4} .

The method is simple because every update reads one old array and writes a second array. Its drawback is slow propagation of long-wavelength error across a large domain. Gauss–Seidel iteration updates one node at a time and uses freshly updated neighbor values as soon as they are available. Ordering the nodes row by row creates a sweep direction, yet the converged solution is independent of that ordering when the iteration converges.

Successive over-relaxation accelerates Gauss–Seidel by moving beyond the local average:

V_{i, j}^{(k + 1)} = (1 - ω) V_{i, j}^{(k)} + ω V_{i, j}^{⋆},

where $V_{i, j}^{⋆}$ is the Gauss–Seidel stencil value. Values $1 < ω < 2$ are over-relaxed. The best $ω$ depends on grid dimensions and geometry. A value too close to 2 can create oscillatory error or divergence. A documented iteration should record the update order, $ω$ , the initial guess, and which nodes were excluded by boundary masks.

An update size measures how much the stored array changes.

δ^{(k)} = i, j max ∣ V_{i, j}^{(k + 1)} - V_{i, j}^{(k)} ∣.

Use update size to detect stalled progress. A slowly varying error can change little from one sweep to the next while still violating Laplace’s equation over a large area. Use the discrete residual

r_{i, j} = \frac{V _{i + 1, j} - 2 V _{i, j} + V _{i - 1, j}}{h _{x}^{2}} + \frac{V _{i, j + 1} - 2 V _{i, j} + V _{i, j - 1}}{h _{y}^{2}} .

The residual has units of $V m^{- 2}$ . A dimensionless criterion divides its maximum by $V_{scale} / h_{min}^{2}$ , where $V_{scale}$ is a representative boundary voltage. The chosen tolerance must lie below the accuracy required from the final potential or field estimate.

A typical relaxation convergence history. The residual drops quickly while short-wavelength error is removed, then trails off slowly as broad modes persist; a small change per sweep can occur before the residual meets tolerance.

Residual maps localize remaining numerical trouble. A large residual band next to an electrode can signal an incorrectly applied boundary value, a coarse stair-step approximation, or an unresolved narrow gap. A checkerboard residual can indicate a parity or indexing error. A low residual everywhere with a solution that disagrees with measurements points toward boundary data, material assumptions, or apparatus geometry.

Convergence has three distinct tests.

Algebraic convergence: the residual and update meet their numerical tolerances.
Discretization convergence: halving the grid spacing changes the requested potential differences and field components by less than the stated spatial-error allowance.
Model convergence: changing the outer boundary, electrode shape representation, or measured boundary data within plausible uncertainty keeps the reported conclusion inside its uncertainty range.

Each test answers a different question. A million iterations can reduce algebraic error on a grid whose electrodes were placed one millimetre from their measured positions. Refining an inaccurate geometric model can produce a stable but physically misplaced solution. Numerical work should report the three checks separately.

Worked finite-difference solve on a square grid

Worked example. A square domain is meshed with four grid intervals per side. The two side boundaries and the bottom are grounded at $0 V$ ; the two noncorner top nodes are driven to $100 V$ (the top corners keep the side value $0 V$ ). Solve for the interior potential and estimate the vertical field.

Four interior nodes reduce to two unknowns by left-right symmetry: let $A$ be either upper interior node (next to the driven edge) and $B$ either lower one. Each node is the average of its four neighbors:

A = \frac{A + B + 100}{4}, B = \frac{A + B}{4} ⟹ 3 A - B = 100, 3 B - A = 0.

Then $A = 3 B$ gives $8 B = 100$ , so

A = 37.5 V, B = 12.5 V .

Both lie inside the boundary range $[0, 100] V$ , as the discrete maximum principle requires. Centered vertical-field estimates use the neighbors directly:

E_{y} (B) \approx - \frac{A - 0}{2 h} = - \frac{18.75}{h} V/m, E_{y} (A) \approx - \frac{100 - 12.5}{2 h} = - \frac{43.75}{h} V/m,

larger near the driven edge, where the local voltage difference is larger.

The direct equations solve the coarse discrete system. A physical square represented with more nodes develops additional spatial structure near the top corners. Grid refinement changes the discrete approximation to the same physical boundary data. It should retain the same outer dimensions and electrode positions in metres.

Jacobi iteration reaches the same pair from the zero interior guess.

A^{(k + 1)} = \frac{A ^{(k)} + B ^{(k)} + 100}{4}, B^{(k + 1)} = \frac{A ^{(k)} + B ^{(k)}}{4} .

The first three updates are

(A^{(1)}, B^{(1)}) = (25, 0),

(A^{(2)}, B^{(2)}) = (31.25, 6.25),

(A^{(3)}, B^{(3)}) = (34.375, 9.375) .

The upper value responds in the first sweep because it touches the driven edge. Lower-node response appears in the next sweep because Jacobi uses old neighbor values. A Gauss–Seidel sweep uses newly written values during the same pass, so its history changes while the converged values remain the same.

Residuals for the two symmetry-reduced equations are

r_{A} = 4 A - (A + B + 100), r_{B} = 4 B - (A + B) .

They vanish at the exact discrete values. At the first Jacobi update, both are $- 25 V$ in the unscaled four-neighbor form. Residual magnitude therefore measures the remaining imbalance of each local average equation. Dividing by $h^{2}$ produces the dimensional discrete Laplacian used on a general grid.

Late Jacobi updates become small, yet a residual tolerance remains necessary. Large domains retain broad error modes that change little from one sweep to the next. A solver can pass an update-size threshold while those modes still affect a requested voltage difference or derived field component. A whole-domain average field would likewise discard the spatial variation already present in the discrete potential, seen above in the larger $E_{y}$ near the driven edge.

Boundary, symmetry, maximum, residual, and refinement checks turn this compact calculation into a reference. A larger apparatus adds curved electrodes, leads, dielectric supports, and three-dimensional fringing. The same checks apply after those physical details are encoded in the boundary model.

Field reconstruction, measurement, and model validation

Laplace solvers return potential values. Electric field is a derived quantity and needs its own numerical treatment. At an interior Cartesian node, second-order centered differences give

E_{x, i, j} \approx - \frac{V _{i + 1, j} - V _{i - 1, j}}{2 h _{x}}, E_{y, i, j} \approx - \frac{V _{i, j + 1} - V _{i, j - 1}}{2 h _{y}} .

The derivative uses a difference of nearby potential values, so it is more sensitive to local noise and grid error than potential itself. A potential array can appear converged and smooth while its reconstructed field remains too noisy for a surface-charge estimate or a trajectory calculation.

Centered formulas require nodes on both sides of the evaluation point. Near a conductor boundary, use a one-sided derivative or an interpolation to the physical surface. Suppose a grounded conductor lies at $x = 0$ , vacuum occupies $x > 0$ , and grid nodes lie at $0$ , $h$ , and $2 h$ . The outward normal field from the metal into vacuum is approximated by

E_{n} \approx - \frac{- 3 V _{0} + 4 V _{1} - V _{2}}{2 h} .

In vacuum, the exterior surface charge density follows

σ = ε_{0} E_{n},

with the normal directed outward from the conductor. The sign convention must be written beside the result. A reversed normal reverses both $E_{n}$ and $σ$ without changing the physical potential.

Field reconstruction has a built-in integral check. For a path from point $A$ to point $B$ , sum the projected field values along short path segments:

ℓ \sum E_{ℓ} \cdot Δ ℓ_{ℓ} \approx - (V_{B} - V_{A}) .

The same potential array provides both sides. Large disagreement indicates an inconsistent derivative stencil, a coordinate-spacing error, or a path interpolation error. Closed-path sums should approach zero:

closed \sum E_{ℓ} \cdot Δ ℓ_{ℓ} \approx 0.

The test is valuable after field interpolation onto particle paths or probe locations, because the interpolation can break the exact discrete gradient relationship that existed at grid nodes.

Path-integral check on a reconstructed field. Summing the field projected along each segment from A to B recovers the potential drop between the two points, and a closed loop must sum to zero.

Electrostatic field must also be irrotational. A discrete curl estimate in two dimensions is

(\nabla \times E)_{z} \approx \frac{E _{y, i + 1, j} - E _{y, i - 1, j}}{2 h _{x}} - \frac{E _{x, i, j + 1} - E _{x, i, j - 1}}{2 h _{y}} .

Computing $E$ from one consistently differentiated potential array makes this quantity small up to truncation and roundoff. Combining field components measured by different probes can produce a nonzero discrete curl through calibration offsets or inconsistent coordinate registration. The curl map provides a diagnostic of the measured field data before fitting a potential surface.

Potential uncertainty propagates strongly through differentiation. Independent voltage uncertainty $u_{V}$ at two equal-spacing nodes gives approximately

u_{E_{x}} \approx \frac{u _{V}}{2 h _{x}} .

Halving probe separation improves spatial resolution and doubles this noise-derived field uncertainty when voltage precision is unchanged. A measurement layout must choose a spacing that resolves the variation of interest while producing voltage differences large enough for the instrument. Repeated readings reduce random voltage noise. A common reference offset or an uncertain physical separation persists across repeated readings.

Potential probes should sample known physical coordinates. A measurement grid drawn on paper becomes unreliable if the probe tip is displaced by a support, if conductor edges are referenced from a different datum, or if the probe height changes over a nominally two-dimensional apparatus. Record each coordinate relative to the same origin used by the numerical model. A high-input-impedance voltmeter reduces loading, but an isolated small electrode can still shift when a probe, cable, or nearby grounded object is introduced.

Model-to-measurement comparison begins with interpolation. Most probes lie between numerical nodes, so evaluate $V_{model} (x_{k}, y_{k})$ with bilinear interpolation or a higher-order method appropriate to grid smoothness. For each measurement, form a standardized residual

z_{k} = \frac{V _{k} - V _{model} ( x _{k} , y _{k} )}{u _{V, k}^{2} + u _{model, k}^{2}} .

The model uncertainty can include grid-refinement change, boundary-voltage uncertainty, coordinate uncertainty, and solver tolerance. A set of standardized residuals scattered around zero with size near one is consistent with the stated uncertainty model. A spatial pattern carries more information than the average: a sign change near one electrode can indicate a misplaced edge, while a common shift can indicate a reference-voltage offset.

Validation needs independent boundary information. Verify the actual voltage on each driven conductor under the same load and measurement arrangement used for the interior map. Record ground connections, supply polarity, electrode dimensions, dielectric pieces, probe height, and outer enclosure geometry. An outer boundary placed near the region of interest can distort a solution even when the local electrode values are exact. Repeat the solve after moving the artificial outer boundary outward; the change in the reported quantity is a boundary-truncation sensitivity estimate.

A complete validation record includes the following items.

Boundary record: measured electrode voltages, conductor locations, material interfaces, ground connection, and outer-domain placement.
Numerical record: grid spacings, mask construction, derivative stencil, iteration method, relaxation parameter, residual tolerance, and grid-refinement comparison.
Measurement record: probe coordinates, voltage calibration, input impedance, reference terminal, repeated readings, and probe-height control.
Comparison record: interpolated model values, residuals, uncertainty components, rejected measurements with reasons, and sensitivity to plausible boundary changes.

The records make a disagreement diagnosable. A potential mismatch that grows near a curved conductor points toward geometric representation. A uniform mismatch over the map points toward a voltage-reference offset. Mismatch that changes after moving a probe points toward probe loading or coordinate error. Mismatch confined to a narrow gap can arise from insufficient grid spacing or unmodeled three-dimensional fringing. Each pattern directs a specific revision to the physical model or measurement procedure.

Practical limits and error diagnosis

Laplace’s equation is appropriate only where free volume charge is negligible and the permittivity model matches the material. A charge cloud, resistive medium carrying current, or strongly varying dielectric changes the governing equation. In a piecewise dielectric with no free volume charge, the more general relation is

\nabla \cdot (ε \nabla V) = 0.

Using $\nabla^{2} V = 0$ across a discontinuous permittivity treats both materials as if they had the same response. The resulting potential can violate the normal displacement condition at the interface even when every grid residual is small.

Coordinate systems impose limits. A two-dimensional Cartesian solve assumes that geometry and boundary values remain unchanged along the omitted direction. Leads, electrode thickness, mounting screws, and finite plate length can break that assumption. A two-dimensional solution remains accurate when measurements show that the omitted variation is below the required accuracy in the study region. Cross-sections near leads and edges commonly require a three-dimensional model or an empirical boundary correction.

Boundary errors often dominate interior numerical accuracy. A potential value that is wrong by one percent along a nearby electrode can alter the entire harmonic solution by a comparable scale. In contrast, a residual tolerance made one hundred times smaller after grid convergence may change the reported voltage difference by far less. Allocate effort according to sensitivity: measure nearby boundary voltages and geometry first, then refine the grid and solver until numerical error falls below those physical uncertainties.

The maximum principle, residual map, path-integral check, grid-refinement comparison, and measurement residuals form a compact diagnostic set. They test boundedness, local equation satisfaction, gradient consistency, discretization dependence, and agreement with the apparatus. Used together, they distinguish a converged numerical array from a validated electrostatic boundary-value solution.

Boundary-condition audit and sensitivity analysis.

Every reported quantity should be linked back to the boundary information that controls it. A potential at a target point may depend strongly on one nearby electrode voltage and weakly on a distant enclosure wall. A normal field near a sharp conductor can depend strongly on the local surface representation and weakly on the far boundary. Sensitivity analysis identifies those dependencies before time is spent refining an insensitive part of the model.

Let $Q$ denote a reported quantity, such as $V$ at a target point, $E_{n}$ at a conductor surface, or a potential difference between probes. For a model parameter $p$ , estimate local sensitivity with paired solves:

S_{p} \approx \frac{Q ( p + Δ p ) - Q ( p - Δ p )}{2Δ p} .

Choose $Δ p$ large enough that the change in $Q$ exceeds solver noise and small enough that the local response remains representative. Parameters can include electrode voltage, electrode position, outer-boundary distance, dielectric permittivity, probe coordinate, and an unknown isolated-conductor potential. The product $∣ S_{p} ∣ u_{p}$ estimates how much uncertainty $u_{p}$ in that parameter contributes to the reported result.

Common sources of correlation deserve explicit treatment. One supply can set several electrode voltages, giving their calibration errors a common sign. One mechanical datum can define several electrode positions. A voltage probe and a supply monitor can share the same reference terminal. Treating such errors as independent and adding them in quadrature understates the uncertainty of a difference or gradient. Sensitivity coefficients make the correlations visible: shared parameters enter the model once, with their combined effect evaluated directly.

Mesh convergence uses a sequence of grids with spacing $2 h$ , $h$ , and $h /2$ for the same physical boundary model. Let $Q_{2 h}$ , $Q_{h}$ , and $Q_{h /2}$ be the reported quantity after algebraic convergence. For a smooth second-order calculation, an observed order estimate is

p_{obs} = log_{2} \frac{Q _{2 h} - Q _{h}}{Q _{h} - Q _{h /2}} .

Values near two support the expected centered-difference behavior in the resolved region. A lower observed order can arise from a curved stair-step boundary, an abrupt boundary function, a one-sided derivative, or a quantity evaluated too close to a corner. The result identifies the feature limiting the claimed accuracy.

When observed order is established, Richardson extrapolation estimates the zero-spacing value:

Q_{⋆} \approx Q_{h /2} + \frac{Q _{h /2} - Q _{h}}{2 ^{p_{obs}} - 1} .

Use the extrapolated value as a diagnostic. Physical resolution still requires a mesh that resolves every feature controlling the reported quantity. A narrow electrode gap represented by one or two cells often lacks a geometry error that follows a smooth power law. Compare the reported quantity at successively smaller spacings.

Boundary audit begins with a node-by-node check. List each boundary segment, its physical meaning, the measured or imposed value, its units, the coordinate range, and the discrete equation used there. Dirichlet segments need voltage and geometric location. Neumann segments need outward-normal convention and derivative value. Isolated conductors need surface-potential constancy and the total-charge constraint. Dielectric interfaces need permittivity values and free surface charge. The list prevents a common failure: correct equations applied to a boundary that was assigned the wrong physical type.

The audit should include corners. A mathematical corner lies on two boundary segments and can receive two incompatible values when data are copied from different measurement records. State which condition owns the corner node or replace the sharp idealization with a small rounded geometry. A high residual or strong mesh dependence near a corner can be physical behavior of the ideal boundary, a discretization artifact, or a transcription error. Comparing nearby probe readings and refined grids separates these cases.

Analytic cases provide strong verification tests. The one-dimensional linear solution checks uniform spacing, Dirichlet handling, and reconstructed constant field. The rectangular single-sine solution checks two-dimensional stencils, boundary-series evaluation, and vertical attenuation. A zero boundary problem checks that iteration preserves zero potential. A mirror-symmetric boundary problem checks coordinate registration and left-right equality. Run these cases before interpreting a complicated electrode geometry.

Measurement validation benefits from a separate prediction set. Use one group of potential readings to establish uncertain boundary quantities or geometry, then reserve additional probe locations for comparison. Reusing every reading to tune the model can hide an incorrect boundary representation. The reserved points should sample the high-gradient region, a central region, and a region near any material interface. Their standardized residuals reveal whether a calibration chosen from one part of the domain transfers across the apparatus.

Sensitivity, grid order, boundary audit, and independent readings complete the verification chain. The final claim should name the solved quantity, physical domain, boundary conditions, mesh spacing, algebraic residual, grid-refinement change, dominant parameter sensitivities, and measurement comparison. A student or reviewer can then distinguish an exact analytic result, a numerically converged approximation, and a model whose boundary assumptions have been tested against the apparatus.

Interpreting systematic comparison patterns.

Residual patterns provide more information than a single root-mean-square number. Plot measured minus modeled potential at the physical probe coordinates, using the same color scale and sign convention across repeated runs. A nearly uniform residual commonly points to a reference-voltage offset or a ground connection difference. A residual that changes approximately linearly across the domain can arise from a probe-coordinate scale error, a tilted apparatus, or a distant boundary whose imposed potential differs from the model.

Local residual concentration near a driven conductor often traces a geometric issue. The conductor edge may be displaced, rounded, thickened, or connected to a lead absent from the two-dimensional boundary. Local residuals near a dielectric interface can indicate an incorrect permittivity, an omitted free surface charge, or a probe height that crosses into a different material region. A sharp residual at one isolated coordinate deserves a repeat measurement before it is treated as a model defect; contact resistance, digitizer noise, and probe placement can create single-point outliers.

Repeat the comparison after changing one controlled parameter. Moving an outer boundary outward tests enclosure sensitivity. Reversing a supply polarity tests sign registration. Translating a probe pattern by a measured small distance tests coordinate alignment. Reducing the grid spacing tests spatial discretization. The observed residual change should follow the predicted direction for the proposed cause. A candidate explanation that fails this controlled test should be removed from the error budget.

Residual analysis also separates uncertainty in potential from uncertainty in derived field. A smooth model can agree with potential probes at their locations while its derivative differs between probes because the sampling gap is too large. Add measurements where the reconstructed field changes quickly, especially near electrode edges and narrow gaps. Conversely, dense potential samples in a low-gradient interior region add little information about boundary geometry. Allocate probe locations according to the sensitivity of the reported quantity and the expected gradient structure.