Chapter 5: Fisher-KPP Traveling Waves

5.1 Urban Sprawl as a Reaction-Diffusion Wave

Urban sprawl does not happen everywhere simultaneously. It propagates outward from city centres as a traveling wave: a front of urbanisation advancing into undeveloped land. The mathematical framework for such waves is the Fisher-KPP equation, which couples logistic growth (reaction) with spatial diffusion:

$$\frac{\partial u}{\partial t} = D \nabla^2 u + r \, u(1 - u)$$

Here $u(\mathbf{x}, t) \in [0, 1]$ represents the fraction of developed land (or normalised population density),$D$ is the spatial diffusion coefficient (capturing outward migration and development spread), and$r$ is the local growth rate.

This equation was independently introduced by R.A. Fisher (1937) in population genetics and by Kolmogorov, Petrovskii & Piskunov (1937) in combustion theory. Its application to urban dynamics captures the essential physics of how cities expand.

5.2 Traveling Wave Reduction

We seek solutions of the form $u(x, t) = U(\xi)$where $\xi = x - ct$ is the co-moving coordinate and $c > 0$ is the wave speed. The boundary conditions are:

$$U(-\infty) = 1 \quad \text{(fully urbanised behind the front)}, \qquad U(+\infty) = 0 \quad \text{(undeveloped ahead)}$$

Substituting into the PDE, using$\partial u/\partial t = -cU'$ and$\partial^2 u/\partial x^2 = U''$:

$$D U'' + c U' + r U(1 - U) = 0$$

This is a second-order autonomous ODE. We convert it to a phase-plane system by defining$V = U'$:

$$\begin{cases} U' = V \\ V' = -\frac{c}{D}V - \frac{r}{D}U(1 - U) \end{cases}$$

5.3 Phase Plane Analysis and Minimum Speed

The fixed points of the phase-plane system are$(U, V) = (0, 0)$ and$(1, 0)$. The traveling wave corresponds to a heteroclinic orbit connecting$(1, 0)$ (urbanised) to$(0, 0)$ (undeveloped).

Linearisation at (0, 0)

The Jacobian at $(0, 0)$ is:

$$J = \begin{pmatrix} 0 & 1 \\ -r/D & -c/D \end{pmatrix}$$

The eigenvalues are:

$$\lambda_{\pm} = \frac{-c/D \pm \sqrt{c^2/D^2 - 4r/D}}{2} = \frac{-c \pm \sqrt{c^2 - 4Dr}}{2D}$$

For a physically meaningful wave (monotone front,$U \geq 0$ everywhere), the eigenvalues must bereal and negative, which requires:

$$c^2 - 4Dr \geq 0 \qquad \Longrightarrow \qquad c \geq 2\sqrt{Dr}$$

The celebrated result:

$$\boxed{c_{\min} = 2\sqrt{Dr}}$$

KPP proved that for compactly supported initial data, the asymptotic speed of the front is exactly$c_{\min}$. This is the speed of urban sprawl.

Urban Interpretation

The sprawl speed $c = 2\sqrt{Dr}$ increases with both diffusivity (ease of migration/development) and growth rate. A city with good transportation (high $D$) and strong economic growth (high$r$) sprawls faster. Doubling$D$ increases sprawl speed by$\sqrt{2} \approx 41\%$.

5.4 Numerical Method: Explicit Finite Differences

We discretise on a uniform grid with spacing$\Delta x$ and time step$\Delta t$. The diffusion term uses the central difference stencil:

$$u_j^{n+1} = u_j^n + \frac{D \Delta t}{\Delta x^2}\left(u_{j+1}^n - 2u_j^n + u_{j-1}^n\right) + \Delta t \, r \, u_j^n (1 - u_j^n)$$

CFL Stability Condition

The explicit scheme is conditionally stable. The CFL condition for the diffusion operator requires:

$$\Delta t \leq \frac{\Delta x^2}{2D}$$

The reaction term imposes an additional constraint$\Delta t \leq 1/r$, but in practice the diffusion CFL is more restrictive.

5.5 Python: Fisher-KPP Wave Propagation

We solve the Fisher-KPP equation using explicit finite differences and track the propagating front.

Fisher-KPP: Urban Sprawl Wave Propagation

Python

script.py86 lines

import numpy as np

# ── Fisher-KPP 1D solver with explicit finite differences ──

# Parameters
D = 1.0        # diffusion coefficient (km^2/year)
r = 0.1        # growth rate (1/year)
L_domain = 100.0  # domain length (km)
nx = 500       # spatial grid points
dx = L_domain / nx
dt_max = dx**2 / (2 * D)  # CFL limit
dt = 0.8 * dt_max  # safety factor
T_final = 80.0     # total simulation time (years)
n_steps = int(T_final / dt)
c_theory = 2 * np.sqrt(D * r)

print("=" * 60)
print("  Fisher-KPP Traveling Wave Simulation")
print("=" * 60)
print(f"  D = {D} km^2/yr,  r = {r} /yr")
print(f"  dx = {dx:.3f} km,  dt = {dt:.4f} yr")
print(f"  CFL limit: dt <= {dt_max:.4f}")
print(f"  Theoretical wave speed: c_min = 2*sqrt(D*r) = {c_theory:.4f} km/yr")
print()

# Initial condition: step function with smooth transition
x = np.linspace(0, L_domain, nx)
u = np.where(x < 10, 1.0, np.exp(-(x - 10)**2 / 4))
u = np.minimum(u, 1.0)

# Store front positions for speed measurement
front_times = []
front_positions = []

# Time integration
snapshot_interval = n_steps // 5
print(f"  {'Time (yr)':>10} {'Front pos (km)':>16} {'Front speed':>14} {'Max u':>8} {'Min u':>8}")
print("-" * 60)

for step in range(n_steps + 1):
    t = step * dt

# Track front position (u = 0.5 contour)
    above = np.where(u > 0.5)[0]
    if len(above) > 0:
        front_pos = x[above[-1]]
    else:
        front_pos = 0.0

if step % snapshot_interval == 0:
        if len(front_times) > 0 and (t - front_times[-1]) > 0:
            speed = (front_pos - front_positions[-1]) / (t - front_times[-1])
        else:
            speed = 0.0
        print(f"  {t:10.1f} {front_pos:16.2f} {speed:14.4f} {np.max(u):8.4f} {np.min(u):8.4f}")
        front_times.append(t)
        front_positions.append(front_pos)

if step < n_steps:
        # Explicit update: diffusion + reaction
        u_new = u.copy()
        u_new[1:-1] = (u[1:-1]
                        + D * dt / dx**2 * (u[2:] - 2*u[1:-1] + u[:-2])
                        + dt * r * u[1:-1] * (1 - u[1:-1]))
        # Boundary conditions
        u_new[0] = 1.0   # urbanised at left
        u_new[-1] = 0.0  # undeveloped at right
        u = u_new

print()

# Measure asymptotic speed from last two front positions
if len(front_positions) >= 2:
    c_measured = (front_positions[-1] - front_positions[-2]) / (front_times[-1] - front_times[-2])
    print(f"  Measured asymptotic speed: {c_measured:.4f} km/yr")
    print(f"  Theoretical c_min:        {c_theory:.4f} km/yr")
    print(f"  Relative error:           {abs(c_measured - c_theory)/c_theory * 100:.2f}%")
print()

# Print wave profile at final time
print("  Final wave profile (every 10th point near front):")
print(f"  {'x (km)':>10} {'u':>10}")
front_region = np.where((x > front_pos - 20) & (x < front_pos + 20))[0]
for idx in front_region[::5]:
    bar = '#' * int(u[idx] * 40)
    print(f"  {x[idx]:10.2f} {u[idx]:10.4f}  {bar}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

5.6 Fortran: Explicit Upwind Scheme

A Fortran implementation with explicit CFL checking and front tracking.

Fortran: Fisher-KPP Explicit Solver

Fortran

program.f9074 lines

program fisher_kpp
  implicit none
  integer, parameter :: nx = 400
  real(8), parameter :: D = 1.0d0, r = 0.1d0
  real(8), parameter :: L_dom = 80.0d0
  real(8) :: dx, dt, dt_max, c_theory
  real(8) :: u(nx), u_new(nx), x(nx)
  real(8) :: t, T_final, front_pos, front_old, speed
  integer :: i, step, n_steps, print_every

dx = L_dom / dble(nx - 1)
  dt_max = dx*dx / (2.0d0 * D)
  dt = 0.8d0 * dt_max
  T_final = 60.0d0
  n_steps = int(T_final / dt)
  c_theory = 2.0d0 * sqrt(D * r)
  print_every = n_steps / 6

write(*,'(A)')       '================================================='
  write(*,'(A)')       '  Fortran Fisher-KPP Solver'
  write(*,'(A)')       '================================================='
  write(*,'(A,F8.4)')  '  dx       = ', dx
  write(*,'(A,F8.5)')  '  dt       = ', dt
  write(*,'(A,F8.5)')  '  CFL max  = ', dt_max
  write(*,'(A,F8.4)')  '  c_theory = ', c_theory
  write(*,'(A,I8)')    '  n_steps  = ', n_steps
  write(*,'(A)')       '-------------------------------------------------'

! Grid and initial condition
  do i = 1, nx
    x(i) = dble(i-1) * dx
    if (x(i) < 8.0d0) then
      u(i) = 1.0d0
    else
      u(i) = exp(-(x(i) - 8.0d0)**2 / 4.0d0)
    end if
    if (u(i) > 1.0d0) u(i) = 1.0d0
  end do

front_old = 8.0d0
  write(*,'(A)') '   Step     Time     FrontPos    Speed'
  write(*,'(A)') '-------------------------------------------------'

do step = 1, n_steps
    t = dble(step) * dt

! Explicit update
    u_new(1) = 1.0d0
    u_new(nx) = 0.0d0
    do i = 2, nx-1
      u_new(i) = u(i) &
        + D*dt/(dx*dx) * (u(i+1) - 2.0d0*u(i) + u(i-1)) &
        + dt * r * u(i) * (1.0d0 - u(i))
    end do
    u = u_new

if (mod(step, print_every) == 0) then
      ! Find front (u = 0.5)
      front_pos = 0.0d0
      do i = 1, nx
        if (u(i) > 0.5d0) front_pos = x(i)
      end do
      speed = (front_pos - front_old) / (dble(print_every) * dt)
      front_old = front_pos
      write(*,'(I8, F10.2, F12.3, F10.4)') step, t, front_pos, speed
    end if
  end do

write(*,'(A)')       '-------------------------------------------------'
  write(*,'(A,F10.4)') '  Theoretical speed: ', c_theory
  write(*,'(A)')       '  Front converges to c_min = 2*sqrt(D*r)'
  write(*,'(A)')       '================================================='

end program fisher_kpp

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

5.7 Summary & Key Takeaways

• Fisher-KPP: $\partial_t u = D\nabla^2 u + ru(1-u)$ — reaction-diffusion model of urban sprawl
• Traveling wave ansatz reduces PDE to ODE: $DU'' + cU' + rU(1-U) = 0$
• Minimum wave speed: $c_{\min} = 2\sqrt{Dr}$
• Phase-plane analysis: monotone front requires real eigenvalues at $(0,0)$
• CFL stability: $\Delta t \leq \Delta x^2 / (2D)$ for the explicit scheme
• Urban interpretation: sprawl speed depends on diffusivity (transportation) and growth rate (economics)

← Logistic ODE Next: Instanton Formalism →