Mathematical foundations of Crank-Nicolson Heat PDE.
The mathematical modeling of diffusive transport—from transient thermal conduction in aerospace heat shields to chemical solute dispersion and financial option pricing via the Black-Scholes PDE—is governed by parabolic partial differential equations. Discretizing the heat equation ∂u/∂t = α ∂²u/∂x² with explicit forward differences leads to severe numerical instability: the time step must satisfy the crippling Courant-type condition r = α Δt / Δx² ≤ 1/2, forcing billions of tiny time steps on fine meshes. In 1947, British mathematicians John Crank and Phyllis Nicolson devised an elegant implicit formulation that averages spatial derivatives at time levels n and n + 1. By centering the stencil at t_{n+1/2}, Crank-Nicolson achieves unconditional von Neumann stability (|G| ≤ 1 for all r) and second-order accuracy O(Δt² + Δx²). In this laboratory, we implement the Crank-Nicolson scheme, invert the tridiagonal system via the Thomas algorithm, and verify the second-order O(h²) error convergence.
Viva voce questions & analytical theory
Why is the explicit Forward Time Central Space (FTCS) scheme problematic for the heat equation?
FTCS has local error O(Δt + Δx²), but von Neumann stability analysis reveals that the amplification factor satisfies |G| ≤ 1 only when the diffusion number r = α Δt / Δx² ≤ 1/2. Halving Δx requires a 4-fold reduction in Δt, leading to astronomical numbers of tiny time steps on fine spatial meshes.
How did John Crank and Phyllis Nicolson achieve unconditional stability in 1947?
Crank and Nicolson averaged the central difference spatial stencils at the current time level t_n and the future time level t_{n+1} (a trapezoidal rule in time). This symmetric time-centering cancels the leading first-order temporal error, achieving second-order accuracy O(Δt² + Δx²) while guaranteeing |G| ≤ 1 for any choice of r = Δt / Δx².
How does von Neumann stability analysis prove Crank-Nicolson is unconditionally stable?
Substituting the discrete Fourier mode u_j^n = G^n e^{i k j Δx} into the Crank-Nicolson stencil yields G = [1 - 2r sin²(θ/2)] / [1 + 2r sin²(θ/2)], where θ = k Δx. Since the denominator is always greater than or equal to the absolute value of the numerator for any r > 0, |G| ≤ 1 holds unconditionally for all wavenumbers k.
What is the computational cost of inverting the Crank-Nicolson matrix at each time step?
In one spatial dimension, the implicit matrix is tridiagonal and strictly diagonally dominant. The Thomas algorithm inverts it in exactly O(N) operations per time step, meaning the implicit step is essentially as fast as an explicit step while permitting time steps orders of magnitude larger.
Can the Crank-Nicolson scheme produce spurious numerical oscillations even though it is stable?
Yes. When r = Δt / Δx² ≫ 1, the amplification factor for the highest spatial frequency (θ = π) approaches G ≈ -1. While bounded (|G| ≤ 1), the negative sign causes high-frequency components to oscillate sign from step to step instead of decaying monotonically, requiring spatial smoothing or L-stable implicit methods (BDF2 or backward Euler).
How is Crank-Nicolson extended to two spatial dimensions (2D heat equation)?
In 2D, the direct implicit matrix is no longer tridiagonal but pentadiagonal with large bandwidth, making direct inversion O(N³) expensive. The Alternating Direction Implicit (ADI) method of Peaceman and Rachford splits each time step into two half-steps, inverting tridiagonal systems along x then y, restoring O(N) efficiency.
What is the physical interpretation of the heat equation as a dissipative system?
The heat equation represents diffusive thermal equilibration driven by Fourier's law of conduction. As t → ∞, all high-frequency spatial temperature gradients decay exponentially as e^{-k² t}, smoothing any initial discontinuity and dissipating total thermodynamic entropy into uniform equilibrium.
Algorithm precautions and stability rules
- Verify that the tridiagonal matrix diagonal 1 + r strictly exceeds off-diagonal sums to preserve strict diagonal dominance in the Thomas algorithm.
- Incorporate Dirichlet boundary values into the right-hand vector to prevent artificial boundary reflection errors.
- Be cautious of very large diffusion numbers r > 10 when initial data has sharp discontinuities, as G ≈ -1 can induce spurious oscillations.
- Ensure floating point precision remains 64-bit IEEE 754 to prevent roundoff accumulation across hundreds of implicit time steps.
- Verify that spatial node count N and time step count M are updated consistently to maintain the chosen scaling relation r = Δt / Δx².
- For non-linear thermal conductivity α(u), solve the resulting non-linear tridiagonal system using Picard or Newton-Raphson iteration.
Theoretical limitations
- While Crank-Nicolson is A-stable, it is not L-stable: it does not damp stiff high-frequency error components as rapidly as backward Euler or BDF2.
- In 2D and 3D, direct Crank-Nicolson produces banded matrices with large bandwidth, requiring Alternating Direction Implicit (ADI) split solvers.
- Crank-Nicolson cannot handle hyperbolic wave propagation without introducing severe dispersion and phase errors.