PaperlayerView exact revision on arXiv

Source-linked AI summary

A discontinuous Galerkin approximation of the Dean--Kawasaki equation

Kamran Arora, Federico Cornalba, Tony Shardlow

arXiv:2609.06681v1math.NAcs.NAmath.APmath.PR

TL;DR

The Dean–Kawasaki SPDE models density fluctuations but is highly singular and formally derived, motivating spatial discretizations in the high-density regime. The paper constructs a cross-variation-preserving finite-element noise discretization and an arbitrary-order dG method; its weak error is O(h^p) with an exponentially small correction when Nh^d >>1.

  • Problem

    The Dean–Kawasaki SPDE models particle-density fluctuations but is highly singular and formally derived, motivating spatial discretizations for high-density regimes.

    Source: §1 · pp. 1–2

  • Method

    The paper constructs an explicit arbitrary-order finite-element discretization of multiplicative divergence-form noise that preserves its cross-variation structure.

    Source: §1 · p. 2

  • Results

    Weak error is O(h^p), with an exponentially small correction for the scheme’s negative part when Nh^d >>1.

    Source: §5.3 · p. 17, §5.6 · p. 20

  • Takeaways & Limitations

    The noise construction applies to continuous and discontinuous finite-element spaces and extends to general mobilities.

    Source: §1 · p. 2

    The analysis covers independent Brownian motions; potentials and zero-flux Neumann boundaries are shown numerically but lie beyond its theory.

    Source: §7 · p. 32

Abstract

from arXiv · show

We introduce and analyse an arbitrary order spatial discontinuous Galerkin (dG) method for the Dean--Kawasaki equation, a highly singular SPDE modelling density fluctuations of $N$ diffusing particles in the regime of large particle number $N \gg 1$. Our starting point is a general procedure for discretising multiplicative, divergence-form noise on finite element spaces whilst preserving its cross-variation structure at the discrete level; the construction is explicit, elementwise, applies to continuous and discontinuous spaces alike, and extends to general mobilities. Using it, we prove weak error estimates of order $O(h^p)$ between fluctuations of the semi-discrete scheme and those of the underlying particle system, together with a correction that is exponentially small in the scaling regime $Nh^d \gg 1$ and arises because the scheme does not preserve positivity. The resulting method is locally and globally conservative and applies on unstructured simplicial meshes. Quantitative numerical experiments support the analysis, and further experiments illustrate the method beyond the scope of the theory: indicator function observables, external and interaction potentials, convection-dominated regimes and reflecting boundary conditions.

1. Introduction

The Dean–Kawasaki equation models density fluctuations in diffusing particle systems, but its singularity and formal derivation make its analysis challenging. This work develops an arbitrary-order dG method and analyzes its weak error against particle-system fluctuations, while situating the method among existing regularizations and discretizations.

1. Introduction

  • The Dean–Kawasaki equation models density fluctuations in diffusing particle systems, yet is highly singular and has only trivial martingale solutions.In the high-density regime Nh^d >>1, discretized simulation costs O(h^-d), compared with at least O(N) for direct particle simulation.Source: §1 · pp. 1–2

1.1. Our contributions.

  • The paper introduces an explicit arbitrary-order procedure that discretizes multiplicative divergence-form noise while preserving its discrete cross-variation structure.The construction is elementwise, works for continuous and discontinuous finite element spaces, and extends to general mobilities.Source: §1 · p. 2
  • Weak error estimates for dG fluctuations have rate O(h^p), with an exponentially small correction associated with negative solution values.The analysis is semi-discrete and assumes a mesh-size restriction because the interior penalty method lacks a discrete maximum principle.Source: §1.1 · pp. 2–3
  • Numerical experiments support the weak error analysis and explore indicator-function observables, potentials, and reflecting boundaries beyond its theoretical scope.The analysis treats independent Brownian motions; external and interaction potentials and zero-flux Neumann conditions are demonstrated numerically.Source: §1.1 · p. 3

1.2. Related literature.

  • Earlier work established weak-error results for finite-difference and piecewise-linear finite-element approximations, while related studies consider regularized models and other discretizations.The cited literature includes analytical regularization, inertial-particle models, spatially correlated noise, finite elements, finite volumes, and hybrid particle-SPDE approaches.Source: §1.2 · pp. 3–4

1.3. Structure of the paper.

  • The paper proceeds from notation and scheme construction through noise analysis, conservation properties, weak-error estimates, and numerical experiments, with extensions discussed last.Appendices provide proofs and auxiliary estimates, including adjoint heat-equation estimates and a method for constructing initial data.Source: §1.3 · p. 4

2. Notation

The notation describes shape-regular, quasi-uniform simplicial meshes and piecewise-polynomial finite element spaces. It also defines the mesh scale, interpolation, and operators used to handle discontinuities.

2. Notation

  • Meshes are shape-regular, quasi-uniform, affine simplicial partitions; h is the maximum element diameter and scales comparably to each element diameter.The facet sizes are also assumed not to degenerate relative to element sizes.Source: Mesh · pp. 4–5
  • The discontinuous Galerkin space consists of functions that are polynomials of degree p on each mesh element.The associated vector-valued polynomial spaces provide the discrete-gradient range.Source: Mesh · p. 5
  • The nodal interpolation operator matches a function at each element’s nodal points, and the paper defines scalar and vector jump conventions on facets.Scalar inputs have vector-valued jumps; vector inputs have scalar-valued jumps, independent of element labeling.Source: Mesh · p. 5, Average and jump operators · p. 5
  • The broken gradient differentiates element interiors, while the discrete gradient additionally lifts facet-jump contributions into an L2 vector field.The paper uses the discrete gradient with index p − 1 unless stated otherwise.Source: Broken and discrete gradients · p. 6

3. The numerical scheme

The scheme discretises the Dean–Kawasaki noise through a variational formulation that preserves its cross-variation structure on the discontinuous space. It combines this noise with a symmetric interior penalty discretisation of the Laplacian.

3.1. Variational formulation of the noise.

  • The noise discretisation uses the discrete gradient to account for facet jumps and seeks a finite-dimensional martingale with the prescribed cross-variation structure.The central numerical flux is symmetric and conservative, and its conservation property yields global mass conservation.Source: §3.1 · p. 7, cross-variation structure · pp. 7–8

3.2. Discretisation of the Laplacian.

  • The Laplacian is discretised with the symmetric interior penalty method, which is symmetric and coercive when its penalty parameter is sufficiently large.The required penalty size depends on dimension, polynomial degree, and mesh shape regularity.Source: §3.2 · p. 8

3.3. Definition of the numerical scheme.

  • The dG Dean–Kawasaki scheme combines the Laplacian discretisation with an explicit Brownian representation of the noise, whose coefficients mimic √ρ_h while recovering the prescribed cross-variation.The coefficients are constructed using the positive part of the numerical solution.Source: cross-variation structure · p. 7, §3.3 · p. 8

4. Construction and properties of the noise and weak existence of a solution to the numerical scheme

The paper constructs an explicit discrete noise that preserves the desired cross-variation and establishes weak existence and moment bounds for the semi-discrete scheme. The construction also extends to general mobilities, and the scheme has global and local conservation properties.

4. Construction and properties of the noise and weak existence of a solution to the numerical scheme

  • The explicit noise construction yields a martingale with the desired cross-variation structure and supports probabilistically weak solutions with moment bounds.The construction applies to arbitrary-order discontinuous spaces and can also be adapted to continuous Galerkin spaces beyond degree one.Source: §4 · pp. 9–10
  • The noise construction is a general procedure for discretising multiplicative, divergence-form noise: a suitable non-negative numerical approximation of prescribed mobility leaves the remainder unchanged and preserves cross-variation structure at the discrete level.The Dean–Kawasaki setting uses m(ρ) = ρ; for the Simple Symmetric Exclusion Process, m(ρ) = ρ(1 −ρ).Source: §4 · pp. 11–12
  • For degree-one elements, the noise has a simpler representation that avoids matrix factorisation by using diagonal mass matrices and projection onto the piecewise-constant space.This reduces the computation to inverting a diagonal matrix.Source: §4 · p. 12
  • The scheme conserves global mass and is locally conservative: an element’s mass changes only through deterministic and stochastic fluxes across its boundary.The local conservation statement is established for every element and time.Source: §4 · pp. 12–13

5. Weak error analysis

This section develops positivity bounds, adjoint-heat-flow tools, and a recursive moment argument to establish weak moment-error estimates comparing particle-system fluctuations with the semi-discrete dG scheme. The analysis separates numerical residuals from effects of negative solutions and notes a mass-preserving post-processing option when element averages are non-negative.

5.1. Mean-field limit.

  • The continuous mean-field PDE and its SIP semi-discretization provide the centering profiles for comparing fluctuations, and the scheme has expectation equal to its mean-field solution.The equality holds for positive times because the drift is linear and the noise is a true martingale, and at time zero because both start from the same initial condition.Source: §5.1 · pp. 13–14

5.2. Assumptions.

  • Under dG1 and dG2, sufficiently small h keeps the semi-discrete mean-field profile between ρmin − δ and ρmax + δ, ensuring the strict positivity needed by the analysis.The SIP-dG heat discretization need not preserve nonnegativity, so the proof obtains positivity through a mesh-size restriction. If the conditional SIP-dG maximum-norm estimate holds, the scaling improves from p + 1 − d/2 to p + 1 up to a logarithmic factor.Source: §5.2 · pp. 15–16

5.3. Statement of the weak error estimate.

  • The weak error compares polynomial moments of tested particle-system and dG observables, splitting the bound into an O(h^p) discretization error and an exponentially small negative-part error when Nh^d >> 1.The comparison is weak because the particle system is a measure while the scheme is a density; both are tested against sufficiently regular functions. The theorem assumes dG1 and dG2 and uses a uniformly positive semi-discrete mean-field profile.Source: §5.3 · pp. 16–17

5.4. Outline of the proof.

  • The proof preserves the noise cross-variation structure in the semi-discrete scheme and compares the resulting noise terms on suitable test functions.Adjoint heat-flow test functions cancel the drift; because dG L∞ estimates are unavailable, the argument uses an L∞-L1 Hölder inequality and L∞(L2) gradient estimates.Source: §5.4 · p. 18, Choice of test functions · p. 18

5.7. Setting up the recursive structure.

  • The proof sets up a recursion for moment errors whose order is linked to test-function regularity; closing the recursion accumulates the error terms and yields the weak estimate.For notational convenience, the proof begins with equal final times.Source: Recursion · p. 18, §5.7 · p. 20

5.5. Backwards (adjoint) heat flows.

  • The analysis introduces continuous and SIP semi-discrete backward heat flows, which evolve final data and provide the adjoint test functions used in the convergence argument.The flows are defined from a final datum at a final time, with the discrete initial datum obtained by interpolation; Appendix C collects their error estimates.Source: §5.5 · p. 19

5.6. L∞(L∞) estimates.

  • The scheme’s expected negative part is exponentially small when Nh^d >> 1; bounding it in L∞ avoids an h^-d loss from an inverse inequality.The preceding L∞(L∞) estimates yield the negative-part bound, while the scaling permits uniform-in-h estimates.Source: §5.6 · p. 20

5.7. Setting up the recursive structure.

  • The dG-DK and Brownian particle system have agreeing first moments, while the second-moment estimate exposes negative-solution and numerical errors, including a numerical error of order O(h^p).The higher-moment estimates are organized through a recursive formula under assumptions dG1 and dG2.Source: §5.7 · pp. 21, 23

5.8. Proof of Theorem 5.7.

  • The proof exhausts a recursion of depth j − 2, controlling negative-solution and numerical residuals combinatorially and requiring test-function regularity W^p+j,∞∩H^p+j+1.Each recursion step introduces a test function requiring one derivative; the base regularity is W^p+2,∞∩H^p+3.Source: §5.8 · p. 24
  • When the solution has non-negative element averages, a Zhang–Shu type limiter can restore non-negativity while preserving elementwise mass, but applicability is only bounded at scaling 1 − E(N,h).The limiter is post-processing and is not fed back into the semi-discrete evolution; the scheme’s negative solutions occur with exponentially small probability.Source: §5.8 · pp. 24–25

6. Numerical experiments

The experiments numerically check weak-error behavior for smooth and indicator observables, then illustrate dG schemes in convection-dominated, interacting-particle, and reflecting-boundary settings.

6.1. Matrix formulation of (dG-DK).

  • The matrix formulation exploits block-diagonal mass matrices and factorization, while stiffness and discrete-gradient matrices retain local sparsity from inter-element flux coupling.Source: §6.1 · p. 25

6.2. Implementation details.

  • The numerical method is implemented in Firedrake, with figure-reproduction code available in an open-source repository.Source: §6.2 · p. 25

6.3. Error with smooth observables.

  • p=1 shows O(h^2), while p=2 initially shows O(h^2) before Monte Carlo noise prevents fine-mesh rate conclusions.For the smooth observable, p=2 appears to have lower overall error than p=1.Source: §6.3 · p. 26, Figure 1 · p. 27

6.4. Error and cost comparison with discontinuous observables.

  • For indicator observables, p=1 converges at O(h^3/4) and p=2 at O(h); at comparable degrees of freedom, p=2 has lower overall error after cost adjustment.Source: §6.4 · pp. 26–27, Figure 2 · p. 28

6.5. Convection-dominated settings.

  • In convection-dominated tests, cG has large spurious oscillations and negative undershoots, while upwind dG remains stable and appears to suppress transport oscillations.The analysis of exponentially rare negativity in the N h^d ≫ 1 regime does not cover this setting.Source: Figure 3 · p. 29, §6.5 · p. 29

6.6. Weakly interacting particle systems.

6.7. Reflecting Brownian motions.

  • On a disk with reflecting boundaries, the dG solution reproduces sampled particle behavior, including an annular boundary layer with similar density values.The authors note that zero-flux conditions extend the method to irregular domains.Source: Figure 5 · p. 31, §6.7 · p. 32

7. Extensions and future work

The analysis is limited to independent Brownian motions, while experiments also consider interaction and external potentials and zero-flux boundaries. Future directions include positivity preservation, fluctuation-dissipation structure, and fully discrete analysis.

7. Extensions and future work

  • The analysis covers independent Brownian motions; extending it to interaction potentials, external potentials, and zero-flux Neumann conditions remains future work.Those settings are explored numerically beyond the theory.Source: §7 · p. 32
  • Positivity preservation remains an open question; proposed work includes improving elementwise post-processing and analyzing its use at every timestep.The latter requires fully discrete analysis and could also address mass loss from the positive part in the noise.Source: §7 · p. 33
  • For SIP dG, the penalty parameter adds dissipation without matching stochasticity, breaking the linearized fluctuation-dissipation structure.A penalty-free form using a higher-order discrete gradient could make the balance structural, but increases computational cost and loses exact consistency.Source: §7 · p. 33
  • Multilevel Monte Carlo could reduce observable-estimation cost, but requires designing coarse–fine noise couplings compatible with the high-order noise discretisation.Such couplings remain a direction for investigation.Source: §7 · p. 33

Appendix A. Proof of weak existence and moment bounds

The appendix recasts the variational scheme as a finite-dimensional SDE and establishes weak existence and moment bounds using growth estimates and stopping-time arguments.

Appendix A. Proof of weak existence and moment bounds

  • Expanding the variational solution in a fixed basis yields a finite-dimensional SDE whose linear growth, with continuous drift and diffusion, gives a non-explosive probabilistically weak solution.The cited finite-dimensional SDE results establish existence on [0,T].Source: Appendix A · p. 36
  • The coefficient representation uses the mass matrix, SIP stiffness matrix, and a diffusion matrix formed from noise functions paired with discrete gradients.These matrices express the variational scheme in coordinates.Source: Appendix A · p. 36
  • Moment bounds follow by applying growth estimates, Hölder and BDG inequalities to a stopped SDE, then removing the stopping via non-explosion and Fatou's lemma.The estimate is independent of the stopping threshold.Source: Appendix A · pp. 36–37

Appendix B. Recursive relations for moments

The appendices derive recursive moment identities for the scheme and gather heat-equation estimates used in the convergence analysis.

Appendix B. Recursive relations for moments

  • Itô differentials for products of test-function pairings yield recursive moment relations, with induction and Hölder estimates extending the bounds to higher moments.The proof uses the scheme itself for the first-order case and adapts earlier arguments to the lack of L∞(L∞) solution estimates.Source: Appendix B · pp. 37–38

Appendix C. Heat estimates

  • Time reversal reduces estimates for the adjoint heat flow to the regular heat equation, whose maximum-principle and energy estimates are recalled for later use.The reversed problem uses the observable as initial data.Source: Appendix C · p. 39, §C.1 · p. 39

C.1. Standard estimates.

  • For initial data in H^(p+2), the standard heat estimates provide the preliminary bounds used in the appendix's later elliptic and parabolic estimates.The cited lemma states these bounds with a constant independent of h.Source: §C.1 · p. 39

C.2. Elliptic estimates.

The appendix develops elliptic estimates for the SIP discretisation and then uses them to establish parabolic heat-equation error estimates under regularity and coercivity assumptions.

C.2. Elliptic estimates.

  • The Ritz projection is defined through the SIP bilinear form; a constraint ensures uniqueness because that form does not distinguish constants.The SIP form is bounded, and is coercive when the penalty parameter is sufficiently large.Source: §C.2 · pp. 39–40

C.3. Parabolic estimates.

  • For φ in H^(p+2) and sufficiently large penalty, the SIP semi-discretisation admits heat-equation estimates with constants independent of h.The parabolic estimates build on the elliptic bounds and assume coercivity of the SIP form.Source: §C.3 · pp. 40, 42

C.4. Gradient estimates.

This section gathers gradient-related estimates for the SIP semi-discretisation of the heat equation. The bounds rely on coercivity and include constants independent of mesh size.

C.4. Gradient estimates.

  • The section recalls a lemma estimating the lifting operators defined in (2.3) and (2.4).The supplied passage does not state the estimates themselves.Source: §C.4 · p. 43
  • The corollary's estimate requires η large enough to ensure coercivity of a_h(·,·).Source: §C.4 · p. 43
  • The SIP semi-discretisation is considered with initial datum I_hφ_i, alongside the continuous heat equation with initial datum φ_i ∈ H^{p+2}.Source: §C.4 · p. 43
  • The stated estimates use constants independent of h.Source: §C.4 · p. 43

Appendix D. Proofs for the convergence analysis

Appendix D provides proofs for the results stated in Section 5.

Appendix D. Proofs for the convergence analysis

  • Appendix D supplies proofs for the various results stated in Section 5.Source: Appendix D · p. 44

D.1. L∞(L∞) estimates.

The appendix develops estimates used in the convergence analysis and proves the recursive formula for higher moments. Its argument relies on mesh-dependent bounds and controls numerical-error terms at order O(h^p).

D.1. L∞(L∞) estimates.

  • Lemma D.1 relates polynomial norms on an element to values at a finite set of points, with a constant depending on p and d.Source: §D.1 · p. 44
  • Lemma D.2 gives estimates for nodal basis functions and the element mass matrix on affine, shape-regular, quasi-uniform simplicial meshes.Source: §D.1 · p. 44
  • Lemma D.3 provides an estimate for the SIP bilinear form using element nodal points and basis functions.Source: §D.1 · p. 45
  • In Step 4, the estimate is extended to all times in [0, T]; for β ≥ 6d + 8 and h sufficiently small, this requires an h^-2 estimate on the discrete Laplacian, supplied by Lemma D.3.Steps 1 and 3 require significant changes from [8, Proposition 13], while the remaining steps follow the same argument with only minor changes to the constants.Source: §D.1 · pp. 46, 48

D.2. Recursive formula for higher moments.

  • The higher-moment proof uses L∞–L1 Hölder estimates and combines numerical-error terms of order O(h^p).The exposition identifies these estimates as part of the proof of Proposition 5.13.Source: §D.2 · pp. 48, 52

Appendix E. Construction of initial data

Appendix E constructs initial densities and deterministic particle placements for the convergence assumptions. The explicit barycentre construction is restricted to polynomial degrees one and two.

Appendix E. Construction of initial data

  • The initial density ρ0,h is bounded strictly above zero and below a finite upper bound for sufficiently small h.Source: Appendix E · p. 52
  • A partition into N axis-parallel rectangles is constructed with a bound whose constant depends only on d and ρmax/ρmin, not N.Source: Appendix E · p. 52
  • The explicit initial particle placement uses each cell's ρ0-weighted barycentre.Source: Appendix E · p. 53
  • For p ∈ {1, 2}, the barycentre construction gives an initial-data error bounded by O(h^p) under the stated scaling N h^d ≳ 1.The proof bounds the interpolation contribution by O(h^{p+1}) and uses the N^{-2/d} term.Source: Appendix E · p. 54
  • The barycentre construction is structurally limited to p ≤ 2 because it reproduces mass and first moments, making its rule exact on P1 but no higher.For general p, the text describes placing q particles per cell for a rule exact on Pp, but says that construction is non-constructive and is not pursued.Source: Appendix E · p. 54
Loading 2609.06681v1…