Fokker-Planck Equation

The Fokker-Planck equation (FPE) governs the time evolution of the probability density of a stochastic process. It is the deterministic counterpart to the Langevin equation -- tracking the distribution rather than individual realisations. We derive the Chang-Cooper scheme for numerically solving the FPE with guaranteed positivity and correct steady state.

1. Derivation from the Langevin Equation

Given the Itô SDE:

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

the probability density $P(x, t)$ satisfies the Fokker-Planck equation:

$$\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^2(x) P\bigr]$$

Derivation via Characteristic Function

Define the characteristic function $\phi(k,t) = \langle e^{ikx_t} \rangle$. By Itô's lemma with $f(x) = e^{ikx}$:

$$d(e^{ikx}) = \left(ika \, e^{ikx} - \frac{k^2 b^2}{2} e^{ikx}\right)dt + ikb \, e^{ikx} \, dW$$

Taking expectations (the $dW$ term vanishes):

$$\frac{\partial \phi}{\partial t} = \langle ik \, a(x) e^{ikx} \rangle - \frac{k^2}{2}\langle b^2(x) e^{ikx} \rangle$$

Inverse Fourier transforming and using integration by parts twice yields the Fokker-Planck equation. The first term $-\partial(aP)/\partial x$ is the drift (probability current), and the second $\frac{1}{2}\partial^2(b^2 P)/\partial x^2$ is the diffusion.

Written in conservation form with probability current $J$:

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

2. Stationary Solution

Setting $\partial P / \partial t = 0$ and $J = 0$ (detailed balance):

$$a(x) P_{\text{st}} = \frac{1}{2}\frac{d}{dx}\bigl[b^2(x) P_{\text{st}}\bigr]$$

This is a first-order ODE. Expanding the right side:

$$a \, P_{\text{st}} = \frac{1}{2}\bigl[2b b' P_{\text{st}} + b^2 P_{\text{st}}'\bigr]$$

Rearranging:

$$\frac{P_{\text{st}}'}{P_{\text{st}}} = \frac{2a - 2bb'}{b^2} = \frac{2a}{b^2} - \frac{2b'}{b}$$

Integrating:

$$\boxed{P_{\text{st}}(x) = \frac{\mathcal{N}}{b^2(x)} \exp\!\left(2\int^x \frac{a(x')}{b^2(x')}\,dx'\right)}$$

where $\mathcal{N}$ is the normalisation constant. For constant diffusion $b = \sigma$ and gradient drift $a(x) = -V'(x)$, this reduces to the Boltzmann distribution $P_{\text{st}} \propto \exp(-2V/\sigma^2)$.

3. Chang-Cooper Numerical Scheme

Standard finite-difference schemes for the FPE can produce negative probabilities or fail to converge to the correct steady state. The Chang-Cooper scheme (1970) cures both issues by interpolating between upwind and downwind differencing.

Discretise the domain $[x_{\min}, x_{\max}]$ into $N_x$ cells with spacing $\Delta x$. The discrete probability current at cell boundary $j+1/2$ is:

$$J_{j+1/2} = a_{j+1/2}\bigl[(1-\delta_j)P_j + \delta_j P_{j+1}\bigr] - \frac{D_{j+1/2}}{\Delta x}(P_{j+1} - P_j)$$

where $D_{j+1/2} = b^2_{j+1/2}/2$ and the Chang-Cooper parameter $\delta_j$ is chosen so that the scheme exactly preserves the stationary distribution. Define:

$$w_j = \frac{a_{j+1/2} \, \Delta x}{D_{j+1/2}}$$

Then the optimal interpolation parameter is:

$$\delta_j = \frac{1}{w_j} - \frac{1}{e^{w_j} - 1}$$

This is the Bernoulli function. When $w \to 0$,$\delta \to 1/2$ (central differencing); when $|w| \gg 1$, it approaches upwind differencing. The time update uses the tridiagonal system:

$$\frac{P_j^{n+1} - P_j^n}{\Delta t} = -\frac{J_{j+1/2} - J_{j-1/2}}{\Delta x}$$

4. Python: Chang-Cooper FPE Solver

We solve the Fokker-Planck equation for the bistable potential and compare the evolving distribution to histograms from Langevin simulations.

Chang-Cooper Fokker-Planck Solver

Python

script.py163 lines

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import solve_banded

# --- Fokker-Planck for bistable: a(x) = x - x^3, b(x) = sigma ---
sigma = 0.5
D_const = sigma**2 / 2.0

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

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

# --- Grid ---
x_min, x_max = -2.5, 2.5
Nx = 200
dx = (x_max - x_min) / Nx
x = np.linspace(x_min + dx/2, x_max - dx/2, Nx)  # cell centres
x_half = np.linspace(x_min, x_max, Nx + 1)         # cell boundaries

# --- Chang-Cooper parameter ---
a_half = drift(x_half)  # drift at boundaries
D_half = D_const * np.ones(Nx + 1)

w = a_half * dx / D_half
# Bernoulli function: delta = 1/w - 1/(exp(w)-1), careful for small w
delta = np.zeros_like(w)
for i in range(len(w)):
    if abs(w[i]) < 1e-6:
        delta[i] = 0.5  # limit
    else:
        delta[i] = 1.0/w[i] - 1.0/(np.exp(w[i]) - 1.0)

# --- Build tridiagonal matrix for implicit Euler ---
dt = 0.005
T_final = 5.0
N_steps = int(T_final / dt)

# Initial condition: Gaussian near x=0
P = np.exp(-x**2 / (2 * 0.1**2))
P /= np.sum(P) * dx

# Store snapshots
snapshots = [(0, P.copy())]
snap_times = [0.5, 1.0, 2.0, 5.0]

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

# Build tridiagonal: A P^{n+1} = P^n
    # J_{j+1/2} = a_{j+1/2}[(1-delta_j)*P_j + delta_j*P_{j+1}] - D_{j+1/2}/dx*(P_{j+1}-P_j)
    # dP_j/dt = -(J_{j+1/2} - J_{j-1/2})/dx

# Coefficients for P_j in J_{j+1/2}: C_j = a_{j+1/2}(1-delta_j) + D_{j+1/2}/dx
    # Coefficients for P_{j+1} in J_{j+1/2}: F_j = a_{j+1/2}*delta_j - D_{j+1/2}/dx

C = a_half * (1 - delta) + D_half / dx   # size Nx+1
    F = a_half * delta - D_half / dx           # size Nx+1

# Sub-diagonal (a_i), diagonal (b_i), super-diagonal (c_i)
    diag = np.ones(Nx)
    lower = np.zeros(Nx)
    upper = np.zeros(Nx)

for j in range(Nx):
        # From J_{j+1/2}: contribution to dP_j/dt is -J_{j+1/2}/dx
        #   P_j term: -C[j+1]/dx (note j+1/2 -> index j+1 in x_half? No: j+1/2 boundary is x_half[j+1])
        # Actually: J at right boundary of cell j is J_{j+1/2} = boundary index j+1 in x_half
        # J at left boundary of cell j is J_{j-1/2} = boundary index j in x_half

# Right boundary (index j+1 in x_half):
        # J_{right} = C[j+1]*P_j + F[j+1]*P_{j+1}  (if j+1 < Nx)
        # Left boundary (index j in x_half):
        # J_{left} = C[j]*P_{j-1} + F[j]*P_j  (if j-1 >= 0)

# dP_j/dt = -(J_right - J_left)/dx

# Implicit: (1 + dt/dx * terms) P^{n+1} = P^n
        coeff_j = 0.0
        coeff_jm1 = 0.0
        coeff_jp1 = 0.0

# Right boundary contribution: -C[j+1]*P_j/dx - F[j+1]*P_{j+1}/dx
        if j < Nx - 1:
            coeff_j += C[j+1] / dx
            coeff_jp1 += F[j+1] / dx
        # Left boundary contribution: +C[j]*P_{j-1}/dx + F[j]*P_j/dx
        if j > 0:
            coeff_j -= F[j] / dx
            coeff_jm1 -= C[j] / dx

diag[j] = 1.0 + dt * coeff_j
        if j > 0:
            lower[j] = dt * coeff_jm1
        if j < Nx - 1:
            upper[j] = dt * coeff_jp1

# Solve tridiagonal system
    ab = np.zeros((3, Nx))
    ab[0, 1:] = upper[:-1]     # super-diagonal
    ab[1, :] = diag             # diagonal
    ab[2, :-1] = lower[1:]     # sub-diagonal

P = solve_banded((1, 1), ab, P)
    P = np.maximum(P, 0)
    P /= np.sum(P) * dx

for st in snap_times:
        if abs(t - st) < dt/2:
            snapshots.append((t, P.copy()))

# --- Also run Langevin for comparison ---
np.random.seed(42)
N_lang = 20000
x_lang = 0.01 * np.random.randn(N_lang)
dt_lang = 0.01
for _ in range(int(T_final / dt_lang)):
    dW = np.sqrt(dt_lang) * np.random.randn(N_lang)
    x_lang += drift(x_lang) * dt_lang + sigma * dW

# --- Plot ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))

# Left: FPE evolution
ax = axes[0]
colors = ['#94a3b8', '#06b6d4', '#10b981', '#f59e0b', '#ef4444']
for i, (t_snap, P_snap) in enumerate(snapshots):
    ax.plot(x, P_snap, color=colors[i % len(colors)], linewidth=2, label=f't = {t_snap:.1f}')
# Theoretical steady state
P_th = np.exp(-2 * potential(x) / sigma**2)
P_th /= np.sum(P_th) * dx
ax.plot(x, P_th, 'w--', linewidth=1.5, label='Steady state (theory)')
ax.set_xlabel('x', fontsize=12, color='white')
ax.set_ylabel('P(x, t)', fontsize=12, color='white')
ax.set_title('Fokker-Planck: Chang-Cooper Evolution', fontsize=13, color='white')
ax.legend(fontsize=9, facecolor='#1e293b', edgecolor='#10b981', labelcolor='white')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2)

# Right: FPE vs Langevin histogram
ax = axes[1]
ax.hist(x_lang, bins=100, density=True, color='#10b981', alpha=0.5, edgecolor='#064e3b', label='Langevin histogram')
ax.plot(x, snapshots[-1][1], 'w-', linewidth=2.5, label='FPE (Chang-Cooper)')
ax.plot(x, P_th, '--', color='#f59e0b', linewidth=2, label='Analytical $P_{st}$')
ax.set_xlabel('x', fontsize=12, color='white')
ax.set_ylabel('Probability density', fontsize=12, color='white')
ax.set_title('FPE vs Langevin at Steady State', fontsize=13, 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 in axes:
    ax.set_facecolor('#0f172a')
plt.tight_layout()
plt.savefig('output.png', dpi=130, bbox_inches='tight', facecolor='#0f172a')
plt.show()

# Verify positivity
print(f"Min P value: {snapshots[-1][1].min():.2e} (should be >= 0)")
print(f"Total probability: {np.sum(snapshots[-1][1]) * dx:.6f} (should be 1.0)")
print(f"L2 error vs theory: {np.sqrt(np.sum((snapshots[-1][1] - P_th)**2) * dx):.4e}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

5. Fortran: Chang-Cooper with Tridiagonal System

For large-scale problems, a Fortran implementation with the Thomas algorithm provides optimal performance. The tridiagonal solver runs in $O(N)$ operations.

Thomas Algorithm (Tridiagonal Solver) for Chang-Cooper FPE

Python

script.py101 lines

# We demonstrate the Thomas algorithm (tridiagonal solver) in Python
# mirroring what the Fortran implementation does

import numpy as np
import matplotlib.pyplot as plt

def thomas_solve(a, b, c, d):
    """Solve tridiagonal system Ax = d using Thomas algorithm.
    a: sub-diagonal (a[0] unused), b: diagonal, c: super-diagonal (c[-1] unused)
    """
    n = len(d)
    c_ = np.zeros(n)
    d_ = np.zeros(n)

c_[0] = c[0] / b[0]
    d_[0] = d[0] / b[0]

for i in range(1, n):
        m = a[i] / (b[i] - a[i] * c_[i-1])
        c_[i] = c[i] / (b[i] - a[i] * c_[i-1]) if i < n-1 else 0
        d_[i] = (d[i] - a[i] * d_[i-1]) / (b[i] - a[i] * c_[i-1])

x = np.zeros(n)
    x[-1] = d_[-1]
    for i in range(n-2, -1, -1):
        x[i] = d_[i] - c_[i] * x[i+1]
    return x

# --- FPE parameters ---
sigma = 0.5
D = sigma**2 / 2
Nx = 300
x_min, x_max = -2.5, 2.5
dx = (x_max - x_min) / Nx
x_c = np.linspace(x_min + dx/2, x_max - dx/2, Nx)
x_b = np.linspace(x_min, x_max, Nx + 1)  # boundaries

drift_b = x_b - x_b**3  # drift at boundaries
D_b = D * np.ones(Nx + 1)

# Chang-Cooper delta
w = drift_b * dx / D_b
delta_cc = np.where(np.abs(w) < 1e-8, 0.5, 1.0/w - 1.0/(np.exp(w) - 1.0))

# Initial: narrow Gaussian
P = np.exp(-x_c**2 / (2 * 0.05**2))
P /= np.sum(P) * dx

dt = 0.002
T_total = 8.0
N_steps = int(T_total / dt)

# Build and solve tridiagonal at each step using Thomas algorithm
for step in range(N_steps):
    sub = np.zeros(Nx)  # a (sub-diagonal)
    dia = np.ones(Nx)   # b (main diagonal)
    sup = np.zeros(Nx)  # c (super-diagonal)

for j in range(Nx):
        C_r = drift_b[j+1] * (1 - delta_cc[j+1]) + D_b[j+1] / dx
        F_r = drift_b[j+1] * delta_cc[j+1] - D_b[j+1] / dx
        C_l = drift_b[j] * (1 - delta_cc[j]) + D_b[j] / dx
        F_l = drift_b[j] * delta_cc[j] - D_b[j] / dx

coeff_j = C_r / dx - F_l / dx
        dia[j] = 1.0 + dt * coeff_j

if j > 0:
            sub[j] = -dt * C_l / dx
        if j < Nx - 1:
            sup[j] = dt * F_r / dx

P = thomas_solve(sub, dia, sup, P)
    P = np.maximum(P, 0)
    P /= np.sum(P) * dx

# Theoretical steady state
V = x_c**4/4 - x_c**2/2
P_th = np.exp(-2 * V / sigma**2)
P_th /= np.sum(P_th) * dx

fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(x_c, P, color='#10b981', linewidth=2.5, label='Chang-Cooper (Thomas alg.)')
ax.plot(x_c, P_th, 'w--', linewidth=2, label='Analytical $P_{st}$')
ax.fill_between(x_c, P, alpha=0.15, color='#10b981')
ax.set_xlabel('x', fontsize=13, color='white')
ax.set_ylabel('P(x)', fontsize=13, color='white')
ax.set_title('Fokker-Planck Steady State via Thomas Algorithm', fontsize=14, color='white')
ax.legend(fontsize=11, facecolor='#1e293b', edgecolor='#10b981', labelcolor='white')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2)
fig.patch.set_facecolor('#0f172a')
ax.set_facecolor('#0f172a')
plt.tight_layout()
plt.savefig('output.png', dpi=130, bbox_inches='tight', facecolor='#0f172a')
plt.show()

L2_err = np.sqrt(np.sum((P - P_th)**2) * dx)
print(f"L2 error after {N_steps} steps: {L2_err:.4e}")
print(f"Min P: {P.min():.2e}, Max P: {P.max():.4f}")
print(f"Thomas algorithm: O(N) = O({Nx}) per timestep")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Langevin Equation Schrödinger Bridge →