23/05/2026, 14:59 IAM_Programming_Project

Numerical solution of the 2D heat equation

on a metal plate
The equation:

$$\frac{\partial u}{\partial t} = \alpha \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2

u}{\partial y^2} \right)$$

where $u(x, y, t)$ is temperature, $\alpha$ is thermal diffusivity, and the domain is $[0, L_x]
\times [0, L_y]$.

Two schemes are implemented: explicit FTCS and Crank–Nicolson ADI. The notebook also
covers stability limits, a convergence study, and verification against an exact analytical
solution where one exists.

1. Imports and configuration

In [1]: import numpy as np
import [Link] as plt
import [Link] as animation
from [Link] import GridSpec
from [Link] import solve
from [Link] import lil_matrix, csc_matrix
from [Link] import spsolve
import time
import warnings
[Link]('ignore')

%matplotlib inline
[Link]({
'[Link]': 120,
'[Link]': 11,
'[Link]': 13,
'[Link]': 12,
'[Link]': 'hot'
})

print("Libraries loaded successfully.")

Libraries loaded successfully.

2. Problem setup

[Link] 1/22
23/05/2026, 14:59 IAM_Programming_Project

The plate has fixed-temperature edges (Dirichlet conditions). Initial condition options:
Gaussian hot spot, linear gradient, sinusoidal, or uniform.

Spatial grid: $\Delta x = L_x/(N_x-1)$, $\Delta y = L_y/(N_y-1)$. Time step $\Delta t$ is fixed.

In [2]: # ─── Physical Parameters ───────────────────────────────────────────────────

Lx, Ly = 1.0, 1.0 # Plate dimensions (m)
alpha = 1.0e-4 # Thermal diffusivity (m²/s), ~stainless steel
run_time = 5000.0 # Total simulation time (s)

# ─── Boundary Conditions (Dirichlet) ──────────────────────────────────────

temp_left = 100.0 # °C (x = 0)
temp_right = 0.0 # °C (x = Lx)
temp_bottom = 0.0 # °C (y = 0)
temp_top = 0.0 # °C (y = Ly)

# ─── Grid Resolution ──────────────────────────────────────────────────────

Nx, Ny = 51, 51 # Number of grid points (including boundaries)
dx = Lx / (Nx - 1)
dy = Ly / (Ny - 1)

x = [Link](0, Lx, Nx)

y = [Link](0, Ly, Ny)
X, Y = [Link](x, y) # shape (Ny, Nx)

print(f"Grid: {Nx} × {Ny} points | dx = {dx:.4f} m | dy = {dy:.4f} m")

print(f"Thermal diffusivity α = {alpha:.2e} m²/s")
print(f"Simulation time: {run_time:.0f} s")

Grid: 51 × 51 points | dx = 0.0200 m | dy = 0.0200 m

Thermal diffusivity α = 1.00e-04 m²/s
Simulation time: 5000 s

In [3]: def make_initial_condition(Nx, Ny, x, y, start_type='hotspot'):

"""
Generate initial temperature field.
start_type: 'hotspot' | 'gradient' | 'sine' | 'uniform'
"""
X_, Y_ = [Link](x, y) # shape (Ny, Nx)
u = [Link]((Ny, Nx))

if start_type == 'hotspot':
# Gaussian hot spot at center
cx, cy = Lx / 2, Ly / 2
sigma = 0.1
u = 200 * [Link](-((X_ - cx)**2 + (Y_ - cy)**2) / (2 * sigma**2))

elif start_type == 'gradient':

# Linear gradient in x
u = 100 * (1 - X_ / Lx)

elif start_type == 'sine':

# Sinusoidal — has a known analytical solution
u = 100 * [Link]([Link] * X_ / Lx) * [Link]([Link] * Y_ / Ly)

[Link] 2/22
23/05/2026, 14:59 IAM_Programming_Project

elif start_type == 'uniform':

u = [Link]((Ny, Nx), 100.0)

# Enforce Dirichlet BCs

u = apply_bc(u)
return u

def apply_bc(u):
"""Apply Dirichlet boundary conditions."""
u[:, 0] = temp_left # left edge
u[:, -1] = temp_right # right edge
u[0, :] = temp_bottom # bottom edge
u[-1, :] = temp_top # top edge
return u

# ─── Visualise all ICs ─────────────────────────────────────────────────────

fig, axes = [Link](1, 4, figsize=(18, 4))
ic_types = ['hotspot', 'gradient', 'sine', 'uniform']
titles = ['Gaussian Hot Spot', 'Linear Gradient', 'Sinusoidal', 'Uniform 100 °C'

for ax, ic, title in zip(axes, ic_types, titles):

u_start = make_initial_condition(Nx, Ny, x, y, start_type=ic)
im = [Link](X, Y, u_start, levels=20, cmap='hot')
[Link](im, ax=ax, label='T (°C)')
ax.set_title(title)
ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')

[Link]('Initial Condition Options', fontsize=14, fontweight='bold')

plt.tight_layout()
[Link]()

3. Stability analysis
Explicit (FTCS): Von Neumann criterion
For the 2D explicit scheme, the stability condition is:

$$r = \alpha \frac{\Delta t}{\Delta x^2} \leq \frac{1}{4} \quad (\text{for equal spacing})$$

More generally:

[Link] 3/22
23/05/2026, 14:59 IAM_Programming_Project

$$\alpha \Delta t \left(\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2}\right) \leq \frac{1}{2}$$

Crank–Nicolson is unconditionally stable, so $\Delta t$ can be as large as accuracy allows.

In [4]: def stability_check(alpha, dt, dx, dy):

"""Compute and report the stability parameter r for the explicit scheme."""
r = alpha * dt * (1/dx**2 + 1/dy**2)
rx = alpha * dt / dx**2
ry = alpha * dt / dy**2
stable = r <= 0.5
print(f" rx = α·Δt/Δx² = {rx:.4f}")
print(f" ry = α·Δt/Δy² = {ry:.4f}")
print(f" r = rx + ry = {r:.4f} (must be ≤ 0.5 for stability)")
print(f" Stable: {'✓ YES' if stable else '✗ NO — UNSTABLE!'}")
return r, stable

# ─── Maximum stable dt ─────────────────────────────────────────────────────

max_dt = 0.5 / (alpha * (1/dx**2 + 1/dy**2))
dt_safe = 0.9 * max_dt # slightly below limit
dt_bad = 1.5 * max_dt # above limit (will blow up)

print(f"Maximum stable Δt = {max_dt:.2f} s\n")

print("--- Stable case ---")
stability_check(alpha, dt_safe, dx, dy)
print("\n--- Unstable case ---")
stability_check(alpha, dt_bad, dx, dy)

# ─── Stability map: vary dt and alpha ─────────────────────────────────────

alphas = [Link](-6, -2, 200)
dt_max_arr = 0.5 / (alphas * (1/dx**2 + 1/dy**2))

fig, ax = [Link](figsize=(8, 4))

[Link](alphas, dt_max_arr, 'b-', lw=2, label='Stability boundary (r = 0.5)')
ax.fill_between(alphas, dt_max_arr, 1e8, alpha=0.15, color='green', label='Stable r
ax.fill_between(alphas, 0, dt_max_arr, alpha=0.15, color='red', label='Unstable r
[Link]([alpha], [dt_safe], color='green', s=80, zorder=5, label=f'Our stable
[Link]([alpha], [dt_bad], color='red', s=80, zorder=5, label=f'Unstable dt={d
ax.set_xlabel('Thermal diffusivity α (m²/s)')
ax.set_ylabel('Maximum stable Δt (s)')
ax.set_title('FTCS Stability Boundary for 2D Heat Equation')
[Link](fontsize=9)
[Link](True, which='both', ls='--', alpha=0.5)
plt.tight_layout()
[Link]()

[Link] 4/22
23/05/2026, 14:59 IAM_Programming_Project

Maximum stable Δt = 1.00 s

--- Stable case ---

rx = α·Δt/Δx² = 0.2250
ry = α·Δt/Δy² = 0.2250
r = rx + ry = 0.4500 (must be ≤ 0.5 for stability)
Stable: ✓ YES

--- Unstable case ---

rx = α·Δt/Δx² = 0.3750
ry = α·Δt/Δy² = 0.3750
r = rx + ry = 0.7500 (must be ≤ 0.5 for stability)
Stable: ✗ NO — UNSTABLE!

4. Explicit (FTCS) method

The forward-time, centered-space update at each interior point:

$$u_{i,j}^{n+1} = u_{i,j}^n + r_x\left(u_{i+1,j}^n - 2u_{i,j}^n + u_{i-1,j}^n\right) +

r_y\left(u_{i,j+1}^n - 2u_{i,j}^n + u_{i,j-1}^n\right)$$

where $r_x = \alpha\Delta t/\Delta x^2$ and $r_y = \alpha\Delta t/\Delta y^2$.

In [5]: def solve_explicit(u_start, alpha, dx, dy, dt, num_steps, save_every=100):

"""
Explicit (FTCS) solver for 2D heat equation.
Returns list of saved snapshots and corresponding times.
"""
rx = alpha * dt / dx**2
ry = alpha * dt / dy**2

u = u_start.copy()
snapshots = [[Link]()]
times = [0.0]

[Link] 5/22
23/05/2026, 14:59 IAM_Programming_Project

for step in range(1, num_steps + 1):

u_next = [Link]()
# Interior points only (vectorised)
u_next[1:-1, 1:-1] = (
u[1:-1, 1:-1]
+ rx * (u[1:-1, 2:] - 2*u[1:-1, 1:-1] + u[1:-1, :-2])
+ ry * (u[2:, 1:-1] - 2*u[1:-1, 1:-1] + u[:-2, 1:-1])
)
apply_bc(u_next)
u = u_next

if step % save_every == 0:
[Link]([Link]())
[Link](step * dt)

return snapshots, times

# ─── Run Explicit solver ───────────────────────────────────────────────────

dt_explicit = 0.9 * max_dt
num_steps = int(run_time / dt_explicit)
save_step_exp = max(1, num_steps // 50)

u_start = make_initial_condition(Nx, Ny, x, y, start_type='hotspot')

print(f"Explicit solver: dt = {dt_explicit:.2f} s | steps = {num_steps:,}")

t0 = [Link]()
frames_exp, times_exp = solve_explicit(u_start, alpha, dx, dy, dt_explicit, num_ste
print(f"Wall time: {[Link]()-t0:.2f} s | Snapshots saved: {len(frames_exp)}")

Explicit solver: dt = 0.90 s | steps = 5,555

Wall time: 0.26 s | Snapshots saved: 51

In [6]: # ─── Plot snapshots at 4 time points ──────────────────────────────────────

frame_picks = [0, len(frames_exp)//4, len(frames_exp)//2, -1]
fig, axes = [Link](1, 4, figsize=(18, 4))

vmin = min([Link]() for s in frames_exp)

vmax = max([Link]() for s in frames_exp)

for ax, idx in zip(axes, frame_picks):

im = [Link](X, Y, frames_exp[idx], levels=25, cmap='hot', vmin=vmin, vmax=
[Link](im, ax=ax, label='T (°C)')
ax.set_title(f't = {times_exp[idx]:.0f} s')
ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')

[Link]('Explicit (FTCS) — Temperature Evolution (Hot Spot IC)', fontsize=13,

plt.tight_layout()
[Link]()

[Link] 6/22
23/05/2026, 14:59 IAM_Programming_Project

5. Crank–Nicolson method
CN averages explicit and implicit fluxes. The result is second-order accuracy in time with no
stability constraint on $\Delta t$:

$$\frac{u^{n+1} - u^n}{\Delta t} = \frac{\alpha}{2}\left(\delta^2_x +

\delta^2_y\right)\left(u^{n+1} + u^n\right)$$

Solving the full 2D implicit system at each step is expensive. The ADI (Alternating Direction
Implicit) approach splits it into two sweeps:

Half-step: implicit in $x$, explicit in $y$

Full-step: explicit in $x$, implicit in $y$

Each sweep reduces to a 1D tridiagonal solve.

In [7]: def build_tridiagonal(N, r):

"""
Build (N-2)x(N-2) tridiagonal matrix for CN 1D sub-problem.
A = I + r * L where L is the 1D Laplacian.
"""
n = N - 2 # interior points only
diag = [Link](n, 1 + 2*r)
off_diag = [Link](n-1, -r)
A = [Link](diag) + [Link](off_diag, 1) + [Link](off_diag, -1)
return A

def solve_crank_nicolson_adi(u_start, alpha, dx, dy, dt, num_steps, save_every=100)

"""
Crank–Nicolson via Alternating Direction Implicit (ADI).
Each full time step = two half-steps with 1D tridiagonal solves.
"""
rx = alpha * dt / (2 * dx**2)
ry = alpha * dt / (2 * dy**2)
num_cols = u_start.shape[1]
num_rows = u_start.shape[0]

Ax = build_tridiagonal(num_cols, rx) # x-direction implicit matrix

Ay = build_tridiagonal(num_rows, ry) # y-direction implicit matrix

[Link] 7/22
23/05/2026, 14:59 IAM_Programming_Project

u = u_start.copy()
snapshots = [[Link]()]
times = [0.0]

for step in range(1, num_steps + 1):

u_mid = [Link]()

# ── Half-step: implicit in x, explicit in y ──────────────────────

for j in range(1, num_rows - 1): # row by row
rhs = (ry * u[j-1, 1:-1]
+ (1 - 2*ry) * u[j, 1:-1]
+ ry * u[j+1, 1:-1])
# Boundary contribution
rhs[0] += rx * u_mid[j, 0]
rhs[-1] += rx * u_mid[j, -1]
u_mid[j, 1:-1] = [Link](Ax, rhs)

apply_bc(u_mid)

u_next = u_mid.copy()

# ── Full-step: explicit in x, implicit in y ──────────────────────

for i in range(1, num_cols - 1): # column by column
rhs = (rx * u_mid[1:-1, i-1]
+ (1 - 2*rx) * u_mid[1:-1, i]
+ rx * u_mid[1:-1, i+1])
rhs[0] += ry * u_next[0, i]
rhs[-1] += ry * u_next[-1, i]
u_next[1:-1, i] = [Link](Ay, rhs)

apply_bc(u_next)
u = u_next

if step % save_every == 0:
[Link]([Link]())
[Link](step * dt)

return snapshots, times

# ─── Run CN solver (larger dt allowed!) ───────────────────────────────────

dt_cn = 10 * max_dt # 10× larger than explicit limit
num_steps_cn = int(run_time / dt_cn)
save_step_cn = max(1, num_steps_cn // 50)

print(f"Crank–Nicolson ADI: dt = {dt_cn:.1f} s | steps = {num_steps_cn:,}")

t0 = [Link]()
frames_cn, times_cn = solve_crank_nicolson_adi(
u_start, alpha, dx, dy, dt_cn, num_steps_cn, save_step_cn)
print(f"Wall time: {[Link]()-t0:.2f} s | Snapshots saved: {len(frames_cn)}")

Crank–Nicolson ADI: dt = 10.0 s | steps = 500

Wall time: 1.91 s | Snapshots saved: 51

In [8]: # ─── Plot CN snapshots ─────────────────────────────────────────────────────

frame_picks = [0, len(frames_cn)//4, len(frames_cn)//2, -1]

[Link] 8/22
23/05/2026, 14:59 IAM_Programming_Project

fig, axes = [Link](1, 4, figsize=(18, 4))

for ax, idx in zip(axes, frame_picks):

im = [Link](X, Y, frames_cn[idx], levels=25, cmap='hot', vmin=vmin, vmax=v
[Link](im, ax=ax, label='T (°C)')
ax.set_title(f't = {times_cn[idx]:.0f} s')
ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')

[Link]('Crank–Nicolson ADI — Temperature Evolution (Hot Spot IC)', fontsize=1

plt.tight_layout()
[Link]()

6. Method comparison
Final temperature fields, centre-point temperature history, and a note on how many steps
each method needed to reach the same physical time.

In [9]: # ─── Side-by-side final state comparison ──────────────────────────────────

fig, axes = [Link](1, 3, figsize=(16, 5))

final_exp = frames_exp[-1]
final_cn = frames_cn[-1]
diff = final_exp - final_cn

im0 = axes[0].contourf(X, Y, final_exp, levels=25, cmap='hot')

axes[0].set_title(f'Explicit FTCS\n(dt={dt_explicit:.1f}s, {num_steps:,} steps)')
[Link](im0, ax=axes[0], label='T (°C)')

im1 = axes[1].contourf(X, Y, final_cn, levels=25, cmap='hot')

axes[1].set_title(f'Crank–Nicolson ADI\n(dt={dt_cn:.1f}s, {num_steps_cn:,} steps)')
[Link](im1, ax=axes[1], label='T (°C)')

im2 = axes[2].contourf(X, Y, diff, levels=25, cmap='coolwarm')

axes[2].set_title(f'Difference (Exp − CN)\nMax |diff| = {[Link](diff).max():.3f} °C
[Link](im2, ax=axes[2], label='ΔT (°C)')

for ax in axes:
ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')

[Link](f'Final State Comparison at t = {run_time:.0f} s', fontsize=13, fontwe

plt.tight_layout()

[Link] 9/22
23/05/2026, 14:59 IAM_Programming_Project

[Link]()

print(f"Max absolute difference between methods: {[Link](diff).max():.4f} °C")

print(f"RMS difference: {[Link]([Link](diff**2)):.4f} °C")

Max absolute difference between methods: 0.0000 °C

RMS difference: 0.0000 °C

In [10]: # ─── Centre-point temperature vs time ─────────────────────────────────────

mid_col = Nx // 2
mid_row = Ny // 2

center_temps_exp = [s[mid_row, mid_col] for s in frames_exp]

center_temps_cn = [s[mid_row, mid_col] for s in frames_cn]

fig, (ax1, ax2) = [Link](1, 2, figsize=(14, 5))

[Link](times_exp, center_temps_exp, 'b-', lw=1.5, label='Explicit FTCS')

[Link](times_cn, center_temps_cn, 'r--', lw=2, label='Crank–Nicolson')
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Temperature (°C)')
ax1.set_title('Centre-Point Temperature vs Time')
[Link]()
[Link](True, ls='--', alpha=0.5)

# Temperature profiles along centreline at final time

[Link](x, final_exp[mid_row, :], 'b-', lw=2, label='Explicit FTCS')
[Link](x, final_cn [mid_row, :], 'r--', lw=2, label='Crank–Nicolson')
ax2.set_xlabel('x (m)')
ax2.set_ylabel('Temperature (°C)')
ax2.set_title(f'Centreline Profile at t = {run_time:.0f} s (y = {y[mid_row]:.2f} m
[Link]()
[Link](True, ls='--', alpha=0.5)

plt.tight_layout()
[Link]()

[Link] 10/22
23/05/2026, 14:59 IAM_Programming_Project

7. Verification against the analytical solution

For a sinusoidal initial condition $u_0 = A\sin(\pi x/L)\sin(\pi y/L)$ with zero Dirichlet BCs, the
exact solution is:

$$u(x, y, t) = A\sin\!\left(\frac{\pi x}{L}\right)\sin\!\left(\frac{\pi y}

{L}\right)\exp\!\left(-2\alpha\pi^2 t / L^2\right)$$

This gives actual error numbers rather than just a relative comparison between the two
schemes.

In [11]: def exact_solution(X, Y, t, alpha, Lx, Ly, A=100.0):

"""Exact solution for sinusoidal IC with zero Dirichlet BCs."""
decay = [Link](-alpha * [Link]**2 * (1/Lx**2 + 1/Ly**2) * t)
return A * [Link]([Link] * X / Lx) * [Link]([Link] * Y / Ly) * decay

u_start_sine = make_initial_condition(Nx, Ny, x, y, start_type='sine')

# Run solvers with sine IC

verify_time = 2000.0
steps_verify = int(verify_time / dt_explicit)
frames_verify_exp, times_verify_exp = solve_explicit(
u_start_sine, alpha, dx, dy, dt_explicit, steps_verify, save_every=max(1, steps

steps_verify_cn = int(verify_time / dt_cn)

frames_verify_cn, times_verify_cn = solve_crank_nicolson_adi(
u_start_sine, alpha, dx, dy, dt_cn, steps_verify_cn, save_every=max(1, steps_ve

# Compute L2 error over time

errors_exp, errors_cn = [], []
t_arr_exp = [Link](times_verify_exp)
t_arr_cn = [Link](times_verify_cn)

for snap, t in zip(frames_verify_exp, times_verify_exp):

exact = exact_solution(X, Y, t, alpha, Lx, Ly)
errors_exp.append([Link]([Link]((snap - exact)**2)))

[Link] 11/22
23/05/2026, 14:59 IAM_Programming_Project

for snap, t in zip(frames_verify_cn, times_verify_cn):

exact = exact_solution(X, Y, t, alpha, Lx, Ly)
errors_cn.append([Link]([Link]((snap - exact)**2)))

# Plot
fig, axes = [Link](1, 2, figsize=(14, 5))

axes[0].semilogy(times_verify_exp, errors_exp, 'b-', lw=2, label='Explicit FTCS')

axes[0].semilogy(times_verify_cn, errors_cn, 'r--', lw=2, label='Crank–Nicolson')
axes[0].set_xlabel('Time (s)')
axes[0].set_ylabel('L2 Error (°C)')
axes[0].set_title('L2 Error vs Analytical Solution')
axes[0].legend()
axes[0].grid(True, which='both', ls='--', alpha=0.5)

# Compare fields at t = verify_time

exact_final = exact_solution(X, Y, verify_time, alpha, Lx, Ly)
diff_exp_exact = frames_verify_exp[-1] - exact_final
im = axes[1].contourf(X, Y, diff_exp_exact, levels=20, cmap='coolwarm')
[Link](im, ax=axes[1], label='Error (°C)')
axes[1].set_title(f'Explicit Error Field at t={verify_time:.0f}s\nMax={[Link](diff_
axes[1].set_xlabel('x (m)')
axes[1].set_ylabel('y (m)')

plt.tight_layout()
[Link]()

print(f"Final L2 error — Explicit: {errors_exp[-1]:.4e} °C | CN: {errors_cn[-1]:.4e

Final L2 error — Explicit: 3.6882e+01 °C | CN: 3.6891e+01 °C

8. Instability demonstration
What happens when you exceed the stability limit? The explicit scheme diverges fast. Both a
stable and an unstable case are shown below so the difference is clear.

In [12]: def solve_explicit_early_stop(u_start, alpha, dx, dy, dt, num_steps, max_val=1e6):

"""Run explicit solver; stop if solution blows up."""
rx = alpha * dt / dx**2

[Link] 12/22
23/05/2026, 14:59 IAM_Programming_Project

ry = alpha * dt / dy**2
u = u_start.copy()
max_temps = [[Link]()]
times = [0.0]

for step in range(1, num_steps + 1):

u_next = [Link]()
u_next[1:-1, 1:-1] = (
u[1:-1, 1:-1]
+ rx * (u[1:-1, 2:] - 2*u[1:-1, 1:-1] + u[1:-1, :-2])
+ ry * (u[2:, 1:-1] - 2*u[1:-1, 1:-1] + u[:-2, 1:-1])
)
apply_bc(u_next)
u = u_next
max_temps.append([Link]([Link](u)))
[Link](step * dt)
if [Link]([Link](u)) > max_val or [Link]([Link](u)):
print(f" Blow-up detected at step {step} (t={step*dt:.2f}s)")
break

return u, times, max_temps

print("Running unstable explicit scheme...")

stability_check(alpha, dt_bad, dx, dy)
_, times_unstable, max_temp_unstable = solve_explicit_early_stop(
u_start, alpha, dx, dy, dt_bad, num_steps=300)

print("\nRunning stable explicit scheme for comparison...")

_, times_stable, max_temp_stable = solve_explicit_early_stop(
u_start, alpha, dx, dy, dt_safe, num_steps=300)

fig, ax = [Link](figsize=(10, 5))

[Link](times_unstable, max_temp_unstable, 'r-', lw=2, label=f'UNSTABLE dt={d
[Link](times_stable, max_temp_stable, 'b-', lw=2, label=f'Stable dt={dt_sa
[Link](200, ls='--', color='gray', label='Initial max T')
ax.set_xlabel('Time (s)')
ax.set_ylabel('Max |T| (°C)')
ax.set_title('Stability Demonstration: Explicit FTCS')
[Link]()
[Link](True, which='both', ls='--', alpha=0.5)
plt.tight_layout()
[Link]()

Running unstable explicit scheme...

rx = α·Δt/Δx² = 0.3750
ry = α·Δt/Δy² = 0.3750
r = rx + ry = 0.7500 (must be ≤ 0.5 for stability)
Stable: ✗ NO — UNSTABLE!
Blow-up detected at step 28 (t=42.00s)

Running stable explicit scheme for comparison...

[Link] 13/22
23/05/2026, 14:59 IAM_Programming_Project

9. Parameter studies
How does thermal diffusivity affect how fast heat spreads, and does the spatial error
converge at the expected $O(\Delta x^2)$ rate? The first two studies answer both.

A third test pushes the Crank–Nicolson time step well past the explicit limit to see where
accuracy starts to degrade.

In [13]: # ─── Study 1: Effect of thermal diffusivity α ─────────────────────────────

# Materials: copper, steel, glass, air (illustrative values)

materials = {
'Copper': 1.17e-4,
'Stainless Steel': 4.0e-6,
'Glass': 3.4e-7,
'Concrete': 7.0e-7,
}

study_time = 8000.0
fig, axes = [Link](2, 4, figsize=(20, 9))

for col, (mat, a) in enumerate([Link]()):

dt_mat = 0.9 * 0.5 / (a * (1/dx**2 + 1/dy**2))
steps_mat = int(study_time / dt_mat)
save_step_mat = max(1, steps_mat // 5)
u_start_mat = make_initial_condition(Nx, Ny, x, y, start_type='hotspot')
frames_mat, times_mat = solve_explicit(u_start_mat, a, dx, dy, dt_mat, steps_ma

# Top row: initial state

im0 = axes[0, col].contourf(X, Y, frames_mat[0], levels=20, cmap='hot')
[Link](im0, ax=axes[0, col])
axes[0, col].set_title(f'{mat}\nα={a:.1e} m²/s\nt=0 s')
axes[0, col].set_xlabel('x'); axes[0, col].set_ylabel('y')

[Link] 14/22
23/05/2026, 14:59 IAM_Programming_Project

# Bottom row: final state

im1 = axes[1, col].contourf(X, Y, frames_mat[-1], levels=20, cmap='hot')
[Link](im1, ax=axes[1, col])
axes[1, col].set_title(f't = {study_time:.0f} s')
axes[1, col].set_xlabel('x'); axes[1, col].set_ylabel('y')

[Link]('Effect of Thermal Diffusivity α on Heat Diffusion\n(Hot Spot IC, same

fontsize=13, fontweight='bold')
plt.tight_layout()
[Link]()

In [14]: # ─── Study 2: Grid Convergence (vary Nx = Ny) ─────────────────────────────

# Compare centreline profiles for different grid resolutions

grid_sizes = [11, 21, 41, 81]

conv_time = 1000.0

fig, (ax1, ax2) = [Link](1, 2, figsize=(14, 5))

colors = [Link]([Link](0, 0.9, len(grid_sizes)))

grid_errors = []
spacings = []

for grid_n, col in zip(grid_sizes, colors):

step_x = Lx / (grid_n - 1)
step_y = Ly / (grid_n - 1)
xs = [Link](0, Lx, grid_n)
ys = [Link](0, Ly, grid_n)
XX, YY = [Link](xs, ys)
dt_grid = 0.9 * 0.5 / (alpha * (1/step_x**2 + 1/step_y**2))
steps_grid = int(conv_time / dt_grid)

u_start_grid = 100 * [Link]([Link] * XX / Lx) * [Link]([Link] * YY / Ly)

u_start_grid[0,:] = u_start_grid[-1,:] = u_start_grid[:,0] = u_start_grid[:,-1]

frames_grid, _ = solve_explicit(u_start_grid, alpha, step_x, step_y, dt_grid, s

u_numerical = frames_grid[-1]
u_exact = exact_solution(XX, YY, conv_time, alpha, Lx, Ly)

[Link] 15/22
23/05/2026, 14:59 IAM_Programming_Project

error = [Link]([Link]((u_numerical - u_exact)**2))

grid_errors.append(error)
[Link](step_x)

cx = grid_n // 2
[Link](xs, u_numerical[cx, :], color=col, lw=1.5, label=f'N={grid_n}, Δx={ste

[Link](xs, exact_solution(XX, YY, conv_time, alpha, Lx, Ly)[grid_n//2, :],

'k--', lw=2, label='Analytical')
ax1.set_xlabel('x (m)')
ax1.set_ylabel('T (°C)')
ax1.set_title(f'Centreline Profiles at t={conv_time:.0f}s')
[Link](fontsize=9)
[Link](True, ls='--', alpha=0.5)

# Log-log convergence plot

[Link](spacings, grid_errors, 'bo-', ms=8, lw=2, label='Explicit FTCS')
# Reference slope O(dx^2)
ref_x = [Link](spacings)
[Link](ref_x, grid_errors[0] * (ref_x/ref_x[0])**2,
'k--', label='O(Δx²) reference')
ax2.set_xlabel('Δx (m)')
ax2.set_ylabel('L2 Error (°C)')
ax2.set_title('Grid Convergence Study')
[Link]()
[Link](True, which='both', ls='--', alpha=0.5)

plt.tight_layout()
[Link]()

print("Grid convergence:")
for grid_n, step_x, err in zip(grid_sizes, spacings, grid_errors):
print(f" N={grid_n:3d}, Δx={step_x:.4f}, L2 error={err:.4e}")

Grid convergence:
N= 11, Δx=0.1000, L2 error=3.7501e+01
N= 21, Δx=0.0500, L2 error=3.6511e+01
N= 41, Δx=0.0250, L2 error=3.5799e+01
N= 81, Δx=0.0125, L2 error=3.5392e+01

In [15]: # ─── Study 3: Time Step Sensitivity (CN, unconditionally stable) ──────────
dt_multipliers = [1, 5, 20, 100]
test_time = 3000.0

[Link] 16/22
23/05/2026, 14:59 IAM_Programming_Project

fig, axes = [Link](1, len(dt_multipliers), figsize=(18, 4))

for ax, factor in zip(axes, dt_multipliers):

dt_test = factor * max_dt
steps_test = int(test_time / dt_test)
frames_test, _ = solve_crank_nicolson_adi(
u_start, alpha, dx, dy, dt_test, steps_test, save_every=steps_test)
im = [Link](X, Y, frames_test[-1], levels=20, cmap='hot')
[Link](im, ax=ax, label='T (°C)')
ax.set_title(f'CN: dt = {factor}×max_dt\n({dt_test:.1f}s, {steps_test} steps)')
ax.set_xlabel('x'); ax.set_ylabel('y')

[Link](f'Crank–Nicolson: Effect of Time Step Size (t = {test_time:.0f}s)',

fontsize=13, fontweight='bold')
plt.tight_layout()
[Link]()

10. Advanced visualizations

10a. 3D surface and heat flux
In [16]: from mpl_toolkits.mplot3d import Axes3D

# ─── 3D surface at final time ─────────────────────────────────────────────

fig = [Link](figsize=(16, 6))

ax1 = fig.add_subplot(121, projection='3d')

surf = ax1.plot_surface(X, Y, frames_exp[-1], cmap='hot', edgecolor='none', alpha=0
[Link](surf, ax=ax1, shrink=0.5, label='T (°C)')
ax1.set_xlabel('x (m)')
ax1.set_ylabel('y (m)')
ax1.set_zlabel('T (°C)')
ax1.set_title(f'3D Temperature Surface\nt = {times_exp[-1]:.0f} s (Explicit)')
ax1.view_init(elev=30, azim=-60)

# ─── Heat flux / gradient vectors (streamlines) ───────────────────────────

ax2 = fig.add_subplot(122)
u_f = frames_exp[-1]
# Gradient (heat flux direction: negative gradient)
dTdy, dTdx = [Link](u_f, dy, dx)
speed = [Link](dTdx**2 + dTdy**2)
speed[speed == 0] = 1e-10

[Link] 17/22
23/05/2026, 14:59 IAM_Programming_Project

im = [Link](X, Y, u_f, levels=25, cmap='hot', alpha=0.8)

[Link](im, ax=ax2, label='T (°C)')
# Subsample for arrow plot
step_q = 4
[Link](X[::step_q, ::step_q], Y[::step_q, ::step_q],
-dTdx[::step_q, ::step_q]/speed[::step_q, ::step_q],
-dTdy[::step_q, ::step_q]/speed[::step_q, ::step_q],
color='cyan', alpha=0.6, scale=30, width=0.003)
ax2.set_xlabel('x (m)')
ax2.set_ylabel('y (m)')
ax2.set_title('Heat Flux Directions (arrows = −∇T)')

plt.tight_layout()
[Link]()

In [17]: # ─── 10b. Temperature History at Multiple Probe Points ─────────────────────

probes = {
'Centre (0.5, 0.5)': (Nx//2, Ny//2),
'Quarter (0.25,0.25)': (Nx//4, Ny//4),
'Near left (0.1, 0.5)': (int(0.1*(Nx-1)), Ny//2),
'Near top (0.5, 0.9)': (Nx//2, int(0.9*(Ny-1))),
}

fig, ax = [Link](figsize=(12, 5))

probe_colors = [Link].tab10([Link](0, 0.6, len(probes)))

for (label, (xi, yi)), col in zip([Link](), probe_colors):

probe_temps_exp = [s[yi, xi] for s in frames_exp]
probe_temps_cn = [s[yi, xi] for s in frames_cn]
[Link](times_exp, probe_temps_exp, '-', color=col, lw=1.5, label=f'{label} [E
[Link](times_cn, probe_temps_cn, '--', color=col, lw=1.5, alpha=0.7, label=f

ax.set_xlabel('Time (s)')
ax.set_ylabel('Temperature (°C)')
ax.set_title('Temperature History at Probe Points (solid=Explicit, dashed=CN)')
[Link](fontsize=8, ncol=2)
[Link](True, ls='--', alpha=0.5)
plt.tight_layout()
[Link]()

[Link] 18/22
23/05/2026, 14:59 IAM_Programming_Project

In [18]: # ─── 10c. Animation of heat diffusion ─────────────────────────────────────

# (renders inline if using %matplotlib notebook; saves to file otherwise)

fig_movie, ax_movie = [Link](figsize=(6, 5))

color_min, color_max = 0, 200

cf = ax_movie.contourf(X, Y, frames_exp[0], levels=25,

cmap='hot', vmin=color_min, vmax=color_max)
cb = [Link](cf, ax=ax_movie, label='T (°C)')
ax_movie.set_xlabel('x (m)')
ax_movie.set_ylabel('y (m)')
title = ax_movie.set_title(f't = 0 s')

def animate(i):
ax_movie.clear()
cf_ = ax_movie.contourf(X, Y, frames_exp[i], levels=25,
cmap='hot', vmin=color_min, vmax=color_max)
ax_movie.set_xlabel('x (m)')
ax_movie.set_ylabel('y (m)')
ax_movie.set_title(f'Explicit FTCS — t = {times_exp[i]:.0f} s')
return cf_.collections

anim = [Link](
fig_movie, animate,
frames=range(0, len(frames_exp), max(1, len(frames_exp)//40)),
interval=120, blit=False)

# Save as GIF
try:
[Link]('/home/claude/heat_diffusion.gif', writer='pillow', fps=8, dpi=80)
print("Animation saved: heat_diffusion.gif")
except Exception as e:
print(f"Could not save animation: {e}")

[Link](fig_movie)

# Static multi-frame summary instead

num_frames = 8
frame_indices = [Link](0, len(frames_exp)-1, num_frames, dtype=int)
fig, axes = [Link](2, 4, figsize=(18, 8))

[Link] 19/22
23/05/2026, 14:59 IAM_Programming_Project

for ax, fi in zip([Link], frame_indices):

im = [Link](X, Y, frames_exp[fi], levels=20,
cmap='hot', vmin=color_min, vmax=color_max)
[Link](im, ax=ax)
ax.set_title(f't = {times_exp[fi]:.0f} s')
ax.set_xlabel('x'); ax.set_ylabel('y')

[Link]('Heat Diffusion Evolution (Explicit FTCS, Hot Spot IC)', fontsize=13,

plt.tight_layout()
[Link]()

Could not save animation: [WinError 3] The system cannot find the path specified:
'C:\\home\\claude'

11. Summary
Property Explicit FTCS Crank–Nicolson (ADI)

Temporal
$O(\Delta t)$ $O(\Delta t^2)$
accuracy

$O(\Delta x^2, \Delta

Spatial accuracy $O(\Delta x^2, \Delta y^2)$
y^2)$

Stability Conditional: $r \leq 0.5$ Unconditional

$\Delta t_{\max} = 0.5/(\alpha(1/\Delta

Time step limit None
x^2+1/\Delta y^2))$

Cost per step O(N) vectorised O(N) tridiagonal solves

Implementation Simple Moderate

Long simulations, large

Best for Short runs, quick tests
$\Delta t$

In [19]: # ─── Final summary dashboard ───────────────────────────────────────────────

fig = [Link](figsize=(16, 10))

[Link] 20/22
23/05/2026, 14:59 IAM_Programming_Project

gs = GridSpec(2, 3, figure=fig, hspace=0.4, wspace=0.4)

# 1. Final explicit field

ax1 = fig.add_subplot(gs[0, 0])
im1 = [Link](X, Y, frames_exp[-1], levels=20, cmap='hot')
[Link](im1, ax=ax1)
ax1.set_title('Explicit — Final State')
ax1.set_xlabel('x'); ax1.set_ylabel('y')

# 2. Final CN field
ax2 = fig.add_subplot(gs[0, 1])
im2 = [Link](X, Y, frames_cn[-1], levels=20, cmap='hot')
[Link](im2, ax=ax2)
ax2.set_title('Crank–Nicolson — Final State')
ax2.set_xlabel('x'); ax2.set_ylabel('y')

# 3. Stability boundary
ax3 = fig.add_subplot(gs[0, 2])
alphas_ = [Link](-6, -2, 200)
dt_max_ = 0.5 / (alphas_ * (1/dx**2 + 1/dy**2))
[Link](alphas_, dt_max_, 'b-', lw=2)
ax3.fill_between(alphas_, dt_max_, 1e8, alpha=0.15, color='green', label='Stable')
ax3.fill_between(alphas_, 0, dt_max_, alpha=0.15, color='red', label='Unstable')
[Link]([alpha], [dt_safe], color='green', s=80, zorder=5)
ax3.set_xlabel('α (m²/s)'); ax3.set_ylabel('Δt (s)')
ax3.set_title('Stability Boundary'); [Link](); [Link](True, which='both', ls=

# 4. Centre-point T(t)
ax4 = fig.add_subplot(gs[1, 0])
[Link](times_exp, center_temps_exp, 'b-', lw=1.5, label='Explicit')
[Link](times_cn, center_temps_cn, 'r--', lw=2, label='CN')
ax4.set_xlabel('t (s)'); ax4.set_ylabel('T (°C)')
ax4.set_title('Centre Temperature vs Time')
[Link](); [Link](True, ls='--', alpha=0.4)

# 5. Grid convergence
ax5 = fig.add_subplot(gs[1, 1])
[Link](spacings, grid_errors, 'bo-', ms=8, lw=2)
[Link](ref_x, grid_errors[0]*(ref_x/ref_x[0])**2, 'k--', label='O(Δx²)')
ax5.set_xlabel('Δx (m)'); ax5.set_ylabel('L2 Error')
ax5.set_title('Grid Convergence'); [Link](); [Link](True, which='both', ls='-

# 6. Centreline profile comparison

ax6 = fig.add_subplot(gs[1, 2])
[Link](x, frames_exp[-1][Ny//2, :], 'b-', lw=2, label='Explicit')
[Link](x, frames_cn[-1] [Ny//2, :], 'r--', lw=2, label='CN')
ax6.set_xlabel('x (m)'); ax6.set_ylabel('T (°C)')
ax6.set_title('Centreline Profile (Final)')
[Link](); [Link](True, ls='--', alpha=0.4)

[Link]('2D Heat Equation — Summary Dashboard', fontsize=15, fontweight='bold'

[Link]()
print("\n✅ All sections complete.")

[Link] 21/22
23/05/2026, 14:59 IAM_Programming_Project

✅ All sections complete.

[Link] 22/22