Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

By the end of this section, you will be able to:

Derive finite difference approximations for the heat equation using Taylor series expansion
Implement the explicit FTCS (Forward Time, Centered Space) method and understand its conditional stability
Apply von Neumann stability analysis to determine the CFL condition $r = \alpha\Delta t/\Delta x^2 \leq 0.5$
Construct the implicit method and solve the resulting tridiagonal system using the Thomas algorithm
Compare explicit, implicit, and Crank-Nicolson methods in terms of stability, accuracy, and computational cost
Analyze truncation error and understand the order of accuracy of different schemes
Connect finite difference methods to neural network architectures and physics-informed learning

The Big Picture: When Analytical Solutions Fail

"The purpose of computing is insight, not numbers." — Richard Hamming

In previous sections, we solved the heat equation analytically using Fourier series. This works beautifully for simple geometries and initial conditions. But what about:

Complex domains: Irregular shapes, holes, or boundaries that don't align with coordinate axes
Variable coefficients: Materials with temperature-dependent conductivity
Nonlinear problems: Radiation boundary conditions, phase changes
High dimensions: 3D problems with complicated geometry

For these problems, we turn to numerical methods. The most intuitive approach is finite differences: replace derivatives with algebraic differences on a discrete grid.

Engineering Simulations

Heat transfer in engines, electronics cooling, building thermal analysis with realistic 3D geometry

Weather & Climate

Atmospheric models, ocean circulation, climate prediction all rely on finite difference schemes

Financial PDEs

Option pricing via Black-Scholes, interest rate models — all solved numerically in practice

Neural Network Layers

Diffusion in latent space, physics-informed neural networks (PINNs) discretize PDEs on grids

Image Processing

Anisotropic diffusion for denoising, inpainting, and edge-preserving smoothing

Material Science

Diffusion in alloys, semiconductor doping, crystal growth simulations

Historical Context

The finite difference method has a rich history stretching back to the earliest days of scientific computing:

Lewis Fry Richardson (1910)

Attempted to predict weather by hand-calculating finite differences. His "forecast factory" vision inspired modern numerical weather prediction. The dimensionless number r in our stability condition is sometimes called the "mesh ratio" or "Fourier number."

John Crank & Phyllis Nicolson (1947)

Developed the Crank-Nicolson scheme while studying heat conduction in uranium during the Manhattan Project. Their method remains the gold standard for parabolic PDEs due to its unconditional stability and second-order accuracy.

Courant-Friedrichs-Lewy (1928)

The famous CFL condition (actually for hyperbolic equations) set the foundation for understanding numerical stability. The idea that information cannot travel faster than one grid cell per time step is central to all time-stepping schemes.

Discretizing Space and Time

To solve the heat equation numerically, we replace the continuous domain with a grid of discrete points:

The Discrete Grid

Spatial Grid

x_i = i \cdot \Delta x, \quad i = 0, 1, \ldots, N_x

where $\Delta x = L/N_x$

Temporal Grid

t^n = n \cdot \Delta t, \quad n = 0, 1, \ldots, N_t

where $\Delta t = T/N_t$

The unknown $u(x, t)$ becomes a grid function $u_i^n \approx u(x_i, t^n)$

Finite Difference Grid and Stencil

Show Stencil

Method:

Explicit (FTCS) Stencil

Uses known values at time n to compute the unknown value at time n+1.

Future (unknown)CenterNeighbors

Grid Parameters

Spatial step: dx = 0.2

Time step: dt = 0.05

Grid ratio: r = dt/dx² = 1.250

Selected point: (i=3, n=2)

Tip: Click on interior grid points to move the stencil and see how the finite difference formula applies at different locations.

Finite Difference Approximations

Using Taylor series, we can approximate derivatives with differences:

Derivative	Approximation	Error Order
Forward time	(u_i^{n+1} - u_i^n) / dt	O(dt)
Backward time	(u_i^n - u_i^{n-1}) / dt	O(dt)
Centered space (1st)	(u_{i+1}^n - u_{i-1}^n) / (2dx)	O(dx^2)
Centered space (2nd)	(u_{i+1}^n - 2u_i^n + u_{i-1}^n) / dx^2	O(dx^2)

The centered second difference is derived from:

u(x+dx) = u + u'dx + (1/2)u''dx² + (1/6)u'''dx³ + ...

u(x-dx) = u - u'dx + (1/2)u''dx² - (1/6)u'''dx³ + ...

Adding: u(x+dx) + u(x-dx) = 2u + u''dx² + O(dx⁴)

⟹ u'' = [u(x+dx) - 2u(x) + u(x-dx)] / dx² + O(dx²)

The FTCS Explicit Method

The simplest finite difference scheme for the heat equation is FTCS: Forward Time, Centered Space.

The FTCS Scheme

\frac{u_i^{n+1} - u_i^n}{\Delta t} = \alpha \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2}

Rearranging for the unknown

u_i^{n+1}

u_i^{n+1} = u_i^n + r(u_{i+1}^n - 2u_i^n + u_{i-1}^n)

where

r = \alpha\Delta t/\Delta x^2

is the mesh ratio (or Fourier number)

Interpretation: Weighted Average

We can rewrite the FTCS formula as:

u_i^{n+1} = r \cdot u_{i-1}^n + (1-2r) \cdot u_i^n + r \cdot u_{i+1}^n

This is a weighted average of the current point and its neighbors! The weights sum to 1 (conservation). But notice: if $r > 0.5$ , the center weight $(1-2r)$ becomes negative, which leads to instability.

Explicit FTCS Method Demonstration

Spatial Step (dx)0.05

Time Step (dt)0.0005

Diffusivity (alpha)1.0

Initial Condition

Stability Parameter: r = alpha * dt / dx² = 0.2000

STABLE (r <= 0.5)

The FTCS Formula

u_iⁿ⁺¹ = u_iⁿ + r(u_i+1ⁿ - 2u_iⁿ + u_i-1ⁿ)

Forward difference in time, centered difference in space. Simple but conditionally stable.

Stability Condition (CFL)

r = alpha * dt / dx² <= 0.5

If violated, numerical errors grow exponentially, causing the "blow-up" you see when r > 0.5.

The Explicit Method is Conditionally Stable

The FTCS method only works when $r = \alpha\Delta t/\Delta x^2 \leq 0.5$ . This means:

\Delta t \leq \frac{\Delta x^2}{2\alpha}

Halving the spatial grid spacing requires quartering the time step! This makes explicit methods expensive for high-resolution simulations.

Stability Analysis: Von Neumann Method

Why does the explicit method blow up when $r > 0.5$ ? The answer comes from von Neumann stability analysis, which examines how individual Fourier modes evolve.

The Key Idea

Any numerical error can be decomposed into Fourier modes. We analyze a single mode:

u_i^n = G^n e^{ik\theta} \quad \text{where } \theta = \text{wave number}

The amplification factor G determines whether the mode grows (|G| > 1, unstable) or decays (|G| < 1, stable).

Amplification Factor for FTCS

Substituting the Fourier mode into the FTCS scheme and simplifying:

G = 1 + r(e^iθ - 2 + e^-iθ)

G = 1 + r(2cos(θ) - 2)

G = 1 - 4r sin²(θ/2)

For stability, we need $|G| \leq 1$ for all wave numbers θ ∈ [0, π]:

At θ = π (highest frequency), G = 1 - 4r. For |G| ≤ 1, need 1 - 4r ≥ -1, giving r ≤ 0.5.
At θ = 0 (lowest frequency), G = 1 (no decay) — always stable.
High-frequency modes (large θ) decay fastest when stable, grow fastest when unstable.

Von Neumann Stability Analysis

The amplification factor G(theta) determines how each Fourier mode evolves. For stability, we need |G| <= 1 for all theta in [0, pi].

Grid Ratio r = alpha * dt / dx²:0.25Stable

Show all methods

Stability Region for Explicit Method

Explicit (FTCS)

G = 1 - 4r sin²(theta/2)

Stable only for r <= 0.5. High frequencies (large theta) decay fastest.

Implicit

G = 1 / (1 + 4r sin²(theta/2))

Unconditionally stable (|G| < 1 for all r). Requires solving linear system.

Crank-Nicolson

G = (1-2r sin²)/(1+2r sin²)

Unconditionally stable AND second-order accurate in time. Best of both worlds!

The Implicit Method

To avoid the stability restriction, we can evaluate the spatial derivative at the new time level n+1:

The Fully Implicit Scheme

\frac{u_i^{n+1} - u_i^n}{\Delta t} = \alpha \frac{u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1}}{\Delta x^2}

Rearranging:

-r u_{i-1}^{n+1} + (1+2r) u_i^{n+1} - r u_{i+1}^{n+1} = u_i^n

The Tridiagonal System

The implicit scheme couples all unknowns at time n+1, creating a tridiagonal linear system:


┌                              ┐ ┌     ┐   ┌     ┐
│ (1+2r)   -r      0    ...   │ │u₁ⁿ⁺¹│   │u₁ⁿ  │
│  -r    (1+2r)   -r    ...   │ │u₂ⁿ⁺¹│ = │u₂ⁿ  │
│   0     -r    (1+2r)  ...   │ │u₃ⁿ⁺¹│   │u₃ⁿ  │
│  ...    ...    ...    ...   │ │ ... │   │ ... │
└                              ┘ └     ┘   └     ┘

This system is solved efficiently using the Thomas algorithm (tridiagonal matrix algorithm), which runs in O(n) operations.

Why Is Implicit Unconditionally Stable?

The amplification factor for the implicit method is:

G = \frac{1}{1 + 4r\sin^2(\theta/2)}

Since the denominator is always > 1 for any r > 0, we have |G| < 1 for all wave numbers and all time steps. The method is unconditionally stable!

Trade-offs

Pro: Can use large time steps (limited only by accuracy, not stability)
Con: Must solve a linear system at each time step
Con: Only first-order accurate in time (same as explicit)

The Crank-Nicolson Method

The Crank-Nicolson method averages the explicit and implicit schemes, using the spatial derivative at both time levels:

The Crank-Nicolson Scheme

\frac{u_i^{n+1} - u_i^n}{\Delta t} = \frac{\alpha}{2}\left[\frac{\delta^2 u_i^{n+1}}{\Delta x^2} + \frac{\delta^2 u_i^n}{\Delta x^2}\right]

where

\delta^2 u_i = u_{i+1} - 2u_i + u_{i-1}

Rearranging into matrix form:

-\frac{r}{2}u_{i-1}^{n+1} + (1+r)u_i^{n+1} - \frac{r}{2}u_{i+1}^{n+1} = \frac{r}{2}u_{i-1}^n + (1-r)u_i^n + \frac{r}{2}u_{i+1}^n

Why Crank-Nicolson Is Special

Property	Explicit	Implicit	Crank-Nicolson
Stability	r <= 0.5	Unconditional	Unconditional
Time accuracy	O(dt)	O(dt)	O(dt^2)
Space accuracy	O(dx^2)	O(dx^2)	O(dx^2)
Linear system?	No	Yes	Yes
Overall accuracy	O(dt, dx^2)	O(dt, dx^2)	O(dt^2, dx^2)

The amplification factor for CN is:

G = \frac{1 - 2r\sin^2(\theta/2)}{1 + 2r\sin^2(\theta/2)}

Since both numerator and denominator are positive for reasonable r, |G| < 1 always. CN is unconditionally stable ANDsecond-order accurate in time.

The Best of Both Worlds

Crank-Nicolson achieves the stability of implicit methods with better accuracy than either explicit or implicit. It's the default choice for parabolic PDEs in most applications.

Method Comparison

Let's compare all three methods on the same problem: the heat equation with initial condition $u(x,0) = \sin(\pi x)$ and homogeneous boundary conditions.

Method Comparison: Explicit vs Implicit vs Crank-Nicolson

explicitimplicitCrank-Nicolson

Parameters: dx = 0.05, dt = 0.001, r = 0.4000, t = 0, step = 0

Explicit Error

0.000e+0

O(dt) + O(dx²)

Implicit Error

0.000e+0

O(dt) + O(dx²)

Crank-Nicolson Error

0.000e+0

O(dt²) + O(dx²)

Observation: Crank-Nicolson typically achieves the lowest error because it's second-order accurate in both space and time, while explicit and implicit methods are only first-order in time.

When to Use Each Method

Explicit (FTCS)

Simple to implement
Good for quick prototypes
Fine spatial resolution needed
Short simulation times

Implicit

Long-time simulations
Stiff problems
When accuracy is less critical
Steady-state convergence

Crank-Nicolson

Most practical applications
Best accuracy-stability trade-off
Standard choice for parabolic PDEs
Option pricing (Black-Scholes)

Truncation Error Analysis

The truncation error is the error introduced by replacing exact derivatives with finite difference approximations. It determines how fast the numerical solution converges to the true solution as we refine the grid.

Truncation Error and Convergence Analysis

The truncation error is the error from approximating derivatives with finite differences. As dx decreases, the error should decrease at a rate determined by the approximation order.

Finite Difference Formulas

Forward: f'(x) ≈ [f(x+dx) - f(x)] / dxO(dx)

Centered 1st: f'(x) ≈ [f(x+dx) - f(x-dx)] / (2dx)O(dx²)

Centered 2nd: f''(x) ≈ [f(x+dx) - 2f(x) + f(x-dx)] / dx²O(dx²)

Key Insight: Order of Accuracy

An O(dx^p) method means the error scales as dx^p:

• Halve dx → O(dx) error halves (1/2)
• Halve dx → O(dx^2) error quarters (1/4)

Higher order = faster convergence = less computation needed for same accuracy.

Why Centered Differences Are Better

Centered differences use symmetric information from both sides of a point. Taylor expansion shows the odd-order error terms cancel, leaving only even-order terms. This is why the centered second derivative formula is the standard choice for the heat equation: it achieves O(dx²) accuracy with minimal extra computation.

Consistency, Stability, and Convergence

The Lax Equivalence Theorem states that for a consistent finite difference scheme, stability is equivalent to convergence:

The Lax Equivalence Theorem

Consistency + Stability = Convergence

Consistency

Truncation error → 0 as dx, dt → 0

Stability

Errors don't grow unboundedly

Convergence

Numerical solution → exact solution

Handling Boundary Conditions

Different boundary conditions require different numerical treatment:

Dirichlet (Prescribed Value)

If $u(0,t) = g_L(t)$ , simply set $u_0^n = g_L(t^n)$ . The boundary value appears on the right-hand side of the matrix system.

Neumann (Prescribed Flux)

If $\partial u/\partial x|_{x=0} = h(t)$ , use a one-sided or centered difference:

Centered (2nd order): (u₁ - u₋₁)/(2dx) = h → u₋₁ = u₁ - 2*dx*h

Introduces a "ghost point" u₋₁ that we eliminate using the PDE at i=0

Robin (Mixed)

A combination: $\alpha u + \beta \partial u/\partial x = g$ . Discretize the derivative and solve for the boundary value or ghost point.

Machine Learning Connections

Finite difference methods have deep connections to modern machine learning:

1. Convolutional Neural Networks

The finite difference stencil $[1, -2, 1]$ for the second derivative is a convolution kernel. CNN filters learn similar spatial operations from data.

The FTCS Update as Convolution

u[n+1] = conv1d(u[n], kernel=[r, 1-2r, r])

One time step of the explicit method is equivalent to applying a 1D convolution with a specific kernel determined by r.

2. Physics-Informed Neural Networks (PINNs)

PINNs embed the PDE into the neural network loss function. The network learns to approximate the solution while automatically satisfying the heat equation:

L_physics = || du/dt - alpha * d²u/dx² ||²

The spatial and temporal derivatives are computed using automatic differentiation, not finite differences!

3. Neural Operators

Modern architectures like Fourier Neural Operators (FNO) learn the solution operator directly, mapping initial conditions to solutions at any time. They generalize the Green's function approach we saw in analytical solutions.

4. Diffusion Models in Generative AI

Image generation models (Stable Diffusion, DALL-E) simulate a discrete diffusion process. Each denoising step is like running the heat equation backward in time — the neural network learns the score function to reverse diffusion.

Python Implementation

Let's implement the explicit and implicit methods in Python:

Explicit and Implicit Finite Difference Methods

🐍heat_finite_differences.py

Explanation(5)

Code(159)

3The Explicit FTCS Function

FTCS = Forward Time, Centered Space. This is the simplest finite difference method: march forward in time using known values to compute new values.

27Stability Check

The stability parameter r = alpha*dt/dx^2 must satisfy r <= 0.5 for the explicit method. Violating this causes exponential error growth.

39The FTCS Update Formula

The core formula: u_new[i] = u_old[i] + r*(u_old[i+1] - 2*u_old[i] + u_old[i-1]). This is a weighted average of the point and its neighbors.

46Implicit Method

The implicit method uses values at the NEW time level on the right-hand side, requiring us to solve a linear system. The payoff: unconditional stability.

72Thomas Algorithm

Tridiagonal systems can be solved in O(n) operations using the Thomas algorithm (a specialized form of Gaussian elimination). This makes implicit methods efficient.

154 lines without explanation

1import numpy as np
2import matplotlib.pyplot as plt
3
4def solve_heat_explicit(L, T, nx, nt, alpha, initial_func, bc_left=0, bc_right=0):
5    """
6    Solve the 1D heat equation using the explicit FTCS method.
7
8    Parameters:
9    -----------
10    L : float - Length of the rod
11    T : float - Total simulation time
12    nx : int - Number of spatial grid points
13    nt : int - Number of time steps
14    alpha : float - Thermal diffusivity
15    initial_func : callable - Initial temperature distribution
16    bc_left, bc_right : float - Boundary conditions (Dirichlet)
17
18    Returns:
19    --------
20    x : array - Spatial grid
21    t : array - Time grid
22    u : 2D array - Solution u[time_index, space_index]
23    """
24    # Create grids
25    dx = L / (nx - 1)
26    dt = T / (nt - 1)
27    x = np.linspace(0, L, nx)
28    t = np.linspace(0, T, nt)
29
30    # Stability parameter
31    r = alpha * dt / dx**2
32    print(f"Grid ratio r = alpha*dt/dx^2 = {r:.4f}")
33    if r > 0.5:
34        print("WARNING: r > 0.5 - method may be unstable!")
35
36    # Initialize solution array
37    u = np.zeros((nt, nx))
38
39    # Set initial condition
40    u[0, :] = initial_func(x)
41
42    # Apply boundary conditions
43    u[:, 0] = bc_left
44    u[:, -1] = bc_right
45
46    # Time-stepping loop (explicit FTCS scheme)
47    for n in range(0, nt - 1):
48        for i in range(1, nx - 1):
49            # The FTCS formula
50            u[n+1, i] = u[n, i] + r * (u[n, i+1] - 2*u[n, i] + u[n, i-1])
51
52    return x, t, u
53
54
55def solve_heat_implicit(L, T, nx, nt, alpha, initial_func, bc_left=0, bc_right=0):
56    """
57    Solve the 1D heat equation using the fully implicit method.
58
59    The implicit method is unconditionally stable, allowing larger time steps.
60    Requires solving a tridiagonal system at each time step.
61    """
62    dx = L / (nx - 1)
63    dt = T / (nt - 1)
64    x = np.linspace(0, L, nx)
65    t = np.linspace(0, T, nt)
66
67    r = alpha * dt / dx**2
68    print(f"Grid ratio r = {r:.4f} (implicit is unconditionally stable)")
69
70    # Initialize
71    u = np.zeros((nt, nx))
72    u[0, :] = initial_func(x)
73    u[:, 0] = bc_left
74    u[:, -1] = bc_right
75
76    # Build the tridiagonal matrix coefficients
77    # System: -r*u_{i-1}^{n+1} + (1+2r)*u_i^{n+1} - r*u_{i+1}^{n+1} = u_i^n
78    n_interior = nx - 2
79
80    # Thomas algorithm coefficients
81    a = np.full(n_interior, -r)           # Lower diagonal
82    b = np.full(n_interior, 1 + 2*r)      # Main diagonal
83    c = np.full(n_interior, -r)           # Upper diagonal
84
85    # Time stepping
86    for n in range(0, nt - 1):
87        # Right-hand side (current time level values)
88        d = u[n, 1:-1].copy()
89        d[0] += r * bc_left        # Boundary contribution
90        d[-1] += r * bc_right
91
92        # Solve tridiagonal system using Thomas algorithm
93        u[n+1, 1:-1] = thomas_algorithm(a, b, c, d)
94
95    return x, t, u
96
97
98def thomas_algorithm(a, b, c, d):
99    """
100    Solve tridiagonal system using the Thomas algorithm (O(n)).
101
102    The system has form:
103    b[0]*x[0] + c[0]*x[1] = d[0]
104    a[i]*x[i-1] + b[i]*x[i] + c[i]*x[i+1] = d[i]  for i = 1..n-2
105    a[n-1]*x[n-2] + b[n-1]*x[n-1] = d[n-1]
106    """
107    n = len(d)
108    c_prime = np.zeros(n)
109    d_prime = np.zeros(n)
110    x = np.zeros(n)
111
112    # Forward sweep
113    c_prime[0] = c[0] / b[0]
114    d_prime[0] = d[0] / b[0]
115
116    for i in range(1, n):
117        denom = b[i] - a[i] * c_prime[i-1]
118        c_prime[i] = c[i] / denom if i < n-1 else 0
119        d_prime[i] = (d[i] - a[i] * d_prime[i-1]) / denom
120
121    # Back substitution
122    x[-1] = d_prime[-1]
123    for i in range(n-2, -1, -1):
124        x[i] = d_prime[i] - c_prime[i] * x[i+1]
125
126    return x
127
128
129# Example: Compare explicit and implicit methods
130L = 1.0          # Length
131T = 0.1          # Total time
132nx = 21          # Spatial points
133alpha = 1.0      # Diffusivity
134
135# Initial condition: sine wave
136initial = lambda x: np.sin(np.pi * x)
137
138# Analytical solution for comparison
139def analytical(x, t, alpha=1.0):
140    return np.exp(-np.pi**2 * alpha * t) * np.sin(np.pi * x)
141
142# Solve with both methods (using different time steps)
143nt_explicit = 501   # Many steps needed for stability
144nt_implicit = 51    # Fewer steps OK (unconditionally stable)
145
146x_exp, t_exp, u_exp = solve_heat_explicit(L, T, nx, nt_explicit, alpha, initial)
147x_imp, t_imp, u_imp = solve_heat_implicit(L, T, nx, nt_implicit, alpha, initial)
148
149# Plot comparison at final time
150plt.figure(figsize=(10, 5))
151plt.plot(x_exp, analytical(x_exp, T), 'k-', linewidth=2, label='Analytical')
152plt.plot(x_exp, u_exp[-1, :], 'b--', linewidth=2, label=f'Explicit (nt={nt_explicit})')
153plt.plot(x_imp, u_imp[-1, :], 'r-.', linewidth=2, label=f'Implicit (nt={nt_implicit})')
154plt.xlabel('Position x')
155plt.ylabel('Temperature u')
156plt.title(f'Heat Equation Solutions at t={T}')
157plt.legend()
158plt.grid(True, alpha=0.3)
159plt.show()

The Crank-Nicolson Implementation

Crank-Nicolson Method and Convergence Study

🐍crank_nicolson.py

Explanation(4)

Code(133)

3Crank-Nicolson Method

CN averages the explicit and implicit schemes, using spatial derivatives at both time n and n+1. This clever averaging achieves second-order accuracy in time.

34The CN Coefficients

The left-hand side has coefficients -r/2, (1+r), -r/2 while the right-hand side has r/2, (1-r), r/2. The half-r terms come from averaging the two time levels.

42Computing the Right-Hand Side

The explicit part: we apply the centered difference operator to the known values at time n, then add boundary condition contributions.

78Convergence Study

We refine the grid and measure how errors decrease. CN shows faster convergence because its O(dt^2) temporal accuracy matches the O(dx^2) spatial accuracy.

129 lines without explanation

1import numpy as np
2import matplotlib.pyplot as plt
3
4def solve_heat_crank_nicolson(L, T, nx, nt, alpha, initial_func, bc_left=0, bc_right=0):
5    """
6    Solve the 1D heat equation using the Crank-Nicolson method.
7
8    Crank-Nicolson averages the spatial derivatives at times n and n+1,
9    achieving:
10    - Unconditional stability (like implicit)
11    - Second-order accuracy in time (better than explicit/implicit)
12    - Second-order accuracy in space
13
14    The scheme:
15    -r/2 * u_{i-1}^{n+1} + (1+r) * u_i^{n+1} - r/2 * u_{i+1}^{n+1}
16    = r/2 * u_{i-1}^n + (1-r) * u_i^n + r/2 * u_{i+1}^n
17    """
18    dx = L / (nx - 1)
19    dt = T / (nt - 1)
20    x = np.linspace(0, L, nx)
21    t = np.linspace(0, T, nt)
22
23    r = alpha * dt / dx**2
24    print(f"Crank-Nicolson: r = {r:.4f}")
25    print("Note: Unconditionally stable AND O(dt^2) accurate!")
26
27    # Initialize
28    u = np.zeros((nt, nx))
29    u[0, :] = initial_func(x)
30    u[:, 0] = bc_left
31    u[:, -1] = bc_right
32
33    n_interior = nx - 2
34    half_r = r / 2
35
36    # Left-hand side matrix coefficients (implicit part)
37    a = np.full(n_interior, -half_r)      # Lower diagonal
38    b = np.full(n_interior, 1 + r)        # Main diagonal
39    c = np.full(n_interior, -half_r)      # Upper diagonal
40
41    # Time stepping
42    for n in range(0, nt - 1):
43        # Right-hand side: explicit part contribution
44        d = np.zeros(n_interior)
45        for i in range(n_interior):
46            j = i + 1  # Index in full array
47            d[i] = half_r * u[n, j-1] + (1 - r) * u[n, j] + half_r * u[n, j+1]
48
49        # Add boundary contributions
50        d[0] += half_r * bc_left
51        d[-1] += half_r * bc_right
52
53        # Solve tridiagonal system
54        u[n+1, 1:-1] = thomas_algorithm(a, b, c, d)
55
56    return x, t, u
57
58
59def thomas_algorithm(a, b, c, d):
60    """Thomas algorithm for tridiagonal systems."""
61    n = len(d)
62    c_prime = np.zeros(n)
63    d_prime = np.zeros(n)
64    x = np.zeros(n)
65
66    c_prime[0] = c[0] / b[0]
67    d_prime[0] = d[0] / b[0]
68
69    for i in range(1, n):
70        denom = b[i] - a[i] * c_prime[i-1]
71        c_prime[i] = c[i] / denom if i < n-1 else 0
72        d_prime[i] = (d[i] - a[i] * d_prime[i-1]) / denom
73
74    x[-1] = d_prime[-1]
75    for i in range(n-2, -1, -1):
76        x[i] = d_prime[i] - c_prime[i] * x[i+1]
77
78    return x
79
80
81# Compare error convergence
82def convergence_study():
83    """Compare truncation errors of all three methods."""
84    L, alpha = 1.0, 1.0
85    T = 0.05
86
87    initial = lambda x: np.sin(np.pi * x)
88    analytical = lambda x, t: np.exp(-np.pi**2 * alpha * t) * np.sin(np.pi * x)
89
90    nx_values = [11, 21, 41, 81, 161]
91    errors = {'explicit': [], 'implicit': [], 'crank-nicolson': []}
92    dx_values = []
93
94    for nx in nx_values:
95        dx = L / (nx - 1)
96        dx_values.append(dx)
97
98        # Choose dt based on stability for explicit method
99        dt_exp = 0.4 * dx**2 / alpha  # r = 0.4 (stable)
100        nt_exp = int(T / dt_exp) + 1
101
102        # For implicit/CN, can use larger dt
103        dt_imp = 2 * dx**2 / alpha  # r = 2 (implicit is stable)
104        nt_imp = int(T / dt_imp) + 1
105
106        x, t, u_exp = solve_heat_explicit(L, T, nx, nt_exp, alpha, initial)
107        _, _, u_imp = solve_heat_implicit(L, T, nx, nt_imp, alpha, initial)
108        _, _, u_cn = solve_heat_crank_nicolson(L, T, nx, nt_imp, alpha, initial)
109
110        u_exact = analytical(x, T)
111
112        errors['explicit'].append(np.max(np.abs(u_exp[-1, :] - u_exact)))
113        errors['implicit'].append(np.max(np.abs(u_imp[-1, :] - u_exact)))
114        errors['crank-nicolson'].append(np.max(np.abs(u_cn[-1, :] - u_exact)))
115
116    # Plot convergence
117    plt.figure(figsize=(8, 6))
118    plt.loglog(dx_values, errors['explicit'], 'b-o', label='Explicit (FTCS)')
119    plt.loglog(dx_values, errors['implicit'], 'g-s', label='Implicit')
120    plt.loglog(dx_values, errors['crank-nicolson'], 'r-^', label='Crank-Nicolson')
121
122    # Reference lines
123    dx_ref = np.array(dx_values)
124    plt.loglog(dx_ref, 0.5*dx_ref**2, 'k--', alpha=0.5, label='O(dx^2)')
125
126    plt.xlabel('Grid spacing dx')
127    plt.ylabel('Max error')
128    plt.title('Error Convergence Comparison')
129    plt.legend()
130    plt.grid(True, alpha=0.3)
131    plt.show()
132
133convergence_study()

Common Pitfalls

Forgetting the Stability Condition

The most common mistake is using an explicit method with too large a time step. Always check $r = \alpha\Delta t/\Delta x^2 \leq 0.5$ before running. Better yet, compute dt from dx automatically.

Wrong Index Handling

Off-by-one errors are extremely common in finite difference codes. Remember:

Boundary values (i = 0 and i = N) are fixed for Dirichlet conditions
Interior points are i = 1 to N-1
The matrix system has size (N-1) × (N-1), not N × N

Numerical Precision Issues

For very long simulations or very fine grids, floating-point errors can accumulate. Use double precision (float64) and consider conservation checks (total energy should be constant for insulated boundaries).

Debugging Strategy

Test with known solutions: Use $u(x,0) = \sin(\pi x)$ which has analytical solution $e^{-\pi^2\alpha t}\sin(\pi x)$
Check conservation: For insulated boundaries, total heat should be constant
Verify convergence rate: Error should decrease as expected when refining the grid
Visualize intermediate steps: Plot the solution at several time points to catch oscillations early

Test Your Understanding

Question 1 of 8

What does the stability condition r = alpha*dt/dx^2 <= 0.5 ensure for the explicit FTCS method?

Summary

Finite difference methods transform the continuous heat equation into a discrete system that can be solved on a computer. The key insight is the trade-off between stability, accuracy, and computational cost.

Key Formulas

Method	Update Formula	Stability	Accuracy
FTCS (Explicit)	u[n+1] = u[n] + r*(u[i+1] - 2u[i] + u[i-1])	r <= 0.5	O(dt, dx^2)
Implicit	Tridiagonal system with r coefficients	Unconditional	O(dt, dx^2)
Crank-Nicolson	Average of explicit and implicit	Unconditional	O(dt^2, dx^2)

Key Takeaways

Discretization replaces continuous derivatives with algebraic differences on a grid of points (x_i, t^n)
The explicit FTCS method is simple but conditionally stable: requires $r = \alpha\Delta t/\Delta x^2 \leq 0.5$
Von Neumann stability analysis examines how Fourier modes evolve; stability requires |G| ≤ 1 for all wave numbers
The implicit method is unconditionally stable but requires solving a tridiagonal system at each time step
Crank-Nicolson achieves the best of both worlds: unconditional stability AND second-order accuracy in time
Truncation error determines convergence rate; O(dx^2) means halving dx quarters the error
Boundary conditions (Dirichlet, Neumann, Robin) require careful discretization at the domain edges
Finite differences connect to machine learning: CNN kernels, physics-informed networks, and diffusion models all use related concepts

The Essence of Finite Differences:

"Replace derivatives with differences, continuous equations with discrete ones — but respect the physics through stability and accuracy."

Coming Next: In the next section, we'll explore Applications: Thermal Analysis — real-world engineering problems where finite difference methods shine.