Crank-Nicolson for 1D Diffusion

Step a 1D heat equation with the unconditionally stable Crank-Nicolson scheme by solving a tridiagonal system each step.

Problem

Crank-Nicolson averages the explicit and implicit updates, giving second-order accuracy in time and unconditional stability. Unlike the explicit scheme it places no upper bound on the time step, at the cost of solving a tridiagonal linear system each step.

Scheme

For u_t = alpha u_xx, define r = alpha dt / (2 dx^2). The update reads (1+2r) u_i^{n+1} - r(u_{i+1}^{n+1}+u_{i-1}^{n+1}) = (1-2r) u_i^n + r(u_{i+1}^n+u_{i-1}^n). The left side is a tridiagonal matrix solved with the Thomas algorithm.

import numpy as np
n=51; dx=1/(n-1); alpha=1.0; dt=0.01; r=alpha*dt/(2*dx*dx)
x=np.linspace(0,1,n); u=np.sin(np.pi*x)
A=(1+2*r)*np.eye(n)-r*np.eye(n,k=1)-r*np.eye(n,k=-1)
B=(1-2*r)*np.eye(n)+r*np.eye(n,k=1)+r*np.eye(n,k=-1)
A[0]=A[-1]=0; A[0,0]=A[-1,-1]=1     # Dirichlet ends
for _ in range(50):
    b=B@u; b[0]=b[-1]=0
    u=np.linalg.solve(A,b)
exact=np.sin(np.pi*x)*np.exp(-np.pi**2*0.5)
print('max error',round(np.max(np.abs(u-exact)),5))

Result

Starting from a sine profile, the analytic solution decays as exp(-pi^2 alpha t). Crank-Nicolson tracks this decay accurately even with a step size that would make the explicit scheme blow up, and the maximum error stays small. The scheme can produce mild oscillations for very large r on sharp initial data, a known trade-off for its high accuracy.

• Unconditional stability lets you choose dt for accuracy, not stability, a major advantage for stiff diffusion.
• Each step is a tridiagonal solve, cheap in 1D and still manageable in higher dimensions with ADI splitting.
• Kronos thermal-transient models of first-wall components use implicit schemes for the same stiffness reasons.