Finite-Differencing the 1D Heat Equation

Solve u_t = alpha u_xx on a rod with an explicit forward-time centered-space scheme, and watch the stability limit bite.

The model

The heat equation u_t = alpha u_xx describes diffusion of temperature along a rod. Given an initial profile and fixed end temperatures, we march it forward in time on a grid.

The FTCS scheme

Forward difference in time, centered second difference in space gives u_new[i] = u[i] + r*(u[i+1] - 2u[i] + u[i-1]), where r = alpha*dt/dx^2. Each new value is a weighted blend of the old neighbours.

Stability

This explicit scheme is stable only when r <= 1/2. Exceed it and errors amplify each step, producing wild oscillations. To halve dx you must quarter dt - a strong constraint that motivates implicit schemes.

import numpy as np
nx=51; L=1.0; dx=L/(nx-1); alpha=0.01
dt=0.4*dx*dx/alpha            # r=0.4 < 0.5, stable
r=alpha*dt/dx**2
u=np.sin(np.pi*np.linspace(0,L,nx)); u[0]=u[-1]=0
for n in range(2000):
    u[1:-1]=u[1:-1]+r*(u[2:]-2*u[1:-1]+u[:-2])
print('peak temp:',round(u.max(),4))  # decays toward 0

What you should see

A single sine hump decays smoothly toward zero, and its analytic decay rate is exp(-alpha pi^2 t). Comparing the numerical peak to that exponential is a good correctness check. Sharp initial features diffuse fastest, exactly as physical heat spreads.

The same stencil, made implicit, becomes the Crank-Nicolson scheme, which removes the stability limit at the cost of solving a tridiagonal system each step.