Equation-Free Methods

Urban models often exist only as microscopic simulations (agent-based, cellular automata) with no closed-form macroscopic equations. The equation-free framework of Kevrekidis and collaborators enables macroscopic analysis—bifurcation diagrams, stability, projective integration—using only short bursts of the microscopic simulator as a “black box.”

1. The Lift-Evolve-Restrict Framework

The central idea is to bridge micro and macro scales through three operators:

Lift $\mathcal{L}$: Create a microscopic state consistent with a given macroscopic observable. For an urban CA, given a target density profile $\mathbf{U}$, generate a CA configuration whose coarse-grained density matches $\mathbf{U}$.
Evolve $\Psi^{\delta t}$: Run the microscopic simulator for a short time $\delta t$. This is a “black box”—we do not need to know the equations, just the code.
Restrict $\mathcal{M}$: Extract the macroscopic observable from the evolved microscopic state. For example, compute the density profile by averaging over CA cells in each coarse bin.

Together, these define a coarse timestepper:

$$\mathbf{U}^{n+1} = \mathcal{M}\circ\Psi^{\delta t}\circ\mathcal{L}(\mathbf{U}^n)$$

This gives us access to the macroscopic time derivative without knowing the macroscopic equation. The numerical derivative is:

$$\frac{d\mathbf{U}}{dt} \approx \frac{\mathcal{M}\circ\Psi^{\delta t}\circ\mathcal{L}(\mathbf{U}) - \mathbf{U}}{\delta t}$$

2. Coarse Projective Euler

The key speedup comes from projective integration: use a short microscopic burst to estimate the macroscopic time derivative, then take a large macroscopic Euler step:

$$\mathbf{U}^{n+1} = \mathbf{U}^n + \Delta T \cdot \frac{\mathcal{M}\circ\Psi^{\delta t}\circ\mathcal{L}(\mathbf{U}^n) - \mathbf{U}^n}{\delta t}$$

where:

$\delta t$ = micro step (short burst of CA simulation, e.g., 3 steps)
$\Delta T$ = macro step (projective leap, e.g., 30 steps)
Speedup factor: $\Delta T / \delta t$ (e.g., 10x)

The method is valid when the microscopic dynamics have a spectral gap: fast microscopic modes decay quickly during $\delta t$, leaving only the slow macroscopic evolution that we project forward.

When Does It Work?

The method succeeds when there is a separation of timescales. For urban CA models, the fast modes are local neighbourhood fluctuations (which equilibrate in a few steps), while the slow modes are large-scale density waves (which evolve over hundreds or thousands of steps). The ratio of these timescales determines the maximum speedup.

3. Stability and Accuracy

The projective Euler method has amplification factor:

$$G(\lambda) = 1 + \frac{\Delta T}{\delta t}\bigl(e^{\lambda\,\delta t} - 1\bigr)$$

where $\lambda$ is the eigenvalue of the macroscopic linear operator. For stability, we need $|G(\lambda)| \leq 1$. The key insight: fast modes (large negative $\lambda$) must have $e^{\lambda\,\delta t} \approx 0$, so $G \approx 1 - \Delta T/\delta t$. For stability:

$$\left|1 - \frac{\Delta T}{\delta t}\right| \leq 1 \qquad \Longrightarrow \qquad \Delta T \leq 2\,\delta t$$

This seems restrictive, but it applies only to fast modes. If we ensure these modes are sufficiently damped ($|\lambda_{\text{fast}}| \cdot \delta t \gg 1$), the stability constraint on the slow modes is simply:

$$\Delta T \leq \frac{2}{|\lambda_{\text{slow}}|}$$

which is exactly the CFL condition for the macroscopic equation we never wrote down. The accuracy is first-order in $\Delta T$ (Euler method), but can be improved to higher order using Runge-Kutta variants of the projective integrator.

4. Speedup Example

Consider an urban CA with:

Micro step: $\delta t = 3$ CA iterations (enough for local equilibration)
Macro step: $\Delta T = 30$ effective CA iterations
Total simulation: 3000 CA-equivalent time steps

Full simulation: 3000 micro steps. Equation-free: 100 macro steps, each requiring 3 micro steps = 300 micro steps total. Speedup: $10\times$.

For more expensive microscopic models (e.g., agent-based models with complex interactions), the speedup can be even more dramatic—each micro step is expensive, so reducing their count by an order of magnitude provides substantial wall-clock savings.

5. Implementation: Lift-Evolve-Restrict for Urban CA

We implement the full equation-free framework for a 1D diffusion CA. The code compares coarse projective integration with full microscopic simulation, measuring both speedup and accuracy.

Equation-Free Coarse Projective Integration for Urban CA

Python

script.py168 lines

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# ====================================================
# Microscopic model: 1D diffusion CA (probabilistic)
# ====================================================
def micro_step(state, p_move=0.3):
    """One CA step: each occupied cell moves left/right with probability p_move."""
    N = len(state)
    new_state = state.copy().astype(float)
    # Diffusive averaging (deterministic version for clarity)
    new_state[1:-1] = (1-2*p_move)*state[1:-1] + p_move*(state[:-2] + state[2:])
    new_state[0] = state[0] + p_move*(state[1] - state[0])
    new_state[-1] = state[-1] + p_move*(state[-2] - state[-1])
    return new_state

def micro_evolve(state, n_steps):
    """Run microscopic simulator for n_steps."""
    s = state.copy()
    for _ in range(n_steps):
        s = micro_step(s)
    return s

# ====================================================
# Equation-Free operators
# ====================================================
N_micro = 200      # microscopic grid
N_macro = 20       # macroscopic bins
bin_size = N_micro // N_macro

def restrict(micro_state):
    """Restrict: micro -> macro (bin averaging)."""
    macro = np.zeros(N_macro)
    for i in range(N_macro):
        macro[i] = np.mean(micro_state[i*bin_size:(i+1)*bin_size])
    return macro

def lift(macro_state):
    """Lift: macro -> micro (piecewise constant interpolation)."""
    micro = np.zeros(N_micro)
    for i in range(N_macro):
        micro[i*bin_size:(i+1)*bin_size] = macro_state[i]
    return micro

# ====================================================
# Coarse Projective Integration
# ====================================================
def coarse_projective_step(U_macro, delta_t, Delta_T):
    """One coarse projective Euler step.
       delta_t: micro healing time (number of micro steps)
       Delta_T: macro projection step (in micro-step units)
    """
    # Lift to micro
    micro = lift(U_macro)
    # Evolve micro for delta_t steps
    micro_evolved = micro_evolve(micro, delta_t)
    # Restrict back to macro
    U_evolved = restrict(micro_evolved)
    # Estimate time derivative
    dU_dt = (U_evolved - U_macro) / delta_t
    # Project forward
    U_new = U_macro + Delta_T * dU_dt
    return U_new

# ====================================================
# Run both methods
# ====================================================
# Initial condition: Gaussian bump
x_micro = np.linspace(0, 1, N_micro)
x_macro = np.linspace(0, 1, N_macro)
sigma0 = 0.08
initial_micro = np.exp(-(x_micro - 0.5)**2 / (2*sigma0**2))
initial_macro = restrict(initial_micro)

# Parameters
delta_t = 3       # micro healing steps
Delta_T = 30      # macro projection step
total_T = 3000    # total simulation time (micro-step units)

# --- Full micro simulation ---
n_full_steps = total_T
state_full = initial_micro.copy()
snapshots_full = {0: restrict(state_full.copy())}
for step in range(1, n_full_steps + 1):
    state_full = micro_step(state_full)
    if step % 300 == 0:
        snapshots_full[step] = restrict(state_full.copy())
micro_count_full = n_full_steps

# --- Equation-free simulation ---
n_macro_steps = total_T // Delta_T
state_ef = initial_macro.copy()
snapshots_ef = {0: state_ef.copy()}
micro_count_ef = 0
for step in range(1, n_macro_steps + 1):
    state_ef = coarse_projective_step(state_ef, delta_t, Delta_T)
    micro_count_ef += delta_t  # only delta_t micro steps per macro step
    t_equiv = step * Delta_T
    if t_equiv % 300 == 0:
        snapshots_ef[t_equiv] = state_ef.copy()

# --- Compute error over time ---
state_full2 = initial_micro.copy()
state_ef2 = initial_macro.copy()
errors = []
times_err = []
for step in range(1, n_full_steps + 1):
    state_full2 = micro_step(state_full2)
    if step % Delta_T == 0:
        state_ef2 = coarse_projective_step(state_ef2, delta_t, Delta_T)
        err = np.sqrt(np.mean((restrict(state_full2) - state_ef2)**2))
        errors.append(err)
        times_err.append(step)

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

# (a) Profiles: full micro
ax = axes[0]
colors = plt.cm.YlGn(np.linspace(0.3, 0.95, len(snapshots_full)))
for idx, (t, prof) in enumerate(sorted(snapshots_full.items())):
    ax.plot(x_macro, prof, color=colors[idx], linewidth=2, label=f't={t}')
ax.set_xlabel('x')
ax.set_ylabel('Density')
ax.set_title('Full Micro Simulation')
ax.legend(fontsize=7, ncol=2)

# (b) Profiles: equation-free
ax = axes[1]
colors2 = plt.cm.YlGn(np.linspace(0.3, 0.95, len(snapshots_ef)))
for idx, (t, prof) in enumerate(sorted(snapshots_ef.items())):
    ax.plot(x_macro, prof, '--', color=colors2[idx], linewidth=2, label=f't={t}')
ax.set_xlabel('x')
ax.set_ylabel('Density')
ax.set_title('Equation-Free (CPI)')
ax.legend(fontsize=7, ncol=2)

# (c) Error and speedup
ax = axes[2]
ax.semilogy(times_err, errors, '-', color='teal', linewidth=2)
ax.set_xlabel('Time (micro-step units)')
ax.set_ylabel('L2 Error (eq-free vs full)')
ax.set_title('Accuracy of CPI')

# Speedup annotation
speedup = micro_count_full / max(micro_count_ef, 1)
ax.text(0.05, 0.95, f'Micro steps (full): {micro_count_full}\n'
        f'Micro steps (CPI): {micro_count_ef}\n'
        f'Speedup: {speedup:.1f}x',
        transform=ax.transAxes, fontsize=10, verticalalignment='top',
        bbox=dict(boxstyle='round', facecolor='black', alpha=0.8, edgecolor='teal'))

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

# --- Summary ---
print("=== Equation-Free Results ===")
print(f"Total simulation time: {total_T} micro-steps")
print(f"Full micro: {micro_count_full} micro-step evaluations")
print(f"Eq-free CPI: {micro_count_ef} micro-step evaluations")
print(f"Speedup: {speedup:.1f}x")
print(f"Final L2 error: {errors[-1]:.6f}")
print(f"Max L2 error: {max(errors):.6f}")
print(f"\nParameters: delta_t={delta_t}, Delta_T={Delta_T}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

6. Beyond Projective Integration

The equation-free framework enables several powerful macroscopic analyses:

Coarse Newton-GMRES: Find steady states of the macroscopic equation by using the coarse timestepper inside a Newton iteration with GMRES as the linear solver. This finds urban equilibria without ever writing down the PDE.
Coarse bifurcation analysis: Track how steady states change as parameters (e.g., growth rate, planning stringency) vary. Detect tipping points and phase transitions in urban dynamics.
Coarse stability analysis: Compute leading eigenvalues of the Jacobian via matrix-free methods (Arnoldi iteration using the coarse timestepper). Determine whether an urban steady state is stable or unstable.
Projective Runge-Kutta: Higher-order projective integrators for improved accuracy at the same speedup cost.

Key Insight

The equation-free framework is the ultimate multi-scale tool: it lets us performmacroscopic mathematical analysis (stability, bifurcation, optimization) on models that exist only as microscopic simulators. For urban science, where agent-based models and CA simulators are often the most natural representation, this bridges the gap between simulation and mathematical understanding—without requiring the modeler to derive the macroscopic equations by hand.

← Renormalization Group Part 10: Network Theory →