Numerical stationary states for nonlocal Fokker-Planck equations via fixed points of consistency maps

We propose a fixed-point-based numerical framework for computing stationary states of nonlocal Fokker-Planck-type equations. Instead of discretising the differential operators directly, we reformulate the stationary problem as a nonlinear fixed-point map built from the original PDE and its nonlocal interaction terms, and solve the resulting finite-dimensional problem with a matrix-free Newton-Krylov method. We compare implementations using the analytic Frechet derivative of this map with a simple central-difference approximation. Because the method does not rely on time evolution, it is agnostic to dynamical stability and can detect both stable and unstable stationary states. Its accuracy is determined mainly by the numerical treatment of convolutions and quadrature, rather than by differentiation stencils. We apply the approach to three model problems with linear diffusion, use existing analytical results to verify the outputs, and reproduce known bifurcation diagrams, as well as new bifurcation behaviour not previously observed in this kind of problem.

💡 Research Summary

This paper introduces a novel numerical framework for computing stationary states of nonlocal Fokker‑Planck‑type equations without discretising differential operators directly. By setting the time derivative to zero in the governing PDE, the authors derive a nonlinear fixed‑point map T that encapsulates the original equation’s diffusion, external potential, and nonlocal interaction terms. Stationary solutions are then identified as fixed points of T, i.e., solutions of F(u)=T(u)−u=0.

To solve the high‑dimensional nonlinear system, the authors employ a Jacobian‑free Newton‑Krylov (JFNK) method. The outer Newton iteration updates the iterate u_k by solving a linear system J(u_k) δu_k = −F(u_k) with an inner Krylov subspace solver (GMRES). Crucially, the Jacobian J = I − DT(u) can be expressed analytically for the considered models because the inverse of the entropy derivative (H′)⁻¹ is explicit (exponential for linear diffusion). This analytic Fréchet derivative yields a well‑conditioned operator that is essentially the identity plus a compact perturbation, allowing rapid convergence without any preconditioning. For comparison, the authors also implement a finite‑difference approximation of the Jacobian, which is shown to be slower and less accurate, especially near boundaries where the nonlocal kernel is poorly resolved.

Three representative models are examined: (i) the McKean‑Vlasov equation with Gaussian interaction kernels, (ii) a Cucker‑Smale flocking model featuring distance‑dependent alignment, and (iii) a neural Fokker‑Planck equation describing population‑level firing‑rate dynamics. All three use linear diffusion (H