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

TL;DR

The paper addresses the challenge of computing stationary states for nonlocal Fokker–Planck-type equations by recasting the stationary problem as a fixed-point map $\mathcal{T}$ and solving $F(u)=\mathcal{T}u-u=0$ with a matrix-free Newton–Krylov method, leveraging $H'$-inversion and mass conservation via $c_0(u)$. The approach relies on quadrature and FFT-based convolutions for nonlocal terms and allows analytic Fréchet derivatives when available, with central-difference approximations as a fallback, enabling efficient, Jacobian-aware solves without time evolution. The method is validated on three linear-diffusion nonlocal models (McKean–Vlasov, Cucker–Smale, neural FP), reproduces known bifurcation diagrams, reveals new bifurcation behavior, and extends naturally to two spatial dimensions, including analysis of the role of the initial iterate. By enabling access to unstable stationary states and providing a robust, model-agnostic framework, the work offers a powerful tool for exploring stationary patterns and bifurcations in nonlocal interacting systems. The practical impact lies in a scalable, accurate numerical pipeline that can handle diverse kernels, domains, and boundary conditions while delivering high-resolution bifurcation information beyond time-stepping limitations.

Abstract

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.

Figures (8)

• Figure 1: Error plots for the numerical solutions obtained to the McKean-Vlasov equation \ref{['eq:McKean_Vlasov_stationary']} with kernel $W = - w_k(x)$ as developed in Section \ref{['sec:kuramoto_model_1']} at the fixed value $\kappa = 3$. For the left panel, we use 2,017 logarithmically spaced (odd) values of $N_n \in [11, 10001]$. Note that due to the logarithmic spacing, the first values of $N_n$ will be consecutive (i.e $N_1 = 11, N_2 = 13, N_3 = 15$). For the centre and right panels, we employ 160 logarithmically spaced (odd) values of $N_n \in [11, 10001]$ and $h \in [10^{-15}, 10^{-1}]$. As an initial guess for Algorithm \ref{['algorithm:jfnk-fixed-points_intro']}, $f(x) = 1/\pi + \cos(2x)$ was used. Any white points in the heatmaps correspond to instances in which Algorithm \ref{['algorithm:jfnk-fixed-points_intro']} failed to converge when using the central-differences approximation of the Fréchet derivative.
• Figure 2: Top-left panel: bifurcation diagram produced using our algorithm. Top-right panel: the two solution profiles (homogeneous and single nontrivial solution) at $\kappa = 3$. Bottom-left panel: the error between the bifurcation diagrams obtained using the analytical reference solution and our numerical approximation. Bottom-right panel: the error between the critical analytical threshold $\kappa^*$ and that predicted by our numerical approximations, across 160 logarithmically spaced (odd) values of $N_n \in [11, 10001]$.
• Figure 3: Bifurcation curves (left panels) and solution profiles (right panels) for the two-mode kernel defined in \ref{['eq:two-mode-kernel']} for problem \ref{['eq:McKean_Vlasov_stationary']}. We consider $2001$ evenly spaced values of $\kappa \in [2,3]$, and $2001$ evenly spaced values of the domain $\mathbb{T} = (-\pi/2, \pi/2]$. We report the average and standard deviation of the run times, across 10 runs.
• Figure 4: Bifurcation curves obtained for three different instances of problem \ref{['eq:McKean_Vlasov_stationary']}. We consider $2001$ evenly spaced values of the domain $\mathbb{T} = (-\pi/2, \pi/2]$ for (A), (B), and $4001$ samples for (C).
• Figure 5: The basins of attraction for initial guesses of the form \ref{['eq:initial-guesses-suite']} when fed into Algorithm \ref{['algorithm:jfnk-fixed-points_intro']}. Note carefully that these are not related to the temporal dynamics of problem \ref{['eq:time_pde_general']}. The suite of initial guesses is parametrised by the wavenumber $k$ and the amplitude $b$ as defined in \ref{['eq:initial-guesses-suite']}. The heat map indicates which branch the corresponding initial iterate converged to, where the basin colours correspond directly with those branches displayed in Figure \ref{['fig:supercritical_bifurcation']} (e.g., red regions correspond to those initial guesses which yield Branch 5, the homogeneous solution, in Figure \ref{['fig:supercritical_bifurcation']}).