43.11.02 · numerical-analysis / 11-finite-difference-pdes

The method of lines and stability for parabolic problems

shipped3 tiersLean: none

Anchor (Master): LeVeque 2007 *Finite Difference Methods for Ordinary and Partial Differential Equations* (SIAM) Ch. 9-10 (the method of lines, the stiff semidiscrete system, the $\theta$-method family, accuracy and absolute-stability analysis, and the $O(k+h^2)$/$O(k^2+h^2)$ convergence orders); Thomée 2006 *Galerkin Finite Element Methods for Parabolic Problems* 2e (Springer) Ch. 1, 7-9 (the semigroup/semidiscrete framework, the backward Euler and Crank-Nicolson time discretisations, and the smoothing/stability estimates); Hairer-Wanner 1996 *Solving Ordinary Differential Equations II* 2e (Springer) §IV.2-IV.3 (stiffness, A- and L-stability, applied to the parabolic semidiscretisation)

Intuition Beginner

Heat spreads. Drop a hot spot into a cold metal bar and the warmth flows toward the cold, smoothing every bump until the whole bar settles. The continuous rule that governs this smoothing is the heat equation from 02.13.03, and like every partial differential equation it asks for a value at every point of the bar and at every instant of time. A computer cannot store infinitely many points, so we make the problem finite in two separate stages, one for space and one for time.

The first stage handles space. Lay a row of grid points along the bar and keep only the temperatures there. The bending of the temperature profile at a point — how far it sits below the average of its left and right neighbours — is exactly what the spatial part of the elliptic unit 43.11.01 already taught us to measure. Replacing the continuous spatial bending by this neighbour comparison turns the single partial differential equation into a large bundle of ordinary differential equations: one equation for how each grid temperature changes in time. That step is called the method of lines, because each grid point traces a line through time.

The second stage handles time. We now have an ordinary-differential-equation system, and the earlier numerical-ODE units give us integrators for exactly that. The simplest one marches forward: the new temperatures are the old ones plus a small time step times the current rate of change. This explicit march is cheap, but it has a catch. If the time step is too large relative to the grid spacing, the computed temperatures do not smooth — they oscillate and explode, the opposite of how heat behaves.

That catch is the heart of this unit. The fine spatial grid makes the ODE system stiff, in the sense of 43.10.05: fast modes that the true solution barely notices still dictate how big a step the explicit method may take. The cure is to integrate in time with a method built for stiffness, one that stays calm at any step size. Those are the implicit schemes, and they cost one tridiagonal solve per step in exchange for removing the restriction entirely.

Visual Beginner

The picture has two axes: space across the bar and time upward. Space is already chopped into grid points by the method of lines; the question is how big a time step the marching scheme may take before it misbehaves.

The table below contrasts the three standard time-marching choices. The explicit one is cheap per step but capped: the time step must stay below half the squared grid spacing. The two implicit ones cost a small solve each step but accept any time step the accuracy allows.

 time
  ^      explicit (forward) step:                implicit step:
  |      uses only known values below            ties new values together,
  |        U^{n+1}  <- U^n                         solve a tridiagonal system
  |          o                                       o---o---o   (all unknown)
  |         /|\                                       |   |   |
  |        o o o   (known)                          o   o   o   (known)
  |
  +-------------------------------------------------------------> space

scheme	per-step cost	time-step rule	order in time
forward Euler (explicit, FTCS)	one multiply	step $\leq$ half the squared spacing	first
backward Euler (implicit)	one tridiagonal solve	any step	first
Crank-Nicolson (implicit)	one tridiagonal solve	any step	second

The takeaway: the method of lines splits the job into a space step and a time step, and the whole stability story lives in the time step. Explicit marching is fast but leashed to a tiny step on a fine grid; implicit marching pays a small solve to run off the leash.

Worked example Beginner

Let us watch the explicit scheme stay calm and then blow up, on a tiny grid, so the stability limit becomes concrete. Take a bar of length $1$ with the two ends held at temperature $0$ . Put down grid points at spacing $h = 1/4$ , so the interior has three unknown temperatures, called $U_{1}, U_{2}, U_{3}$ from left to right.

The spatial rule says each grid temperature changes at a rate equal to its left neighbour plus its right neighbour minus twice itself, all divided by the squared spacing $h^{2} = 1/16$ . The explicit march updates each temperature by adding the time step $k$ times that rate. Collecting the factor, the new middle temperature is $U_{2} + r (U_{1} - 2 U_{2} + U_{3})$ , where $r = k / h^{2}$ is the single number that controls everything. The stability limit, derived in the formal part, is $r \leq 1/2$ .

Step 1. Start from a simple spike: $U_{1} = 0$ , $U_{2} = 1$ , $U_{3} = 0$ , with the ends fixed at $0$ .

Step 2. Take a safe step, $r = 1/2$ (so $k = h^{2} /2 = 1/32$ ). The new middle value is $1 + \frac{1}{2} (0 - 2 + 0) = 1 - 1 = 0$ . The new left value is $0 + \frac{1}{2} (0 - 0 + 1) = 1/2$ , and by symmetry the new right value is $1/2$ . The spike of height $1$ has spread into two half-height shoulders. The values stay between $0$ and $1$ — smoothing, as heat should.

Step 3. Take an unsafe step, $r = 1$ (so $k = h^{2} = 1/16$ , twice the limit). The new middle value is $1 + 1 (0 - 2 + 0) = 1 - 2 = - 1$ . The new shoulders are $0 + 1 (0 - 0 + 1) = 1$ . The profile is now $1, - 1, 1$ : it did not smooth, it flipped sign and kept its size.

Step 4. Repeat the unsafe step once more from $1, - 1, 1$ . The new middle is $- 1 + 1 (1 + 1 - 2 (- 1)) = - 1 + 4 = 3$ . The middle value has grown from $1$ to $- 1$ to $3$ , swinging wider each step. The computed temperatures explode while real heat would have faded.

What this tells us: the same explicit scheme is trustworthy at $r = 1/2$ and ruinous at $r = 1$ . The dividing line is the parabolic step restriction $k \leq h^{2} /2$ . Because the limit involves the square of the spacing, halving the grid forces the time step down by a factor of four — the practical pain that pushes people toward the implicit schemes.

Check your understanding Beginner

Exercise (easy, multiple choice).

What does an implicit time-stepping scheme such as backward Euler buy you, compared with the explicit march, for the heat equation?

A. A larger error at the same step size. B. Freedom from the time-step cap, at the cost of one small solve per step. C. No need for a spatial grid. D. Exact answers with no error at all.

Hint

The implicit schemes are the ones built for stiffness; they stay calm at any step size but each step ties the new values together.

Answer

B. Feedback-correct: correct; backward Euler and Crank-Nicolson are stable at any time step, paying one tridiagonal solve per step to escape the $k \leq h^{2} /2$ cap. Feedback-wrong: they still have spatial error and still need a grid (ruling out C and D), and their stability is the gain, not a larger error (ruling out A).

Formal definition Intermediate+

Consider the one-dimensional heat (diffusion) equation 02.13.03 on the interval $(0, 1)$ with homogeneous Dirichlet data, $$ u_t = u_{xx} \ \text{ in } (0,1)\times(0,T], \qquad u(0,t) = u(1,t) = 0, \qquad u(x,0) = u_0(x), $$ whose classical solvability and smoothing are supplied by the continuous theory and are not reproved here. Fix an integer $m \geq 1$ , set the spatial mesh width $h = 1/ (m + 1)$ and grid points $x_{j} = j h$ for $0 \leq j \leq m + 1$ , and let $U_{j} (t) \approx u (x_{j}, t)$ for the interior indices $1 \leq j \leq m$ , with $U_{0} = U_{m + 1} = 0$ .

Semidiscretisation: the method of lines. Replace the spatial derivative $u_{xx}$ by the second-difference operator $D^{2}$ of 43.11.01 at each interior node, leaving time continuous: $$ \frac{dU_j}{dt} = \frac{U_{j-1}(t) - 2U_j(t) + U_{j+1}(t)}{h^2}, \qquad 1 \le j \le m. $$ Collecting the interior values into a vector $U (t) \in R^{m}$ , this is the linear, constant-coefficient, stiff ODE system $$ U'(t) = A_h,U(t), \qquad A_h = \frac{1}{h^2},\mathrm{tridiag}(1,-2,1) \in \mathbb{R}^{m\times m}, $$ the semidiscrete system. Here $A_{h}$ is the negative of the symmetric positive-definite operator $T = \frac{1}{h ^{2}} tridiag (- 1, 2, - 1)$ of 43.11.01; it is symmetric negative definite, with eigenvalues $$ \mu_p = -\frac{4}{h^2}\sin^2!\Big(\frac{p\pi h}{2}\Big), \qquad 1 \le p \le m, $$ spreading from $μ_{1} \to - π^{2}$ (the slow mode) down to $μ_{m} \to - 4/ h^{2}$ (the fast mode) as $h \to 0$ . The huge spread $∣ μ_{m} ∣/∣ μ_{1} ∣ = O (h^{- 2})$ is the stiffness of the system in the sense of 43.10.05.

Full discretisation: the $θ$ -method family. Apply a one-step time integrator with step $k$ and time levels $t_{n} = nk$ , writing $U^{n} \approx U (t_{n})$ . The $θ$ -method takes the spatial operator at a convex combination of the new and old levels, $$ \frac{U^{n+1} - U^n}{k} = A_h\big((1-\theta)U^n + \theta,U^{n+1}\big), \qquad \theta \in [0,1]. $$ Three choices are standard. With $θ = 0$ this is forward Euler, the forward-time centred-space (FTCS) scheme $$ U^{n+1} = (I + kA_h),U^n, $$ explicit: each new value is an affine combination of known values. With $θ = 1$ it is backward Euler, $$ (I - kA_h),U^{n+1} = U^n, $$ implicit: one tridiagonal solve per step. With $θ = \frac{1}{2}$ it is Crank-Nicolson (the trapezoidal rule applied to the semidiscrete system), $$ \Big(I - \tfrac{k}{2}A_h\Big)U^{n+1} = \Big(I + \tfrac{k}{2}A_h\Big)U^n, $$ also one tridiagonal solve per step. In every case $U^{n + 1} = r (k A_{h}) U^{n}$ for the method's stability function $r$ of 43.10.04: $r (z) = 1 + z$ (forward Euler), $r (z) = 1/ (1 - z)$ (backward Euler), $r (z) = (1 + z /2) / (1 - z /2)$ (Crank-Nicolson), evaluated at the matrix argument $k A_{h}$ .

Local truncation error. The local truncation error $τ^{n}$ is the residual obtained by inserting the exact solution into the scheme and dividing by $k$ . For FTCS, $$ \tau^n_j = \frac{u(x_j,t_{n+1}) - u(x_j,t_n)}{k} - \frac{u(x_{j-1},t_n) - 2u(x_j,t_n) + u(x_{j+1},t_n)}{h^2}. $$ A scheme is consistent of order $(p, q)$ if $∥ τ ∥_{\infty} = O (k^{p} + h^{q})$ . The discrete maximum norm is $∥ V ∥_{\infty} = max_{j} ∣ V_{j} ∣$ , the spatial mesh ratio is $r = k / h^{2}$ , and the symbols $A_{h}$ , $D^{2}$ , $U^{n}$ , $r (z)$ , $μ_{p}$ , $θ$ , $τ^{n}$ , and $O (k + h^{2})$ are recorded in _meta/NOTATION.md.

Counterexamples to common slips

Semidiscretising is not yet a numerical method. The method of lines produces an exact ODE system $U^{'} = A_{h} U$ ; the approximation in time only enters when an integrator is chosen. Conflating the two hides where each error term — $O (h^{2})$ from space, $O (k)$ or $O (k^{2})$ from time — actually comes from.
The restriction is $k \leq h^{2} /2$ , not $k \leq h /2$ . The fast eigenvalue is $μ_{m} \approx - 4/ h^{2}$ , and forward Euler needs $k ∣ μ_{m} ∣ \leq 2$ , giving $k \leq h^{2} /2$ . Reading the cap as linear in $h$ underestimates the cost of refinement by a full power of $h$ .
Unconditional stability is not unconditional accuracy. Backward Euler and Crank-Nicolson place no cap on $k$ for stability, but a large $k$ still destroys accuracy; the schemes merely refuse to blow up, they do not promise to resolve the solution.
Crank-Nicolson is A-stable but not L-stable. Its stability function tends to $- 1$ as $z \to - \infty$ along the real axis, so the stiffest unresolved spatial modes are damped only marginally and can ring; backward Euler, with $r (\infty) = 0$ , annihilates them. This is the 43.10.05 L-stability distinction surfacing in the PDE.

Key theorem with proof Intermediate+

The signature result is the parabolic stability restriction: the explicit FTCS scheme for the heat equation is stable in the maximum and $ℓ^{2}$ norms exactly when the mesh ratio satisfies $r = k / h^{2} \leq 1/2$ , while the backward-Euler and Crank-Nicolson schemes are stable for every $k$ . The mechanism is the 43.10.04 eigen-decoupling: $A_{h}$ is symmetric, so the one-step operator $r (k A_{h})$ is diagonalised by the same sine modes, and stability reduces to forcing each scaled eigenvalue $k μ_{p}$ into the method's absolute-stability region ^{[LeVeque, R. J. — Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems; Hairer, E. & Wanner, G. — Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems (2nd ed.)]}.

Theorem (parabolic CFL and unconditional implicit stability). Let $A_{h} = \frac{1}{h ^{2}} tridiag (1, - 2, 1)$ with eigenvalues $μ_{p} = - \frac{4}{h ^{2}} sin^{2} (p π h /2)$ , and let the time-stepping operator be $U^{n + 1} = r (k A_{h}) U^{n}$ . Then, in the $ℓ^{2}$ norm:

(i) The FTCS scheme $r (z) = 1 + z$ satisfies $∥ r (k A_{h}) ∥_{2} \leq 1$ if and only if $k \leq h^{2} /2$ .

(ii) The backward-Euler scheme $r (z) = 1/ (1 - z)$ and the Crank-Nicolson scheme $r (z) = (1 + z /2) / (1 - z /2)$ satisfy $∥ r (k A_{h}) ∥_{2} \leq 1$ for every $k > 0$ .

Proof. Since $A_{h}$ is real symmetric, the spectral theorem gives an orthonormal eigenbasis $v_{j}^{(p)} = 2 h sin (p π x_{j})$ with $A_{h} v^{(p)} = μ_{p} v^{(p)}$ (the discrete sine modes of 43.11.01). For any rational $r$ analytic on the spectrum, $r (k A_{h})$ has the same eigenvectors with eigenvalues $r (k μ_{p})$ , and because the eigenbasis is orthonormal the spectral norm is the largest eigenvalue modulus: $$ |r(kA_h)|2 = \max{1\le p\le m} |r(k\mu_p)|. $$ Thus stability $∥ r (k A_{h}) ∥_{2} \leq 1$ holds exactly when $k μ_{p}$ lies in the absolute-stability region $S = {z : ∣ r (z) ∣ \leq 1}$ for every $p$ .

(i) Forward Euler. Here $r (z) = 1 + z$ and $S = [- 2, 0]$ on the real axis 43.10.04. The eigenvalues $k μ_{p} = - 4 r sin^{2} (p π h /2)$ are real and negative, where $r = k / h^{2}$ ; the most negative is at $p = m$ , where $sin^{2} (mπ h /2) = sin^{2} (\frac{π}{2} - \frac{π h}{2}) = cos^{2} (\frac{π h}{2}) \to 1$ as $h \to 0$ , so $k μ_{m} \to - 4 r$ . The condition $k μ_{p} \geq - 2$ for all $p$ is $4 r sin^{2} (p π h /2) \leq 2$ , that is $r sin^{2} (p π h /2) \leq 1/2$ ; the binding case is the largest $sin^{2}$ , giving $r \leq \frac{1}{2 s i n ^{2} ( mπ h /2 )}$ . Since $sin^{2} (mπ h /2) < 1$ for every finite $m$ but tends to $1$ , the uniform condition that holds for all $h$ is $r \leq 1/2$ , i.e. $k \leq h^{2} /2$ ; this is also sufficient, as then $k μ_{p} \in [- 2, 0]$ for all $p$ . Hence $∥ r (k A_{h}) ∥_{2} \leq 1 ⟺ k \leq h^{2} /2$ .

(ii) Backward Euler and Crank-Nicolson. For backward Euler $r (z) = 1/ (1 - z)$ ; since every $k μ_{p} < 0$ , the denominator $1 - k μ_{p} > 1$ , so $∣ r (k μ_{p}) ∣ = 1/ (1 - k μ_{p}) < 1$ for all $p$ and all $k > 0$ . For Crank-Nicolson $r (z) = (1 + z /2) / (1 - z /2)$ ; for real $z = k μ_{p} \leq 0$ , $$ |r(z)|^2 = \frac{(1 + z/2)^2}{(1 - z/2)^2}, \qquad (1-z/2)^2 - (1+z/2)^2 = -2z \ge 0, $$ so $∣ r (z) ∣ \leq 1$ for every $z \leq 0$ , with equality only at $z = 0$ . Both stability functions are bounded by one on the entire negative real axis where the spectrum lives — both methods are A-stable 43.10.05 — so $∥ r (k A_{h}) ∥_{2} \leq 1$ holds for every $k > 0$ with no restriction. $□$

Bridge. This theorem is the foundational reason the parabolic time step splits into two regimes: the explicit scheme inherits the small forward-Euler stability disc, and feeding the discrete Laplacian's largest eigenvalue $\sim - 4/ h^{2}$ into that disc is exactly what forces $k \leq h^{2} /2$ , while the implicit schemes inherit A-stable regions that swallow the whole negative axis. This is exactly the eigen-decoupling of 43.10.04 carried from a generic linear system to the specific matrix $A_{h}$ built in 43.11.01, and it builds toward the von Neumann analysis of 43.11.03, where the same eigenvalues reappear as Fourier amplification factors $g (ξ)$ and the criterion $∣ g ∣ \leq 1$ recovers $r \leq 1/2$ by a different route. The stability half here is dual to the consistency half computed below — one bounds how the scheme amplifies modes, the other bounds the residual the exact solution leaves — and putting these together is the central insight that a consistent, stable parabolic scheme converges at its consistency order. The same pattern appears again in the Lax-Richtmyer equivalence theorem 43.11.05, where the power-bounded operator $r (k A_{h})^{n}$ is the abstract object whose uniform boundedness is the stability hypothesis.

Exercises Intermediate+

Exercise 3 (medium, symbolic).

Show that the FTCS update for the middle value, $U_{j}^{n + 1} = U_{j}^{n} + r (U_{j - 1}^{n} - 2 U_{j}^{n} + U_{j + 1}^{n})$ with $r = k / h^{2}$ , is a convex combination of the three old values exactly when $r \leq 1/2$ , and explain why this gives the discrete maximum principle (no overshoot).

Hint

Rewrite as $U_{j}^{n + 1} = r U_{j - 1}^{n} + (1 - 2 r) U_{j}^{n} + r U_{j + 1}^{n}$ and ask when all three coefficients are nonnegative.

Answer

Rearranging, $U_{j}^{n + 1} = r U_{j - 1}^{n} + (1 - 2 r) U_{j}^{n} + r U_{j + 1}^{n}$ . The three coefficients sum to $r + (1 - 2 r) + r = 1$ , so this is an affine combination; it is convex (all coefficients nonnegative) precisely when $r \geq 0$ and $1 - 2 r \geq 0$ , that is $r \leq 1/2$ . When $r \leq 1/2$ each new value is a weighted average of old values, so $min_{j} U_{j}^{n} \leq U_{j}^{n + 1} \leq max_{j} U_{j}^{n}$ : the scheme cannot overshoot the data, the discrete maximum principle. For $r > 1/2$ the middle coefficient $1 - 2 r$ turns negative and averaging is lost, which is the same threshold as the $ℓ^{2}$ restriction. Rubric: full credit for the coefficient sum, the convexity threshold $r \leq 1/2$ , and the no-overshoot conclusion.

Exercise 4 (medium, symbolic).

Compute the local truncation error of the FTCS scheme by Taylor expansion and show it is $O (k + h^{2})$ . Identify the leading time and space error terms.

Hint

Expand $u (x_{j}, t_{n + 1})$ in $k$ and the centred space difference in $h$ ; use $u_{t} = u_{xx}$ to simplify.

Answer

The time difference gives $\frac{u ( x _{j} , t _{n + 1} ) - u ( x _{j} , t _{n} )}{k} = u_{t} + \frac{k}{2} u_{tt} + O (k^{2})$ . The centred space difference gives $\frac{u _{j - 1} - 2 u _{j} + u _{j + 1}}{h ^{2}} = u_{xx} + \frac{h ^{2}}{12} u_{xxxx} + O (h^{4})$ , all evaluated at $(x_{j}, t_{n})$ . Subtracting and using $u_{t} = u_{xx}$ , $$ \tau^n_j = \big(u_t - u_{xx}\big) + \tfrac{k}{2}u_{tt} - \tfrac{h^2}{12}u_{xxxx} + O(k^2 + h^4) = \tfrac{k}{2}u_{tt} - \tfrac{h^2}{12}u_{xxxx} + O(k^2 + h^4). $$ So $∥ τ ∥_{\infty} = O (k + h^{2})$ : first order in time (the $\frac{k}{2} u_{tt}$ term) and second order in space (the $\frac{h ^{2}}{12} u_{xxxx}$ term). Rubric: full credit for both expansions, the cancellation of $u_{t} - u_{xx}$ , and the $O (k + h^{2})$ order.

Exercise 5 (medium, symbolic).

Show that the Crank-Nicolson local truncation error is $O (k^{2} + h^{2})$ , gaining a power in time over FTCS. Explain which feature of the trapezoidal rule causes the time order to rise.

Hint

Crank-Nicolson centres the time difference at $t_{n + 1/2}$ ; the truncation error of the trapezoidal rule is one order higher than Euler's.

Answer

Crank-Nicolson is the trapezoidal rule applied to $U^{'} = A_{h} U$ : it evaluates the time difference $\frac{u ^{n + 1} - u ^{n}}{k}$ and the spatial operator both effectively at the midpoint $t_{n + 1/2}$ . The time difference is a centred (midpoint) approximation to $u_{t} (t_{n + 1/2})$ , whose error is $O (k^{2})$ rather than the $O (k)$ of a one-sided Euler difference, because the symmetric expansion cancels the $u_{tt}$ term. The spatial second difference contributes its usual $O (h^{2})$ . Hence $τ = O (k^{2}) + O (h^{2}) = O (k^{2} + h^{2})$ . The order gain is exactly the symmetry of the trapezoidal/midpoint evaluation, the same mechanism that makes the trapezoidal rule second-order for ODEs 43.10.04. Rubric: full credit for identifying the midpoint centring, the $O (k^{2})$ time error, and the combined order.

Exercise 7 (hard, symbolic).

Prove that the Crank-Nicolson scheme is unconditionally stable in the $ℓ^{2}$ norm by diagonalising $A_{h}$ , and show that its amplification factor for the fastest mode tends to $- 1$ as $h \to 0$ at fixed mesh ratio $r$ . Interpret this last fact.

Hint

Use $∥ r (k A_{h}) ∥_{2} = max_{p} ∣ r (k μ_{p}) ∣$ with $r (z) = (1 + z /2) / (1 - z /2)$ , and study $r (k μ_{m})$ as $k μ_{m} \to - \infty$ .

Answer

Since $A_{h}$ is symmetric with orthonormal sine-mode eigenbasis, $∥ r (k A_{h}) ∥_{2} = max_{p} ∣ r (k μ_{p}) ∣$ where $r (z) = (1 + z /2) / (1 - z /2)$ . For each real $z = k μ_{p} \leq 0$ , $∣ r (z) ∣^{2} = (1 + z /2)^{2} / (1 - z /2)^{2}$ and $(1 - z /2)^{2} - (1 + z /2)^{2} = - 2 z \geq 0$ , so $∣ r (z) ∣ \leq 1$ ; thus $∥ r (k A_{h}) ∥_{2} \leq 1$ for every $k > 0$ , unconditional stability. For the fastest mode $k μ_{m} \approx - 4 r$ at fixed $r$ , but as $h \to 0$ with $k = r h^{2}$ the product $k μ_{m} \to - 4 r$ stays finite; to see the $- 1$ limit fix instead a coarse temporal resolution of the fast mode, i.e. let $k μ_{m} \to - \infty$ : then $r (z) = (2/ z + 1) / (2/ z - 1) \to - 1$ . So the stiffest unresolved spatial transient is multiplied by a factor of size approaching $1$ each step and alternates sign — it rings rather than decaying. This is why Crank-Nicolson, though A-stable, is not L-stable 43.10.05; backward Euler with $r (\infty) = 0$ damps such transients. Rubric: full credit for the spectral-norm reduction, the $∣ r ∣ \leq 1$ proof, the $r (\infty) = - 1$ limit, and the ringing interpretation.

Exercise 8 (hard, symbolic).

Using the eigenvalues $μ_{p} = - \frac{4}{h ^{2}} sin^{2} (p π h /2)$ , show that the exact semidiscrete solution $U (t) = e^{t A_{h}} U (0)$ decays mode by mode at rate $e^{μ_{p} t}$ , and compare the slowest discrete decay rate with the slowest decay rate $e^{- π^{2} t}$ of the continuous heat equation.

Hint

In the eigenbasis, $U (t) = \sum_{p} c_{p} e^{μ_{p} t} v^{(p)}$ . The slowest mode is $p = 1$ ; compare $μ_{1}$ with the continuous first eigenvalue $- π^{2}$ .

Answer

Expand the initial data in the eigenbasis, $U (0) = \sum_{p} c_{p} v^{(p)}$ . Since $A_{h} v^{(p)} = μ_{p} v^{(p)}$ , the semidiscrete solution is $U (t) = e^{t A_{h}} U (0) = \sum_{p} c_{p} e^{μ_{p} t} v^{(p)}$ : each mode decays at its own rate $e^{μ_{p} t}$ , with $μ_{p} < 0$ . The slowest is $p = 1$ : $μ_{1} = - \frac{4}{h ^{2}} sin^{2} (π h /2)$ . Using $sin (π h /2) \approx π h /2$ for small $h$ , $μ_{1} \approx - \frac{4}{h ^{2}} \cdot \frac{π ^{2} h ^{2}}{4} = - π^{2} + O (h^{2})$ . The continuous heat equation has first Dirichlet eigenvalue $π^{2}$ , so its slowest mode decays as $e^{- π^{2} t}$ ; the discrete slowest decay $e^{μ_{1} t}$ matches this to $O (h^{2})$ , confirming that the semidiscretisation reproduces the dominant long-time behaviour to second order in space. The fast modes $μ_{m} \approx - 4/ h^{2}$ decay almost instantly and carry the stiffness. Rubric: full credit for the modal decay formula, the $μ_{1} \to - π^{2}$ limit, and the comparison with the continuous rate.

Advanced results Master

Theorem 1 (stability function, A-stability, and the parabolic restriction). The full discretisation is $U^{n + 1} = r (k A_{h}) U^{n}$ with $r$ the stability function of the time integrator 43.10.04, and because $A_{h}$ is symmetric the spectral norm is $∥ r (k A_{h}) ∥_{2} = max_{p} ∣ r (k μ_{p}) ∣$ . The explicit scheme $r (z) = 1 + z$ has stability region $S = {∣1 + z ∣ \leq 1}$ meeting the negative real axis in $[- 2, 0]$ , so stability demands $k ∣ μ_{m} ∣ \leq 2$ , the parabolic restriction $k \leq h^{2} /2$ as $∣ μ_{m} ∣ \to 4/ h^{2}$ . The implicit schemes are A-stable — $S \supseteq {Re z \leq 0}$ , and the spectrum of $A_{h}$ lies on the negative real axis — so $∥ r (k A_{h}) ∥_{2} \leq 1$ for all $k$ ; backward Euler is moreover L-stable ( $r (\infty) = 0$ ) while Crank-Nicolson has $r (\infty) = - 1$ ^{[Hairer, E. & Wanner, G. — Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems (2nd ed.)]}.

Theorem 2 (consistency orders of the $θ$ -method). The $θ$ -method has local truncation error $τ = (θ - \frac{1}{2}) k u_{tt} + \frac{k ^{2}}{\dots} + \frac{h ^{2}}{12} u_{xxxx} + O (k^{2} + h^{2} \cdot k)$ , so it is $O (k + h^{2})$ for $θ \neq = \frac{1}{2}$ and $O (k^{2} + h^{2})$ for $θ = \frac{1}{2}$ , the Crank-Nicolson value, where the $k u_{tt}$ term cancels. Forward and backward Euler ( $θ = 0, 1$ ) are first order in time; Crank-Nicolson is second. The space order is $O (h^{2})$ throughout, inherited from the second-difference operator $D^{2}$ of 43.11.01 ^{[LeVeque, R. J. — Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems]}.

Theorem 3 (convergence by the discrete Duhamel/Lax synthesis). For $u \in C^{4, 2}$ the global error $E^{n} = U^{n} - u (\cdot, t_{n})$ satisfies the recursion $E^{n + 1} = r (k A_{h}) E^{n} - k τ^{n}$ with $E^{0} = 0$ , whence $E^{n} = - k \sum_{j = 0}^{n - 1} r (k A_{h})^{n - 1 - j} τ^{j}$ . Power-boundedness $∥ r (k A_{h})^{ℓ} ∥_{2} \leq 1$ (Theorem 1, in the stable regime) gives $∥ E^{n} ∥_{2} \leq k \sum_{j = 0}^{n - 1} ∥ τ^{j} ∥_{2} \leq t_{n} max_{j} ∥ τ^{j} ∥_{2}$ , so the global error is bounded by $T$ times the truncation error: $O (k + h^{2})$ for the Euler schemes and $O (k^{2} + h^{2})$ for Crank-Nicolson. This is the discrete-Duhamel form of the continuous variation-of-constants representation of 02.13.03, and the abstract version of the Lax-Richtmyer theorem 43.11.05 ^{[Thomée, V. — Galerkin Finite Element Methods for Parabolic Problems (2nd ed.)]}.

Theorem 4 (smoothing and the L-stability advantage for rough data). The continuous heat semigroup $e^{t A_{h}}$ smooths: $∥ A_{h} e^{t A_{h}} ∥_{2} \leq C / t$ , so high modes are annihilated instantly. A fully discrete scheme inherits a discrete smoothing estimate $∥ A_{h} r (k A_{h})^{n} ∥_{2} \leq C / t_{n}$ only when $r$ is L-stable, $r (\infty) = 0$ . Backward Euler enjoys this and converges at its full rate for non-smooth initial data; Crank-Nicolson, with $r (\infty) = - 1$ , fails the discrete smoothing estimate and exhibits persistent oscillations from under-resolved high modes unless the data are smooth or the first steps are damped (Rannacher time-stepping: a few backward-Euler steps before switching to Crank-Nicolson) ^{[Thomée, V. — Galerkin Finite Element Methods for Parabolic Problems (2nd ed.)]}.

Synthesis. The method of lines is the device that makes a parabolic PDE into an object the entire numerical-ODE apparatus already controls, and the foundational reason the time step behaves as it does is that the spatial operator $A_{h}$ of 43.11.01 is symmetric negative definite with spectrum stretched from $- π^{2}$ to $- 4/ h^{2}$ : this single spectral fact carries the stiffness, the parabolic restriction $k \leq h^{2} /2$ , the unconditional stability of the A-stable implicit schemes, and the smoothing that separates L-stable backward Euler from merely A-stable Crank-Nicolson. This is exactly the eigen-decoupling of 43.10.04 specialised to $A_{h}$ — each discrete eigenvalue $μ_{p}$ becomes one scalar test equation, and $k μ_{p}$ must land in the integrator's stability region — and it is dual to the elliptic theory of 43.11.01, where the same matrix's smallest eigenvalue set the conditioning while here its largest sets the time-step budget.

The central insight is that convergence is again the product of a consistency order, computed from a Taylor expansion separating the $O (k^{p})$ temporal and $O (h^{2})$ spatial errors, with a stability bound, the power-boundedness $∥ r (k A_{h})^{n} ∥_{2} \leq 1$ ; putting these together yields the $O (k + h^{2})$ and $O (k^{2} + h^{2})$ rates by a discrete Duhamel sum that generalises the continuous variation-of-constants formula of 02.13.03. The bridge to the wider FD-PDE theory is that this concrete eigenvalue computation is the special case in which the power-boundedness hypothesis of the Lax-Richtmyer equivalence theorem 43.11.05 is verified by an explicit spectral decomposition, the same hypothesis that the Fourier amplification analysis of 43.11.03 verifies mode by mode for constant-coefficient schemes on the whole line.

Full proof set Master

Proposition 1 (spectral characterisation of $ℓ^{2}$ stability). Let $A_{h}$ be real symmetric with orthonormal eigenbasis ${v^{(p)}}$ and eigenvalues ${μ_{p}}$ , and let $r$ be rational and analytic on $σ (A_{h})$ . Then $∥ r (k A_{h}) ∥_{2} = max_{p} ∣ r (k μ_{p}) ∣$ , and the scheme $U^{n + 1} = r (k A_{h}) U^{n}$ satisfies $∥ U^{n} ∥_{2} \leq ∥ U^{0} ∥_{2}$ for all $n$ if and only if $k μ_{p} \in S = {z : ∣ r (z) ∣ \leq 1}$ for every $p$ .

Proof. Diagonalise $A_{h} = Q Λ Q^{⊤}$ with $Q$ orthogonal ( $Q^{⊤} Q = I$ ) and $Λ = diag (μ_{p})$ . For a rational $r$ analytic on the spectrum, $r (k A_{h}) = Q r (k Λ) Q^{⊤}$ with $r (k Λ) = diag (r (k μ_{p}))$ . An orthogonal similarity preserves the spectral norm, so $∥ r (k A_{h}) ∥_{2} = ∥ r (k Λ) ∥_{2} = max_{p} ∣ r (k μ_{p}) ∣$ , the largest diagonal modulus of a diagonal matrix. Iterating, $r (k A_{h})^{n} = Q r (k Λ)^{n} Q^{⊤}$ , so $∥ r (k A_{h})^{n} ∥_{2} = max_{p} ∣ r (k μ_{p}) ∣^{n}$ . This is $\leq 1$ for all $n$ iff $max_{p} ∣ r (k μ_{p}) ∣ \leq 1$ , i.e. $∣ r (k μ_{p}) ∣ \leq 1$ for every $p$ , which is $k μ_{p} \in S$ for every $p$ ; and then $∥ U^{n} ∥_{2} = ∥ r (k A_{h})^{n} U^{0} ∥_{2} \leq ∥ U^{0} ∥_{2}$ . Conversely if some $k μ_{p_{0}} \in / S$ then $∣ r (k μ_{p_{0}}) ∣ > 1$ and the data $U^{0} = v^{(p_{0})}$ grows as $∣ r (k μ_{p_{0}}) ∣^{n} \to \infty$ . $□$

Proposition 2 (the explicit restriction is exactly $k \leq h^{2} /2$ ). For FTCS, $r (z) = 1 + z$ , applied to the heat semidiscretisation, $∥ r (k A_{h}) ∥_{2} \leq 1$ holds for all admissible $h = 1/ (m + 1)$ if and only if $k \leq h^{2} /2$ .

Proof. By Proposition 1, stability is $∣1 + k μ_{p} ∣ \leq 1$ for all $p$ , i.e. $- 2 \leq k μ_{p} \leq 0$ . Since $μ_{p} = - \frac{4}{h ^{2}} sin^{2} (p π h /2) < 0$ , the upper bound $k μ_{p} \leq 0$ is automatic, and the lower bound is $k \cdot \frac{4}{h ^{2}} sin^{2} (p π h /2) \leq 2$ , that is $r sin^{2} (p π h /2) \leq 1/2$ with $r = k / h^{2}$ . The strongest constraint is at the largest $sin^{2}$ , namely $p = m$ , where $sin^{2} (mπ h /2) = cos^{2} (π h /2) = 1 - sin^{2} (π h /2)$ . For a fixed grid the sharp bound is $r \leq \frac{1}{2 c o s ^{2} ( π h /2 )}$ , but since $cos^{2} (π h /2) ↑ 1$ as $h ↓ 0$ , the condition that is sufficient on every grid and necessary in the limit is $r \leq 1/2$ , i.e. $k \leq h^{2} /2$ . With $r \leq 1/2$ one has $r sin^{2} \leq 1/2$ for all $p$ , so every $k μ_{p} \in [- 2, 0]$ and $∥ r (k A_{h}) ∥_{2} \leq 1$ . $□$

Proposition 3 (unconditional stability of the implicit schemes). Backward Euler ( $r (z) = 1/ (1 - z)$ ) and Crank-Nicolson ( $r (z) = (1 + z /2) / (1 - z /2)$ ) satisfy $∥ r (k A_{h}) ∥_{2} \leq 1$ for every $k > 0$ and every $h$ .

Proof. By Proposition 1 it suffices to show $∣ r (z) ∣ \leq 1$ for every $z$ in the spectrum, and $σ (k A_{h}) \subset (- \infty, 0]$ since each $k μ_{p} < 0$ . Backward Euler: for $z < 0$ , $1 - z > 1$ , so $∣ r (z) ∣ = 1/∣1 - z ∣ = 1/ (1 - z) < 1$ . Crank-Nicolson: for $z \leq 0$ , $$ |r(z)|^2 = \frac{(1 + z/2)^2}{(1 - z/2)^2}, \qquad (1 - z/2)^2 - (1 + z/2)^2 = -2z \ge 0, $$ so the denominator dominates the numerator and $∣ r (z) ∣ \leq 1$ , with equality only at $z = 0$ . In both cases $∣ r (z) ∣ \leq 1$ throughout the spectrum, so $max_{p} ∣ r (k μ_{p}) ∣ \leq 1$ and $∥ r (k A_{h}) ∥_{2} \leq 1$ for every $k > 0$ . The pole of $r$ sits at $z = 1$ (backward Euler) or $z = 2$ (Crank-Nicolson), both off the negative axis, so $r (k A_{h})$ is well defined: $I - k A_{h}$ and $I - \frac{k}{2} A_{h}$ are positive definite, hence invertible, for every $k > 0$ . $□$

Proposition 4 (global error bound: consistency + stability ⇒ convergence). Let $u$ be smooth enough that $∥ τ^{n} ∥_{2} \leq C (k^{p} + h^{2})$ for the $θ$ -method ( $p = 1$ for $θ \neq = \frac{1}{2}$ , $p = 2$ for $θ = \frac{1}{2}$ ), and suppose the scheme is stable, $∥ r (k A_{h})^{ℓ} ∥_{2} \leq 1$ for $0 \leq ℓ k \leq T$ . Then $max_{0 \leq t_{n} \leq T} ∥ U^{n} - u (\cdot, t_{n}) ∥_{2} \leq C T (k^{p} + h^{2})$ .

Proof. The exact solution satisfies $u^{n + 1} = r (k A_{h}) u^{n} + k τ^{n}$ by the definition of $τ^{n}$ (insert $u$ into the one-step form and collect the residual), while $U^{n + 1} = r (k A_{h}) U^{n}$ . Subtract: the error $E^{n} = U^{n} - u^{n}$ obeys $E^{n + 1} = r (k A_{h}) E^{n} - k τ^{n}$ with $E^{0} = 0$ . Unrolling the recursion, $$ E^n = -k\sum_{j=0}^{n-1} r(kA_h)^{,n-1-j},\tau^j . $$ Take norms and use stability $∥ r (k A_{h})^{n - 1 - j} ∥_{2} \leq 1$ : $$ |E^n|2 \le k\sum{j=0}^{n-1}|\tau^j|_2 \le n k,\max_j|\tau^j|_2 = t_n \max_j|\tau^j|_2 \le T,C(k^p + h^2). $$ Hence the global error is $O (k^{p} + h^{2})$ , uniformly for $t_{n} \leq T$ , which is the stated convergence rate. $□$

Connections Master

The elliptic finite-difference unit 43.11.01 builds the spatial operator this unit time-steps: the matrix $A_{h} = \frac{1}{h ^{2}} tridiag (1, - 2, 1)$ is the negative of the discrete Laplacian assembled there, its sine-mode eigenvectors and eigenvalues $μ_{p} = - \frac{4}{h ^{2}} sin^{2} (p π h /2)$ are taken directly from that unit's spectral theorem, and where the elliptic unit used the smallest eigenvalue ( $\to - π^{2}$ ) to fix the conditioning $κ = O (h^{- 2})$ , this unit uses the largest ( $\to - 4/ h^{2}$ ) to fix the parabolic time-step restriction. The two units are the steady-state and time-dependent faces of the same matrix.
The stiffness and A-stability theory of 43.10.05 is the lens through which the semidiscrete system is read: $U^{'} = A_{h} U$ is the canonical stiff system, its eigenvalue spread $O (h^{- 2})$ is the textbook source of stiffness, and the A-stability that makes backward Euler and Crank-Nicolson unconditionally usable here, together with the L-stability that distinguishes them on rough data, are exactly the properties classified in that unit. The parabolic method of lines is the principal PDE application of stiff-ODE stability theory.
The absolute-stability and stability-region machinery of 43.10.04 supplies the eigen-decoupling that turns matrix stability into a scalar condition: each $k μ_{p}$ must lie in the integrator's region $S$ , the forward-Euler interval $[- 2, 0]$ yielding $k \leq h^{2} /2$ and the A-stable left half-plane yielding no restriction. The stability functions $r (z) = 1 + z$ , $1/ (1 - z)$ , $(1 + z /2) / (1 - z /2)$ used throughout this unit are the very objects defined there.
The continuous heat equation 02.13.03 is the problem being discretised: its heat kernel, its Duhamel/variation-of-constants representation, and its eigenfunction decay $e^{- π^{2} t}$ are the continuous facts the semidiscrete decay $e^{μ_{p} t}$ and the discrete Duhamel error sum of Proposition 4 approximate. The smoothing estimate that separates L-stable from A-stable schemes is the discrete shadow of the parabolic smoothing of the continuous semigroup.
The von Neumann stability analysis of 43.11.03 is the constant-coefficient Fourier counterpart of the eigenvalue argument here: where this unit diagonalises the finite matrix $A_{h}$ by sine modes on a bounded interval, the next unit substitutes a Fourier mode $e^{i ξ j h}$ into the constant-coefficient scheme on the whole line and reads stability from the amplification factor $g (ξ)$ , recovering the same restriction $r \leq 1/2$ from $∣ g (ξ) ∣ \leq 1$ . Both feed the Lax-Richtmyer equivalence theorem 43.11.05, whose power-boundedness hypothesis this unit verifies by an explicit spectral decomposition.

Historical & philosophical context Master

The idea of discretising only the space variables of a partial differential equation, leaving a system of ordinary differential equations in time, was used by Erich Rothe in 1930 for parabolic problems — the Rothe method of horizontal lines — and the complementary method of lines that discretises space and integrates the resulting ODE system in time became a standard computational technique through the 1960s, formalised in the monograph literature and in W. E. Schiesser's later systematic treatment. The explicit forward-time centred-space scheme and its stability limit $k \leq h^{2} /2$ for the heat equation were already implicit in the foundational 1928 analysis of Richard Courant, Kurt Friedrichs, and Hans Lewy ^{[Courant, Friedrichs & Lewy 1928]}, whose work on difference equations of mathematical physics identified the constraint linking the time and space steps that the next unit recasts as the CFL condition.

The unconditionally stable implicit scheme that bears their names was introduced by John Crank and Phyllis Nicolson in their 1947 paper A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type (Proceedings of the Cambridge Philosophical Society 43, 50–67) ^{[Crank & Nicolson 1947]}, motivated by the need to integrate heat-conduction problems without the punishing explicit step restriction. The systematic stability theory for these schemes — eigenvalue decoupling, A-stability, and the smoothing estimates that distinguish backward Euler from Crank-Nicolson on non-smooth data — was consolidated in the numerical-analysis literature by Vidar Thomée in Galerkin Finite Element Methods for Parabolic Problems (Springer, 1984; 2nd ed. 2006) ^{[Thomée, V. — Galerkin Finite Element Methods for Parabolic Problems (2nd ed.)]} and, on the stiff-ODE side, by Ernst Hairer and Gerhard Wanner ^{[Hairer, E. & Wanner, G. — Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems (2nd ed.)]}.

Bibliography Master

@book{leveque2007fdmparabolic,
  author    = {LeVeque, Randall J.},
  title     = {Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems},
  publisher = {Society for Industrial and Applied Mathematics (SIAM)},
  year      = {2007}
}

@article{cranknicolson1947,
  author  = {Crank, John and Nicolson, Phyllis},
  title   = {A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type},
  journal = {Proceedings of the Cambridge Philosophical Society},
  volume  = {43},
  number  = {1},
  year    = {1947},
  pages   = {50--67}
}

@article{courantfriedrichslewy1928parabolic,
  author  = {Courant, Richard and Friedrichs, Kurt and Lewy, Hans},
  title   = {{\"U}ber die partiellen Differenzengleichungen der mathematischen Physik},
  journal = {Mathematische Annalen},
  volume  = {100},
  number  = {1},
  year    = {1928},
  pages   = {32--74}
}

@book{thomee2006galerkin,
  author    = {Thom\'{e}e, Vidar},
  title     = {Galerkin Finite Element Methods for Parabolic Problems},
  edition   = {2},
  series    = {Springer Series in Computational Mathematics},
  volume    = {25},
  publisher = {Springer-Verlag},
  year      = {2006}
}

@book{hairerwanner1996parabolic,
  author    = {Hairer, Ernst and Wanner, Gerhard},
  title     = {Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems},
  edition   = {2},
  series    = {Springer Series in Computational Mathematics},
  volume    = {14},
  publisher = {Springer-Verlag},
  year      = {1996}
}

@book{iserles2009firstparabolic,
  author    = {Iserles, Arieh},
  title     = {A First Course in the Numerical Analysis of Differential Equations},
  edition   = {2},
  publisher = {Cambridge University Press},
  year      = {2009}
}

@book{strang2007cseparabolic,
  author    = {Strang, Gilbert},
  title     = {Computational Science and Engineering},
  publisher = {Wellesley-Cambridge Press},
  year      = {2007}
}

Prerequisites

43.11.01
43.10.05
02.13.03

Tier anchors

beginner: LeVeque 2007 *Finite Difference Methods for Ordinary and Partial Differential Equations* (SIAM) §9.1-9.4 (the method of lines for the heat equation, the explicit forward-time scheme, and the step restriction $k\le h^2/2$); Strang 2007 *Computational Science and Engineering* (Wellesley-Cambridge) §6.4 (the heat equation discretised, explicit vs implicit time-stepping, and the stability limit)
intermediate: LeVeque 2007 *Finite Difference Methods for Ordinary and Partial Differential Equations* (SIAM) §9.1-9.5 (the semidiscrete system $U'=A_hU$, forward Euler/FTCS with $k\le h^2/2$, backward Euler, Crank-Nicolson, and the eigenvalue-based absolute-stability argument); Iserles 2009 *A First Course in the Numerical Analysis of Differential Equations* 2e (Cambridge) §16.1-16.4 (semidiscretisation, the heat-equation schemes, and their orders and stability)
master: LeVeque 2007 *Finite Difference Methods for Ordinary and Partial Differential Equations* (SIAM) Ch. 9-10 (the method of lines, the stiff semidiscrete system, the $\theta$-method family, accuracy and absolute-stability analysis, and the $O(k+h^2)$/$O(k^2+h^2)$ convergence orders); Thomée 2006 *Galerkin Finite Element Methods for Parabolic Problems* 2e (Springer) Ch. 1, 7-9 (the semigroup/semidiscrete framework, the backward Euler and Crank-Nicolson time discretisations, and the smoothing/stability estimates); Hairer-Wanner 1996 *Solving Ordinary Differential Equations II* 2e (Springer) §IV.2-IV.3 (stiffness, A- and L-stability, applied to the parabolic semidiscretisation)

References

LeVeque, R. J. — Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems · SIAM 2007. Chapters 9-10 develop the time-dependent theory. §9.1 introduces the method of lines: discretise $u_t=u_{xx}$ in space by the second-difference operator $D^2$ to obtain the stiff ODE system $U'(t)=A_hU(t)$ with $A_h=\tfrac1{h^2}\mathrm{tridiag}(1,-2,1)$, then apply an ODE integrator in time. §9.2 derives the explicit forward-time centred-space (FTCS) scheme $U^{n+1}=U^n+kA_hU^n$ and its stability restriction $k\le h^2/2$ by requiring $k\lambda_p(A_h)$ to lie in the forward-Euler absolute-stability interval $[-2,0]$, where $\lambda_p=-\tfrac4{h^2}\sin^2(p\pi h/2)\to-4/h^2$. §9.3 treats the implicit backward-Euler scheme $(I-kA_h)U^{n+1}=U^n$ (unconditionally stable, one tridiagonal solve per step) and the Crank-Nicolson scheme $(I-\tfrac k2A_h)U^{n+1}=(I+\tfrac k2A_h)U^n$ (trapezoidal, second-order in time, unconditionally stable). §9.5 assembles the local truncation orders into the global convergence rates $O(k+h^2)$ for the Euler schemes and $O(k^2+h^2)$ for Crank-Nicolson, and discusses the stiffness of $A_h$ (eigenvalues spread from $O(1)$ to $-4/h^2$).
Thomée, V. — Galerkin Finite Element Methods for Parabolic Problems (2nd ed.) · Springer 2006. Chapters 1 and 7-10 give the semidiscrete/fully discrete framework for parabolic problems: the spatially semidiscrete equation $U'+A_hU=0$ as a stiff system whose solution operator is $e^{-tA_h}$; the backward Euler discretisation $U^{n}=(I+kA_h)^{-1}U^{n-1}=r(kA_h)U^{n-1}$ with $r(z)=1/(1+z)$, and the Crank-Nicolson choice $r(z)=(1-z/2)/(1+z/2)$ (here with $A_h$ positive, the discrete $-\Delta_h$); the unconditional stability $\|r(kA_h)\|\le1$ following from $|r(z)|\le1$ for $z\ge0$ and the spectral theorem on the symmetric $A_h$; smoothing estimates and the $O(k^2+h^2)$ error bound for Crank-Nicolson under suitable regularity.
Hairer, E. & Wanner, G. — Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems (2nd ed.) · Springer 1996. §IV.2-IV.3 develop stiffness and A-/L-stability: the method-of-lines semidiscretisation of a parabolic PDE is the canonical example of a stiff system, with eigenvalues spreading to $-4/h^2$ as the mesh refines; A-stability (the stability region contains the whole left half-plane) makes backward Euler and Crank-Nicolson unconditionally stable on it, while the explicit forward-Euler region forces $k\lambda\in[-2,0]$ and hence $k\le h^2/2$; L-stability ($R(\infty)=0$) of backward Euler versus the trapezoidal $R(\infty)=-1$ of Crank-Nicolson, explaining the damping of unresolved high-frequency parabolic transients.

Estimated time

beginner: 20m
intermediate: 50m
master: 90m