Solution of Partial Differential Equations

Heat flowing through a metal plate, the electric potential around a charged conductor, the pressure of groundwater under a d a m: each is described by an equation in two or three space variables.

This final unit shows how to turn such a partial differential equation (PDE) into a system of ordinary algebraic equations on a grid, and solve it with the tools of Unit 4.

6.1 Introduction and classification

A second-order linear PDE in two variables, A uₓₓ + B uₓᵧ + C uᵧᵧ + (lower terms) = 0, is classified by B^{2} − 4AC:

This unit treats the elliptic equations of Laplace and Poisson.

For them the solution is determined by the values of u on the boundary of the region (a Dirichlet problem), and the interior values depend on all boundary values at once, so the whole grid is solved together.

6.2 Deriving difference equations

Cover the region with a grid of spacing h in both directions and write u_{ij} = u (x_{i} ,\ y_{j}). Adding the Taylor series of u (x \pm h) gives the central second difference

uₓₓ ≈ (uᵢ₊₁,ⱼ − 2uᵢⱼ + uᵢ₋₁,ⱼ)/h², and likewise uᵧᵧ ≈ (uᵢ,ⱼ₊₁ − 2uᵢⱼ + uᵢ,ⱼ₋₁)/h²,

each with error O (h^{2}). Substituting into the PDE:

• Laplace (u_{xx} + uᵧᵧ = 0): uᵢ₊₁,ⱼ + uᵢ₋₁,ⱼ + uᵢ,ⱼ₊₁ + uᵢ,ⱼ₋₁ − 4uᵢⱼ = 0, i.e. u_{ij} = \frac{1}{4} ( \mathrm{sum} of the four neighbours): each interior value is the average of its neighbours (the five-point formula).
• Poisson (u_{xx} + uᵧᵧ = f): uᵢ₊₁,ⱼ + uᵢ₋₁,ⱼ + uᵢ,ⱼ₊₁ + uᵢ,ⱼ₋₁ − 4uᵢⱼ = h^{2} f_{ij}.

Writing one equation for each interior grid point gives a linear system, with 4 on the diagonal and −1 for each neighbour: diagonally dominant, so the Gauss–Seidel method of Unit 4 converges.

Applied to this system, Gauss–Seidel is called Liebmann's iterative method: sweep over the grid and replace each value by the average of its neighbours until nothing changes.

6.3 Laplace's equation on a square plate

A finer grid has 3 \times 3 interior points, nine unknowns u_{ij} ( i = row from the top, 1 to 3 ;\ j = column, 1 to 3), with the top boundary at 100 and the others 0.

Other boundary conditions. If part of the boundary is insulated (\frac{\partial u}{\partial n} = 0, a Neumann condition), a ghost point outside the grid is introduced with the value of its mirror image inside, so the neighbour sum doubles one interior term.

On irregular boundaries, special difference formulas with unequal spacing are used.

6.4 Poisson's equation

6.5 Accuracy and practical notes

• The difference equations are accurate to O (h^{2}); halving h quarters the truncation error but quadruples the number of unknowns.
• Direct solution (Gauss elimination, banded LU) costs more memory than iteration for large grids; iterative methods exploit the sparsity (only five nonzeros per row).
• For parabolic (heat) problems, the explicit scheme u_{i}^{j+1} = r u_{i-1}^{j} + (1 - 2 r) u_{i}^{j} + r u_{i+1}^{j} with r = \frac{c^{2} \Delta t}{h^{2}} is stable only if r \leq \frac{1}{2}; the implicit Crank–Nicolson scheme is stable for all r. Hyperbolic problems need the Courant condition \frac{c \Delta t}{h} \leq 1. These are outside the syllabus but explain why time-dependent problems need more care than the steady-state ones here.

6.6 Seeing the ideas: the five-point stencil, the solved plate and Gauss-Seidel on the grid

A cup of coffee loses heat to its surroundings; the temperature inside any object settles to a smooth pattern once its boundaries are held steady. That settled pattern obeys Laplace's equation, and the grid methods of this unit compute it.

(a) Heat flowing out of a cup

(b) Five-point stencil

6.6.1 Eight uses of Laplace and Poisson equations

1. Steady heat in a plate or wall. After enough time the temperature in a wall stops changing and obeys Laplace's equation. Engineers use grid methods to estimate heat loss through brick walls and cooling fins of electronics. This helps in designing warm houses for the Nepali hills.
2. Electric potential. The voltage in the empty space between conductors satisfies Laplace's equation. The grid method gives the potential around chips, cables and insulators. Designers find where the field is strongest and sparks may occur.
3. Groundwater and seepage. Water moving slowly through soil under a d a m follows Laplace's equation. The solution shows the pressure and the flow paths. Dam designers use it to avoid leaks and failure.
4. Fluid flow. Ideal flow of a fluid around a pipe bend or a bridge pier can be described by a potential that satisfies Laplace's equation. The grid solution shows speed and pressure. It is the first step before more advanced models.
5. Poisson's equation and charge. If charges or heat sources are present, the equation becomes Poisson's equation. The source term is added to the stencil. Computing the field of a charged body or heat inside an electronic part uses it.
6. Deflection of membranes and plates. The shape of a stretched drum skin under pressure satisfies Poisson's equation. A civil engineer can estimate how a thin slab bends under a load. The same equation appears in many different fields.
7. Image processing. Filling a hole in a photograph can be done by solving Laplace's equation with the edge of the hole as the boundary. The result blends smoothly with its surroundings. Photo-editing tools use this idea.
8. Teaching and checking bigger programs. A simple plate with known boundary values has a solution one can verify by symmetry, as in the exercise here. It is the standard test for larger simulation software. Learning it builds the habit of checking results.

(a) A thermal image: a temperature field

(b) Centre-point error of the five-point scheme

Formula and fact sheet

• uₓₓ ≈ (uᵢ₊₁,ⱼ − 2uᵢⱼ + uᵢ₋₁,ⱼ)/h².
• Laplace: uᵢⱼ = ¼(uᵢ₊₁,ⱼ + uᵢ₋₁,ⱼ + uᵢ,ⱼ₊₁ + uᵢ,ⱼ₋₁).
• Poisson: uᵢ₊₁,ⱼ + uᵢ₋₁,ⱼ + uᵢ,ⱼ₊₁ + uᵢ,ⱼ₋₁ − 4uᵢⱼ = h² fᵢⱼ.
• Classification by B^{2} − 4AC: < 0 elliptic, = 0 parabolic, > 0 hyperbolic.
• Liebmann = Gauss–Seidel on the grid.

Practice

1. Classify u_{xx} + uᵧᵧ = 0 ,\ u_{t} = u_{xx} and u_{tt} = u_{xx}.
2. Derive the five-point difference formula for Laplace's equation.
3. A square plate has its top edge at 100 ^{\circ} and the other edges at 0 ^{\circ}. Using a mesh with four interior points, find the temperatures.
4. Do two Liebmann sweeps for the 3 \times 3 interior problem starting from zero, in row order, and compare with the exact solution.
5. Write the difference equations for \nabla^{2} u = - 10 (x^{2} + y^{2} + 10) with h = 1 on [0 ,\ 3] \times [0 ,\ 3] ,\ u = 0 on the boundary, and solve them.
6. Why is the Gauss–Seidel iteration guaranteed to converge for the five-point system?
7. Explain what changes in the difference equation at an insulated edge.
8. What is the effect of halving h on the accuracy and on the size of the system?

Model answers: medium (5 marks) and long (10 marks) questions

Q1 (5 marks). A square plate has its top edge at 100 degrees and the other three edges at 0 degrees. Find the temperatures at the nine interior points of a 4 x 4 mesh using the five-point formula and symmetry.

Formula. For Laplace's equation, u (i ,\ j) = \frac{u (i + 1 ,\ j) + u (i - 1 ,\ j) + u (i ,\ j + 1) + u (i ,\ j - 1)}{4}. By left-right symmetry, name the unknowns: top row a ,\ b ,\ a; middle row c ,\ d ,\ c; bottom row e ,\ f ,\ e.

Equations. a = \frac{100 + b + c}{4} ;\ b = \frac{100 + 2 a + d}{4} ;\ c = \frac{a + e + d}{4} ;\ d = \frac{b + f + 2 c}{4} ;\ e = \frac{c + f}{4} ;\ f = \frac{d + 2 e}{4}.

Solving gives a = 42.857 ,\ b = 52.679 ,\ c = 18.750 ,\ d = 25.000 ,\ e = 7.143 ,\ f = 9.821.

Check. a: \frac{100 + 52.679 + 18.75}{4} = 42.857. d: \frac{52.679 + 9.821 + 37.5}{4} = 25. The centre is one quarter of 100, as superposition of the four edges (each contributing 25 at the centre by symmetry) requires (Figure 6.1).

Conclusion temperature falls smoothly from the hot edge towards the cold edges.

Q2 (5 marks). Explain how the Gauss-Seidel iteration is used for the plate problem. Show its early progress.

Method. Start with a guess (all zeros) and sweep through the grid, replacing each interior value by the average of its four neighbours using the newest values; repeat until changes are below a tolerance.

The system is diagonally dominant (each unknown has coefficient 4 against four neighbours of 1), so convergence is guaranteed.

The error shrinks by a factor of about two each sweep here (Figure 6.2).

Speed-up successive over-relaxation (SOR) with a factor between 1 and 2 converges much faster on fine grids.

Conclusion iterative solution is the practical way to solve large grids.

Q3 (10 marks). Classify second-order partial differential equations. Derive the difference equations for Laplace's and Poisson's equations and explain how a Dirichlet problem is solved numerically.

Classification. For A u_{x} x + B u_{x} y + C u_{y} y + lower terms = 0, the sign of B^{2} - 4AC decides the type: elliptic if negative (Laplace u_{x} x + u_{y} y = 0, Poisson u_{x} x + u_{y} y = f; steady state), parabolic if zero (heat u_{t} = c^{2} u_{x} x; diffusion in time), hyperbolic if positive (wave u_{t} t = c^{2} u_{x} x; vibration). This unit treats the elliptic type.

Difference equations. Cover the region with a grid of spacing h. Taylor series give u_{x} x is about \frac{u (i + 1 ,\ j) - 2 u (i ,\ j) + u (i - 1 ,\ j)}{h^{2}} and likewise for u_{y} y, with error O (h^{2}).

Adding them gives for Laplace: u (i ,\ j) = \frac{u (i + 1 ,\ j) + u (i - 1 ,\ j) + u (i ,\ j + 1) + u (i ,\ j - 1)}{4}, the average of the four neighbours (the stencil of Plate 6.1).

For Poisson: u (i ,\ j) = [sum of four neighbours - h^{2} f (i ,\ j) ] / 4.

Dirichlet problem. The values of u on the boundary are given. Each interior point gives one linear equation, so an n x n interior grid gives n^{2} equations (nine for the 3 x 3 plate of Figure 6.1).

They are solved by Gauss elimination or by iteration (Gauss-Seidel, SOR; Figure 6.2), the latter being preferred for fine grids.

Worked Poisson example. For u_{x} x + u_{y} y = - 10 on a square with u = 0 on all edges and h = 1, the single interior point of a 3 x 3 grid has 4u = 0 + 0 + 0 + 0 - h^{2} (- 10), so u = 2.5. Halving h gives a more accurate answer with more points.

Accuracy. The error is O (h^{2}); halving h cuts it to about a quarter but multiplies the number of unknowns by four. Corner and curved boundaries need special treatment.

Conclusion Elliptic problems reduce to large sparse linear systems that links this unit to Unit 4; symmetry checks and a convergence test on the sweeps make the answers trustworthy.