Diagonal Frog: High-order positivity-preserving FD schemes for anisotropic Fokker-Planck equations

2026-06-22

AI Summary

This work addresses the challenge of non-physical negative probabilities arising in conventional discretizations of multidimensional Fokker–Planck equations under strong anisotropic diffusion and jump processes. To overcome this, the authors propose a novel operator-splitting “diagonal leapfrog” scheme that employs matrix exponential Krylov approximations for directional sub-operators and a factored resolvent solver to handle mixed derivative terms. By integrating a Zeno-type infinite subdivision strategy, the method guarantees non-negativity and mass conservation without requiring flux limiters. Although the resulting discrete operator is not a local M-matrix, it exhibits eventual M-matrix properties. The scheme achieves second-order accuracy in both space and time, with computational complexity O(m²N + m³), and demonstrates robustness and efficiency in high-Péclet-number regimes and scenarios featuring strong convection coupled with cross-diffusion.

Abstract

The Fokker-Planck equation is fundamental to statistical mechanics, yet in settings with multiple state variables, anisotropic (cross-) diffusion, and jumps, conventional discretizations frequently produce non-physical negative probability densities. Building on the operator approach of "A. Itkin, Pricing derivatives under Levy models. Modern finite difference and pseudo-differential operators approach, Springer, 2017, ISBN 978-1-4939-6792-6", we introduce a family of "Diagonal Frog" discretizations whose spatial operators are eventually M-matrices (EM-matrices). Although these operators lack a local M-matrix structure, positivity of the directional sub-operators emerges in the spirit of Zeno's paradox: the matrix exponential, assembled as the limit of infinitely many ever-smaller substeps, is provably nonnegative after a short transient even though no single substep is. For the mixed-derivative block, whose generator is not eventually nonnegative, positivity instead rests on a factorized resolvent solver and holds conditionally, on an explicit step-size window; discrete mass is conserved exactly by the splitting for every step size. The resulting schemes are second-order accurate in time and space and require O(m 2 N + m 3) operations per time step, where m is the dimension of the Krylov subspace used to apply the exponential. As stress tests, we solve a two-dimensional anisotropic Fokker-Planck equation in the strong cross-diffusion regime against an exact Gaussian reference, a Kramers escape problem in a double-well potential, and an advection-dominated problem, and observe that the schemes remain stable, nonnegative, and mass-conservative for a wide range of Pécklet numbers (so, don't need any flux limiter). Finally, we extend the construction to multidimensional processes and to the backward Kolmogorov equation with jumps.

Problem

• Fokker-Planck equation
• anisotropic diffusion
• positivity preservation
• negative probability density
• cross-diffusion

Innovation

• eventually M-matrices
• positivity-preserving
• anisotropic Fokker-Planck equation
• matrix exponential splitting
• mass conservation