pith. sign in

arxiv: 2606.23980 · v1 · pith:7IA6TWMZnew · submitted 2026-06-22 · 🧮 math.NA · cs.NA· physics.bio-ph· physics.comp-ph· q-fin.CP· q-fin.PR

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

Andrey Itkin This is my paper

Pith reviewed 2026-06-26 07:14 UTC · model grok-4.3

classification 🧮 math.NA cs.NAphysics.bio-phphysics.comp-phq-fin.CPq-fin.PR
keywords Fokker-Planck equationpositivity-preserving schemesfinite difference methodsanisotropic diffusionM-matricesmatrix exponentialnumerical stability
0
0 comments X

The pith

Diagonal Frog finite-difference schemes keep solutions to anisotropic Fokker-Planck equations nonnegative and mass-conserving without flux limiters.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper develops a family of finite-difference schemes called Diagonal Frog for Fokker-Planck equations that feature anisotropic cross-diffusion and jumps. These schemes construct spatial operators that are eventually M-matrices so that directional parts become nonnegative through the matrix exponential taken as the limit of infinitely many substeps. A sympathetic reader would care because standard discretizations often generate unphysical negative probability densities in multi-variable settings, while the new methods remain stable and exactly mass-conserving for wide ranges of Péclet numbers.

Core claim

The Diagonal Frog discretizations build spatial operators that are eventually M-matrices. Positivity of the directional sub-operators emerges via the matrix exponential assembled as the limit of infinitely many ever-smaller substeps, even though no single substep is nonnegative. For the mixed-derivative block positivity instead rests on a factorized resolvent solver and holds conditionally on an explicit step-size window. The resulting schemes are second-order accurate in time and space, conserve mass exactly by the splitting for every step size, and require O(m² N + m³) operations per time step where m is the Krylov subspace dimension.

What carries the argument

The eventually M-matrix (EM-matrix) spatial operators, whose positivity for directional sub-operators emerges from the matrix exponential limit of infinitely many substeps.

If this is right

  • The schemes remain stable, nonnegative, and mass-conservative for a wide range of Péclet numbers without any flux limiter.
  • Discrete mass is conserved exactly by the splitting for every step size.
  • The construction extends directly to multidimensional processes and to the backward Kolmogorov equation with jumps.
  • Computational cost per time step scales as O(m² N + m³) with m the dimension of the Krylov subspace.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The same splitting and exponential-limit idea could be tested on other transport equations that lose positivity under strong cross terms.
  • Adaptive time-step control might be used to stay inside the conditional window for the mixed-derivative block.
  • Extension to processes with jumps raises the question of whether the same EM-matrix property survives after adding jump operators.

Load-bearing premise

Positivity of the directional sub-operators emerges via the matrix exponential as the limit of infinitely many ever-smaller substeps, while mixed-derivative positivity holds only inside an explicit step-size window.

What would settle it

A numerical test on the two-dimensional anisotropic Fokker-Planck equation against the exact Gaussian reference that produces negative density values for any step size inside the reported stability window.

Figures

Figures reproduced from arXiv: 2606.23980 by Andrey Itkin.

Figure 1
Figure 1. Figure 1: Minimum of e ∆tLp0 versus ∆t (κ = 6, σ = 0.5, max Pe ≈ 7.7). The DF propagator (solid) is nonnegative for ∆t ≥ τ0 ≈ 10−4 – an inverted CFL condition (Corollary 2); the centred scheme (dashed) recovers positivity only near ∆t ≈ 0.2 [PITH_FULL_IMAGE:figures/full_fig_p031_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: OU density after T = 0.3 in the advection-dominated regime (max Pe ≈ 7.7, ∆t = 0.01). At this fixed mesh the centred scheme (dashed) has a small negative undershoot on the leading flank; DF (solid) remains nonnegative. The dotted curve is the initial datum. Relaxation to the stationary law. Finally we verify long-time behaviour. Starting from an off-centre narrow Gaussian, [PITH_FULL_IMAGE:figures/full_fi… view at source ↗
Figure 3
Figure 3. Figure 3: Relaxation of the DF solution to the stationary OU density (κ = 2, σ = 1, ∆t = 0.05). The ℓ1 error to p∞ follows the envelope e −κt (dashed); positivity and mass are preserved throughout. Page 32 of 64 [PITH_FULL_IMAGE:figures/full_fig_p032_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: shows the result. No linear second-order scheme can be monotone here – this is exactly Godunov’s barrier (Remark 2) – so the Diagonal-Frog scheme is not positive in this extreme: it carries a small Godunov ripple, minx p ≈ −7 × 10−3 . The distinction is in the character of the failure. The centred scheme (and Crank–Nicolson on it) produce classical dispersive ringing: a deep negative sink minx p ≈ −0.34 – … view at source ↗
Figure 5
Figure 5. Figure 5: Density p(x, T) for the Kramers escape problem with double-well potential Eq. (51): κ = 8, m = 1, σ = 0.2, x0 = −1.9, T = 0.04. The cubic drift produces a local cell Péclet number Pemax = µ(−1.9)h/D ≈ 198 on the outer flank traversed by the pulse. The centred scheme (dashed) and Crank–Nicolson (green) develop non-physical oscillations there (minx p ≈ −3.67); the Diagonal-Frog scheme (blue) remains nonnegat… view at source ↗
Figure 6
Figure 6. Figure 6: Positivity window Θ of the trapezoidal central substep Φxy(∆t) of Eq. (35) on a resolved Gaussian (ρ = 0.8), versus mesh n. The window grows as the datum is resolved; the dashed line is the reference slope Θ ∝ h 2 . Positivity is therefore conditional but benign for resolved solutions, consistent with Proposition 11 (a) [PITH_FULL_IMAGE:figures/full_fig_p038_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: Most negative value | minn p n| over a run from a steep Gaussian (n = 64, ρ = 0.8, T = 0.4), versus ∆t, for the central factor realised three ways. The bare exponential e ∆tAxy is catastrophic and ∆t-independent (∼ 104 ), reflecting the ill-posed sub-flow of Remark 8. The implicit trapezoidal factor Φxy on its own is far better and improves steeply with ∆t. The full Strang step Eq. (34) stays nonnegative t… view at source ↗
Figure 8
Figure 8. Figure 8: Time-dependent anisotropic Fokker–Planck equation Eq. (56). Left: terminal density p(x, y, T) in the advection-dominated Regime II. Centre: error under joint refinement ∆t ∼ h in the strong-coupling Regime I; the Diagonal Frog scheme and the unsplit BDF2 reference both approach slope ≈ 2, with comparable constants ( [PITH_FULL_IMAGE:figures/full_fig_p043_8.png] view at source ↗
read the original 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\'ecklet numbers (so, don't need any flux limiter). Finally, we extend the construction to multidimensional processes and to the backward Kolmogorov equation with jumps.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

2 major / 1 minor

Summary. The manuscript proposes 'Diagonal Frog' finite-difference schemes for anisotropic Fokker-Planck equations in multiple dimensions. Directional sub-operators are constructed as eventually M-matrices whose positivity is obtained in the limit of the matrix exponential (Zeno construction of infinitely many infinitesimal substeps). Mixed-derivative blocks use a factorized resolvent solver whose positivity holds inside an explicit step-size window. The schemes are stated to be second-order accurate in space and time, exactly mass-conserving for any step size, and positivity-preserving without flux limiters across a wide range of Péclet numbers; this is supported by tests against exact Gaussian solutions, a Kramers escape problem, and an advection-dominated case. The construction is extended to multidimensional processes and to the backward Kolmogorov equation with jumps.

Significance. If the positivity claims can be shown to hold without hidden step-size restrictions that tighten with Péclet number, the work would provide a useful addition to structure-preserving discretizations for Fokker-Planck and related kinetic equations. Exact mass conservation by the splitting and the numerical comparisons to closed-form Gaussian solutions are concrete strengths. The approach builds on the author's earlier operator framework but introduces new splitting and eventual-nonnegativity arguments.

major comments (2)
  1. [Abstract and mixed-derivative block] Abstract and mixed-derivative block: the central claim that the schemes remain nonnegative for a wide range of Péclet numbers without flux limiters rests on the factorized resolvent solver for the mixed-derivative term, which is positivity-preserving only inside an explicit step-size window. No derivation is supplied showing that this window remains non-restrictive when cross-diffusion or advection strength increases, and the manuscript does not verify whether the time steps chosen for the strong-cross-diffusion and advection-dominated numerical examples lie inside the window. This directly affects the 'no limiter needed' assertion.
  2. [Directional sub-operators] Directional sub-operators and implementation: the nonnegativity argument via the matrix-exponential limit of infinitely many substeps is presented, yet the manuscript supplies neither a complete derivation of the transient behavior for finite Krylov-subspace approximations nor an error analysis that quantifies how many substeps are required in practice to reach the nonnegative regime.
minor comments (1)
  1. [Abstract] Abstract: 'P\'ecklet' is a typographical error for 'Péclet'.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for the careful reading and constructive comments. We address each major comment below, indicating revisions where appropriate to strengthen the presentation of the positivity results.

read point-by-point responses

  1. Referee: [Abstract and mixed-derivative block] Abstract and mixed-derivative block: the central claim that the schemes remain nonnegative for a wide range of Péclet numbers without flux limiters rests on the factorized resolvent solver for the mixed-derivative term, which is positivity-preserving only inside an explicit step-size window. No derivation is supplied showing that this window remains non-restrictive when cross-diffusion or advection strength increases, and the manuscript does not verify whether the time steps chosen for the strong-cross-diffusion and advection-dominated numerical examples lie inside the window. This directly affects the 'no limiter needed' assertion.

    Authors: We agree that the manuscript would benefit from an explicit derivation bounding the step-size window in terms of the cross-diffusion and advection coefficients, together with verification that the time steps used in the numerical examples satisfy the positivity condition. The abstract already notes the conditional nature of the result, but additional analysis will be added in the revision to support the claim of applicability across a wide range of Péclet numbers without limiters. revision: yes

  2. Referee: [Directional sub-operators] Directional sub-operators and implementation: the nonnegativity argument via the matrix-exponential limit of infinitely many substeps is presented, yet the manuscript supplies neither a complete derivation of the transient behavior for finite Krylov-subspace approximations nor an error analysis that quantifies how many substeps are required in practice.

    Authors: The nonnegativity proof applies to the exact matrix exponential in the infinite-substep (Zeno) limit. The Krylov method with dimension m is used only as a practical means to evaluate the exponential action, with accuracy controlled by standard Krylov error bounds. We acknowledge that a dedicated discussion of transient behavior under finite-m approximations and practical selection of m would improve the implementation section. This material will be incorporated in the revision. revision: yes

Circularity Check

0 steps flagged

Minor self-citation to 2017 book; new splitting and positivity arguments do not reduce to prior results

full rationale

The derivation relies on the 2017 book only for the general operator approach. The core constructions—Diagonal Frog discretizations, eventually M-matrices, Zeno-paradox limit for directional sub-operators, and factorized resolvent solver for the mixed block—are presented as new mathematical arguments inside the paper. No step equates a claimed prediction or positivity result to a fitted parameter or to the cited book by construction. Numerical experiments against an exact Gaussian reference supply independent external checks. This is the normal, non-circular case of building on prior technique while adding self-contained analysis.

Axiom & Free-Parameter Ledger

0 free parameters · 2 axioms · 0 invented entities

Based on abstract only; relies on standard properties of matrix exponentials and M-matrices for eventual nonnegativity.

axioms (2)
  • standard math Matrix exponential of eventually nonnegative operators yields nonnegative result after transient
    Invoked to establish positivity via infinite substeps limit.
  • domain assumption Factorized resolvent solver preserves positivity inside explicit step-size window
    Stated for the mixed-derivative block.

pith-pipeline@v0.9.1-grok · 5876 in / 1310 out tokens · 24248 ms · 2026-06-26T07:14:16.867411+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

21 extracted references

  1. [1]

    Comment on “On the numerical simulation of particle dynamics in the radiation belt. Part I: Implicit and semi-implicit schemes

    Albert, J. M. (2013). “Comment on “On the numerical simulation of particle dynamics in the radiation belt. Part I: Implicit and semi-implicit schemes” and “On the numerical simulation of particle dynamics in the radiation belt. Part II: Procedure based on the diagonalization of the diffusion tensor” by E. Camporeale et al.” In:Journal of Geophysical Resea...

    2013

  2. [2]

    The Magnus Expansion and Some of Its Applications

    Classics in Applied Mathematics. Philadelphia, PA: Society for Industrial and Applied Mathematics.isbn: 978-0-89871-321-3. Blanes, S. et al. (2009). “The Magnus Expansion and Some of Its Applications”. In:Physics Reports470, 5–6, pp. 151–238. Bolley, C. and M. Crouzeix (1978). “Conservation de la positivité lors de la discrétisation des problèmes d’évolut...

    2009

  3. [3]

    Berlin: Springer

    Springer Series in Computational Mathematics. Berlin: Springer. Hairer, E., S. P. Nørsett, and G. Wanner (1993).Solving Ordinary Differential Equations I: Nonstiff Problems. 2nd. Vol

    1993

  4. [4]

    High resolution schemes for hyperbolic conservation laws

    Springer Series in Computational Mathematics. Berlin: Springer. Harten, A. (1983). “High resolution schemes for hyperbolic conservation laws”. In:Journal of Computational Physics49, 3, pp. 357–393. Higham, N. J. (2009). “The Scaling and Squaring Method for the Matrix Exponential Revisited”. In:SIAM Review51, 4, pp. 747–764. Hirsa, A. (2013).Computational ...

    1983

  5. [5]

    (2011).Options, Futures, and Other Derivatives

    Hull, J. (2011).Options, Futures, and Other Derivatives. 8th. Upper Saddle River, NJ: Prentice Hall. Hundsdorfer, W. and J. G. Verwer (2003).Numerical Solution of Time-Dependent Advection-Diffusion- Reaction Equations. Vol

    2011

  6. [6]

    Componentwise splitting methods for pricing American options under stochastic volatility

    Springer Series in Computational Mathematics. Berlin: Springer. Ikonen, S. and J. Toivanen (2007). “Componentwise splitting methods for pricing American options under stochastic volatility”. In:Int. J. Theor. Appl. Finance10, pp. 331–361. Ikonen, S. and J. Toivanen (2008). “Efficient numerical methods for pricing American options under stochastic volatili...

    2007

  7. [7]

    Jumps without Tears: A New Splitting Technology for Barrier Options

    Basel: Birkhauser. Itkin, A. and P. Carr (2011). “Jumps without Tears: A New Splitting Technology for Barrier Options”. In: International Journal of Numerical Analysis and Modeling8, 4, pp. 667–704. Itkin, A. and A. Lipton (2015). “Efficient solution of structural default models with correlated jumps and mutual obligations”. In:International Journal of Co...

    2011

  8. [8]

    On the Exponential Solution of Differential Equations for a Linear Operator

    Magnus, W. (1954). “On the Exponential Solution of Differential Equations for a Linear Operator”. In: Communications on Pure and Applied Mathematics7, 4, pp. 649–673. McFarland, D. M. et al. (2025). “Efficient solution of Fokker–Planck equations in two dimensions”. In: Mathematics13, 3, p

    1954

  9. [9]

    On Perron–Frobenius property of matrices having some negative entries

    Noutsos, D. (2006). “On Perron–Frobenius property of matrices having some negative entries”. In:Linear Algebra and its Applications412, 2–3, pp. 132–153. Page 49 of 64 Positivity-preserving FD schemes Noutsos, D. and M. J. Tsatsomeros (2008). “Reachability and holdability of nonnegative states”. In:SIAM Journal on Matrix Analysis and Applications30, 2, pp...

    2006

  10. [10]

    On the construction and comparison of difference schemes

    Springer Series in Synergetics. Berlin, Heidelberg: Springer-Verlag. Strang, G. (1968). “On the construction and comparison of difference schemes”. In:SIAM J. Numerical Analysis5, pp. 509–517. Thomas, J. W. (1995).Numerical Partial Differential Equations: Finite Difference Methods. Vol

    1968

  11. [11]

    A componentwise splitting method for pricing American options under the Bates model

    Texts in Applied Mathematics. New York: Springer.isbn: 978-0-387-97999-1. Toivanen, J. (2010). “A componentwise splitting method for pricing American options under the Bates model”. In:Computational Methods in Applied Sciences15, pp. 213–227. Varah, J. M. (1975). “A lower bound for the smallest singular value of a matrix”. In:Linear Algebra and its Applic...

    2010

  12. [12]

    Appendices A M-matrices, EM-matrices and positivity preserving solutions To define non-negativity of the solution of some problem (which is equivalent oto the solution of the corresponding PDE or PIDE), we introduce the following definitions [Kuzmin, Löhner, and Turek, 2012]. Let( M,µ)be a σ-finite measure space withµbeing a positive measure defined on aσ...

    2012

  13. [13]

    Using the above definitions and properties of M-matrices, we can now address our main goal

    All diagonal elements ofA are positive, and there exists a positive diagonal matrixD such thatAD is strictly diagonally dominant, or equivalently D−1ADis strictly diagonally dominant; in particular,ADhas all positive row sums. Using the above definitions and properties of M-matrices, we can now address our main goal. As noted in the Introduction, we requi...

    2016

  14. [14]

    That is, a Metzler matrix is any matrixA which satisfiesA= (a ij);a ij≥0, i̸=j

    To avoid working directly with matrices of the form−B, we introduce the following definition Definition A.8.AMetzleror essentially nonnegative is a matrix if all of its elements are non-negative except for those on the main diagonal, which are unconstrained. That is, a Metzler matrix is any matrixA which satisfiesA= (a ij);a ij≥0, i̸=j. It can also be see...

    1994

  15. [15]

    We also denote a unit matrix asI

    Definition A.10.The second order approximations of∇≡∂x and∇2≡∂x,x are defined, respectively, by forward(F),backward(B) and central (C) discretizations as FB 2 C(x) = 3C(x)−4C(x−h) +C(x−2h) 2h ,F C 2 =FF 1 +F B 1 2 , FF 2 C(x) =−3C(x) + 4C(x+h)−C(x+ 2h) 2h ,S B 2 C(x) = 2C(x)−5C(x−h) + 4C(x−2h)−C(x−3h) h2 , SC 2 C(x) = C(x+h)−2C(x) +C(x−h) h2 ,S F 2 C(x) =...

    2011

  16. [16]

    Lemma 17.LetA∈RN×N

    We also require the following lemma from [Noutsos and Tsatsomeros, 2008]. Lemma 17.LetA∈RN×N. The following conditions are equivalent: 1.Ais eventually exponentially nonnegative. 2.A+bIis eventually nonnegative for someb≥0. 3.A T +bIis eventually nonnegative for someb≥0. For a proof, see [Itkin, 2017b]. Definition A.12.An N×Nmatrix A = [aij]is called anEM...

    2008

  17. [17]

    B Conservativeness of Strang splitting for the FPE We start with some definitions and assumptions. We equip Rnxny with the scaled Euclidean norm ∥v∥2 2 = hxhy ∑ i,jv2 ij, the discrete analogue of theL2(Ω)norm, and denote byp( t)the solution of the semi-discrete (method-of-lines) system ˙p(t) =A(t)p(t), A(t) =A x(t) +Ay(t) +Axy(t),p(0) =p 0,(B.1) with Aα= ...

    2007

  18. [18]

    , withC(T)independent ofh,∆tandn. Proof. The argument is the standard Lady Windermere fan, carried out inℓ1 on the positivity regime; the only changes from theℓ2 version are that the uniform stability bound is theℓ1 Markov bound of Assumption 1 and that the local defects, bounded inℓ2 by Assumption 2, are converted toℓ1 through ∥v∥1≤|Ω|1/2∥v∥2 (a conseque...

    1954

  19. [19]

    are unrestricted. Our numerical experiments in Section 5 useρ= 0.8, comfortably away from the degeneracy; the degenerate limit¯ρ= 1, which would require a hypoelliptic regularity argument, is left for future work. △ C Proof of Proposition 9 Let us assumeρ>0since the proof forρ≤0can be done in the same way. (i)Consider Ty :=QI+ρ √ ∆tA B 2,y. Withineach1Dbl...

    1975

  20. [20]

    Positivity of TR-BDF2

    For the classicalγ= 2− √ 2 one hasω= 1/ √ 2, a = 1+ √ 2 2 , b = √ 2−1 2 , and the special property α=γ/2: both stages involve the same matrixI−γk 2 A [Bank et al., 1985; Hosea and Shampine, 1996]. Positivity of TR-BDF2. Although the right-hand side of Eq. (E.10) carries a negative weight−b<0onp n, so that no stage-wise positivity argument applies, the neg...

    1985

  21. [21]

    The window Eq.(E.12)exceeds the CN window Eq.(E.5)by the factor(1 + √ 2)/2≈1.21. Proof. In Eq. (E.11) the factor(I−γk 2 A)−2≥0for all k >0by Theorem 7; the factorI + ( √ 2−1)kA is Metzler off the diagonal and has diagonal1 + ( √ 2−1)kbi≥0iff Eq. (E.12) holds. Page 63 of 64 Positivity-preserving FD schemes No second-order one-step rational scheme can do be...

    1978