Langevin Equation

The Langevin equation provides the foundational stochastic differential equation (SDE) for modelling noisy dynamics in urban systems. We derive it from microscopic jump processes, explore the Itô and Stratonovich interpretations, and simulate noise-induced transitions in bistable urban potentials.

1. From Master Equation to Langevin

Consider a continuous-time random walk on a lattice with spacing $\Delta x$. The probability $P(n, t)$ of being at site $n$ obeys the master equation:

$$\frac{\partial P(n,t)}{\partial t} = W_{n-1 \to n} P(n-1,t) + W_{n+1 \to n} P(n+1,t) - (W_{n \to n-1} + W_{n \to n+1}) P(n,t)$$

where $W_{m \to n}$ are transition rates. Setting $x = n \Delta x$ and performing a Kramers-Moyal expansion (Taylor expand the transition probabilities to second order in $\Delta x$):

$$\frac{\partial P(x,t)}{\partial t} = -\frac{\partial}{\partial x}\bigl[a_1(x) P\bigr] + \frac{1}{2}\frac{\partial^2}{\partial x^2}\bigl[a_2(x) P\bigr] + \mathcal{O}(\Delta x^3)$$

where the jump moments are:

$$a_k(x) = \lim_{\Delta t \to 0} \frac{1}{\Delta t} \langle (\Delta x)^k \rangle_x = \int (x' - x)^k \, W(x'|x) \, dx'$$

Truncating at second order gives the Fokker-Planck equation. The equivalent stochastic process at the individual trajectory level is the Langevin equation:

$$dx = a(x)\,dt + b(x)\,dW_t$$

where $a(x) = a_1(x)$ is the drift coefficient,$b(x) = \sqrt{a_2(x)}$ is the diffusion coefficient, and $W_t$ is the standard Wiener process with $\langle dW_t \rangle = 0$ and$\langle dW_t^2 \rangle = dt$.

2. Properties of the Wiener Process

The Wiener process $W_t$ (Brownian motion) has the following properties:

$W_0 = 0$
Independent increments: $W_t - W_s$ is independent of $\{W_u : u \leq s\}$ for $t > s$
Gaussian increments: $W_t - W_s \sim \mathcal{N}(0, t-s)$
Continuous but nowhere differentiable paths

The crucial property for stochastic calculus is the quadratic variation:

$$(dW_t)^2 = dt, \qquad dW_t \cdot dt = 0, \qquad (dt)^2 = 0$$

These "Itô calculus rules" are the foundation of all subsequent manipulations. The fact that $(dW)^2 = dt$ rather than zero (as in ordinary calculus) is what makes stochastic calculus fundamentally different.

3. Itô vs Stratonovich Interpretation

The SDE $dx = a\,dt + b\,dW$ is ambiguous when $b$ depends on $x$. The two standard interpretations differ in how the stochastic integral is evaluated:

Itô Interpretation

Evaluate the integrand at the left endpoint:

$$\int_0^T b(x_t)\,dW_t \approx \sum_k b(x_{t_k})(W_{t_{k+1}} - W_{t_k})$$

This gives a martingale: $\langle \int b \, dW \rangle = 0$.

Stratonovich Interpretation

Evaluate at the midpoint:

$$\int_0^T b(x_t) \circ dW_t \approx \sum_k b\!\left(\frac{x_{t_k} + x_{t_{k+1}}}{2}\right)(W_{t_{k+1}} - W_{t_k})$$

The ordinary chain rule applies, but the integral is no longer a martingale.

The two are related by a drift correction:

$$dx = a_S \, dt + b \circ dW_t = \underbrace{\left(a_S - \frac{1}{2}b\frac{db}{dx}\right)}_{a_I} dt + b \, dW_t$$

where $a_S$ is the Stratonovich drift and $a_I$ is the Itô drift. The noise-induced drift $-\frac{1}{2}b b'$ is called the spurious drift or noise-induced drift.

4. Itô's Lemma

The chain rule for stochastic calculus. If $x_t$ satisfies the Itô SDE $dx = a\,dt + b\,dW$ and $f(x,t)$ is a twice-differentiable function, then:

$$df = \left(\frac{\partial f}{\partial t} + a\frac{\partial f}{\partial x} + \frac{1}{2}b^2 \frac{\partial^2 f}{\partial x^2}\right)dt + b\frac{\partial f}{\partial x}\,dW_t$$

Derivation

Taylor expand $f(x + dx, t + dt)$ to second order:

$$df = f_t \, dt + f_x \, dx + \frac{1}{2}f_{xx}(dx)^2 + f_{xt}\,dx\,dt + \cdots$$

Substitute $dx = a\,dt + b\,dW$:

$$(dx)^2 = (a\,dt + b\,dW)^2 = a^2(dt)^2 + 2ab\,dt\,dW + b^2(dW)^2$$

Apply the Itô rules: $(dt)^2 = 0$, $dt \cdot dW = 0$,$(dW)^2 = dt$:

$$(dx)^2 = b^2 \, dt$$

Collecting terms of order $dt$ and $dW$ yields the lemma. The $\frac{1}{2}b^2 f_{xx}$ term -- absent in ordinary calculus -- is the hallmark of stochastic calculus and arises directly from $(dW)^2 = dt \neq 0$.

5. Euler-Maruyama Simulation

The simplest numerical scheme for an Itô SDE is the Euler-Maruyama method:

$$x_{n+1} = x_n + a(x_n)\,\Delta t + b(x_n)\,\Delta W_n, \qquad \Delta W_n \sim \mathcal{N}(0, \Delta t)$$

The strong order of convergence is $1/2$ and the weak order is $1$. For the Stratonovich interpretation, we use the Heun method(predictor-corrector):

$$\tilde{x}_{n+1} = x_n + b(x_n)\,\Delta W_n$$

$$x_{n+1} = x_n + a(x_n)\,\Delta t + \frac{1}{2}\bigl[b(x_n) + b(\tilde{x}_{n+1})\bigr]\,\Delta W_n$$

We now simulate Langevin dynamics in a bistable urban potential, modelling a city that can exist in two stable states (e.g., high-density vs low-density equilibrium).

Langevin Dynamics: Euler-Maruyama and Ito vs Stratonovich

Python

script.py135 lines

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# Bistable potential: V(x) = x^4/4 - x^2/2  ->  a(x) = -V'(x) = x - x^3
# Two stable points at x = +/-1, unstable at x = 0

def drift(x):
    return x - x**3

def diffusion(x):
    return 0.5  # constant noise

def potential(x):
    return x**4/4 - x**2/2

# --- Euler-Maruyama (Ito) ---
dt = 0.01
T = 50.0
N_steps = int(T / dt)
N_paths = 5

fig, axes = plt.subplots(2, 2, figsize=(14, 9))

# Top left: sample paths
ax = axes[0, 0]
for i in range(N_paths):
    x = np.zeros(N_steps + 1)
    x[0] = 0.1 * np.random.randn()
    for n in range(N_steps):
        dW = np.sqrt(dt) * np.random.randn()
        x[n+1] = x[n] + drift(x[n]) * dt + diffusion(x[n]) * dW
    t_arr = np.linspace(0, T, N_steps + 1)
    ax.plot(t_arr, x, linewidth=0.8, alpha=0.8)
ax.set_xlabel('Time', fontsize=12, color='white')
ax.set_ylabel('x(t)', fontsize=12, color='white')
ax.set_title('Euler-Maruyama: Sample Paths (Bistable)', fontsize=12, color='white')
ax.axhline(y=1, color='#10b981', linestyle='--', alpha=0.5, label='Stable equilibria')
ax.axhline(y=-1, color='#10b981', linestyle='--', alpha=0.5)
ax.axhline(y=0, color='#ef4444', linestyle=':', alpha=0.5, label='Unstable')
ax.legend(fontsize=9, facecolor='#1e293b', edgecolor='#10b981', labelcolor='white')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2)

# Top right: potential
ax = axes[0, 1]
xx = np.linspace(-1.8, 1.8, 200)
ax.plot(xx, potential(xx), color='#10b981', linewidth=2.5)
ax.fill_between(xx, potential(xx), potential(xx).max(), alpha=0.1, color='#10b981')
ax.set_xlabel('x', fontsize=12, color='white')
ax.set_ylabel('V(x)', fontsize=12, color='white')
ax.set_title('Bistable Potential $V(x) = x^4/4 - x^2/2$', fontsize=12, color='white')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2)

# Bottom left: Ito vs Stratonovich comparison with multiplicative noise
# SDE: dx = -x dt + sigma*x dW  (geometric Brownian motion)
sigma = 0.8
dt2 = 0.001
T2 = 5.0
N2 = int(T2 / dt2)
N_ensemble = 5000

# Ito: E[x(t)] = x0 * exp(-t)
# Stratonovich: E[x(t)] = x0 * exp((-1 + sigma^2/2)*t)
x0 = 1.0

x_ito = np.zeros((N_ensemble, N2+1))
x_strat = np.zeros((N_ensemble, N2+1))
x_ito[:, 0] = x0
x_strat[:, 0] = x0

for n in range(N2):
    dW = np.sqrt(dt2) * np.random.randn(N_ensemble)
    # Euler-Maruyama (Ito)
    x_ito[:, n+1] = x_ito[:, n] + (-x_ito[:, n]) * dt2 + sigma * x_ito[:, n] * dW
    # Heun (Stratonovich)
    x_pred = x_strat[:, n] + sigma * x_strat[:, n] * dW
    x_strat[:, n+1] = x_strat[:, n] + (-x_strat[:, n]) * dt2 + 0.5 * (sigma * x_strat[:, n] + sigma * x_pred) * dW

t2 = np.linspace(0, T2, N2+1)
mean_ito = np.mean(x_ito, axis=0)
mean_strat = np.mean(x_strat, axis=0)

ax = axes[1, 0]
ax.plot(t2, mean_ito, color='#10b981', linewidth=2, label='Ito (Euler-Maruyama)')
ax.plot(t2, mean_strat, color='#f59e0b', linewidth=2, label='Stratonovich (Heun)')
ax.plot(t2, x0 * np.exp(-t2), 'w--', linewidth=1.5, label=r'Theory Ito: $e^{-t}$')
ax.plot(t2, x0 * np.exp((-1 + sigma**2/2)*t2), '--', color='#f59e0b', linewidth=1.5, alpha=0.7,
        label=r'Theory Strat: $e^{(-1+\sigma^2/2)t}$')
ax.set_xlabel('Time', fontsize=12, color='white')
ax.set_ylabel(r'$\langle x(t) \rangle$', fontsize=12, color='white')
ax.set_title('Ito vs Stratonovich: Multiplicative Noise', fontsize=12, color='white')
ax.legend(fontsize=9, facecolor='#1e293b', edgecolor='#10b981', labelcolor='white')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2)

# Bottom right: histogram of final positions (bistable)
ax = axes[1, 1]
N_hist = 3000
x_final = np.zeros(N_hist)
T_hist = 20.0
N_hist_steps = int(T_hist / dt)
for i in range(N_hist):
    x = 0.01 * np.random.randn()
    for n in range(N_hist_steps):
        dW = np.sqrt(dt) * np.random.randn()
        x = x + drift(x) * dt + diffusion(x) * dW
    x_final[i] = x

ax.hist(x_final, bins=80, density=True, color='#10b981', alpha=0.6, edgecolor='#064e3b')
# Theoretical stationary: P(x) ~ exp(-2V(x)/b^2)
b_val = 0.5
xx_hist = np.linspace(-2, 2, 300)
P_theory = np.exp(-2 * potential(xx_hist) / b_val**2)
P_theory /= np.trapezoid(P_theory, xx_hist)
ax.plot(xx_hist, P_theory, 'w-', linewidth=2, label='Theory $P_{st}(x)$')
ax.set_xlabel('x', fontsize=12, color='white')
ax.set_ylabel('Probability density', fontsize=12, color='white')
ax.set_title('Stationary Distribution (Bistable)', fontsize=12, color='white')
ax.legend(fontsize=10, facecolor='#1e293b', edgecolor='#10b981', labelcolor='white')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2)

fig.patch.set_facecolor('#0f172a')
for ax_row in axes:
    for ax in ax_row:
        ax.set_facecolor('#0f172a')
plt.tight_layout()
plt.savefig('output.png', dpi=130, bbox_inches='tight', facecolor='#0f172a')
plt.show()
print("Simulations complete.")
print(f"Fraction ending near x=+1: {np.mean(x_final > 0.5):.3f}")
print(f"Fraction ending near x=-1: {np.mean(x_final < -0.5):.3f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

6. Noise-Induced Transitions

In the bistable potential, noise enables transitions between the two stable states. The mean transition time is given by Kramers' escape rate:

$$\tau_{\text{escape}} \sim \frac{2\pi}{\sqrt{|V''(x_{\text{unstable}})| \cdot V''(x_{\text{stable}})}} \exp\!\left(\frac{2\Delta V}{\sigma^2}\right)$$

where $\Delta V = V(x_{\text{unstable}}) - V(x_{\text{stable}})$ is the barrier height. For our potential $V(x) = x^4/4 - x^2/2$:

Stable points: $x_{\pm} = \pm 1$, with $V(\pm 1) = -1/4$
Unstable point: $x_0 = 0$, with $V(0) = 0$
Barrier height: $\Delta V = 1/4$
$V''(\pm 1) = 2$, $V''(0) = -1$

$$\tau_{\text{escape}} \sim \frac{2\pi}{\sqrt{2}} \exp\!\left(\frac{1}{2\sigma^2}\right)$$

For urban systems, this means that a city in one equilibrium (e.g., sprawl) will spontaneously transition to another (e.g., compact) only on exponentially long timescales unless external shocks push it over the barrier. The noise amplitude $\sigma$ represents the intensity of random perturbations (economic shocks, policy changes, migration waves).

← Part V Overview Fokker-Planck Equation →