Solving the 2D Heat Equation

March a discretized 2D diffusion problem forward in time with an explicit five-point stencil, and check the stability limit that constrains the time step.

Problem

The heat equation u_t = alpha (u_xx + u_yy) describes temperature spreading across a plate. We solve it on a unit square with a hot spot in the center and zero temperature on all edges, using a uniform grid and explicit time stepping.

Discretization

Replace second derivatives with the five-point Laplacian: (u_xx+u_yy) at (i,j) is approximately (u[i+1,j]+u[i-1,j]+u[i,j+1]+u[i,j-1]-4 u[i,j]) / dx^2. An explicit Euler step then updates every interior node from its four neighbors.

import numpy as np
n=41; dx=1.0/(n-1); alpha=1.0
dt=0.2*dx*dx/alpha         # inside stability limit
u=np.zeros((n,n)); u[n//2,n//2]=1.0/dx/dx
lap=lambda u:(np.roll(u,1,0)+np.roll(u,-1,0)+np.roll(u,1,1)+np.roll(u,-1,1)-4*u)/dx/dx
for _ in range(2000):
    u[1:-1,1:-1]+=dt*alpha*lap(u)[1:-1,1:-1]
    u[0,:]=u[-1,:]=u[:,0]=u[:,-1]=0.0
print(u.max(), u.sum()*dx*dx)  # peak decays, integral conserved until it reaches edges

Stability

The explicit scheme is stable only when dt <= dx^2 / (4 alpha) in two dimensions. Above that limit small errors grow geometrically and the solution explodes. We chose dt = 0.2 dx^2 / alpha, comfortably inside the bound. Implicit schemes like Crank-Nicolson remove this restriction at the cost of solving a linear system each step.

• Total heat is conserved while the profile stays away from the boundaries, then decays as it leaks out the fixed edges.
• Halving dx quarters the allowed dt, so explicit 2D diffusion becomes expensive on fine grids.
• The same stencil models transient thermal loads on burner first-wall panels; production runs use implicit solvers for the stiff time scales.