Homogenization: From CA to PDE

This chapter builds the rigorous bridge between the cellular automaton models of Module 3 and the PDE continuum of Module 2. Through Taylor expansion and two-scale homogenization, we show that a discrete CA neighbourhood average converges to the diffusion equation—and we quantify how heterogeneity in urban fabric affects the effective diffusion coefficient.

1. Taylor Expansion of the CA Update

Consider a CA on a square lattice with spacing $\varepsilon$. The state at cell $\mathbf{x}$ is updated by a weighted average over neighbours:

$$\rho^{n+1}(\mathbf{x}) = \frac{\sum_k w_k\,\rho^n(\mathbf{x} + \varepsilon\,\boldsymbol{\delta}_k)}{\sum_k w_k}$$

where $\boldsymbol{\delta}_k$ are the neighbour offsets (including the cell itself with $\boldsymbol{\delta}_0 = \mathbf{0}$) and $w_k$ are the weights. Now Taylor-expand each term:

$$\rho(\mathbf{x} + \varepsilon\boldsymbol{\delta}_k) = \rho + \varepsilon\,\boldsymbol{\delta}_k \cdot \nabla\rho + \frac{\varepsilon^2}{2}\,\boldsymbol{\delta}_k^T\,\nabla^2\!\rho\;\boldsymbol{\delta}_k + O(\varepsilon^3)$$

Summing over all neighbours with weights, the first-order terms vanish by symmetry (equal weights in opposite directions). The leading correction is second-order:

$$\rho^{n+1} = \rho^n + \frac{\varepsilon^2}{2}\,\frac{\sum_k w_k\,|\boldsymbol{\delta}_k|^2}{\sum_k w_k}\,\nabla^2\rho + O(\varepsilon^4)$$

Identifying the time step $\Delta t = 1$ and defining the effective diffusion coefficient:

$$D = \frac{\varepsilon^2}{2}\,\frac{\sum_k w_k\,|\boldsymbol{\delta}_k|^2}{\sum_k w_k}$$

For a standard 5-point stencil (von Neumann neighbourhood with uniform weights $w_k = 1$ for the 4 neighbours and $w_0 = 0$ for the centre), each neighbour has $|\boldsymbol{\delta}_k|^2 = 1$, so:

$$D = \frac{\varepsilon^2}{2}\cdot\frac{4 \times 1}{4} = \frac{\varepsilon^2}{2}$$

For the “spread” CA where $w_0 = 1$ (include self) with 4 neighbours: $D = \varepsilon^2 \cdot 4/(2 \cdot 5) = 2\varepsilon^2/5$. More generally, $D = \varepsilon^2 \cdot \text{spread} / 4$ for the standard neighbourhood with equal weights, recovering the diffusion equation $\partial_t \rho = D\nabla^2\rho$.

2. Two-Scale Expansion and the Cell Problem

When the diffusion coefficient varies periodically at the fine scale—representing heterogeneous urban fabric (residential vs. commercial vs. industrial zones)—we introduce two spatial variables:

Slow variable: $\mathbf{x}$ (city scale)
Fast variable: $\mathbf{y} = \mathbf{x}/\varepsilon$ (neighbourhood/block scale)

The diffusion equation with periodic coefficients is:

$$\frac{\partial \rho}{\partial t} = \nabla \cdot \bigl[D(\mathbf{x}/\varepsilon)\,\nabla\rho\bigr]$$

Expanding $\rho = \rho_0(\mathbf{x},t) + \varepsilon\,\rho_1(\mathbf{x},\mathbf{y},t) + \varepsilon^2\,\rho_2 + \cdots$ and collecting powers of $\varepsilon$, the leading-order problem yields the cell problem:

$$\nabla_{\mathbf{y}} \cdot \bigl[D(\mathbf{y})\,(\nabla_{\mathbf{y}}\chi_j + \mathbf{e}_j)\bigr] = 0$$

where $\chi_j(\mathbf{y})$ is the corrector(periodic in $\mathbf{y}$). The effective (homogenised) diffusion tensor is:

$$D_{ij}^* = \frac{1}{|Y|}\int_Y D(\mathbf{y})\,\bigl(\delta_{ij} + \partial_{y_i}\chi_j\bigr)\,d\mathbf{y}$$

A fundamental bound from homogenization theory:

$$D^* \leq \langle D \rangle = \frac{1}{|Y|}\int_Y D(\mathbf{y})\,d\mathbf{y}$$

Heterogeneity always retards effective diffusion (and hence urban sprawl). The effective diffusivity is bounded above by the arithmetic mean and below by the harmonic mean. For a city with alternating high- and low-permeability zones, sprawl is slower than a naive average would predict.

3. Convergence: CA Resolution Study

We verify that the CA converges to the PDE as $\varepsilon \to 0$ by running simulations at 4 resolutions and measuring the error against the analytic diffusion solution. The convergence rate should be $O(\varepsilon^2)$.

CA-to-PDE Convergence and Homogenization

Python

script.py151 lines

import numpy as np
import matplotlib.pyplot as plt

# --- Analytic solution for 1D diffusion ---
def analytic_gaussian(x, t, D, sigma0=1.0):
    """Gaussian spreading: P(x,t) = (1/sqrt(2pi*s2)) exp(-x^2/(2*s2))
       where s2 = sigma0^2 + 2Dt."""
    s2 = sigma0**2 + 2*D*t
    return np.exp(-x**2 / (2*s2)) / np.sqrt(2*np.pi*s2)

# --- CA simulation (1D, weighted average) ---
def run_ca(N, n_steps, L=20.0, sigma0=1.0):
    """Run 1D CA diffusion and return final profile."""
    eps = 2*L / N
    x = np.linspace(-L, L, N)

# Initial Gaussian
    rho = np.exp(-x**2 / (2*sigma0**2))
    rho /= np.sum(rho) * eps  # normalise

# Weights: self=0, two neighbours=1 each (simple average)
    # D_eff = eps^2/2 * (2*1)/(2) = eps^2/2
    D_eff = eps**2 / 2.0

for step in range(n_steps):
        rho_new = np.zeros_like(rho)
        # Interior: average of two neighbours
        rho_new[1:-1] = 0.5 * (rho[:-2] + rho[2:])
        # Boundaries: zero flux
        rho_new[0] = rho[1]
        rho_new[-1] = rho[-2]
        rho = rho_new

t_final = n_steps  # dt = 1
    return x, rho, D_eff, t_final

# --- Resolution study ---
resolutions = [64, 128, 256, 512]
L = 20.0
sigma0 = 1.0

# We want the same physical time for all resolutions
# D_eff = eps^2/2, so for N cells: eps = 2L/N, D = 2L^2/N^2
# To reach same physical diffusion time T_phys, need n_steps = T_phys / dt = T_phys
# But D differs, so we fix T_phys via D*t = const
# Actually: use enough steps so the Gaussian spreads noticeably

# Fix: use the same *physical* end time by adjusting steps
# t_phys = n_steps * 1 (dt=1 for CA)
# D_eff * t_phys should be the same for fair comparison
# Let's fix D*t = 5 for all resolutions

target_Dt = 5.0  # D_eff * t_final = 5

errors = []
eps_vals = []

fig, axes = plt.subplots(1, 3, figsize=(15, 5))

colors = ['#065f46', '#059669', '#34d399', '#a7f3d0']

for idx, N in enumerate(resolutions):
    eps = 2*L / N
    D_eff = eps**2 / 2.0
    n_steps = int(target_Dt / D_eff)

x, rho_ca, _, t_final = run_ca(N, n_steps, L, sigma0)
    rho_exact = analytic_gaussian(x, D_eff * n_steps, D_eff, sigma0)

# L2 error
    err = np.sqrt(np.sum((rho_ca - rho_exact)**2) * eps)
    errors.append(err)
    eps_vals.append(eps)

# Plot profiles for selected resolutions
    if N in [64, 512]:
        ax = axes[0]
        ax.plot(x, rho_ca, color=colors[idx], linewidth=1.5,
                label=f'CA N={N}')
        if N == 512:
            ax.plot(x, rho_exact, 'k--', linewidth=1, label='PDE exact')

eps_vals = np.array(eps_vals)
errors = np.array(errors)

# Fit convergence rate
log_eps = np.log(eps_vals)
log_err = np.log(errors)
slope, intercept = np.polyfit(log_eps, log_err, 1)

# (a) Profile comparison
axes[0].set_xlabel('x')
axes[0].set_ylabel('rho(x)')
axes[0].set_title('CA vs PDE Solution')
axes[0].legend(fontsize=9)

# (b) Convergence plot
ax = axes[1]
ax.loglog(eps_vals, errors, 'o-', color='teal', markersize=8, linewidth=2)
eps_fit = np.linspace(eps_vals.min()*0.8, eps_vals.max()*1.2, 50)
ax.loglog(eps_fit, np.exp(intercept)*eps_fit**slope, '--', color='crimson',
          linewidth=2, label=f'slope = {slope:.2f}')
ax.loglog(eps_fit, 0.5*eps_fit**2, ':', color='gray', alpha=0.5, label='O(eps^2) ref')
ax.set_xlabel('epsilon (grid spacing)')
ax.set_ylabel('L2 error')
ax.set_title('Convergence Rate')
ax.legend(fontsize=10)

# (c) Effective D vs arithmetic mean for heterogeneous case
# Demonstrate D* < <D> with alternating coefficients
n_cells = 200
y = np.linspace(0, 1, n_cells)
D_local = np.where(np.sin(2*np.pi*4*y) > 0, 2.0, 0.5)  # alternating zones

# Solve cell problem: d/dy [D(y) (d chi/dy + 1)] = 0
# => D(y)(d chi/dy + 1) = const = C
# => d chi/dy = C/D(y) - 1
# Periodicity: integral of d chi/dy over period = 0
# => C * integral(1/D) = 1  => C = 1/<1/D> = D_harmonic

D_harmonic = 1.0 / np.mean(1.0 / D_local)
D_arithmetic = np.mean(D_local)

D_star_values = np.linspace(D_harmonic * 0.9, D_arithmetic * 1.1, 5)
labels_d = ['D_harm', 'D*', 'D_arith']
vals_d = [D_harmonic, D_harmonic, D_arithmetic]  # D* = D_harmonic in 1D

ax = axes[2]
ax.fill_between(y, D_local*0.3, color='teal', alpha=0.3)
ax.plot(y, D_local*0.3, color='teal', linewidth=0.5)
ax.axhline(D_harmonic*0.3, color='crimson', linestyle='--', linewidth=2,
           label=f'D* = D_harm = {D_harmonic:.3f}')
ax.axhline(D_arithmetic*0.3, color='gold', linestyle='--', linewidth=2,
           label=f'<D> = {D_arithmetic:.3f}')
ax.set_xlabel('y (fast variable)')
ax.set_ylabel('D(y) (scaled)')
ax.set_title('Heterogeneous D: D* < <D>')
ax.legend(fontsize=9)

plt.tight_layout()
plt.savefig('output.png', dpi=130, bbox_inches='tight')
plt.show()

print(f"Convergence rate: O(eps^{slope:.2f})")
print(f"Expected: O(eps^2)")
print(f"\nHeterogeneous medium (1D):")
print(f"  Arithmetic mean <D> = {D_arithmetic:.4f}")
print(f"  Harmonic mean (= D* in 1D) = {D_harmonic:.4f}")
print(f"  Ratio D*/<D> = {D_harmonic/D_arithmetic:.4f}")
print(f"  Heterogeneity retards sprawl by {(1-D_harmonic/D_arithmetic)*100:.1f}%")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

4. Fortran: Gauss-Seidel for the Cell Problem

In 2D, the cell problem $\nabla_y \cdot [D(y)(\nabla_y \chi + \mathbf{e}_j)] = 0$ must be solved numerically. Gauss-Seidel iteration is the standard approach for this elliptic PDE.

program cell_problem
  ! Solve 2D cell problem for homogenization
  ! nabla_y . [D(y) (nabla_y chi + e_1)] = 0
  ! on unit cell Y = [0,1]^2 with periodic BC
  implicit none
  integer, parameter :: N = 64
  real(8) :: chi(N,N), D_coeff(N,N), rhs(N,N)
  real(8) :: dy, residual, tol, Dstar
  integer :: i, j, ip, im, jp, jm, iter, max_iter

  dy = 1.d0 / N
  tol = 1.d-8
  max_iter = 50000

  ! Heterogeneous D(y): checkerboard pattern
  do i = 1, N
    do j = 1, N
      if (mod((i-1)/(N/4) + (j-1)/(N/4), 2) == 0) then
        D_coeff(i,j) = 2.0d0   ! high-permeability zone
      else
        D_coeff(i,j) = 0.5d0   ! low-permeability zone
      end if
    end do
  end do

  ! RHS = -d/dy_1 [D(y)]  (for e_1 direction)
  do i = 1, N
    do j = 1, N
      ip = mod(i, N) + 1;  im = mod(i-2+N, N) + 1
      rhs(i,j) = -(D_coeff(ip,j) - D_coeff(im,j)) / (2.d0*dy)
    end do
  end do

  ! Gauss-Seidel iteration
  chi = 0.d0
  do iter = 1, max_iter
    residual = 0.d0
    do i = 1, N
      do j = 1, N
        ip = mod(i, N) + 1;  im = mod(i-2+N, N) + 1
        jp = mod(j, N) + 1;  jm = mod(j-2+N, N) + 1

        ! D averaged at half-points
        real(8) :: De, Dw, Dn, Ds, chi_new
        De = 0.5d0*(D_coeff(i,j) + D_coeff(ip,j))
        Dw = 0.5d0*(D_coeff(i,j) + D_coeff(im,j))
        Dn = 0.5d0*(D_coeff(i,j) + D_coeff(i,jp))
        Ds = 0.5d0*(D_coeff(i,j) + D_coeff(i,jm))

        chi_new = (De*chi(ip,j) + Dw*chi(im,j) + &
                   Dn*chi(i,jp) + Ds*chi(i,jm) + &
                   dy*dy*rhs(i,j)) / (De + Dw + Dn + Ds)

        residual = residual + (chi_new - chi(i,j))**2
        chi(i,j) = chi_new
      end do
    end do
    residual = sqrt(residual / (N*N))
    if (residual < tol) exit
  end do

  ! Compute D* = <D(1 + d chi/dy_1)>
  Dstar = 0.d0
  do i = 1, N
    do j = 1, N
      ip = mod(i, N) + 1; im = mod(i-2+N, N) + 1
      Dstar = Dstar + D_coeff(i,j) * &
              (1.d0 + (chi(ip,j) - chi(im,j))/(2.d0*dy))
    end do
  end do
  Dstar = Dstar / (N*N)

  write(*,'(A,I6,A,E12.4)') 'Converged in ', iter, &
        ' iterations, residual = ', residual
  write(*,'(A,F10.6)') 'Effective D* = ', Dstar
  write(*,'(A,F10.6)') 'Arithmetic <D> = ', sum(D_coeff)/(N*N)
end program cell_problem

The Bridge

This chapter completes the loop: Module 3's cellular automaton is shown to be a discrete approximation of Module 2's diffusion PDE, with an effective diffusion coefficient that depends on the urban microstructure. Homogenization theory provides the exact relationship, and the convergence rate $O(\varepsilon^2)$ confirms the accuracy of the continuum approximation even at coarse resolutions.

← Part 9 Overview Renormalization Group →