Numerical solutions of partial differential equations

Start with four spatial values

At 0 s, a supplied temperature-excess profile has values 0 K, 2 K, 2 K and 0 K at positions 0 m, 1 m, 2 m and 3 m. The two endpoint values stay at 0 K. What does a numerical diffusion update give at 1 s? Temperature excess means temperature minus a fixed 280 K reference: a zero excess here corresponds to 280 K, not to absolute zero.

Use a spatial spacing of 1 m, a time step of 1 s and a supplied diffusivity of 0.25 m²/s. The update below gives each old neighbour weight 0.25 and the old center weight 0.5. At 1 m it gives (0.25)(0 K) + (0.5)(2 K) + (0.25)(2 K) = 1.5 K. The value at 2 m is also 1.5 K. These are numerical values for this chosen update; we will check what they establish about the continuous problem.

To read the underlying equation, you need partial derivatives and the distinction between initial and boundary data. The published lessons Derivatives, Taylor series, and local approximation and Partial differential equations and Green functions provide that preparation. You can first follow the finite calculations and graph, then return to the derivative-based error bound. The coupled step later can be read as two simultaneous equations; Matrices, linear algebra, and eigenvalue problems provides further preparation.

State the continuous problem

Let x be position along a fixed straight axis in metres, t elapsed time in seconds, and θ(x,t) temperature excess in kelvin. On 0 < x < 3 m and 0 < t < 6 s, supply the constant-coefficient diffusion equation θₜ = κθₓₓ, with κ = 0.25 m²/s. The subscripts mean a partial time derivative and a second partial spatial derivative. Both sides have K/s units: κ multiplies a K/m² derivative. This is a stipulated heat-equation model, without an internal source term.

Prescribe θ(0,t) = θ(3 m,t) = 0 K for 0 ≤ t ≤ 6 s. These are Dirichlet conditions: they supply field values at fixed spatial endpoints. They do not describe insulated endpoints, which would instead constrain a spatial derivative or flux. Changing the boundary condition changes the problem and may change the numerical boundary rule.

For the whole initial profile, define k = π/(3 m) and prescribe θ(x,0) = ((4/√3) K) sin(kx) on 0 ≤ x ≤ 3 m. The angle kx is dimensionless. The profile vanishes at both endpoints and gives 2 K at x = 1 m and x = 2 m because sin(π/3) = sin(2π/3) = √3/2. A finite list of values alone would not specify this entire initial profile. If the sine notation is unfamiliar, use those four supplied starting values for the grid calculation and return to the exact check in the final section.

Update the grid from one old time

Choose a positive spatial spacing h that divides 3 m into M equal intervals, where M is an integer with M ≥ 2. The nodes are xⱼ = jh for j = 0, …, M. Choose a positive time step Δt with NΔt = 6 s for a positive integer N, and grid times tₙ = nΔt for n = 0, …, N. Write Θⱼⁿ for the numerical excess at node j and step n. The superscript counts steps; it is not a power of Θ.

Approximate the time derivative by a forward difference and the second spatial derivative by a centered difference, both using the old time level. With dimensionless μ = κΔt/h², the explicit rule is Θ_j^(n+1) = μΘ_(j−1)^n + (1−2μ)Θ_j^n + μΘ_(j+1)^n for interior nodes j = 1, …, M − 1. Keep the two endpoint values at their prescribed zero excess. The name forward time, centered space describes this stencil: one old center and two old spatial neighbours determine one new value.

For h = 1 m and Δt = 1 s, M = 3 and μ = 0.25. Starting from [0 K, 2 K, 2 K, 0 K], the first two updates give [0 K, 1.5 K, 1.5 K, 0 K] and [0 K, 1.125 K, 1.125 K, 0 K]. Store all new interior values separately before replacing the old row. Reusing the newly computed left neighbour while updating the right node would be a different algorithm.

The graph shows only these four spatial nodes at three times. Different point shapes distinguish the time levels; the shared zero endpoint values are marked once. No curve is drawn between the points, because the update itself supplies no continuous spatial profile there. The table gives every plotted value.

Check a sufficient step restriction

If 0 ≤ μ ≤ 1/2, the three weights μ, 1 − 2μ and μ are nonnegative and sum to one. Each new value is therefore a weighted average of the three old values. With the fixed zero endpoints, nonnegative starting grid values stay nonnegative and cannot exceed the old grid maximum. This is a useful sufficient condition for this uniform one-dimensional explicit diffusion stencil.

The same argument bounds differences between two computations with the same endpoint data. If every old interior difference has absolute size at most E, each new difference has size at most E. The zero endpoint differences complete the stencil. Thus repeated steps do not amplify an initial grid perturbation in the maximum norm, the largest absolute node difference. This statement assumes exact arithmetic and compares this recurrence with itself; it does not yet bound its discrepancy from the PDE.

For κ > 0, the sufficient restriction is Δt ≤ h²/(2κ). Here h = 1 m permits Δt ≤ 2 s. Halving h quarters that permitted time step: at h = 0.5 m the bound is 0.5 s. The relation is specific to this diffusion update. It is not a universal step rule for waves, transport, other dimensions or every PDE solver, and it does not say that every larger μ is unstable on every fixed finite grid.

Find a hidden amplifying mode

Keep h = 1 m but choose Δt = 3 s, so μ = 0.75. There are still only two interior nodes, and the calculation reaches 6 s after two steps. The clean equal interior pair [2 K, 2 K] is multiplied by 1 − μ = 0.25 per step, giving [0.5 K, 0.5 K] and [0.125 K, 0.125 K]. A decreasing clean sequence can conceal an amplified difference.

Compare a second run whose initial interior values differ by [+ε, −ε], with ε = 0.1 K; all endpoint values are identical. The two-node explicit rule multiplies this alternating pair by 1 − 3μ = −1.25 per step. Its differences become [−0.125 K, +0.125 K] at 3 s and [+0.15625 K, −0.15625 K] at 6 s. The signs alternate while the magnitudes grow.

This is an exact amplified mode of the stated discrete recurrence, not a measured fluctuation or the exact PDE error. Continuing its algebraic recurrence would keep multiplying this difference by −1.25; the declared model calculation here stops at 6 s. The example shows why a plausible-looking clean output does not establish numerical stability. A negative temperature excess would mean a value below the reference, not automatically a negative absolute temperature.

Separate consistency from accumulated error

For a local check, start an explicit step from exact PDE values at its three old nodes. Define dⱼⁿ as the exact excess at the next grid time minus that one-step result. It has K units. The rate defect dⱼⁿ/Δt has K/s units; sources use different truncation-error conventions, so keep the formula and units explicit. This local check removes inherited numerical errors.

Suppose the exact solution has continuous second time and fourth spatial derivatives throughout the required time steps and spatial stencils, with uniform bounds |θₜₜ| ≤ Mₜ and |θₓₓₓₓ| ≤ Mₓ. Taylor remainder bounds give |d_j^n| ≤ (Δt²/2)M_t + (κΔt h²/12)M_x. Mₜ uses K/s² and Mₓ uses K/m⁴. Dividing by Δt gives the rate bound B = (Δt/2)Mₜ + (κh²/12)Mₓ. These are sufficient smoothness assumptions, not quantities inferred from a small value table.

When the convex-weight condition holds, combine that local bound with the perturbation argument. If Eₙ is the largest absolute exact-minus-numerical interior error, the boundary data are exact and the same B bounds every step, then Eₙ₊₁ ≤ Eₙ + Δt B. Over N steps, E_N ≤ E₀ + NΔt B. In particular, at a fixed elapsed time T = NΔt this is E₀ + TB. It excludes additional rounding, boundary, parameter and model errors unless they are included separately.

With derivative bounds that remain valid as the mesh is refined, this estimate decreases with Δt and h², provided the step restriction is retained and the initial grid error also decreases. That is a conditional route to convergence, meaning approach to the exact PDE values at the grid points. Agreement between two grids alone is not this certificate. On a first reading, keep local error and accumulated error distinct; return to the Taylor argument after the named derivative lesson.

Solve a coupled implicit step

An implicit alternative evaluates the centered spatial difference at the new time level. Its interior equations are (1+2μ)Θ_j^(n+1) − μΘ_(j−1)^(n+1) − μΘ_(j+1)^(n+1) = Θ_j^n. Several unknown next values are now coupled. Prescribed endpoint values enter these equations at the new time; they are not additional free unknowns.

For our two interior nodes and μ = 0.75, call their next values a and b. From the initial pair [2 K, 2 K], solve 2.5a − 0.75b = 2 K and −0.75a + 2.5b = 2 K. Subtracting gives 3.25(a − b) = 0 K, so a = b. Each equation then gives 1.75a = 2 K and a = b = 8/7 K at 3 s. The coefficient determinant is 2.5² − 0.75² = 91/16, nonzero, so this pair is unique.

For this two-node problem, an equal perturbation pair is multiplied by 1/(1 + μ), and an alternating pair by 1/(1 + 3μ) in an implicit step. Both factors lie between zero and one for μ > 0; at μ = 0.75 they are 4/7 and 4/13. This removes the particular amplification found in the explicit counterexample. It does not establish accuracy, solve every nonlinear implicit equation or make an arbitrary PDE method physically faithful.

Compare with an exact check and retain the model

For our supplied continuous problem, there is a useful exact check: θ(x,t) = ((4/√3) K) sin(kx) exp(−γt), where k = π/(3 m) and γ = κk² = π²/(36 s). Differentiating gives θₜ = −γθ and θₓₓ = −k²θ, so θₜ = κθₓₓ. At t = 0 it retains the full initial profile, and the two endpoint values remain zero. The exponent γt is dimensionless. The published Complex numbers and complex exponentials lesson supplies preparation for the exponential notation if needed.

At either interior node, this exact check gives 2 exp(−π²/36) K at 1 s and 2 exp(−π²/18) K at 2 s. The explicit grid gives 1.5 K and 1.125 K instead. Likewise the implicit value 8/7 K at 3 s differs from the exact 2 exp(−π²/12) K. All these calculations can follow their formulas correctly while answering different questions: an exact PDE solution and two discrete approximations.

Record the equation, diffusivity, domain, full initial profile, boundary data, spacings, update rule and arithmetic precision. Check boundary enforcement, simultaneous old-level updates, the relevant stability argument and errors at matching positions and times. Validate an actual thermal application separately: neither this supplied coefficient nor the exact sine check verifies a material, apparatus or boundary-control assumption. Other PDEs and geometries require appropriate data, stencils and error analysis.

Worked examples

Try it yourself

1. Starting from the displayed 2 s row, take one more explicit step with h = 1 m and Δt = 1 s. Find both interior values and the endpoint values at 3 s.

   Show the answer and reasoning

   Both interior values are 0.84375 K; both endpoint values remain 0 K.

   Each interior has one zero endpoint neighbour and the other interior value 1.125 K. Thus (0.5 + 0.25)(1.125 K) = 27/32 K = 0.84375 K. Use the complete old row, retain the prescribed endpoints and remain within the 0–6 s interval.

2. With the same κ, halve h to 0.5 m but first leave Δt at 1 s. Find μ. What is the largest step allowed by the convex-weight bound, and which step retains μ = 0.25?

   Show the answer and reasoning

   μ = 1, so the convex-weight bound fails. Its largest permitted step is 0.5 s; Δt = 0.25 s retains μ = 0.25.

   The new h² is 0.25 m². Hence μ = (0.25 m²/s)(1 s)/(0.25 m²) = 1, and h²/(2κ) = 0.5 s. Quartering the original time step as well as h² preserves the original ratio. A finer spatial grid alone does not preserve the sufficient stability condition.

3. For the original two interior nodes with μ = 0.25, start two runs with difference [+0.1 K, −0.1 K] and identical endpoints. Find the interior difference after two steps. Does that establish the exact PDE error?

   Show the answer and reasoning

   The difference is [+0.00625 K, −0.00625 K]. It does not establish the PDE error.

   The alternating multiplier is 1 − 3μ = 1/4. Two steps multiply the initial pair by 1/16, giving ±1/160 K. This compares two discrete runs; their shared discretization error relative to the continuous equation is a separate quantity.

4. For a separate smooth solution of the same diffusion equation, suppose |θₜₜ| ≤ 0.08 K/s² and |θₓₓₓₓ| ≤ 0.48 K/m⁴ over every required stencil and step. Use h = 1 m, Δt = 1 s, exact boundaries and E₀ ≤ 0.01 K in exact arithmetic. Bound the rate defect and the accumulated error after two steps.

   Show the answer and reasoning

   B ≤ 0.05 K/s and E₂ ≤ 0.11 K, under all the supplied hypotheses.

   The time term is (1 s)(0.08 K/s²)/2 = 0.04 K/s. The spatial term is (0.25 m²/s)(1 m)²(0.48 K/m⁴)/12 = 0.01 K/s. Thus B ≤ 0.05 K/s. μ = 0.25 retains the convex bound, and E₂ ≤ 0.01 K + (2 s)(0.05 K/s) = 0.11 K. These derivative bounds are additional assumptions for this exercise, not values inferred from the main table or asserted for its sine profile.

5. As a separate grid-data case, take one implicit step with μ = 0.75 from interior values [2 K, 1 K], keeping both endpoints zero. Solve 2.5a − 0.75b = 2 K and −0.75a + 2.5b = 1 K. Is assigning a from its equation while leaving b unknown a completed step?

   Show the answer and reasoning

   a = 92/91 K and b = 64/91 K. Both equations must be satisfied together.

   Multiply both equations by 4 to get 10a − 3b = 8 K and −3a + 10b = 4 K. Their determinant is 91. Elimination gives 91a = 92 K and 91b = 64 K. Substitution verifies both right sides; neither new value can be determined independently by an old-neighbour explicit rule. This separate pair changes the initial grid data.

6. For the final exact check θ(x,t), write θₜ and θₓₓ in terms of θ, verify the equation and data, and give θ(1 m, 2 s). Would adding a constant 1 K retain the complete original problem?

   Show the answer and reasoning

   θₜ = −γθ and θₓₓ = −k²θ; γ = κk² verifies the PDE. The original data are retained and θ(1 m, 2 s) = 2 exp(−π²/18) K. Adding 1 K would change the prescribed data.

   The exponential differentiates to −γ times itself and the sine's second x derivative is −k² times itself. At t = 0 the exponential is one; at the two spatial endpoints the sine is zero. At 1 m the sine factor gives 2 K before decay. A constant addition leaves these PDE derivatives unchanged but raises every initial and boundary value by 1 K, so it fails the complete original problem.

What to carry forward

• A numerical PDE problem keeps the equation, field meaning, domain, initial profile and boundary data together.
• The explicit stencil uses one complete old time level; discrete values alone do not define a continuous profile.
• Nonnegative weights give a useful sufficient stability bound for this particular diffusion update.
• A decreasing clean sequence can conceal an amplified alternating perturbation.
• Consistency starts from exact old values; accumulated error also includes inherited discrepancies and requires its conditions.
• A coupled implicit solve can damp a test mode without certifying accuracy or validating a physical model.