Poisson's Equation in Electrostatics

The Big Picture

Take a region of empty space and sprinkle some electric charge into it. Two questions immediately arise: where is the voltage high, where is it low, and which way will a test charge accelerate? Both questions are answered by a single scalar field, the electric potential $V(\mathbf{r})$. And that potential is locked to the charge distribution by exactly one equation, Poisson's equation:

Here $\rho(\mathbf{r})$ is the charge per unit volume, $\varepsilon_0 \approx 8.854 \times 10^{-12}\,\mathrm{F/m}$ is the permittivity of free space, and $\nabla^2$ is the Laplacian. Read it left-to-right: tell me how the potential curves at a point, and I'll tell you how much charge is sitting there. Read it right-to-left: tell me where the charges are, and I'll tell you how the potential bends.

Why this equation is so powerful. All of static electricity — every capacitor, every semiconductor device, every lightning rod, every MRI gradient coil — is governed by this single PDE. Solve it once for your geometry, and you know the voltage everywhere, the field everywhere ($\mathbf{E} = -\nabla V$), and the energy stored in the system. There is no "step 2."

The Curvature-Equals-Source Intuition

Forget the Greek letters for a moment. The Laplacian of a function at a point compares the function's value there to the average of its immediate neighbours:

If a point is exactly at the average of its surroundings, $\nabla^2 V = 0$. If it's a peak (higher than average), the Laplacian is negative. If it's a valley, the Laplacian is positive. Now look at Poisson:

Wherever there is positive charge ($\rho > 0$), the right-hand side is negative, so the left-hand side is negative — the potential at that point sits above the average of its neighbours. Positive charge pushes the voltage up locally; negative charge pulls it down. No charge? Then the potential is exactly the average of its surroundings (this is Laplace's equation, and the "mean value property" of harmonic functions).

Analogy 1 — the rubber sheet. Imagine a tight rubber sheet. Now press your finger up from below at some location. The sheet bulges up there and gently flattens far away. That bulge is the potential $V$; your finger is the positive charge. Press from above and you get a dip — that's negative charge. Multiple fingers? Their bulges and dips superpose linearly (because Poisson is linear).

Analogy 2 — temperature with heaters. Replace voltage by temperature and charge by heat sources. In steady state, the temperature distribution in a slab satisfies $\nabla^2 T = -q/k$: heat sources push the temperature above the local average; cold spots are sinks. Exactly the same equation, different letters. This is why the same numerical solver works for both physics problems.

From Gauss's Law to Poisson's Equation

Poisson's equation isn't pulled out of a hat — it is a one-line consequence of Gauss's law plus the definition of the potential. The chain is short enough to memorize.

1. Gauss's law (differential form). The divergence of the electric field equals the charge density divided by $\varepsilon_0$: $\nabla \cdot \mathbf{E} = \rho/\varepsilon_0$.
2. Electric field from a potential. In electrostatics there are no time-varying magnetic fields, so the field is curl-free, which means it can be written as the gradient of a scalar: $\mathbf{E} = -\nabla V$. The minus sign is convention so that $\mathbf{E}$ points from high voltage to low voltage (positive test charges roll downhill in $V$).
3. Substitute step 2 into step 1. $\nabla \cdot \mathbf{E} = \nabla \cdot (-\nabla V) = -\nabla^2 V$. Setting that equal to $\rho/\varepsilon_0$ and flipping signs gives the result:

Notice what got hidden: the operator $\nabla \cdot \nabla$ is the Laplacian. In Cartesian coordinates it expands to second partial derivatives:

The reason every interesting electrostatics problem boils down to second derivatives is precisely this substitution chain. The field is a gradient; the source rule is on the divergence; chained together, you get a second-order operator on the potential.

Boundary Conditions Make It Concrete

Poisson's equation by itself has infinitely many solutions: you can add any harmonic function (a solution to $\nabla^2 U = 0$) to a particular solution and get another valid one. What pins down the unique physical answer is the boundary data:

Practical rule of thumb. Inside any solid conductor in electrostatic equilibrium, $V$ is constant. This means the surface of every conductor in your problem becomes a Dirichlet boundary at the conductor's voltage. Most engineering problems reduce to: "solve Poisson in the air/insulator region, with the metal boundaries held at known voltages."

Worked Example: Uniformly Charged Sphere

This single example unlocks an entire class of problems. Take a ball of radius $R$ with total charge $Q$ spread uniformly through its volume. By spherical symmetry, the potential $V$ depends only on the distance $r$ from the centre, so the Laplacian in spherical coordinates collapses to:

Inside the sphere the charge density is the constant $\rho_0 = Q/(\tfrac{4}{3}\pi R^3)$. Outside there is no charge, so Poisson reduces to Laplace. We solve the two regions separately and glue them together at $r = R$ by demanding $V$ and $dV/dr$ match (no surface charge at the boundary, no sudden jump in field).

After integrating twice in each region and applying the matching conditions plus the physical requirement $V(\infty) = 0$, you get the famous piecewise answer:

Three things to notice in the visualization above:

• Outside the sphere, the answer is identical to a point charge sitting at the centre. The sphere's extended structure is invisible from outside — Newton's "shell theorem," rediscovered electrically.
• Inside, the potential is a downward-opening parabola in $r$, peaking at the centre. The electric field $E = -dV/dr = (Q/4\pi\varepsilon_0)\,r/R^3$ grows linearly from zero at the centre to its maximum at the surface.
• At the surface ($r = R$), both pieces give $V(R) = Q/(4\pi\varepsilon_0 R)$ and both derivatives match — that's the gluing condition working.

Interactive: Charges → Potential Field

Reading the equation is one thing; seeing the potential field rearrange itself as you drop charges is something else. The widget below solves the 2D Poisson equation on a 64×64 grid in real time using Jacobi relaxation (the same algorithm you'll see in Python next). Click anywhere to drop a charge; toggle between positive and negative; watch the equipotential contours and the field arrows reorganize themselves around every new source.

What to play with. Place two charges far apart — see how their influences add (superposition). Place +/- close together — that's a dipole, and the field becomes strongly directional. Place a ring of positive charges — the interior becomes nearly equipotential (a model for a Faraday cage). Crank iterations up: more iterations means closer to the true solution. Drop it way down to see how Jacobi propagates information one cell per sweep (look at the early-iteration speckle pattern near isolated charges).

Parallel-Plate Capacitor with Bulk Charge

Here is a problem you can do entirely by hand, and which models a real device. Take two large flat conductors separated by a gap $d$. Hold the left plate at $V = 0$ and the right plate at $V = V_0$. Suppose the space between them is filled with a uniform space charge of density $\rho$ (think: electron cloud in a vacuum tube, or doped semiconductor depletion region). What does the potential look like in between?

Far from the edges of the plates, the field only depends on $x$ (the direction across the gap). Poisson reduces to a 1D ODE:

Integrate twice: $V(x) = -\rho x^2/(2\varepsilon_0) + A x + B$. The boundary conditions $V(0) = 0$ and $V(d) = V_0$ fix $B = 0$ and $A = V_0/d + \rho d/(2\varepsilon_0)$:

Two limiting cases make this physical:

Numbers for the second case. Take $V_0 = 100\,\text{V}$, $d = 1\,\text{cm}$, $\rho = 10^{-6}\,\text{C/m}^3$. Then at the midplane $V(d/2) \approx 51.4\,\text{V}$ (not 50 V — the parabola bulges upward). The field at the negatively-biased plate is $E(0) \approx -1.06\times 10^{4}\,\text{V/m}$ and at the positive plate $E(d) \approx -9.44\times 10^{3}\,\text{V/m}$. The asymmetry is the signature of bulk charge — and it's why semiconductor depletion-region calculations all start with Poisson.

Computing It: Finite-Difference Jacobi Solver

For anything more complicated than a sphere or parallel plates, you go numerical. The most transparent approach is finite differences + Jacobi relaxation:

1. Replace the continuous Laplacian by its three-point (1D) or five-point (2D) stencil.
2. Solve algebraically for the centre cell: $V_i = \tfrac{1}{2}(V_{i-1} + V_{i+1} + h^2 f_i)$ in 1D, $V_{ij} = \tfrac{1}{4}(\text{4 neighbours} + h^2 f_{ij})$ in 2D.
3. Sweep that update across the grid over and over until nothing changes. Because the discrete Laplacian is a strict contraction (after fixing the boundary), this loop converges from any starting guess.

Below is a complete, runnable Python implementation. It includes a manufactured-solution test: pick $V_{\text{exact}}(x) = \sin(\pi x)$, compute $-V''_{\text{exact}} = \pi^2 \sin(\pi x)$, feed that into the solver, and verify it recovers $\sin(\pi x)$ back. The max error of $\sim 3 \times 10^{-4}$ is the discretization error from the finite stencil (it shrinks as $h^2$ as you refine the grid).

np.array([0,1,4,1,0])  →  shape (5,)

Continuous:  -V''(x) = f(x)
Discrete:    -(V_{i-1} - 2V_i + V_{i+1})/h² = f_i

rho_over_eps.shape = (50,)  →  N = 50

V = [0., 0., 0., ..., 0.]   (length 50)

V[-1] = 0.0  →  V = [0., 0., ..., 0.]

V       = [0., 3., 5., 3., 0.]
V_new = V.copy()  →  separate array

If V = [0,0,0,0,0] and h²f = [_,1,4,1,_] (interior),
after ONE sweep V_new = [0, 0.5, 2.0, 0.5, 0].
Converged answer: V = [0, 3, 5, 3, 0].
Check:  -(V[1] - 2V[2] + V[3])/1² = -(3 - 10 + 3) = 4  ✓ matches f[2].

err history: 0.45, 0.31, 0.22, ..., 9.2e-10

V_exact(0.5) = sin(π/2) = 1.0

np.linspace(0, 1, 5) → [0., 0.25, 0.5, 0.75, 1.0]

f at x=0.5: π² × sin(π/2) = π² × 1 ≈ 9.8696

returns V ≈ [0, 0.064, 0.127, ..., 0.064, 0]

exact[25] = sin(π·0.510) ≈ 0.9995

max(|V − exact|) ≈ 3.42e-04

1import numpy as np
2
3def solve_poisson_1d(rho_over_eps, h, V_left=0.0, V_right=0.0,
4                     max_iter=20000, tol=1e-9):
5    """Solve -V''(x) = rho(x)/eps0 on a 1D grid using Jacobi relaxation.
6
7    Args:
8        rho_over_eps: array of source values f_i = rho_i / eps0, shape (N,)
9        h:            grid spacing
10        V_left:       Dirichlet boundary V(0)
11        V_right:      Dirichlet boundary V(L)
12    Returns:
13        V: array of potential values, shape (N,)
14    """
15    N = rho_over_eps.shape[0]
16    V = np.zeros(N)
17    V[0]  = V_left
18    V[-1] = V_right
19
20    for iteration in range(max_iter):
21        V_new = V.copy()
22        # Jacobi update: V_i = 0.5 * (V_{i-1} + V_{i+1} + h^2 * f_i)
23        V_new[1:-1] = 0.5 * (V[:-2] + V[2:] + h * h * rho_over_eps[1:-1])
24        err = np.max(np.abs(V_new - V))
25        V = V_new
26        if err < tol:
27            break
28    return V, iteration + 1
29
30# --- Test: source with exact answer V(x) = sin(pi x) on [0, 1] ---
31N = 50
32h = 1.0 / (N - 1)
33x = np.linspace(0, 1, N)
34f = np.pi**2 * np.sin(np.pi * x)         # -V'' = pi^2 sin(pi x)
35
36V, n_iter = solve_poisson_1d(f, h)
37exact     = np.sin(np.pi * x)
38max_err   = np.max(np.abs(V - exact))
39print(f"Converged in {n_iter} iterations, max error = {max_err:.2e}")

Why Jacobi? It is the cleanest possible expression of "curvature equals source." Each grid cell literally averages its neighbours and adds a charge bump. The downside is it is slow: information about boundary conditions diffuses by one cell per sweep, so a grid of $N$ points takes $O(N^2)$ iterations to converge in 1D. Real-world solvers use multigrid (which fixes this to $O(N)$) or conjugate-gradient methods. But every one of those advanced methods is built on top of the relaxation idea you see here.

PyTorch on the GPU

The exact same Jacobi update vectorizes trivially to 2D, and PyTorch lets you run it on a GPU with one line change. Why bother? Two reasons. First, modern problems with millions of grid points are GPU-natural — the inner loop is just neighbour-averaging, which is exactly what convolutional kernels do. Second, the iteration is differentiable: PyTorch's autograd lets you backpropagate through it, which means you can use Poisson's equation as a layer inside a neural network (physics-informed networks, learned PDE solvers, neural-style Poisson image editing).

x = torch.zeros(32, 32, device='cuda')

rho_over_eps.device  →  device(type='cuda', index=0)

rho_over_eps.shape = torch.Size([32, 32])  →  N = 32

V = tensor([[0,0,0,...], [0,0,0,...], ...], device='cuda:0')

V_new shares no memory with V; modifying V_new[i,j] does not change V.

Centre cell (16,16) with neighbours all 0 and h²f = 0 → V_new[16,16] = 0.
Charge cell (10,12) with h²·500 source → V_new[10,12] grows over time toward ≈ 0.33.

diff = 0.00000732  →  still iterating

rho_eps.device  →  device(type='cuda', index=0)   # if a GPU is present

rho_eps[10,12] = 500.0  →  cell (10,12) of the source becomes +500

rho_eps[22,20] = -500.0

(V, n_iter) ≈ (32×32 tensor, 1932)

1import torch
2
3def solve_poisson_2d_gpu(rho_over_eps, h, max_iter=5000, tol=1e-8):
4    """Solve nabla^2 V = -rho/eps0 on a 2D grid using vectorized Jacobi on GPU.
5
6    Args:
7        rho_over_eps: (N, N) tensor of source values
8        h:            grid spacing
9    Returns:
10        V: (N, N) tensor with Dirichlet V=0 on the boundary
11    """
12    device = rho_over_eps.device
13    N      = rho_over_eps.shape[0]
14    V      = torch.zeros(N, N, device=device)
15
16    for it in range(max_iter):
17        V_new           = V.clone()
18        V_new[1:-1, 1:-1] = 0.25 * (
19            V[:-2, 1:-1] + V[2:, 1:-1]      # vertical neighbours
20          + V[1:-1, :-2] + V[1:-1, 2:]      # horizontal neighbours
21          + h * h * rho_over_eps[1:-1, 1:-1]
22        )
23        diff = (V_new - V).abs().max().item()
24        V    = V_new
25        if diff < tol:
26            break
27    return V, it + 1
28
29# --- Two opposite point charges ---
30torch.manual_seed(0)
31N        = 32
32h        = 1.0 / (N - 1)
33rho_eps  = torch.zeros(N, N, device="cuda" if torch.cuda.is_available() else "cpu")
34rho_eps[10, 12] = +500.0
35rho_eps[22, 20] = -500.0
36
37V, n_iter = solve_poisson_2d_gpu(rho_eps, h)
38print(f"iters={n_iter}, V[10,12]={V[10,12]:.4f}, V[22,20]={V[22,20]:.4f}, V_center={V[16,16]:.4e}")

Running this on the two-charge dipole above gives $V[10,12] \approx +0.328$, $V[22,20] \approx -0.322$, and the potential at the geometric centre is nearly zero ($\approx 2.7\times 10^{-3}$). The small asymmetry between the two extrema comes from the fact that the two charges aren't equally far from the boundary — the Dirichlet edge pulls the potential slightly more strongly toward zero on whichever side is closer.

Method of Images: A Trick from Symmetry

Before computers, a generation of physicists used a beautiful trick to solve otherwise-hard Poisson problems: the method of images. Suppose you have a single positive charge $+q$ sitting at height $h$ above an infinite grounded conducting plane ($V = 0$ on the plane). The potential in the upper half-space is hard — there are induced surface charges on the plane that you don't know in advance.

The trick: replace the conductor by an imaginary charge $-q$ at depth $-h$. The plane becomes the perpendicular bisector of the line joining $+q$ and $-q$, and by symmetry the potential along that bisector is exactly zero — exactly the boundary condition we need! Now there's no conductor anywhere: just two free-space charges, and the potential in the upper half-plane is simply the sum of their Coulomb potentials.

Why this works mathematically: Poisson's equation has unique solutions given the sources and the boundary conditions. The image construction reproduces both the correct sources (only the real charge is in the physical region) and the correct boundary ($V = 0$ on the plane). Uniqueness then guarantees it's the right answer.

Try this in the interactive widget above. Place a positive charge somewhere in the upper half of the grid and a negative charge of the same magnitude directly below it, mirrored across the horizontal midline. Crank up the iterations. You'll see the equipotential line along the midline is essentially zero — the widget reproduces the method-of-images solution from first principles.

Connections to Machine Learning

Poisson's equation isn't just a physics curiosity — it appears in modern machine learning in surprisingly direct ways.

• Poisson image editing. When you paste a region from one photo into another, the seam usually shows. The Perez/Gangnet 2003 algorithm solves a Poisson equation where the source $f$ is the gradient of the source patch and the Dirichlet boundary is the destination image at the seam. The result blends seamlessly — and the solver is exactly the 2D Jacobi (or multigrid) you saw above. Modern photo editors (including the "Healing Brush") use this.
• Physics-informed neural networks (PINNs). A PINN trains a neural network $V_\theta(x,y)$ by minimizing the residual $\|\nabla^2 V_\theta + \rho/\varepsilon_0\|^2$ at sampled points, plus a boundary loss. The Laplacian inside the loss is computed by autograd — exactly the same tool that backpropagates through every other neural net. The PyTorch code above is the "classical" teacher; a PINN is the neural student.
• Graph Laplacian regularization. In semi-supervised learning, you spread labels across a graph by solving $L f = 0$ where $L$ is the graph Laplacian. This is exactly Laplace's equation on a discrete graph — Poisson's equation with zero source. The labelled nodes act as Dirichlet boundary conditions, and the propagation rule is the discrete "average of neighbours" you saw at the start of this section.
• Diffusion models. Score-based generative models can be reframed as solving a partial differential equation related to the heat (Fokker-Planck) equation, and the steady-state of that equation is Poisson's equation. The same mathematics that places voltage around charges places probability mass around training data points.

Summary

The takeaway in one sentence. Once you can read $\nabla^2 V = -\rho/\varepsilon_0$ as "the potential at a point exceeds the average of its neighbours in proportion to the local charge," every numerical method, every analytic trick, and every machine-learning application of this equation becomes a variation on a single, intuitive theme.