Approximating phenomena that change across space and time — grids, finite differences, the heat equation, the wave equation, the Laplace equation, and the question of stability
This chapter introduces numerical methods for partial differential equations as an accessible, carefully paced preview. Partial differential equations describe phenomena that change across both space and time simultaneously — heat flowing through a rod, waves traveling through water, temperature settling across a metal plate. These equations are among the most important mathematical tools in physics, engineering, biology, climate science, and computational modeling. They are also, in most cases, impossible to solve exactly. Numerical methods turn them into something computable.
Chapter 12 showed how to approximate solutions to ordinary differential equations step by step in time. This chapter extends that idea into two dimensions: we will approximate solutions that depend on both position and time, or on two spatial coordinates at once. The central tool is the finite difference method, which replaces continuous derivatives with differences computed on a grid of discrete points.
This chapter is intentionally a preview, not a full PDE course. The goal is to give students a clear and honest picture of how numerical PDE methods work, what they can do, what can go wrong, and where the subject leads. Students who continue into numerical analysis, scientific computing, computational physics, or engineering will develop these ideas in much greater depth. Here, the emphasis is on understanding the structure of the problem, the geometry of the grid, and the logic of the approximation.
A metal rod is heated at one end and held at a fixed temperature at the other. Heat flows along the rod over time. How does the temperature at each interior point change from one moment to the next? You know from calculus that temperature satisfies a partial differential equation, the heat equation. But that equation has no simple closed-form solution for a general initial temperature profile. How could you compute approximate temperatures at each point of the rod at each moment in time? This is the kind of question numerical PDE methods answer.
Every applied field that studies continuous physical systems eventually encounters partial differential equations. Structural engineers model stress in bridges using elliptic PDEs. Meteorologists simulate atmospheric temperature and pressure using hyperbolic PDEs. Climate scientists track heat flow through ocean layers using parabolic PDEs. Medical imaging reconstructs tissue density by solving inverse PDE problems numerically. Financial derivatives are priced using the Black-Scholes PDE.
None of these problems can be solved analytically for realistic domains and initial conditions. Numerical PDE methods make them computable. A student who understands how a finite difference method discretizes a PDE, what stability means in this context, and how grid spacing and time steps control accuracy is prepared to read, use, and critically evaluate the computational models that underlie modern science and engineering.
By the end of this chapter, you should be able to:
partial differential equation (PDE), ordinary differential equation (ODE), grid, mesh, node, grid spacing, time step, discretization, finite difference method, explicit method, implicit method, boundary condition, initial condition, Dirichlet boundary condition, Neumann boundary condition, heat equation, parabolic PDE, wave equation, hyperbolic PDE, Laplace equation, elliptic PDE, stability, CFL condition, numerical diffusion, convergence, truncation error, iterative grid method, Gauss-Seidel iteration on a grid
An ordinary differential equation involves a function of a single variable and its derivatives with respect to that variable. We write things like
\[ \frac{dy}{dt} = f(t, y) \]
and seek a function $y(t)$ that depends only on time. Chapter 12 showed how Euler's method and Runge-Kutta methods approximate solutions to such equations.
A partial differential equation involves a function of two or more variables and its partial derivatives with respect to those variables. The simplest examples involve a function $u$ that depends on both position $x$ and time $t$, and we write expressions like
\[ \frac{\partial u}{\partial t} = k \frac{\partial^2 u}{\partial x^2} \]
This is the one-dimensional heat equation. The function $u(x,t)$ represents temperature at position $x$ and time $t$, and the equation says that the rate at which temperature changes in time is proportional to the second derivative of temperature in space. Heat diffuses from regions of high temperature to regions of low temperature, and the equation captures that process exactly.
The difficulty is that even this simple-looking equation is rarely solvable in closed form except for very special initial and boundary conditions. For realistic problems — irregular domains, complex boundary conditions, nonlinear terms, multiple spatial dimensions — analytical solutions are essentially unavailable. This is the fundamental reason numerical PDE methods exist.
An ODE asks for a function of one variable. A PDE asks for a function of two, three, or four variables simultaneously. The solution surface is far more complex. The space of possible initial and boundary conditions is far richer. The interactions between spatial structure and temporal evolution require that we track not just one evolving quantity but an entire field of quantities at every point in space at every moment in time.
Numerical methods handle this by replacing the continuous field $u(x,t)$ with a finite array of values $u_i^n$, where $i$ indexes the position in space and $n$ indexes the moment in time. The PDE becomes a system of algebraic equations relating these values, and we solve that system step by step.
This chapter focuses on three classical PDEs, each representing a different physical and mathematical structure.
Each type requires a different numerical strategy. Understanding why each strategy fits its PDE type is a key goal of this chapter.
The heat, wave, and Laplace equations are the simplest representatives of their PDE families. They are chosen for this chapter precisely because their behavior is well understood, which makes it possible to see clearly how each numerical method matches the structure of its equation.
The first step in any numerical PDE method is to lay down a grid — a discrete set of points in the domain where we will compute the approximate solution.
Consider a one-dimensional spatial domain: a rod of length $L$ running from $x = 0$ to $x = L$. We divide this interval into $N$ equal subintervals of width
\[ h = \frac{L}{N} \]
The grid points are the positions $x_i = i \cdot h$ for $i = 0, 1, 2, \ldots, N$. At each of these positions we will track an approximate temperature value.
Now add time. Suppose we want to compute the temperature from time $t = 0$ to time $t = T$. We divide the time interval into $M$ equal steps of width
\[ k = \frac{T}{M} \]
The time levels are $t_n = n \cdot k$ for $n = 0, 1, 2, \ldots, M$.
The result is a two-dimensional grid of points $(x_i, t_n)$, with $i$ ranging over space positions and $n$ ranging over time levels. At each grid point, we compute an approximation $u_i^n \approx u(x_i, t_n)$.
Draw a rectangle with the horizontal axis labeled $x$ (running from $0$ to $L$) and the vertical axis labeled $t$ (running from $0$ to $T$). Mark vertical grid lines at $x_0 = 0$, $x_1 = h$, $x_2 = 2h$, and so on. Mark horizontal grid lines at $t_0 = 0$, $t_1 = k$, $t_2 = 2k$, and so on. At each intersection, mark a node. Shade the bottom row (the initial condition $t = 0$) and the left and right columns (the boundary conditions). Label an interior node $(x_i, t_n)$ and its neighbors $(x_{i-1}, t_n)$, $(x_{i+1}, t_n)$, and $(x_i, t_{n+1})$.
When the problem involves two spatial dimensions — a flat plate, for example — the grid becomes a mesh of points $(x_i, y_j)$ covering a rectangular domain. The solution $u_{i,j}^n$ now has three indices: two for space and one for time. The ideas are the same, but the arrays are larger and the stencils involve more neighbors.
Smaller grid spacing $h$ and smaller time step $k$ generally produce more accurate approximations but require more computation. Choosing appropriate grid spacing involves a tradeoff between accuracy and cost. But there is also a deeper constraint: for some methods, the ratio $k/h^2$ must stay below a certain threshold to prevent the computation from becoming unstable. This is the stability issue we will examine carefully in Section 13.7.
The points on a grid are called nodes or grid points. The values computed at those points are called nodal values. The spacing $h$ is sometimes called the mesh size or grid spacing. The time step is usually called $k$ or $\Delta t$. The spatial step is usually called $h$ or $\Delta x$.
Once the grid is established, the next step is to replace the continuous partial derivatives in the PDE with finite difference approximations computed at grid points. This process is called discretization.
The forward difference approximation for the first partial derivative of $u$ with respect to $t$ at the point $(x_i, t_n)$ is
\[ \frac{\partial u}{\partial t}\bigg|_{(x_i, t_n)} \approx \frac{u_i^{n+1} - u_i^n}{k} \]
This is a forward difference in time. It is first-order accurate: the truncation error is $O(k)$.
The central difference approximation for the second partial derivative of $u$ with respect to $x$ at the point $(x_i, t_n)$ is
\[ \frac{\partial^2 u}{\partial x^2}\bigg|_{(x_i, t_n)} \approx \frac{u_{i-1}^n - 2u_i^n + u_{i+1}^n}{h^2} \]
This is a centered difference in space. It uses the values at the left neighbor, the current node, and the right neighbor, all at the same time level $t_n$. It is second-order accurate: the truncation error is $O(h^2)$.
Both formulas follow directly from Taylor series, exactly as in Chapter 6. Expanding $u(x_i, t_{n+1})$ and $u(x_i, t_{n-1})$ in Taylor series centered at $(x_i, t_n)$, and expanding $u(x_{i-1}, t_n)$ and $u(x_{i+1}, t_n)$ in Taylor series centered at $(x_i, t_n)$, gives the approximations above along with explicit error terms. The details are the same as the numerical differentiation analysis in Chapter 6, extended to partial derivatives.
The set of grid points used in a finite difference approximation is called a stencil. The central difference second derivative uses a three-point stencil: the node itself and its two immediate spatial neighbors. For time-stepping schemes, the stencil specifies which time levels appear in the formula.
By replacing partial derivatives with finite differences, we transform the PDE — an equation relating continuous functions — into a system of algebraic equations relating discrete nodal values. Instead of asking for a function $u(x,t)$, we ask for a large array of numbers $\{u_i^n\}$. The PDE becomes a rule for computing each new row of this array from the previous rows.
Before applying finite differences to specific PDEs, we collect the key formulas and their accuracy.
| Approximation | Formula | Order of Accuracy |
|---|---|---|
| Forward difference in time | $\dfrac{\partial u}{\partial t}\Big|_i^n \approx \dfrac{u_i^{n+1} - u_i^n}{k}$ | $O(k)$ |
| Backward difference in time | $\dfrac{\partial u}{\partial t}\Big|_i^n \approx \dfrac{u_i^n - u_i^{n-1}}{k}$ | $O(k)$ |
| Centered difference in time | $\dfrac{\partial u}{\partial t}\Big|_i^n \approx \dfrac{u_i^{n+1} - u_i^{n-1}}{2k}$ | $O(k^2)$ |
| Centered second difference in space | $\dfrac{\partial^2 u}{\partial x^2}\Big|_i^n \approx \dfrac{u_{i-1}^n - 2u_i^n + u_{i+1}^n}{h^2}$ | $O(h^2)$ |
| Centered first difference in space | $\dfrac{\partial u}{\partial x}\Big|_i^n \approx \dfrac{u_{i+1}^n - u_{i-1}^n}{2h}$ | $O(h^2)$ |
These formulas are the building blocks for all the numerical PDE methods in this chapter. They should feel familiar from Chapter 6. The difference here is that $u$ has two indices — one for space, one for time — and we apply difference quotients to one index while holding the other fixed.
When we write $u_i^n$, we mean the approximate value of the solution at the $i$-th spatial grid point at the $n$-th time step. The superscript is a time index, not an exponent. This notation takes some getting used to but is standard in numerical PDE literature.
A PDE alone does not have a unique solution. It requires boundary conditions that specify what the solution does at the edges of the domain, and for time-dependent problems, initial conditions that specify the solution at time $t = 0$.
For a time-dependent PDE on the interval $[0, L]$, the initial condition specifies the value of the solution at every spatial point at time zero:
\[ u(x, 0) = f(x) \quad \text{for } 0 \leq x \leq L \]
In the discrete setting, this means we know all values $u_i^0 = f(x_i)$ for $i = 0, 1, \ldots, N$. These values fill the bottom row of our grid and serve as the starting point for all subsequent time steps.
A Dirichlet boundary condition specifies the value of the solution at the boundary. For a rod with fixed temperatures at both ends:
\[ u(0, t) = \alpha \quad \text{and} \quad u(L, t) = \beta \quad \text{for all } t \geq 0 \]
In the discrete setting, this means $u_0^n = \alpha$ and $u_N^n = \beta$ for all time steps $n$. These values are fixed throughout the computation; we never update the boundary nodes using the interior scheme.
A Neumann boundary condition specifies the derivative of the solution at the boundary. For a rod with an insulated right end (no heat flux), we have
\[ \frac{\partial u}{\partial x}(L, t) = 0 \quad \text{for all } t \geq 0 \]
In the discrete setting, this is implemented using a one-sided finite difference at the boundary. A common approach is to introduce a ghost point just outside the domain and use the zero-derivative condition to relate the ghost point value to its interior neighbor.
Mixed boundary conditions specify the value at one end and the derivative at the other. Many physical problems have mixed boundary conditions.
Boundary conditions are not merely mathematical formalities. They constrain the solution and, in the numerical setting, they determine the values held fixed throughout the computation. An incorrect or inconsistent boundary condition will corrupt the entire computation, not just the boundary rows. Checking that boundary conditions are correctly implemented is one of the first things to verify when a numerical PDE solution looks wrong.
Revisit the space-time grid from Section 13.2. Highlight the bottom row in blue and label it "initial condition: $u_i^0 = f(x_i)$." Highlight the left column in red and label it "left boundary: $u_0^n = \alpha$." Highlight the right column in red and label it "right boundary: $u_N^n = \beta$." Shade the interior nodes in gray and note that these are computed by the finite difference scheme.
The one-dimensional heat equation is
\[ \frac{\partial u}{\partial t} = \kappa \frac{\partial^2 u}{\partial x^2} \]
where $u(x,t)$ is temperature at position $x$ and time $t$, and $\kappa > 0$ is the thermal diffusivity, a positive constant that describes how quickly heat flows through the material. The equation says: the rate of temperature change at a point equals $\kappa$ times the concavity of the temperature distribution at that point. Where the temperature profile is concave down, the temperature will fall; where it is concave up, the temperature will rise.
Substituting the forward difference in time and the centered difference in space into the heat equation gives
\[ \frac{u_i^{n+1} - u_i^n}{k} = \kappa \frac{u_{i-1}^n - 2u_i^n + u_{i+1}^n}{h^2} \]
Solving for the unknown $u_i^{n+1}$:
\[ u_i^{n+1} = u_i^n + r\left(u_{i-1}^n - 2u_i^n + u_{i+1}^n\right) \]
where we define the dimensionless parameter
\[ r = \frac{\kappa k}{h^2} \]
This formula is explicit: every quantity on the right side is known (it comes from the current or previous time level), so $u_i^{n+1}$ can be computed directly without solving any system of equations. We apply this formula for every interior node $i = 1, 2, \ldots, N-1$, holding the boundary values $u_0^n$ and $u_N^n$ fixed, and advance one time step at a time.
The update formula says that the new temperature at node $i$ is the old temperature at $i$ plus $r$ times the second difference of the current temperatures. The second difference $u_{i-1}^n - 2u_i^n + u_{i+1}^n$ is the discrete analog of the second derivative. It is positive when the profile is concave up (the node is below its neighbors), pulling the temperature up. It is negative when the profile is concave down, pulling the temperature down. This is exactly the physical idea of heat diffusion.
Algorithm: Explicit Finite Difference for the Heat Equation
Purpose:
Approximate u(x,t) satisfying u_t = kappa * u_xx
with given initial and Dirichlet boundary conditions.
Inputs:
kappa — thermal diffusivity (positive constant)
L — length of spatial domain [0, L]
T — total time to simulate
N — number of spatial subintervals
M — number of time steps
f(x) — initial condition function, u(x, 0) = f(x)
alpha — left boundary value, u(0, t) = alpha
beta — right boundary value, u(L, t) = beta
Compute:
h = L / N
k = T / M
r = kappa * k / h^2
Check stability:
If r > 0.5, issue a warning: method may be unstable.
Initialize:
For i = 0 to N:
u[i] = f(i * h)
u[0] = alpha
u[N] = beta
Time-stepping loop:
For n = 0 to M-1:
For i = 1 to N-1:
u_new[i] = u[i] + r * (u[i-1] - 2*u[i] + u[i+1])
u_new[0] = alpha
u_new[N] = beta
Replace u with u_new
Output:
The array u, representing the approximate solution at time T.
Reliability notes:
The method is first-order accurate in time and second-order accurate in space.
Stability requires r <= 0.5.
Smaller h and k improve accuracy but increase computational cost.The stability condition $r \leq \tfrac{1}{2}$ is not a suggestion — it is a requirement. If this condition is violated, the explicit method produces oscillations that grow without bound, completely corrupting the computed solution even though the true solution is smooth and well-behaved.
Suppose we use $h = 0.25$ but choose a time step $k = 0.05$, giving $r = 1 \cdot 0.05 / 0.0625 = 0.8 > 0.5$. Even starting from the smooth sine initial condition, the computed temperatures will oscillate and blow up within a few steps. The update formula amplifies small errors rather than damping them. The computed solution has no physical meaning.
This condition can be derived rigorously using Fourier analysis of the error. Each spatial Fourier mode of the error is multiplied by a factor called the amplification factor at each time step. For the explicit heat equation method, the amplification factor for mode $m$ is
\[ G_m = 1 - 4r\sin^2\left(\frac{m\pi h}{2L}\right) \]
For stability, we require $|G_m| \leq 1$ for all modes $m$. The most restrictive mode gives the condition $|1 - 4r| \leq 1$, which simplifies to
\[ r \leq \frac{1}{2} \]
This is called a von Neumann stability condition after the mathematician John von Neumann, who developed this analysis method in the 1940s.
If you halve the spatial step $h$, you must reduce $k$ by a factor of four to maintain $r \leq 0.5$, since $r = \kappa k / h^2$. This means that achieving twice the spatial resolution requires four times as many time steps. The cost of explicit methods for the heat equation grows rapidly as the grid is refined.
One way to escape the stability constraint is to use an implicit method, where the spatial differences are evaluated at the new time level $t_{n+1}$ rather than the old time level $t_n$. The Crank-Nicolson method, for example, averages the spatial differences at both time levels:
\[ \frac{u_i^{n+1} - u_i^n}{k} = \frac{\kappa}{2}\left(\frac{u_{i-1}^n - 2u_i^n + u_{i+1}^n}{h^2} + \frac{u_{i-1}^{n+1} - 2u_i^{n+1} + u_{i+1}^{n+1}}{h^2}\right) \]
This method is unconditionally stable for any $r > 0$, at the cost of requiring the solution of a tridiagonal linear system at each time step. The trade-off is algorithmic complexity for freedom from the time step constraint. Implicit methods are central to professional PDE solvers but require linear algebra tools beyond what this introductory chapter develops in detail.
Never choose a time step for the heat equation without checking the stability condition $r = \kappa k / h^2 \leq 1/2$. A numerical solution that looks smooth in the early steps but develops wild oscillations later is a classic sign of instability. Stability must be checked analytically before any computation begins, not discovered after the computation fails.
The one-dimensional wave equation is
\[ \frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2} \]
where $u(x,t)$ represents the displacement of a vibrating string at position $x$ and time $t$, and $c > 0$ is the wave speed. Unlike the heat equation, which has a single first derivative in time, the wave equation has a second derivative in time. This means we need two initial conditions: the initial displacement $u(x, 0) = f(x)$ and the initial velocity $\frac{\partial u}{\partial t}(x, 0) = g(x)$.
Using the centered second difference in both time and space:
\[ \frac{u_i^{n+1} - 2u_i^n + u_i^{n-1}}{k^2} = c^2 \frac{u_{i-1}^n - 2u_i^n + u_{i+1}^n}{h^2} \]
Solving for the new time level $u_i^{n+1}$:
\[ u_i^{n+1} = 2u_i^n - u_i^{n-1} + \lambda^2\left(u_{i-1}^n - 2u_i^n + u_{i+1}^n\right) \]
where $\lambda = \dfrac{ck}{h}$ is the Courant number.
The formula for $u_i^{n+1}$ requires both $u_i^n$ (the current level) and $u_i^{n-1}$ (the previous level). At the first time step, $n = 1$, we need $u_i^0 = f(x_i)$ and $u_i^{-1}$, which represents the solution before the start of the simulation. We handle this using the initial velocity condition: a first-order approximation gives
\[ u_i^{-1} \approx u_i^0 - k \cdot g(x_i) \]
Substituting into the update formula at the first step produces a workable starting rule. Higher-order starting rules are available and used in more careful implementations.
For the explicit wave equation method, stability requires
\[ \lambda = \frac{ck}{h} \leq 1 \]
This is the Courant-Friedrichs-Lewy condition, usually called the CFL condition after its three discoverers. It has a beautiful physical interpretation: the numerical domain of dependence must contain the physical domain of dependence. The value $ck$ is the distance a wave travels in one time step, and $h$ is the grid spacing. If the wave travels more than one grid spacing per time step, information propagates faster than the grid can track it, and the method becomes unstable.
Unlike the heat equation, which smooths out initial profiles over time, the wave equation preserves sharp features. A triangular initial displacement propagates as two triangular waves traveling in opposite directions. If the wave equation method is working correctly, these features should be maintained without artificial smoothing or oscillation. Numerical diffusion — the artificial spreading of sharp features — is a concern in wave equation computation, and higher-order methods are used in practice to minimize it.
Draw a time sequence of displacement profiles for a string with a triangular initial displacement. Show three snapshots: the initial triangle at $t = 0$; two smaller triangles moving in opposite directions at an intermediate time $t_1$; the two triangles reaching the boundaries and reflecting at a later time $t_2$. Label the wave speed $c$ and indicate the direction of propagation.
The Laplace equation is
\[ \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = 0 \]
This equation governs steady-state problems: temperature in a plate after it has reached thermal equilibrium, electrostatic potential, pressure in an incompressible fluid. Unlike the heat and wave equations, the Laplace equation has no time variable. The solution $u(x,y)$ depends only on position in a two-dimensional domain.
A function satisfying the Laplace equation is called harmonic. Harmonic functions have the mean-value property: the value at any interior point equals the average of the values at all surrounding points. This is the mathematical expression of equilibrium. Each interior point takes the average of its neighbors.
On a uniform grid with spacing $h$ in both the $x$ and $y$ directions, the centered second differences in $x$ and $y$ give
\[ \frac{u_{i-1,j} - 2u_{i,j} + u_{i+1,j}}{h^2} + \frac{u_{i,j-1} - 2u_{i,j} + u_{i,j+1}}{h^2} = 0 \]
Multiplying through by $h^2$ and rearranging:
\[ u_{i,j} = \frac{1}{4}\left(u_{i-1,j} + u_{i+1,j} + u_{i,j-1} + u_{i,j+1}\right) \]
This is the five-point stencil for the discrete Laplace equation. It says: the value at each interior node equals the average of its four neighbors (left, right, below, above). This discrete mean-value property perfectly mirrors the continuous mean-value property.
Draw a small portion of a two-dimensional grid. Highlight the five-point stencil: a center node $(i,j)$ and its four immediate neighbors $(i-1,j)$, $(i+1,j)$, $(i,j-1)$, $(i,j+1)$. Label all five nodes. Draw arrows from the four neighbors pointing toward the center, indicating that the center value is the average of its neighbors.
The Laplace equation requires boundary conditions on all four sides of the rectangular domain. Dirichlet conditions fix the value of $u$ on each side. For example, a metal plate might be held at $u = 100$ on one side and $u = 0$ on the other three sides. The numerical solution finds the steady-state temperature distribution in the interior.
The discrete Laplace equation gives one equation for each interior node: the value at that node equals the average of its four neighbors. If the domain has $M \times N$ interior nodes, the system has $MN$ equations and $MN$ unknowns. The system is large but sparse: each equation involves only five of the $MN$ unknowns. This structure can be exploited by both direct and iterative solvers.
For the Laplace equation on a large grid, the system of $MN$ equations can be solved iteratively. The discrete mean-value property suggests an elegant iterative approach: start with an initial guess for the interior values, then update each node to be the average of its current four neighbors, and repeat until the values converge.
In the Jacobi method, all interior nodes are updated simultaneously using the current values of all neighbors:
\[ u_{i,j}^{(k+1)} = \frac{1}{4}\left(u_{i-1,j}^{(k)} + u_{i+1,j}^{(k)} + u_{i,j-1}^{(k)} + u_{i,j+1}^{(k)}\right) \]
The superscript $(k)$ denotes the iteration number. We start with an initial guess — often all interior values set to zero or to the average of the boundary values — and iterate until
\[ \max_{i,j} \left|u_{i,j}^{(k+1)} - u_{i,j}^{(k)}\right| < \varepsilon \]
for a chosen tolerance $\varepsilon$.
In the Gauss-Seidel method, updated values are used immediately as they become available during the sweep:
\[ u_{i,j}^{(k+1)} = \frac{1}{4}\left(u_{i-1,j}^{(k+1)} + u_{i+1,j}^{(k)} + u_{i,j-1}^{(k+1)} + u_{i,j+1}^{(k)}\right) \]
Because the updated left and lower neighbors are available immediately, Gauss-Seidel typically converges roughly twice as fast as Jacobi for this problem. The method is directly analogous to the Gauss-Seidel iteration for linear systems studied in Chapter 9.
Algorithm: Gauss-Seidel Iteration for the Laplace Equation
Purpose:
Approximate the solution to u_xx + u_yy = 0 on a rectangular domain
with Dirichlet boundary conditions.
Inputs:
N, M — number of grid points in x and y directions
h — uniform grid spacing (same in both directions)
BC values — boundary values on all four sides
tolerance — convergence threshold epsilon
max_iter — maximum number of iterations
Initialize:
Set all interior values u[i][j] = 0 (or boundary average)
Fix boundary values (do not update them)
Iteration loop:
For iteration = 1 to max_iter:
max_change = 0
For i = 1 to N-1:
For j = 1 to M-1:
u_old = u[i][j]
u[i][j] = 0.25 * (u[i-1][j] + u[i+1][j] + u[i][j-1] + u[i][j+1])
max_change = max(max_change, |u[i][j] - u_old|)
If max_change < tolerance:
Stop; convergence achieved.
Output:
The array u[i][j], representing the approximate steady-state solution.
Reliability notes:
Convergence is guaranteed for the Laplace equation with Dirichlet BCs.
Convergence can be slow for large grids; many professional codes use
acceleration techniques (SOR, multigrid) not covered here.One of the distinguishing features of PDE computation is that the solution is a field, not a single number. For a one-dimensional time-dependent problem, the solution is an array of temperatures at every grid point at every time step. For a two-dimensional steady-state problem, the solution is a matrix of values covering the spatial domain.
A space-time plot displays the solution as a surface or contour plot over the space-time grid. The horizontal axis represents position $x$, the vertical axis represents time $t$, and color or height represents the value of $u(x,t)$. For the heat equation, this plot shows the temperature profile smoothing and decaying toward the steady state. For the wave equation, it shows characteristic lines propagating from the initial profile.
A simpler alternative is to plot the solution $u(x, t_n)$ as a function of $x$ at several chosen time levels $t_n$. Each snapshot is a curve showing the spatial profile at one moment. Overlaying several snapshots on the same axes gives a clear picture of how the solution evolves.
For two-dimensional steady-state problems, the solution $u(x,y)$ is naturally displayed as a heatmap or contour plot over the spatial domain. Contour lines connect points of equal solution value, analogous to elevation contours on a topographic map. For the Laplace equation on a plate with one hot side, the contour lines bow toward the cold boundaries, showing how heat spreads through the interior.
For the heat equation example: draw a series of six curves on the same axes, each showing the temperature profile $u(x, t_n)$ as a function of $x$ at a different time $t_n$. The first curve is a sine arch (the initial condition). Each subsequent curve is a flatter, shorter arch, showing the temperature decaying toward zero as heat flows to the boundaries. Label the time level on each curve.
Visualization is not merely decorative — it is a diagnostic tool. A numerical PDE solution should be checked visually against physical expectations: Does the heat equation solution decay smoothly? Does the wave equation solution preserve wave shape? Does the Laplace equation solution satisfy the mean-value property qualitatively? Solutions that look physically wrong almost certainly have an error somewhere in the numerical implementation.
Partial differential equations are the mathematical language of the physical sciences. The applications of numerical PDE methods are vast, and understanding even the basics of this chapter connects students to major areas of scientific computing.
Students working with numerical PDE methods encounter several recurring errors. Recognizing them early prevents wasted computation and misinterpreted results.
This chapter has introduced the major ideas of numerical PDE computation: grids, discretization, finite differences, boundary conditions, stability, the explicit heat equation method, the wave equation method, the Laplace equation, and iterative grid methods. The examples have been deliberately simple — one-dimensional domains, uniform grids, Dirichlet boundary conditions — so that the core ideas are clear.
Chapter 14 brings together the entire numerical methods course in a capstone context. Students will see how to choose among methods, manage error, document computations, visualize results, and communicate numerical findings. The numerical PDE methods of this chapter will appear as one important part of that broader computational toolkit.
Students who continue into advanced numerical analysis will encounter:
The foundation built in this chapter — understanding what a PDE is, why numerical methods are needed, how a grid works, how finite differences discretize derivatives, and what stability means — supports all of these advanced topics.
Chapter 13 introduced numerical PDE methods as an accessible preview of a vast subject. Chapter 14 closes the textbook by bringing together error analysis, algorithms, root-finding, interpolation, approximation, differentiation, integration, linear systems, eigenvalues, optimization, ODEs, and PDEs into a single computational practice. You will learn to choose a numerical method deliberately, estimate and test error, document and reproduce a computation, visualize and communicate results, and complete a capstone simulation project that draws on the entire course.