Implicit schemes and Crank–Nicolson

The implicit scheme

Lesson 6.5.1 replaced $u_{xx}$ in $u_t = u_{xx}$ by the second difference at the old time $t_n$. Take it at the new time $t_{n + 1}$ instead:

Multiply by $\Delta t$, write $r = \frac{\Delta t}{\Delta x^2}$ as before, and gather the new values on the left:

Seen from the new time, the time difference is now a backward difference (lesson 2.2.2). Its error is about $-\frac{1}{2}u_{tt}\,\Delta t$ (lesson 2.8.1), the same size as the forward difference's, so the implicit scheme is still first order in time.

One step by hand

Take $M = 4$ cells with both ends held at 0. Start with a spike, 1 at the middle point and 0 elsewhere, and take $r = 1$.

Call the three new inside values $a$, $b$ and $c$. Then the rule gives

Elimination (lesson 3.3.1) clears the matrix below its diagonal in two steps:

Back substitution gives $c = \frac{3}{8} \div \frac{21}{8} = \frac{1}{7}$, then $b = (1 + c) \div \frac{8}{3} = \frac{3}{7}$, then $a = \frac{b}{3} = \frac{1}{7}$.

At $r = 1$ the explicit scheme is unstable, and its middle value has gone negative, which heat never does. The implicit step spreads the spike smoothly. Its values add up to $\frac{5}{7}$, less than 1, because heat leaks out through the cold ends.

Boundary values

The first and last equations contain an end value. In the first, $-r\,u_0^{n+1}$ is known from the boundary condition, so it moves to the right side: $(1 + 2r)\,u_1^{n+1} - r\,u_2^{n+1} = u_1^n + r\,u_0^{n+1}$.

The last equation gains $r\,u_M^{n+1}$ the same way. With both ends at 0, as here, nothing is added.

Tridiagonal systems

With more points the matrix keeps the same pattern: $1 + 2r$ down the diagonal, $-r$ just above and just below it, and 0 everywhere else. For four inside points:

A matrix whose only nonzero entries lie on the diagonal and the two lines beside it is tridiagonal. Each step of the implicit scheme solves $A\mathbf{u}^{n+1} = \mathbf{b}$, where $\mathbf{u}^{n+1}$ holds the new inside values and $\mathbf{b}$ the old ones plus the boundary terms.

Elimination is cheap here. Each pivot has just one nonzero entry below it, and clearing it adds nothing new to the rows below. So the work grows like the number of points, not like its cube, as it does for a full matrix.

It is also safe. Each diagonal entry, $1 + 2r$, is larger than the two others in its row put together, $2r$. So elimination never meets a small pivot and needs no row swaps: the pivots above were 3, $\frac{8}{3}$ and $\frac{21}{8}$, all above 2.

The condition number of lesson 3.3.2 is small too: 2.78 for the matrix above, and never more than $1 + 4r$.

For the example, ab is $\begin{pmatrix} \ast & -1 & -1 \\ 3 & 3 & 3 \\ -1 & -1 & \ast \end{pmatrix}$, where $\ast$ marks the unused corners.

A full matrix wastes room. With 10,000 points it holds $10^8$ numbers, 800 MB, almost all of them 0. The three lines hold 30,000 numbers.

Stable for every r

Lesson 6.5.1 measured stability with a zig-zag, whose second difference is $-4$ times itself. Put it into the implicit rule: $(1 + 4r)\,u^{n+1} = u^n$, so each step multiplies the zig-zag by $\frac{1}{1 + 4r}$.

That factor lies between 0 and 1 for every $r > 0$. The zig-zag shrinks without flipping sign, and so does every smoother wave. The implicit scheme is unconditionally stable: no limit on $r$.

Lesson 6.5.1's tent, $\min(x, 1 - x)$ on 20 cells, shows it. At $r = 5$, ten times the explicit limit, the largest value falls to 0.3909, 0.3339, 0.2284 and 0.1272 after 1, 2, 5 and 10 steps. No value ever goes negative.

Now accuracy limits the step.

Take lesson 6.5.1's rod with its ends at 0 and 1. In 20 implicit steps of $\Delta t = 0.05$ ($r = 5$) it reaches the line $u = x$, to within $2.19 \times 10^{-4}$ at $t = 1$. The explicit scheme took 250 steps but got within $2.91 \times 10^{-5}$: bigger steps, bigger errors.

Crank–Nicolson

The explicit scheme takes $u_{xx}$ at the old time, the implicit one at the new time. The Crank–Nicolson scheme takes the average of the two. Write $\delta^2 u_i = u_{i+1} - 2u_i + u_{i-1}$ for the second difference at point $i$:

Gather the new values on the left and the old ones on the right:

Why the average helps: $\frac{u_i^{n+1} - u_i^n}{\Delta t}$ is a central difference about the half step $t_n + \frac{1}{2}\Delta t$. The average of the two second differences sits at that half step too. Lesson 2.8.1's central difference has an error in $h^2$, not $h$, and so does this one: Crank–Nicolson's time error shrinks like $\Delta t^2$.

The zig-zag factor is $\frac{1 - 2r}{1 + 2r}$, which lies between $-1$ and 1 for every $r$. So Crank–Nicolson is unconditionally stable too.

Back to the spike with $r = 1$. The right side is $(0.5, 0, 0.5)$, and the system $2a - 0.5b = 0.5$, $-0.5a + 2b - 0.5c = 0$, $-0.5b + 2c = 0.5$ gives $a = c = \frac{2}{7} \approx 0.2857$ and $b = \frac{1}{7} \approx 0.1429$.

The middle has fallen below its neighbours: a dip the true solution never has. The last section of this lesson comes back to it.

Three schemes compared

Take lesson 6.5.1's rod with the profile $\sin(\pi x)$, whose exact solution is $e^{-\pi^2 t}\sin(\pi x)$, and run to $t = 0.1$. Refine space and time together, keeping $\Delta t = \frac{\Delta x}{10}$, so $r$ grows from 1 to 8: far too large for the explicit scheme. It runs instead at $r = 0.4$, with the many more steps that needs.

Each halving of both steps halves the implicit error: its $\Delta t$ term is the larger one. It divides Crank–Nicolson's error by 4, since both its terms are squares. On 80 cells, Crank–Nicolson takes 80 steps to beat what the explicit scheme gets in 1,600.

Oscillations near a kink

Crank–Nicolson has one weakness. Start from the kink $\max(x, 0)$ on $[-3, 3]$, with the ends held, and take large steps: $\Delta x = 0.1$ and $\Delta t = 0.1$, so $r = 10$. After 5 steps, at $t = 0.5$, compare with the exact solution of lesson 6.4.2, $u = x\,N\!\left(\frac{x}{\sqrt{2t}}\right) + \sqrt{2t}\,\varphi\!\left(\frac{x}{\sqrt{2t}}\right)$.

Here $N$ and $\varphi$ are the standard normal distribution function and density. The values look fine: for $|x| \le 1.5$ they are within $1.0 \times 10^{-2}$ of the exact solution for the implicit scheme, $1.2 \times 10^{-2}$ for Crank–Nicolson. The second differences tell another story:

These are spurious oscillations: a zig-zag that the true solution does not have. The kink holds a little of every wave, the zig-zag included. Crank–Nicolson multiplies the zig-zag by $\frac{1 - 2r}{1 + 2r} = -0.9048$ per step, so it flips sign each step and hardly shrinks.

The oscillation hides in the values and shows in the curvature. Case study 6.5.4 reads an option's gamma, a second derivative, off such a grid, so it matters there.

Two remedies are standard. Keep $r$ moderate near the kink: with $\Delta t = 0.01$ ($r = 1$) the factor is $-\frac{1}{3}$ and the zig-zag dies quickly.

Or start with the implicit scheme, which removes the zig-zag at once. Replace the first step by two implicit steps of half the length, then carry on with Crank–Nicolson. This is Rannacher's start, after the mathematician who studied it.

With Rannacher's start the values are within $1.2 \times 10^{-3}$, ten times closer than either scheme alone. The second differences are 0.383, 0.398, 0.424, 0.398 and 0.383, with a largest error of 0.025 near the kink.

Common mistakes

1. Forgetting the boundary terms. A known end value moves to the right side of the first or last equation. Leaving it out holds that end at 0.
2. Building the full matrix. np.linalg.solve gives the right answer, but its memory grows like the square of the points and its work like the cube. Use solve_banded.
3. Misplacing the lines in ab. The line above the diagonal is shifted right, ab[0, 1:]; the line below is shifted left, ab[2, :-1].
4. Taking stable for accurate. The implicit scheme survives any $\Delta t$, but its error grows with $\Delta t$.
5. Copying the explicit weights into Crank–Nicolson. Its right side has $1 - r$ in the middle, not $1 - 2r$: each side carries half the second difference.
6. Trusting Crank–Nicolson's curvature near a kink. With large $r$ the values look smooth while the second differences zig-zag. Take smaller steps there, or start with implicit steps.