Crank-Nicolson for the Heat Equation

Trade an explicit stability limit for an unconditionally stable, second-order-in-time implicit scheme that solves a tridiagonal system each step.

Motivation

The explicit FTCS scheme demands dt <= dx^2/(2 alpha). Crank-Nicolson averages the spatial operator between the old and new time levels, giving a method that is stable for any dt and second-order accurate in both space and time.

The scheme

With r = alpha dt / dx^2, the update is -r/2 u_new[i-1] + (1+r) u_new[i] - r/2 u_new[i+1] = r/2 u[i-1] + (1-r) u[i] + r/2 u[i+1]. The left side is a tridiagonal matrix applied to the unknowns; the right side is known.

import numpy as np
nx=51; dx=1/(nx-1); alpha=0.01; dt=0.5; r=alpha*dt/dx**2
u=np.sin(np.pi*np.linspace(0,1,nx)); u[0]=u[-1]=0
n=nx-2
A=np.diag((1+r)*np.ones(n))+np.diag(-r/2*np.ones(n-1),1)+np.diag(-r/2*np.ones(n-1),-1)
B=np.diag((1-r)*np.ones(n))+np.diag(r/2*np.ones(n-1),1)+np.diag(r/2*np.ones(n-1),-1)
for step in range(200):
    u[1:-1]=np.linalg.solve(A,B@u[1:-1])
print(round(u.max(),5))

Cost versus benefit

Each step now solves a linear system, but a tridiagonal solve is order n via the Thomas algorithm - cheap. In return dt is limited only by accuracy, not stability, so stiff diffusion problems run in far fewer, larger steps.

A caution

Crank-Nicolson is stable but not strongly damping: with very large dt it can produce slowly decaying oscillations on sharp initial data. When robustness matters more than order, the fully implicit backward-Euler scheme (r on the new level only) damps those artifacts at the price of first-order time accuracy.