Fractional Fokker-Planck Equation

Cities have memory. Past development patterns constrain future growth through infrastructure lock-in, legal precedent, and cultural inertia. The fractional Fokker-Planck equation replaces the ordinary time derivative with a fractional one, introducing power-law memory kernels that capture subdiffusive urban dynamics.

1. The Caputo Fractional Derivative

The Caputo fractional derivative of order $0 < \alpha < 1$ is defined as:

$$\,^C\!D_t^{\alpha}\,f(t) = \frac{1}{\Gamma(1-\alpha)}\int_0^t \frac{f'(\tau)}{(t - \tau)^{\alpha}}\,d\tau$$

This is a convolution of the ordinary derivative with a memory kernel:

$$K(t) = \frac{t^{-\alpha}}{\Gamma(1-\alpha)}, \qquad 0 < \alpha < 1$$

Key properties of the Caputo derivative:

When $\alpha = 1$, it reduces to the ordinary derivative $f'(t)$.
It weights recent changes more heavily, but all past history contributes (power-law memory).
Initial conditions take the familiar form $f(0) = f_0$, unlike the Riemann-Liouville version.
The kernel $K(t) \sim t^{-\alpha}$ is non-integrable at $t = 0$ for $\alpha > 0$, giving infinite weight to the infinitesimal past.

2. The Fractional Fokker-Planck Equation

The standard Fokker-Planck equation governs the probability density $P(x,t)$ of a diffusing particle subject to drift $a(x)$ and diffusion $b(x)$:

$$\frac{\partial P}{\partial t} = -\frac{\partial}{\partial x}\bigl[a(x)\,P\bigr] + \frac{1}{2}\frac{\partial^2}{\partial x^2}\bigl[b(x)^2\,P\bigr]$$

The fractional version replaces the left side with a Caputo derivative:

$$\,^C\!D_t^{\alpha}\,P = -\frac{\partial}{\partial x}\bigl[a(x)\,P\bigr] + \frac{1}{2}\frac{\partial^2}{\partial x^2}\bigl[b(x)^2\,P\bigr]$$

The fundamental consequence is anomalous diffusion. For free diffusion ($a = 0$, $b = \text{const}$), the mean squared displacement becomes:

$$\langle x^2(t) \rangle \sim t^{\alpha}$$

$\alpha = 1$: Normal diffusion. $\langle x^2 \rangle \sim t$.
$\alpha < 1$: Subdiffusion. Spreading is slower than normal. The particle gets “trapped” by memory.
$\alpha > 1$: Superdiffusion (requires fractional spatial derivatives instead).

3. Deriving the Memory Effect

The fractional FP can be rewritten as an integro-differential equation. Using the definition of the Caputo derivative and the Laplace transform relationship $\mathcal{L}\{{}^C\!D_t^{\alpha}f\}(s) = s^{\alpha}\hat{f}(s) - s^{\alpha-1}f(0)$, we can write the solution in terms of the Mittag-Leffler function:

$$E_{\alpha}(z) = \sum_{k=0}^{\infty}\frac{z^k}{\Gamma(\alpha k + 1)}$$

For the free-diffusion case, the fundamental solution (Green's function) is:

$$P(x,t) = \frac{1}{\sqrt{4\pi D_{\alpha}\,t^{\alpha}}}\,\exp\!\left(-\frac{x^2}{4D_{\alpha}\,t^{\alpha}}\right)$$

This is a Gaussian whose width grows as $t^{\alpha/2}$ rather than $t^{1/2}$. The physical picture is a continuous-time random walk (CTRW) where waiting times between jumps follow a power-law distribution:

$$\psi(\tau) \sim \tau^{-(1+\alpha)}, \qquad \langle\tau\rangle = \infty \text{ for } \alpha < 1$$

The infinite mean waiting time is the microscopic origin of subdiffusion: occasional very long trapping events dominate the transport process.

4. Urban Memory and Subdiffusion

Why would urban dynamics be subdiffusive? Several mechanisms create long-memory trapping:

Infrastructure lock-in: Once roads and sewers are built, development patterns are frozen for decades. The “waiting time” before a neighbourhood can be redeveloped follows a heavy-tailed distribution.
Legal and regulatory inertia: Zoning changes require lengthy approval processes with unpredictable delays.
Path dependence: Historical street layouts (some dating to medieval or even Roman times) constrain modern development, creating long-range temporal correlations.
Heterogeneous land tenure: Different ownership structures create a distribution of “resistance to change” times.

The fractional order $\alpha$ becomes a single number that quantifies the “stickiness” of urban history. Values measured from real cities typically fall in the range $0.6 \leq \alpha \leq 0.9$.

5. Numerical Solution via Grünwald-Letnikov

We discretise the fractional time derivative using the Grünwald-Letnikov scheme. The$\alpha$-th order derivative at time step $n$ is approximated by:

$$\,^{GL}\!D_t^{\alpha}\,P_i^n \approx \frac{1}{\Delta t^{\alpha}}\sum_{j=0}^{n} w_j^{(\alpha)}\,P_i^{n-j}$$

where the Grünwald-Letnikov weights satisfy the recursion:

$$w_0^{(\alpha)} = 1, \qquad w_j^{(\alpha)} = \left(1 - \frac{1+\alpha}{j}\right)w_{j-1}^{(\alpha)}$$

Fractional Fokker-Planck: Normal vs Subdiffusive Urban Spreading

Python

script.py121 lines

import numpy as np
import matplotlib.pyplot as plt

# --- Grunwald-Letnikov weights ---
def gl_weights(alpha, n):
    """Compute first n Grunwald-Letnikov weights."""
    w = np.zeros(n)
    w[0] = 1.0
    for j in range(1, n):
        w[j] = w[j-1] * (1.0 - (1.0 + alpha) / j)
    return w

# --- Fractional Fokker-Planck solver ---
def solve_ffp(alpha, Nx=201, Nt=600, L=10.0, T=50.0, D=0.2):
    dx = 2*L / (Nx - 1)
    dt = T / Nt
    x = np.linspace(-L, L, Nx)

# Initial condition: Gaussian
    P = np.exp(-x**2 / 0.5) / np.sqrt(0.5 * np.pi)
    P /= np.trapezoid(P, x)  # normalise

# Storage for history (needed for fractional derivative)
    P_history = [P.copy()]
    weights = gl_weights(alpha, Nt + 1)

snapshots = {}
    snap_times = [0, Nt//6, Nt//3, Nt//2, 2*Nt//3, Nt-1]
    msd = [np.trapezoid(x**2 * P, x)]
    times_out = [0.0]

for n in range(1, Nt):
        # Compute memory sum: sum_{j=1}^{n} w_j * P^{n-j}
        mem_sum = np.zeros(Nx)
        for j in range(1, min(n+1, len(P_history))):
            mem_sum += weights[j] * P_history[n - j]

# Spatial Laplacian of P (central difference)
        Pold = P_history[-1]
        laplacian = np.zeros(Nx)
        laplacian[1:-1] = (Pold[2:] - 2*Pold[1:-1] + Pold[:-2]) / dx**2

# Update: w_0 * P^n = dt^alpha * D * Lap(P) - mem_sum
        # => P^n = dt^alpha * D * Lap(P) - mem_sum  (since w_0 = 1)
        P_new = dt**alpha * D * laplacian - mem_sum
        # Ensure non-negative
        P_new = np.maximum(P_new, 0)
        # Renormalise (fractional schemes can drift)
        norm = np.trapezoid(P_new, x)
        if norm > 1e-10:
            P_new /= norm

P_history.append(P_new.copy())

# Record MSD
        msd.append(np.trapezoid(x**2 * P_new, x))
        times_out.append(n * dt)

if n in snap_times:
            snapshots[n] = (P_new.copy(), n*dt)

return x, np.array(times_out), np.array(msd), snapshots, P_history[0]

# --- Run for two values of alpha ---
x1, t1, msd1, snaps1, P0_1 = solve_ffp(alpha=1.0, Nt=400, T=30.0)
x2, t2, msd2, snaps2, P0_2 = solve_ffp(alpha=0.7, Nt=400, T=30.0)

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

# (a) PDF evolution: alpha=1
ax = axes[0]
ax.plot(x1, P0_1, 'k--', linewidth=1, label='t=0')
colors = plt.cm.YlGn(np.linspace(0.3, 0.9, len(snaps1)))
for idx, (step, (P, t)) in enumerate(sorted(snaps1.items())):
    ax.plot(x1, P, color=colors[idx], linewidth=1.5, label=f't={t:.1f}')
ax.set_xlabel('x')
ax.set_ylabel('P(x,t)')
ax.set_title('Normal diffusion (alpha=1.0)')
ax.legend(fontsize=7)
ax.set_ylim(bottom=0)

# (b) PDF evolution: alpha=0.7
ax = axes[1]
ax.plot(x2, P0_2, 'k--', linewidth=1, label='t=0')
colors2 = plt.cm.YlGn(np.linspace(0.3, 0.9, len(snaps2)))
for idx, (step, (P, t)) in enumerate(sorted(snaps2.items())):
    ax.plot(x2, P, color=colors2[idx], linewidth=1.5, label=f't={t:.1f}')
ax.set_xlabel('x')
ax.set_ylabel('P(x,t)')
ax.set_title('Subdiffusion (alpha=0.7)')
ax.legend(fontsize=7)
ax.set_ylim(bottom=0)

# (c) MSD comparison
ax = axes[2]
mask1 = t1 > 0.1
mask2 = t2 > 0.1
ax.loglog(t1[mask1], msd1[mask1], 'o-', color='teal', markersize=2,
          linewidth=1, label='alpha=1.0')
ax.loglog(t2[mask2], msd2[mask2], 's-', color='#059669', markersize=2,
          linewidth=1, label='alpha=0.7')

# Reference slopes
t_ref = np.linspace(1, 20, 50)
ax.loglog(t_ref, 0.3*t_ref**1.0, '--', color='gray', alpha=0.5, label='slope=1')
ax.loglog(t_ref, 0.3*t_ref**0.7, ':', color='gray', alpha=0.5, label='slope=0.7')
ax.set_xlabel('t')
ax.set_ylabel('<x^2>')
ax.set_title('Mean Squared Displacement')
ax.legend(fontsize=9)

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

print("Fractional Fokker-Planck: Normal vs Subdiffusive spreading")
print(f"  alpha=1.0: final MSD = {msd1[-1]:.4f}")
print(f"  alpha=0.7: final MSD = {msd2[-1]:.4f}")
print(f"  Ratio: {msd1[-1]/max(msd2[-1],1e-10):.2f}x")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

6. Summary and Connections

The fractional Fokker-Planck equation provides a minimal model for urban dynamics with memory. The single parameter $\alpha$ bridges two limits:

$\alpha = 1$: Memoryless (Markovian) diffusion—the standard Fokker-Planck of Module 2.
$\alpha < 1$: Subdiffusive spreading with power-law memory—historically constrained cities.

The framework connects to the CTRW (continuous-time random walk) picture: development events occur at random times with a heavy-tailed waiting-time distribution, leading naturally to fractional calculus as the macroscopic description.

← DLA Fractals Quantum Game Theory →