Prev
Unit 10 · Boundary Value Problems and Partial Differential Equations
Elliptic Equations: Laplace and Poisson
60m100 XP
Next
🌡️Boundary Value Problems and Partial Differential Equations

Elliptic Equations: Laplace and Poisson

The five-point stencil, Dirichlet and Neumann boundaries, sparse assembly and solution, seepage under a dam

Idaho, 5 June 1976
The Teton Dam, a 93 m earth-fill dam on the Teton River, has just been completed by the US Bureau of Reclamation, and its reservoir is filling for the first time. On the morning of Saturday 5 June, workers notice muddy water seeping from the downstream face — water carrying soil with it. Within hours a hole appears and widens; bulldozers sent to plug it are swallowed. Shortly before noon the embankment breaches, and the reservoir empties into the valley below. Eleven people are killed and thousands of homes destroyed. Investigations found that water had flowed through the wind-blown silt used in the dam's core and through fissures in the rock of the abutment, eroding a pipe through the embankment from the inside — internal erosion, or piping. Where water flows through soil, and how steeply its pressure drops along the way, is governed by one of the simplest equations in mathematical physics. The total head hh in steady groundwater flow satisfies Laplace's equation, ∇2h=0\nabla^2h = 0. Its gradient at the downstream ground surface — the exit gradient — is what decides whether soil grains are lifted out and piping begins.
Elliptic equations describe equilibrium. Laplace's equation ∇2u=0\nabla^2u = 0 and Poisson's equation −∇2u=f-\nabla^2u = f govern steady heat conduction, groundwater flow, electrostatic potential, potential flow around a wing, the Prandtl stress function in torsion and the deflection of a membrane. There is no time and no direction of propagation: the value at every point depends on all the boundary data at once, and the solution is as smooth as the data allow. Their discretisations give one large sparse linear system — exactly the problems Units 3 and 4 were built to solve.
💡The five-point stencil
On a square grid of spacing hh, central differences in xx and yy give −∇2u≈4ui,j−ui−1,j−ui+1,j−ui,j−1−ui,j+1h2=fi,j.-\nabla^2u \approx \dfrac{4u_{i,j} - u_{i-1,j} - u_{i+1,j} - u_{i,j-1} - u_{i,j+1}}{h^2} = f_{i,j}. For Laplace's equation (f=0f = 0) every value is the average of its four neighbours — the discrete version of the mean-value property. With m2m^2 interior unknowns the matrix is m2×m2m^2 \times m^2, symmetric positive definite, with five nonzeros per row: kron(I, T) + kron(T, I) from lesson 3-4. The error is O(h2)O(h^2): for −∇2u=2π2sin⁡πxsin⁡πy-\nabla^2u = 2\pi^2\sin\pi x\sin\pi y on the unit square, halving hh divides the error by 4.01, 4.00, 4.00 from 225 to 16,129 unknowns.
Boundary conditions. A Dirichlet condition prescribes the value — a face held at 100 °C, the water head along a reservoir bed — and moves known values to the right-hand side. A Neumann condition prescribes the normal derivative — an insulated face, an impervious boundary, a line of symmetry, all with ∂u/∂n=0\partial u/\partial n = 0. Implement it with a ghost point outside the domain mirrored from inside, which makes the centred difference of the normal derivative zero. In a plate with 100 °C and 0 °C on opposite edges and insulated top and bottom, the computed temperature is linear in xx to rounding error and the heat flux is exactly kΔT/Lk\Delta T/L — a simple check that the Neumann treatment is right.
Derivation
Why a mirrored ghost point imposes zero flux
Suppose the boundary row j=0j = 0 is insulated, ∂u/∂y=0\partial u/\partial y = 0 there. Keep the boundary values as unknowns and add a ghost row j=−1j = -1 outside the domain. The centred difference for the normal derivative at the boundary is second-order accurate (lesson 1-3):
∂u∂y∣i,0≈ui,1−ui,−12h=0⟹ui,−1=ui,1.\frac{\partial u}{\partial y}\Big|_{i,0} \approx \frac{u_{i,1} - u_{i,-1}}{2h} = 0 \qquad\Longrightarrow\qquad u_{i,-1} = u_{i,1}.
The ghost value is simply the mirror image of the first interior value. Substitute it into the ordinary five-point equation at the boundary point:
4ui,0−ui−1,0−ui+1,0−2ui,1=h2fi,0.4u_{i,0} - u_{i-1,0} - u_{i+1,0} - 2u_{i,1} = h^2f_{i,0}.
The ghost has disappeared and the interior neighbour counts twice. (Dividing that row by 2 restores the symmetry of the matrix.) A line of symmetry in a problem is treated the same way, which halves the size of many models.
Seepage from a flow net. In the flow nets of the callout below, each 'square' has equal width and length, so Darcy's law (lesson 4-4) gives it a flow of k×(head drop H/Nd)k\times(\text{head drop } H/N_d) per unit thickness, the same in every channel. NfN_f channels in parallel carry q=kH Nf/Ndq = kH\,N_f/N_d.
When a theory is sufficiently general to cover many fields of application, it acquires some "truth" from each of them.
– Richard Hamming, Methods of Mathematics Applied to Calculus, Probability, and Statistics (1985)
📐Seepage under a concrete dam
A dam with an impervious base 24 m wide sits on a 12 m permeable layer over rock; 10 m of water upstream, none downstream; soil permeability 10−510^{-5} m/s. The head field on a 0.25 m grid — 17,000 unknowns, solved by spsolve — gives:
  • seepage under the dam of 2.98 m³ per day per metre of dam;
  • a head of exactly 5.000 m under the middle of the base, as symmetry requires;
  • an exit gradient of 0.55 one metre downstream of the toe. With a critical gradient of about 1 for typical soils, the factor of safety against piping is 1.81 — below the 3 to 4 that designers usually require. A sheet-pile cutoff or a longer base would be needed.
The exit gradient is a derivative of the solution, and right at the toe it is theoretically infinite; engineers evaluate it over a small finite distance, which is why its value depends on resolution.
🕰Flow nets
Before computers, seepage was solved graphically. Engineers sketched a flow net of equipotentials and flow lines crossing at right angles into curvilinear squares — a graphical solution of Laplace's equation popularised by Arthur Casagrande in the 1930s. The seepage per metre is q=kH Nf/Ndq = kH\,N_f/N_d, the ratio of flow channels to potential drops, and the exit gradient is read off the last square. A good sketch is within a few percent of the grid solution; drawing one is still the best way to develop intuition for where water goes.
💡Where you will use this
Geotechnical Engineering: seepage under dams and into excavations, uplift pressures, exit gradients and piping. Heat transfer: steady conduction in heat sinks, turbine blades, electronics. Electromagnetics: potentials and fields in insulators, motors and transformers. Structures: torsion of non-circular sections (Prandtl), membrane deflection. Aerodynamics: potential flow, the starting point of panel methods.
📚Further reading
  • Randall LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations (SIAM, 2007), chapter 3, 'Elliptic Equations' — the five-point stencil, its accuracy and boundary conditions.
  • Harry Cedergren, Seepage, Drainage, and Flow Nets (3rd ed., Wiley, 1989) — the engineer's book on seepage, flow nets and exit gradients.
  • Milton Harr, Groundwater and Seepage (McGraw-Hill, 1962; Dover reprint 1991) — exact conformal-mapping solutions, including the sheet-pile problem of lesson 4-4.
↓But there is a problem
Square grids are wonderful on rectangles. Real domains have curved boundaries — the arched face of a dam, a fillet in a machine part, a tunnel in rock — and mixed materials with interfaces that do not line up with grid lines. Forcing them onto a grid gives staircase boundaries and smeared interfaces. There is another way to compute a field, and it is older than the computer. Instead of asking the equation to hold at every grid point, ask which shape stores the least energy — the question a bent strip of wood answers by itself (lesson 7-3). In 1909 a 31-year-old Swiss physicist in Göttingen, gravely ill with tuberculosis, used that question to compute the patterns that sand forms on a vibrating square plate — a problem that had defeated mathematicians for a century.