In the previous post I built the Ornstein-Uhlenbeck (OU) model of a moose with a home range, and ended with a list of nine defects and a promise: rebuild the model in the general Fokker-Planck form, repairing the defects one term at a time. That is this post. The OU model was the special case in which the drift is a spring, the diffusion is a constant, and everything stays Gaussian and closed-form. The general equation is the framework in which habitat, barriers, interactions, behavioural states, and jumps can all be written down explicitly. The price is that almost nothing in it is solvable by hand: every result below is computed numerically.
The general equation
The Fokker-Planck equation for the density \(p(x,t)\) of an animal’s position, in the form that can carry all of the repairs at once, is
\[\frac{\partial p}{\partial t} = -\nabla \cdot \bigl(b[x,p]\, p\bigr) + \frac{1}{2}\,\partial_i \partial_j \bigl(D_{ij}(x)\, p\bigr) - \lambda(x)\,p.\]
Compared to the equation in the previous post, three things have become general. The drift \(b\) may depend on space through a habitat field and on the density itself through a nonlocal coupling. The diffusion \(D_{ij}\) is a field rather than a constant, so movement can be easier in some terrain than others. The sink term \(\lambda\) removes mass, modelling mortality or animals leaving the study area. Each change answers a specific defect from the OU list:
- Habitat instead of a spring (defects 1 and 2). Write the drift as \(b = -\nabla U\) for a habitat potential \(U\), large where habitat is poor. The potential-function approach to movement is standard in ecology (Preisler et al. 2013); nothing forces \(U\) to be quadratic or its level sets to be circles, so the home range can be multi-modal and shaped like the landscape.
- A diffusion field (defect 3). \(D(x)\) varies across snow, muskeg, forest, and open ground, and along corridors.
- Barriers as boundary conditions (defect 4). Rivers, lakes, and highways become walls on which the probability current is set to zero, with gaps where crossing is possible.
- Interactions as nonlocal coupling (defect 5). Each animal climbs the local density of conspecifics, \(b[x,p] = -\nabla U + \chi \, \nabla \log (G * p)\), where \(G\) is a smoothing kernel and \(\chi\) the attraction strength. This is a mean-field stand-in for the local attraction and alignment rules that organize real groups (Couzin et al. 2002).
- Behaviour as multiple states (defects 6 and 7). Run several densities \(p_i\) at once, each with its own potential and diffusivity, coupled by switching rates \(q_{ij}\). A moose alternates between a tight core range and a wider exploratory mode; the observable density is the mixture.
- Jumps as fractional diffusion (defect 9). Replace the Laplacian by the fractional Laplacian of order \(2\alpha\), so the diffusion term becomes \(D\,(-\Delta)^{\alpha} p\) with \(\alpha = 1\) the Brownian case and \(\alpha < 1\) a process that makes occasional long relocations (Metzler et al. 1994; Xie et al. 2015).
One property of the OU model survives generalization and becomes the workhorse check for everything below. With a constant scalar \(D\) and no-flux walls, the current \(J = -\nabla U\, p - D \nabla p\) vanishes identically for the Boltzmann equilibrium
\[p_\infty(x) \propto \exp\!\left(-\frac{U(x)}{D}\right),\]
which reduces to exactly the Gaussian of the previous post when \(U\) is the spring \(\kappa(x-c)^2/2\) and \(D = \sigma^2/2\). For constant \(D\), habitat alone determines where the animal concentrates; the numerical question is whether the scheme reproduces it.
Experiment 1: habitat, a river, and a bridge
The Python script below puts the general equation on a two-dimensional landscape: a \(4 \times 4\) km domain on an \(81 \times 81\) grid with two Gaussian habitat wells (two home ranges), a river along \(x = 2\) km with a single bridge near \(y = 0.7\) km, and constant \(D = 0.5\ \text{km}^2/\text{h}\). All density starts in the left well. Three experiments are run: relaxation to equilibrium with and without the river, and a final run in which \(D(x)\) has a hot spot of easier movement near the right well, where the Boltzmann formula no longer applies.
"""A general Fokker-Planck model of moose movement in two dimensions.
This solves the Fokker-Planck equation
dp/dt = -div[ b(x) p ] + d_i d_j [ D_ij(x) p ],
with a habitat potential U so that the drift is b = -grad U (plus an optional
nonlocal interaction term), a spatially varying diffusion D(x), and no-flux
walls around the domain and around a barrier (a river with a crossing).
Three experiments are reported:
1. equilibrium for a constant diffusion: p ~ exp(-U / D), including the effect
of the barrier on how fast that equilibrium is approached;
2. the transient mixing of the two home-range lobes with and without the
barrier;
3. a spatially varying diffusion, for which exp(-U / D) is no longer the
stationary solution.
"""
from __future__ import annotations
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
from moose_fp.numerics import total_variation
LENGTH = 4.0 # Domain side length, km
GRID = 81 # Cells per side
D0 = 0.5 # Base diffusion, km^2/h
DEPTH = 2.0 # Habitat well depth
WIDTH = 0.5 # Habitat well width, km
def build_domain(
length: float = LENGTH,
grid: int = GRID,
depth: float = DEPTH,
well_width: float = WIDTH,
):
"""Return ``(x, y, xx, yy, potential, solid)`` for the toy landscape.
The landscape has two Gaussian habitat wells and a river at ``x = 2`` km
with a single bridge near ``y = 0.7`` km. ``solid`` marks the cells that
are walled off (the river outside the bridge).
"""
x = np.linspace(0.0, length, grid)
y = np.linspace(0.0, length, grid)
xx, yy = np.meshgrid(x, y, indexing="ij")
# Two home-range wells, one on each side of the river.
wells = [(1.0, 2.0), (3.0, 2.0)]
potential = np.zeros_like(xx)
for wx, wy in wells:
potential -= depth * np.exp(
-((xx - wx) ** 2 + (yy - wy) ** 2) / (2 * well_width**2)
)
# A river at x = 2 with a single bridge at y ~ 0.7.
solid = np.zeros_like(xx, dtype=bool)
wall = np.isclose(xx, 2.0, atol=0.75 * (x[1] - x[0]))
bridge = np.abs(yy - 0.7) < 0.2
solid[wall & ~bridge] = True
return x, y, xx, yy, potential, solid
def steady_state(potential, solid, diffusion, weights=None):
"""Boltzmann prediction p ~ exp(-U / D) on the open cells."""
if weights is None:
weights = np.ones_like(potential)
p = np.where(solid, 0.0, weights * np.exp(-potential / diffusion))
p /= p.sum()
return p
def gradient(field, spacing):
"""Central gradient with second-order edges."""
return np.gradient(field, spacing, edge_order=2)
def solve(
potential,
solid,
diffusion,
p0,
t_end,
interaction=0.0,
record_times=(),
):
"""Explicit finite-volume integration of the general Fokker-Planck equation.
``diffusion`` may be a scalar or a two-dimensional array. ``interaction``
is the strength of a nonlocal attraction towards the mean of the
distribution (McKean-Vlasov term); zero means no interaction.
Returns ``(p, dt, saved, open_)`` where ``saved`` maps each requested
record time to a copy of the density at that time.
"""
nx, ny = potential.shape
dx = LENGTH / (nx - 1)
dy = LENGTH / (ny - 1)
if np.isscalar(diffusion):
diffusion = np.full_like(potential, diffusion)
bx, by = gradient(-potential, dx)
d_max = float(np.max(diffusion[~solid]))
b_max = max(float(np.max(np.abs(bx[~solid]))), float(np.max(np.abs(by[~solid]))))
# 2D explicit stability: diffusion number D dt / dx^2 <= 1/4 and
# advective CFL |b| dt / dx <= 1, with a safety factor.
dt = 0.2 * min(dx**2 / d_max, dx / max(b_max, 1e-12))
n_steps = round(t_end / dt)
p = np.zeros_like(potential)
p[:] = p0
p[solid] = 0.0
open_ = ~solid
saved = {}
next_record = 0
for step in range(n_steps):
drift_x, drift_y = bx, by
if interaction != 0.0:
smoothed = smooth(p)
gx, gy = gradient(smoothed, dx)
drift_x = bx - interaction * gx
drift_y = by - interaction * gy
flux_x = face_flux(p, drift_x, diffusion, solid, axis=0, spacing=dx)
flux_y = face_flux(p, drift_y, diffusion, solid, axis=1, spacing=dy)
divergence = divergence_from_flux(flux_x, flux_y, spacing=dx)
p = p - dt * divergence
p[solid] = 0.0
p[p < 0.0] = 0.0
if next_record < len(record_times) and (step + 1) * dt >= record_times[next_record]:
saved[record_times[next_record]] = p.copy()
next_record += 1
return p, dt, saved, open_
def face_flux(p, drift, diffusion, solid, axis, spacing):
"""Advective (upwinded) plus diffusive flux on the cell faces."""
if axis == 0:
solid_face = solid[:-1, :] | solid[1:, :]
drift_face = 0.5 * (drift[:-1, :] + drift[1:, :])
diff_face = 0.5 * (diffusion[:-1, :] + diffusion[1:, :])
left, right = p[:-1, :], p[1:, :]
else:
solid_face = solid[:, :-1] | solid[:, 1:]
drift_face = 0.5 * (drift[:, :-1] + drift[:, 1:])
diff_face = 0.5 * (diffusion[:, :-1] + diffusion[:, 1:])
left, right = p[:, :-1], p[:, 1:]
upwind = np.where(drift_face > 0, left, right)
flux = drift_face * upwind - diff_face * (right - left) / spacing
return np.where(solid_face, 0.0, flux)
def divergence_from_flux(flux_x, flux_y, spacing):
"""Divergence of a face-centred flux on a cell-centred grid."""
nx, ny = flux_x.shape[0] + 1, flux_x.shape[1]
divergence = np.zeros((nx, ny))
face_x = np.zeros((nx + 1, ny))
face_x[1:-1] = flux_x
divergence += (face_x[1:] - face_x[:-1]) / spacing
face_y = np.zeros((nx, ny + 1))
face_y[:, 1:-1] = flux_y
divergence += (face_y[:, 1:] - face_y[:, :-1]) / spacing
return divergence
def smooth(field):
"""Small Gaussian smoothing used for the nonlocal interaction term."""
from scipy.ndimage import gaussian_filter
return gaussian_filter(field, sigma=2.0, mode="nearest")
def centre_delta(x, y, cx, cy, solid):
"""A point mass at the grid node nearest ``(cx, cy)``, normalised."""
p = np.zeros((len(x), len(y)))
i = int(np.argmin(np.abs(x - cx)))
j = int(np.argmin(np.abs(y - cy)))
p[i, j] = 1.0
p[solid] = 0.0
return p / p.sum()
def _plot(xx, yy, potential, solid, equilibrium, record_times, error_barrier,
error_free, open_cells, p_variable, diffusion_field):
fig, axes = plt.subplots(2, 2, figsize=(12, 10.5))
extent = [0, LENGTH, 0, LENGTH]
ax = axes[0, 0]
image = ax.imshow(potential.T, origin="lower", extent=extent, cmap="terrain")
ax.contour(xx.T, yy.T, potential.T, levels=8, colors="black", linewidths=0.4)
ax.imshow(
np.where(solid, 1.0, np.nan).T,
origin="lower",
extent=extent,
cmap="Greys",
vmin=0,
vmax=1,
alpha=0.9,
)
ax.set_title("Habitat potential U (dark = good) and river with bridge")
ax.set_xlabel("x (km)")
ax.set_ylabel("y (km)")
fig.colorbar(image, ax=ax, label="U")
ax = axes[0, 1]
image = ax.imshow(equilibrium.T, origin="lower", extent=extent, cmap="viridis")
ax.set_title(r"Equilibrium for constant $D$: $p \propto e^{-U/D}$")
ax.set_xlabel("x (km)")
ax.set_ylabel("y (km)")
fig.colorbar(image, ax=ax, label="density")
ax = axes[1, 0]
ax.plot(record_times, error_barrier, color="darkred", linewidth=2, label="with river")
ax.plot(record_times, error_free, color="steelblue", linewidth=2, label="no river")
ax.set_yscale("log")
ax.set_xlabel("t (h)")
ax.set_ylabel("total-variation distance to equilibrium")
ax.set_title("The river slows mixing but not the equilibrium")
ax.legend()
ax = axes[1, 1]
scatter = ax.scatter(
-potential[open_cells] / D0,
np.log(p_variable[open_cells]),
c=np.log10(diffusion_field[open_cells]),
s=12,
cmap="plasma",
)
fit = np.polyfit(-potential[open_cells] / D0, np.log(p_variable[open_cells]), 1)
xs = np.linspace(
(-potential[open_cells] / D0).min(), (-potential[open_cells] / D0).max(), 10
)
ax.plot(xs, np.polyval(fit, xs), color="black", linewidth=2, label=f"slope {fit[0]:.2f}")
ax.set_xlabel(r"$-U/D_0$")
ax.set_ylabel(r"$\log p$ (steady state)")
ax.set_title("Variable diffusion breaks the Boltzmann law")
ax.legend()
fig.colorbar(scatter, ax=ax, label=r"$\log_{10} D$")
fig.suptitle("General Fokker-Planck model of moose movement")
fig.tight_layout()
return fig
def run(output: str | Path = "figures/general_fp_2d.png") -> Path:
"""Run the three experiments, print diagnostics, and write the figure."""
x, y, xx, yy, potential, solid = build_domain()
equilibrium = steady_state(potential, solid, D0)
# --- Experiment 1 and 2: barrier versus no barrier -------------------
record_times = np.linspace(0.0, 60.0, 121)
start = centre_delta(x, y, 1.0, 2.0, solid)
_p_barrier, _dt, saved_barrier, _ = solve(
potential, solid, D0, start, 60.0, record_times=tuple(record_times)
)
p_free, _, saved_free, _ = solve(
potential, np.zeros_like(solid), D0, start, 60.0, record_times=tuple(record_times)
)
error_barrier = np.array(
[total_variation(saved_barrier[t], equilibrium) for t in record_times]
)
error_free = np.array([total_variation(saved_free[t], equilibrium) for t in record_times])
# --- Experiment 3: spatially varying diffusion -----------------------
diffusion_field = D0 * (1.0 + 0.9 * np.exp(-((xx - 3.0) ** 2 + (yy - 1.0) ** 2) / 1.5))
p_variable, _, _, _ = solve(potential, solid, diffusion_field, start, 120.0)
prediction = steady_state(potential, solid, D0)
open_cells = ~solid & (prediction > 1e-5)
ratios = p_variable[open_cells] / prediction[open_cells]
print("\nBarrier delays mixing:")
for t in (2.0, 5.0, 10.0, 20.0, 40.0, 60.0):
k = int(np.argmin(np.abs(record_times - t)))
print(
f" t = {record_times[k]:5.1f} h"
f" TV to equilibrium: barrier {error_barrier[k]:.4f}"
f" no barrier {error_free[k]:.4f}"
)
max_rel_err = float(
np.max(
np.abs(
p_free[open_cells] / p_free.sum() - equilibrium[open_cells] / equilibrium.sum()
)
/ equilibrium[open_cells]
)
)
corr = float(
np.corrcoef(np.log(p_variable[open_cells]), -potential[open_cells] / D0)[0, 1]
)
print("\nSolver validation (constant D, no barrier):")
print(f" max relative error versus exp(-U/D): {max_rel_err:.4f}")
print(f" TV distance at t = 60 h: {error_free[-1]:.5f}")
print("\nVariable diffusion versus the constant-D Boltzmann prediction:")
print(f" median ratio p_steady / (exp(-U/D0) normalised): {np.median(ratios):.3f}")
print(
f" interquartile range of ratio: "
f"{np.subtract(*np.percentile(ratios, [75, 25])):.3f}"
)
print(f" correlation of log p with -U/D0: {corr:.3f}")
fig = _plot(
xx,
yy,
potential,
solid,
equilibrium,
record_times,
error_barrier,
error_free,
open_cells,
p_variable,
diffusion_field,
)
output_path = Path(output)
output_path.parent.mkdir(parents=True, exist_ok=True)
fig.savefig(output_path, dpi=300, bbox_inches="tight", facecolor="white")
plt.close(fig)
return output_path
def main() -> None:
"""Console-script entry point."""
run()
if __name__ == "__main__":
main()
Three things are visible in the figure and confirmed by the printed diagnostics.
The equilibrium is bimodal. With two wells, \(e^{-U/D}\) has two humps. Defect 1 is repaired by construction: multi-modal utilization distributions are ordinary output, not something the modeler has to disbelieve.
The river slows mixing but does not change the equilibrium. No-flux walls cannot change a zero-current solution, so both landscapes relax towards the same \(p_\infty\); what the river changes is how long redistribution takes. Measured in total-variation distance to equilibrium:
Barrier delays mixing:
t = 2.0 h TV to equilibrium: barrier 0.4884 no barrier 0.3573
t = 5.0 h TV to equilibrium: barrier 0.4625 no barrier 0.2031
t = 10.0 h TV to equilibrium: barrier 0.4218 no barrier 0.0819
t = 20.0 h TV to equilibrium: barrier 0.3508 no barrier 0.0485
t = 40.0 h TV to equilibrium: barrier 0.2427 no barrier 0.0479
t = 60.0 h TV to equilibrium: barrier 0.1679 no barrier 0.0479After sixty hours the free landscape has essentially settled while the riverine one is still far from it: the bridge is the only channel, and mixing is gated by its width. For management this is the difference between where animals will eventually be and where they can be reached in a season.
The solver has its own physics, again. With the river removed, the numerical steady state is compared against \(e^{-U/D}\) directly:
Solver validation (constant D, no barrier):
max relative error versus exp(-U/D): 0.1683
TV distance at t = 60 h: 0.04789The error is not concentrated in the tails: the median relative error over open cells is about \(0.155\), because first-order upwinding of the drift adds an artificial diffusion of size roughly \(|b|\,\Delta x/2 \approx 0.06\ \text{km}^2/\text{h}\) against a physical \(D = 0.5\) — a twelve percent inflation of the effective diffusivity. Refining the grid from \(81^2\) to \(161^2\) halves both the median relative error (\(0.155 \to 0.078\)) and the total-variation distance (\(0.048 \to 0.024\)), the first-order convergence the scheme owes. This is the same lesson as in the OU post, now in two dimensions: solving the Fokker-Planck equation always means solving a discretization of it.
Finally, the variable-diffusion run. Where \(D\) varies, \(e^{-U/D_0}\) is no longer the stationary density:
Variable diffusion versus the constant-D Boltzmann prediction:
median ratio p_steady / (exp(-U/D0) normalised): 1.461
interquartile range of ratio: 0.721
correlation of log p with -U/D0: 0.943The potential still organizes the density — the correlation stays high — but no longer determines it: the equilibrium is shifted by almost fifty percent at the median, with a spread of the same order. There is no universal closed-form equilibrium left to check a solver against; where movement ease varies across a landscape, the pattern of \(D\) matters as much as the pattern of \(U\), and solving the PDE is not optional.
Experiment 2: three repairs in one dimension
The second script isolates three of the general mechanisms in one dimension, where each can be checked against something.
"""Three ways of repairing the Ornstein-Uhlenbeck model of moose movement.
1. Fractional diffusion: replace the Laplacian by the fractional Laplacian of
order 2 alpha, so the Fokker-Planck operator becomes
dp/dt = -d/dx (b p) + D (-d^2/dx^2)^alpha p ,
and the stationary density acquires power-law tails instead of decaying
like a Gaussian. Solved with a split scheme: the exact Ornstein-Uhlenbeck
contraction, remapped conservatively cell by cell, then exact fractional
diffusion applied in Fourier space.
2. Nonlocal interaction: add a McKean-Vlasov drift that pulls each animal
towards the local density of animals, in the form
b(x) = -kappa x + chi d/dx log(G * p) .
Solved with the explicit finite-volume scheme until it stops changing.
3. Behavioural states: two coupled Fokker-Planck equations for a moose that
switches between a tight core range and a wider exploratory mode.
dp_i/dt = -d/dx (b_i p_i) + D_i d^2p_i/dx^2 + q_ji p_j - q_ij p_i
Solved in time from a point mass, then compared against a single-state
Ornstein-Uhlenbeck fit.
"""
from __future__ import annotations
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
from scipy.ndimage import gaussian_filter1d
from moose_fp.numerics import total_variation
ALPHA = 0.5 # Fractional order of the diffusion operator / 2
KAPPA = 1.0 # OU fidelity, 1/h
DIFFUSION = 1.0 # Diffusion coefficient, km^2/h
def flow_contract(p, x, a):
"""Push cell masses through the exact OU contraction x -> x/a.
The Ornstein-Uhlenbeck drift -kappa x points inward everywhere, so the
deterministic flow of one step maps the domain into itself without any
wraparound: mass in the outer strips moves in and leaves a hole near the
seam, which the diffusion step then refills. The remap reads mass off the
cumulative distribution and conserves it exactly.
"""
nodes = np.append(x, x[-1] + (x[1] - x[0]))
cumulative = np.concatenate(([0.0], np.cumsum(p)))
preimages = np.clip(a * nodes, nodes[0], nodes[-1])
moved = np.interp(preimages, nodes, cumulative)
return moved[1:] - moved[:-1]
def fractional_ou(alpha, n_steps, dt, grid, length):
"""Stationary density of the fractional OU process on a periodic domain.
Returns ``(x, p, spacing, mass_drift)`` where ``p`` is normalised as a
density and ``mass_drift`` is the largest deviation from total mass one
seen over the run.
"""
x = -0.5 * length + np.arange(grid) * (length / grid)
spacing = length / grid
k = 2.0 * np.pi * np.fft.fftfreq(grid, d=spacing)
p = np.exp(-0.5 * KAPPA * x**2)
p /= p.sum()
contraction = np.exp(KAPPA * dt)
evolved = np.exp(-DIFFUSION * np.abs(k) ** (2.0 * alpha) * dt)
low = high = p.sum()
for _ in range(n_steps):
p = flow_contract(p, x, contraction)
p = np.real(np.fft.ifft(np.fft.fft(p) * evolved))
low = min(low, p.sum())
high = max(high, p.sum())
p = p / (p.sum() * spacing)
return x, p, spacing, high - low
def advance_periodic(p, drift, spacing, diffusion, dt):
"""Conservative explicit step for dp/dt = -d/dx (b p) + d^2 (D p)/dx^2."""
b_face = 0.5 * (drift + np.roll(drift, -1))
upwind = np.where(b_face > 0.0, p, np.roll(p, -1))
flux = upwind * b_face - diffusion * (np.roll(p, -1) - p) / spacing
p = p - dt * (flux - np.roll(flux, 1)) / spacing
return np.maximum(p, 0.0)
def drift_habitat(x):
"""The Ornstein-Uhlenbeck restoring drift -kappa x."""
return -KAPPA * x
def interaction_run(chi, grid, length, dt, t_end, width, tolerance=1e-10):
"""Equilibrium with a nonlocal attraction towards the animal density.
Returns ``(x, p, elapsed_time, final_change)``.
"""
x = -0.5 * length + np.arange(grid) * (length / grid)
spacing = length / grid
p = np.exp(-0.5 * KAPPA * x**2)
p /= p.sum()
n_steps = round(t_end / dt)
previous = p.copy()
change = float("nan")
for step in range(n_steps):
drift = drift_habitat(x)
if chi != 0.0:
smoothed = gaussian_filter1d(p, sigma=width / spacing, mode="nearest")
smoothed = np.maximum(smoothed, 1e-14 * smoothed.max())
drift = drift + chi * np.gradient(np.log(smoothed), spacing)
p = advance_periodic(p, drift, spacing, DIFFUSION, dt)
if step % 200 == 0 and step > 0:
change = 0.5 * np.sum(np.abs(p - previous))
if change < tolerance:
return x, p, (step + 1) * dt, change
previous = p.copy()
return x, p, n_steps * dt, change
def self_consistent_variance(chi, width, iterations=400):
"""Fixed point of v = D / (kappa + chi / (v + width^2)) for a Gaussian."""
variance = DIFFUSION / KAPPA
for _ in range(iterations):
variance = DIFFUSION / (KAPPA + chi / (variance + width**2))
return variance
def face_flux(p, drift, spacing, diffusion):
"""Advective (upwinded) plus diffusive flux on the cell faces."""
flux = np.zeros(len(p) + 1)
interior_b = 0.5 * (drift[:-1] + drift[1:])
interior = (
np.where(interior_b > 0.0, p[:-1], p[1:]) * interior_b
- diffusion * (p[1:] - p[:-1]) / spacing
)
flux[1:-1] = interior
return flux
def switch_run(x, spacing, dt, t_end, snapshot_at):
"""Time-dependent solution of the two coupled state equations.
Returns a dict mapping each requested snapshot time to ``(p_a, p_b)``.
"""
state_a = {"kappa": 1.0, "centre": -3.0, "diffusion": 0.25, "out": 0.10}
state_b = {"kappa": 0.5, "centre": 2.5, "diffusion": 0.60, "out": 0.10}
p_a = np.where(np.isclose(x, -3.0), 1.0, 0.0)
p_b = np.zeros_like(p_a)
p_a /= p_a.sum()
n_steps = round(t_end / dt)
snapshots = {}
for step in range(n_steps):
flux_a = face_flux(
p_a, -state_a["kappa"] * (x - state_a["centre"]), spacing, state_a["diffusion"]
)
flux_b = face_flux(
p_b, -state_b["kappa"] * (x - state_b["centre"]), spacing, state_b["diffusion"]
)
p_a_next = p_a - dt * (flux_a[1:] - flux_a[:-1]) / spacing
p_b_next = p_b - dt * (flux_b[1:] - flux_b[:-1]) / spacing
p_a_next += dt * (-state_a["out"] * p_a + state_b["out"] * p_b)
p_b_next += dt * (-state_b["out"] * p_b + state_a["out"] * p_a)
p_a = np.maximum(p_a_next, 0.0)
p_b = np.maximum(p_b_next, 0.0)
for when in snapshot_at:
if abs((step + 1) * dt - when) < 0.5 * dt:
snapshots[when] = (p_a.copy(), p_b.copy())
return snapshots
def sample_path(grid_dt, t_end, seed=20261016):
"""One representative animal path with switching behaviour.
Returns ``(times, position, states)``.
"""
rng = np.random.default_rng(seed)
n = round(t_end / grid_dt)
times = np.arange(n + 1) * grid_dt
states = np.zeros(n + 1, dtype=int)
position = np.zeros(n + 1)
for i in range(n):
states[i + 1] = states[i]
if rng.random() < 0.10 * grid_dt:
states[i + 1] = 1 - states[i]
current = states[i + 1]
kappa, centre, diffusion = [
(1.0, -3.0, 0.25),
(0.5, 2.5, 0.60),
][current]
drift = -kappa * (position[i] - centre)
position[i + 1] = (
position[i] + drift * grid_dt + np.sqrt(2.0 * diffusion * grid_dt) * rng.normal()
)
return times, position, states
def _plot(x_f, p_alpha1, p_alphahalf, gaussian, slope, intercept, x_i, p_free,
p_inter, x_s, mixture, p_a, p_b, gaussian_fit, times, path):
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
ax = axes[0, 0]
core = (x_f > 0.4) & (x_f < 10.0) & (p_alpha1 > 1e-12)
ax.loglog(x_f[core], p_alpha1[core], color="steelblue", linewidth=2,
label=r"$\alpha = 1$ (Gaussian)")
ax.loglog(x_f[core], gaussian[core], color="black", linestyle="--", linewidth=1.5,
label="exact Gaussian")
tail_alpha = (x_f > 0.4) & (x_f < 10.0) & (p_alphahalf > 1e-14)
ax.loglog(x_f[tail_alpha], p_alphahalf[tail_alpha], color="darkred", linewidth=2,
label=r"$\alpha = 0.5$")
reference = np.exp(intercept) * (x_f[tail_alpha] ** slope)
ax.loglog(x_f[tail_alpha], reference, color="darkorange", linestyle=":", linewidth=1.5,
label=f"tail slope {slope:.2f}")
ax.set_xlabel("distance from home range centre (km)")
ax.set_ylabel("stationary density")
ax.set_title("Fractional diffusion gives power-law tails")
ax.legend(fontsize=9)
ax = axes[0, 1]
ax.plot(x_i, p_free, color="steelblue", linewidth=2,
label=r"$\chi = 0$ (independent moose)")
ax.plot(x_i, p_inter, color="darkred", linewidth=2,
label=r"$\chi = 0.6$ (nonlocal attraction)")
ax.set_xlabel("distance from home range centre (km)")
ax.set_ylabel("stationary density")
ax.set_title("Interaction makes the herd bunch together")
ax.legend(fontsize=9)
ax = axes[1, 0]
ax.fill_between(x_s, mixture, color="grey", alpha=0.35, label="total density")
ax.plot(x_s, p_a, color="steelblue", linewidth=2, label="state A: core range")
ax.plot(x_s, p_b, color="darkred", linewidth=2, label="state B: exploratory")
ax.plot(x_s, gaussian_fit, color="black", linestyle="--", linewidth=1.5,
label="single-state fit")
ax.set_xlabel("distance from home range centre (km)")
ax.set_ylabel("stationary density")
ax.set_title("Behavioural states turn one Gaussian into two lobes")
ax.legend(fontsize=9)
ax = axes[1, 1]
for _mode, colour, label in ((0, "steelblue", "state A"), (1, "darkred", "state B")):
ax.scatter(times, path, s=2.5, color=colour, label=label, rasterized=True)
ax.set_xlabel("t (h)")
ax.set_ylabel("position (km)")
ax.set_title("One animal switching between behaviours")
ax.legend(fontsize=9, markerscale=3)
fig.suptitle("Beyond the Ornstein-Uhlenbeck model")
fig.tight_layout()
return fig
def run(output: str | Path = "figures/beyond_gaussian_fp.png") -> Path:
"""Run the three one-dimensional repairs, print diagnostics, write figure."""
# --- Fractional tails ------------------------------------------------
x_f, p_alpha1, spacing_f, _mass_alpha1 = fractional_ou(
alpha=1.0, n_steps=20000, dt=0.005, grid=2048, length=32.0
)
x_f, p_alphahalf, _, _mass_alphahalf = fractional_ou(
alpha=ALPHA, n_steps=20000, dt=0.005, grid=2048, length=32.0
)
gaussian = KAPPA * np.exp(-0.5 * KAPPA * x_f**2) / np.sqrt(2.0 * np.pi * DIFFUSION)
tv_alpha1 = total_variation(p_alpha1, gaussian, spacing_f)
tail = (x_f > 2.5) & (x_f < 6.0) & (p_alphahalf > 1e-14)
slope, intercept = np.polyfit(np.log(x_f[tail]), np.log(p_alphahalf[tail]), 1)
print("Fractional OU, 20000 split steps of 0.005 h, 2048 points on [-16, 16) km")
print(f" alpha = 1: mass drift over the run {_mass_alpha1:.3e}")
print(f" alpha = 0.5: mass drift over the run {_mass_alphahalf:.3e}")
print(f" alpha = 1 TV distance to the exact Gaussian: {tv_alpha1:.4f}")
print(f" alpha = 0.5 tail fit on 2.5-6 km: slope {slope:.3f} (predicted -2.000)")
# --- Nonlocal interaction -------------------------------------------
width = 0.5
interaction_grid = 256
interaction_spacing = 16.0 / interaction_grid
dt = 0.2 * interaction_spacing**2 / DIFFUSION
x_i, p_free, t_free, change_free = interaction_run(
0.0, interaction_grid, 16.0, dt, 400.0, width
)
x_i, p_inter, t_inter, change_inter = interaction_run(
0.6, interaction_grid, 16.0, dt, 400.0, width
)
def variance(p):
mean = float(np.sum(p * x_i))
return float(np.sum(p * x_i**2)) - mean**2
v_free, v_inter = variance(p_free), variance(p_inter)
v_predicted = self_consistent_variance(0.6, width)
gaussian_predicted = np.exp(-0.5 * x_i**2 / v_predicted)
gaussian_predicted /= gaussian_predicted.sum()
tv_prediction = total_variation(p_inter, gaussian_predicted)
tv_equilibria = total_variation(p_inter, p_free)
print(f"\nNonlocal attraction chi = 0.6 km^2/h, interaction width {width} km")
print(f" converged after {t_free:.1f} h (chi = 0) and {t_inter:.1f} h (chi = 0.6)")
print(f" final TV change: {change_free:.2e} and {change_inter:.2e}")
print(f" variance: no interaction {v_free:.3f}, with interaction {v_inter:.3f}")
print(f" self-consistent Gaussian prediction: {v_predicted:.3f}")
print(f" TV distance from that prediction: {tv_prediction:.3f}")
print(f" TV distance between the two equilibria: {tv_equilibria:.3f}")
# --- Behavioural states ----------------------------------------------
grid_s, length_s = 384, 16.0
spacing_s = length_s / grid_s
x_s = -0.5 * length_s + np.arange(grid_s) * spacing_s
dt_s = 0.25 * spacing_s**2 / 0.60
snapshots = switch_run(x_s, spacing_s, dt_s, 100.0, snapshot_at=(60.0, 100.0))
p_a, p_b = snapshots[100.0]
mixture = p_a + p_b
mean = float(np.sum(mixture * x_s))
sd = float(np.sqrt(np.sum(mixture * (x_s - mean) ** 2)))
gaussian_fit = np.exp(-0.5 * ((x_s - mean) / sd) ** 2)
gaussian_fit /= gaussian_fit.sum()
tv_gaussian = total_variation(mixture, gaussian_fit)
tv_stationary = total_variation(
snapshots[60.0][0] + snapshots[60.0][1], mixture
)
midpoint_ratio = float(mixture[np.isclose(x_s, 0.0)][0]) / float(np.max(mixture))
print("\nTwo behavioural states, 100 h of switching (10 h typical per state)")
print(f" mass in state A {np.sum(p_a):.6f}, state B {np.sum(p_b):.6f}")
print(f" change between t = 60 h and t = 100 h: TV {tv_stationary:.2e}")
print(f" mixture mean {mean:.3f} km, sd {sd:.3f} km")
print(f" TV distance from a single Gaussian with the same mean and sd: {tv_gaussian:.3f}")
print(f" density at the midpoint versus the highest mode: {midpoint_ratio:.3f}")
times, path, _states = sample_path(0.005, 100.0)
fraction_wide = float(np.mean(_states == 1))
print(f" example path: {fraction_wide:.3f} of the time in the wide mode")
# --- Figure ----------------------------------------------------------
fig = _plot(
x_f,
p_alpha1,
p_alphahalf,
gaussian,
slope,
intercept,
x_i,
p_free,
p_inter,
x_s,
mixture,
p_a,
p_b,
gaussian_fit,
times,
path,
)
output_path = Path(output)
output_path.parent.mkdir(parents=True, exist_ok=True)
fig.savefig(output_path, dpi=300, bbox_inches="tight", facecolor="white")
plt.close(fig)
return output_path
def main() -> None:
"""Console-script entry point."""
run()
if __name__ == "__main__":
main()
Fractional diffusion gives the moose heavy tails. The stationary density of the fractional OU flow is computed by a split scheme: each step applies the exact Ornstein-Uhlenbeck contraction, remapped conservatively from cell to cell, and then the exact fractional diffusion in Fourier space. As a check, \(\alpha = 1\) must reproduce the ordinary Gaussian, and does, to a total-variation distance of \(0.0037\) over twenty thousand steps with mass drift at the \(10^{-12}\) level. For \(\alpha = 0.5\) the stationary density is no longer Gaussian in the tails:
Fractional OU, 20000 split steps of 0.005 h, 2048 points on [-16, 16) km
alpha = 1: mass drift over the run 9.237e-12
alpha = 0.5: mass drift over the run 2.734e-12
alpha = 1 TV distance to the exact Gaussian: 0.0037
alpha = 0.5 tail fit on 2.5-6 km: slope -1.985 (predicted -2.000)A fitted tail slope of \(-1.985\) against the predicted \(-(1+2\alpha) = -2\): the fractional moose still concentrates near its centre, but its rare multi-kilometre relocations occur at a power-law rate rather than the astronomical improbability the Gaussian assigns them. The numerical point worth keeping is that the flow step must be conservative. On the periodic domain used here the OU contraction points strictly inward, so the exact one-step map folds nothing across the seam; treating it as a wrap-around contraction would manufacture mass on every step, and the scheme would explode.
Nonlocal attraction makes the herd bunch together. With the interaction drift \(\chi \, \partial_x \log(G*p)\), the equilibrium can still be predicted: for a Gaussian \(p\) of variance \(v\), the smoothed density \(G*p\) is Gaussian of variance \(v + w^2\), so \(\log(G*p)\) is exactly quadratic and the self-consistency condition closes,
\[v = \frac{D}{\kappa + \chi/(v + w^2)}.\]
For \(\chi = 0.6\ \text{km}^2/\text{h}\) and interaction width \(w = 0.5\) km its fixed point is \(v = 0.581\), down from the non-interacting value of \(1\). The simulation:
Nonlocal attraction chi = 0.6 km^2/h, interaction width 0.5 km
converged after 9.1 h (chi = 0) and 8.8 h (chi = 0.6)
final TV change: 8.31e-11 and 7.78e-11
variance: no interaction 1.050, with interaction 0.632
self-consistent Gaussian prediction: 0.581
TV distance from that prediction: 0.018
TV distance between the two equilibria: 0.122The interacting equilibrium sits within \(0.018\) of the predicted Gaussian; the residual, and the gap between \(0.632\) and \(0.581\), come from the same upwind numerical diffusion that inflated the free variance to \(1.050\) in the first place. The two equilibria are \(0.122\) apart in total variation: animals that pay attention to each other occupy a visibly tighter range than the same animals alone.
Behavioural states turn one Gaussian into two lobes. The last experiment couples two state equations: state A is a core range (\(\kappa = 1\ \text{h}^{-1}\), \(D = 0.25\ \text{km}^2/\text{h}\), centre \(-3\) km), state B an exploratory mode (\(\kappa = 0.5\ \text{h}^{-1}\), \(D = 0.6\ \text{km}^2/\text{h}\), centre \(+2.5\) km), with symmetric switching at \(0.1\ \text{h}^{-1}\) — a typical ten hours per state:
Two behavioural states, 100 h of switching (10 h typical per state)
mass in state A 0.500000, state B 0.500000
change between t = 60 h and t = 100 h: TV 4.14e-06
mixture mean -0.462 km, sd 2.553 km
TV distance from a single Gaussian with the same mean and sd: 0.340
density at the midpoint versus the highest mode: 0.184
example path: 0.261 of the time in the wide modeBy \(t = 60\) h the mixture is stationary. Its total density has two lobes with a dip at the midpoint to \(0.184\) of the peak height, and it is \(0.340\) away from the best single Gaussian with the same mean and variance — far beyond any discretization error. This is the quantitative form of the argument that made switching state-space models standard in movement ecology (Jonsen et al. 2007; Patterson et al. 2009): a single-state model fitted to such an animal invents a home range centred in the gap between the two lobes, habitat the moose in fact avoids.
What generality does not buy
Repairs have prices, and honesty requires listing them.
- None of this is fitted. Every parameter above was chosen to make the mechanism legible. The OU model at least admitted likelihood fitting through its Kalman filter; the general equation loses the linear-Gaussian structure that made that cheap (Hooten et al. 2017). Model inference, not model construction, is now the bottleneck.
- Identifiability is the next wall. The equilibrium is shaped jointly by \(U\) and \(D\), and from positions alone good habitat is hard to distinguish from easy movement. Multi-lobed densities are produced by multi-modal potentials, by behavioural switching, by heavy tails, and by interactions alike; distinguishing mechanisms will need transients or multi-animal data, not just stationary utilizations.
- Infinite speed is not repaired, only reinterpreted. Fractional diffusion replaces Gaussian steps with explicit jumps, but the operator is nonlocal in space, so influence still propagates instantly. The hard maximum-speed bounds from the first post live outside any Fokker-Planck model and remain complementary to it.
- It is still a stationary framework. Season, snow, and rut enter here only through the choice of behavioural states, not as explicit time-dependence; defect 7 is partially, not fully, answered.
- It costs compute. Two-dimensional grids, nonlocal couplings, and fractional operators turn closed forms into matrix and FFT work. Minutes on a laptop for the toy domains here; realistic landscapes demand more.
Conclusion
The defect list from the Ornstein-Uhlenbeck post served as a specification, and term by term the general Fokker-Planck equation carries what the OU model could not: a habitat potential that makes the home range multi-modal, a diffusion field that makes terrain matter, barriers that gate mixing through single crossings, conspecific attraction that tightens ranges, behavioural states that split one utilization distribution into two, and fractional diffusion that gives rare relocations a power-law rate. The OU model is not discarded — it survives as the special case these terms are switched off into, and as the null model every term is measured against. What is given up is everything that made the OU model computable: the Gaussian, the closed form, the filter. The next serious step is turning from model construction to model fitting, putting real collar fixes on the likelihood of some chosen subset of these terms.