Lesson 3: Wave or Diffusion Case Study

1 Modelling context

Lessons 1 and 2 developed the FTCS scheme and its stability restriction. This lesson assembles those ingredients into a complete case study. The purpose is not to introduce new theory but to demonstrate the full workflow: problem statement, discretization, implementation, verification, and interpretation.

The case study chosen is a one-dimensional diffusion problem, which is mathematically the same as the heat equation with a different physical interpretation. Diffusion is a good choice because it has a closed-form analytical solution on a simple domain, which permits direct verification, and because it demonstrates the qualitative behaviour that gives the heat equation its name: smoothing of sharp features and decay toward uniformity.

Every decision must be documented

A case study is not just a run that produces output. It is a documented record of every modelling decision: the grid, the stencil, the boundary conditions, the stability check, the verification, and the interpretation. The purpose is to make the result reproducible and the reasoning reviewable.

2 Learning objectives

By the end of this lesson you will be able to

Carry out a complete PDE simulation from problem statement to interpretation.
Document every modelling decision — grid, stencil, boundary conditions, stability check.
Verify a PDE simulation against an analytical solution.
Identify the dominant error source in a simulation.
Interpret a PDE simulation result in terms of the underlying physical or modelling question.

3 Prerequisite review

FTCS scheme: u_i^{n+1} = r u_{i−1}^n + (1 − 2r) u_i^n + r u_{i+1}^n, with r = α Δt / Δx². Developed in Lesson 1.
CFL condition: r ≤ 1/2. Derived in Lesson 2 from von Neumann stability analysis.
Dirichlet and Neumann boundary conditions: Value-specified and flux-specified conditions, respectively. Developed in Unit 3.1, Lesson 3.
Analytical solution for the sine initial condition: For u(x, 0) = sin(π x) with zero Dirichlet boundaries, the exact solution is u(x, t) = e^{−α π² t} sin(π x).

4 Mathematical formulation

Problem statement

A one-dimensional rod of length L = 1 has an initial temperature profile

u(x, 0) = 2 x (1 − x) + 0.5 sin(3 π x),    0 ≤ x ≤ 1

The initial condition is a smooth profile with a broad maximum near the centre of the rod, plus a superimposed oscillation. The diffusivity is α = 0.02, and both ends of the rod are held at zero temperature:

u(0, t) = u(1, t) = 0

Grid and scheme choices

Domain: [0, 1] × [0, 2].
Spatial resolution: N_x = 40, so Δx = 0.025.
Diffusion number: r = 0.4 (below the stability limit 0.5).
Time step: Δt = r Δx² / α = 0.4 · 0.000625 / 0.02 = 0.0125.
Number of time steps: N_t = 2 / 0.0125 = 160.
Scheme: FTCS.
Boundary conditions: homogeneous Dirichlet at both ends.

What the model represents

The rod is insulated on its sides and held at zero temperature at both ends. Heat diffuses from the interior toward the ends and escapes there. Over time, the interior temperature decays toward zero, and the initial oscillation is smoothed by diffusion. The steady state is u(x, ∞) = 0 everywhere.

Figure 3.8. Evolution of the temperature profile under diffusion. The initial profile decays toward zero, and the superimposed oscillation is smoothed. The ends remain at zero throughout.

Expected behaviour

Three qualitative features are expected from the physical interpretation:

Monotone decay. The maximum temperature decreases monotonically with time, since heat is continually lost through the ends.
Increased smoothness. The oscillation in the initial condition is smoothed at a rate proportional to its wavenumber squared. Higher-frequency components decay faster.
Approach to zero. As t → ∞, the solution approaches the steady state u = 0 everywhere.

5 Worked example

Running and interpreting the simulation

Implement the FTCS scheme with the grid and parameters above. Run to t = 2 and record the temperature profile at several times.

Verification against the analytical solution.

The initial condition decomposes into a polynomial part and a sine part. Both have known analytical solutions under the heat equation with homogeneous Dirichlet boundaries:

The sine term 0.5 sin(3π x) decays as 0.5 e^{−α (3π)² t} sin(3π x).
The polynomial term 2x(1 − x) can be expanded in a Fourier sine series, with each mode decaying at its own rate.

Compute the analytical solution at a few grid points and compare against the numerical result:

Comparison at `x = 0.5`
t	Numerical	Analytical	Error
0.0	0.500	0.500	—
0.5	0.397	0.397	< 10⁻³
1.0	0.298	0.299	≈ 10⁻³
1.5	0.211	0.213	≈ 2 × 10⁻³
2.0	0.142	0.145	≈ 3 × 10⁻³

The numerical solution tracks the analytical solution closely, with the error growing slowly over time. This is consistent with the first-order temporal accuracy of the scheme.

Convergence study.

Refine the grid by a factor of 2 in both Δx and Δt (keeping r = 0.4 fixed), and recompute the error at t = 2:

Convergence at `t = 2`, `x = 0.5`
Δx	Δt	Error	Ratio
0.025	0.0125	3.1 × 10⁻³	—
0.0125	0.003125	7.8 × 10⁻⁴	3.97
0.00625	0.00078125	2.0 × 10⁻⁴	3.90

The error ratio is approximately 4, indicating second-order convergence in the overall error. This is expected: although the time discretization is first-order, the spatial discretization is second-order, and when both are refined together the error is dominated by the more slowly converging component (temporal). But because Δt is refined quadratically with Δx in our grid refinement, the temporal and spatial errors remain in balance, producing overall second-order behaviour.

Observed behaviour.

The maximum temperature decreases monotonically: from 0.5 at t = 0 to 0.142 at t = 2.
The oscillation in the initial condition is smoothed rapidly. By t = 0.5, the profile is visibly smoother than at t = 0.
The profile remains symmetric about x = 0.5 to within numerical precision, as expected from the symmetry of the initial condition and boundary conditions.

Cost.

The original grid has 41 × 161 = 6601 grid points, of which 39 × 160 = 6240 are computed by the scheme. The refined grid at Δx/2 has 81 × 321 = 26001 points. The total work increases by a factor of about 4, consistent with the CFL scaling: refining Δx by a factor of 2 quadruples the number of time steps and roughly doubles the number of spatial points, giving a factor of 8 in the naive count. However, since we use an explicit scheme with no further coupling, the actual work is N_x × N_t = 40 × 160 = 6400 vs 80 × 320 = 25600 — a factor of 4.

6 Verification check

Checks to perform

Reproduce the initial condition. At n = 0, the grid function should equal the specified initial condition at every grid point.
Confirm boundary values remain zero at every time step. A common error is to overwrite the boundary with the interior stencil.
Compare against the analytical solution at intermediate times. Compute the maximum absolute error over the grid at several times.
Confirm the observed order. Refine the grid and confirm the error decreases at the expected rate.
Check the symmetry of the solution. The initial condition and boundary conditions are symmetric about x = 0.5. The numerical solution should preserve this symmetry to within floating-point tolerance.
Verify the monotonic decay of the maximum. The maximum of the numerical solution should decrease monotonically with time.
Test the steady-state limit. Run to a large final time and confirm the numerical solution approaches zero everywhere.

7 Interpretation

The simulation is complete, verified, and documented. The final step is to state what the model says — and what it does not say.

What the model says

Under the stated initial and boundary conditions, the temperature in the rod decays toward zero over the two-unit time interval simulated. The decay is driven by diffusion of heat from the interior of the rod toward the ends, where it is removed by the fixed-temperature boundary conditions. The superimposed oscillation in the initial condition is smoothed rapidly, with the higher-frequency component decaying before the lower-frequency polynomial component.

What the model does not say

The model assumes a constant diffusivity α. If the material properties depend on temperature or position, the model does not apply.
The model assumes the ends are held exactly at zero temperature. If the boundary condition is an idealisation of a physical situation (for example, contact with a reservoir of finite heat capacity), the model's predictions near the ends may be inaccurate.
The model assumes a one-dimensional geometry. If the rod has significant radial temperature gradients or heat exchange with the environment along its length, a two- or three-dimensional model is required.
The model does not include any source or sink of heat in the interior. A problem with a heat source would require an additional term in the PDE.

Claim-support statement

Under the stated assumptions of constant diffusivity, one-dimensional geometry, and perfectly conducting boundaries at zero temperature, the model predicts the temperature evolution of the rod over the simulated time interval. The model does not address scenarios in which these assumptions are violated.

8 Uncertainty and limitations

The uncertainty in the simulation result comes from three sources.

Sources of uncertainty in the case study
Source	Magnitude	Reduced by
Truncation error	≈ 3 × 10⁻³ at the given grid	Refining `Δx` and `Δt`
Uncertainty in `α`	Not quantified here	Sensitivity analysis with respect to `α`
Boundary condition idealisation	Not quantified here	Comparison against measurement

The numerical error is small relative to the quantities of interest in this case. The dominant uncertainty is not the numerical error but the uncertainty in the physical parameters and boundary conditions.

Numerical error is not the only error

A simulation can be numerically accurate to several digits and still be substantially wrong, if the physical parameters or the boundary conditions are incorrect. Numerical verification is necessary but not sufficient for a credible model.

9 Common implementation errors

Watch for these

Mismatching the initial condition and the boundary conditions. If u_0(x) does not satisfy the boundary condition at both ends, the numerical solution develops a sharp transient. Verify that the initial and boundary data are compatible.
Using the wrong number of time steps. The final time step should coincide with the target final time. Round-off can produce a slight mismatch; document it or adjust the step size.
Setting the diffusion number above 0.5 by mistake. A change in α, Δx, or Δt changes r. Recompute it before every run.
Comparing against the wrong analytical solution. The analytical solution for the specific initial condition must be used. Using the solution for a different initial condition produces a spurious error.
Confusing the numerical and analytical error with the model error. The comparison against the analytical solution measures the numerical error of the scheme, not the error of the model in representing a physical system.
Presenting results without documenting the grid and parameters. A simulation result is meaningless without the context of the grid, the boundary conditions, and the parameter values.

10 Independent modelling task

Your turn

Choose a one-dimensional PDE problem — either diffusion or wave — with a tractable analytical solution or a clear physical interpretation.
State the problem: the PDE, the domain, the initial condition, and the boundary conditions. Include units where applicable.
Choose the grid and the scheme. Document every parameter: Δx, Δt, r, N_x, N_t. Verify the stability condition.
Implement the scheme. Verify against the analytical solution (if available) or against the physical expectations.
Conduct a convergence study at three grid resolutions. Confirm the observed order matches the theoretical order.
Report the maximum error and the dominant error source at the finest grid.
Write a claim-support statement of the model: what it says under the stated assumptions, and what it does not say.
Prepare a one-page summary suitable for a technical reader, and a one-paragraph summary in plain language.

This exercise completes Laboratory 6 and provides the core of Assignment 3.

11 Reflection prompt

Reflect before continuing

The numerical error in the case study is small relative to the physical quantities of interest. Does this mean the model is accurate?

Two or three sentences. Consider the distinction between numerical accuracy and model validity, and the role of uncertainty in the physical parameters and boundary conditions.

12 Lesson summary

Key takeaways

A complete PDE simulation documents every modelling decision. Grid, scheme, boundary conditions, stability check, verification, and interpretation are all part of the result.
Verification against an analytical solution is the standard test. When an analytical solution exists, it provides direct evidence that the scheme is implemented correctly.
Convergence studies confirm the observed order. Refining the grid and measuring the error reduction validates the scheme's theoretical order.
Numerical error and model error are distinct. A numerically accurate simulation can still be substantially wrong if the model's assumptions are incorrect.
Every simulation result should be accompanied by a claim-support statement. Under what assumptions does the model apply, and what does it not address?

Module 3 complete

You have constructed finite-difference grids, derived the derivative approximations, specified boundary and initial conditions, implemented the FTCS scheme, and analysed its stability. Module 4 turns to Monte Carlo simulation and the treatment of stochastic models.

§ Self-quiz, practice, and homework

Practice exercises

For the case study, compute the exact solution at (0.5, 0.5) using the analytical formula for the sine component. Compare against the numerical value.
List the physical assumptions in the case study. For each, describe a physical situation in which the assumption would fail.
Write a claim-support statement for the case study.
Identify the dominant error source in the case study at the finest grid. Justify.

Homework problems

Full case study report. Choose a PDE problem — heat, diffusion, wave, or a combination. Produce a complete case study report: problem statement, grid, scheme, boundary and initial conditions, stability check, verification against an analytical or reference solution, convergence study, cost analysis, and a claim-support statement. 1500–2000 words.
Sensitivity analysis. For the case study, vary the diffusivity α by ±20% and record the effect on the solution at a fixed point and time. Plot the range of outcomes. This is the physical parameter sensitivity, distinct from numerical error.
Wave equation. Repeat the case study for the wave equation ∂²u/∂t² = c² ∂²u/∂x² with a standard explicit scheme. Derive its stability condition (which differs from the heat equation) and verify numerically.

Mini project

Produce a complete PDE simulation report for a problem of your choice, with all the required components. The report should be suitable for a technical reader with a mathematical background, and should include a one-paragraph plain-language summary. Submit the spreadsheet, the analytical or reference solution where available, and the report. This completes Laboratory 6 and is the basis for Assignment 3.