PaperlayerView exact revision on arXiv

Source-linked AI summary

A weighted polygonal discontinuous Galerkin method for hierarchically coupled reaction-diffusion systems with cubic interactions

Mattia Corti

arXiv:2609.34625v1math.NAcs.NA

TL;DR

Pairwise models may not capture higher-order population effects, while cubic reactions make analysis delicate. The paper analyzes a hierarchical SWIP-PolyDG method with dissipative diagonal cubic self-interactions. It proves local-in-time stability and an a priori error estimate, and reports numerical convergence and robustness results.

  • Problem

    Higher-order interactions modify pairwise effects, while cubic reaction terms complicate analysis and require control of higher Lebesgue norms.

    Source: §1 · pp. 1–2

  • Method

    The paper analyzes a SWIP-PolyDG method for hierarchically coupled systems with heterogeneous diffusion, using dissipative diagonal cubic self-interactions to control reactions recursively.

    Source: §1 · p. 2

  • Results

    The analysis establishes recursive local-in-time stability estimates and an a priori combined L2(Ω)–dG error estimate; numerical experiments support expected convergence and diffusion robustness.

    Source: §1 · p. 2, §8 · p. 18

  • Takeaways & Limitations

    The framework accommodates up to cubic reactions and heterogeneous, possibly discontinuous or anisotropic diffusion tensors on general polygonal meshes.

    Source: §8 · p. 18

    The stability analysis is developed in two dimensions; extending it to three dimensions remains open because the available embedding may not provide the higher integrability required.

    Source: §1 · p. 2, §4 · p. 9

Abstract

from arXiv · show

We develop and analyse a symmetric weighted interior penalty polytopal discontinuous Galerkin (SWIP-PolyDG) method for multi-species reaction--diffusion systems with nonlinear reaction terms of up to cubic order. Such terms arise naturally when higher-order interactions are incorporated into population models, allowing non-additive effects among multiple species to influence local growth and conversion mechanisms. The proposed framework accommodates heterogeneous and possibly anisotropic diffusion tensors on general polygonal meshes. To control the nonlinear coupling, we consider a hierarchical block structure in the reaction operator, whereby each population block depends only on its own variables and on those associated with preceding blocks. In addition, the cubic self-interactions within each block are assumed to have a dissipative diagonal structure. Under these hypotheses and proceeding recursively over the hierarchy of population blocks, we derive local-in-time stability estimates in two spatial dimensions for the semi-discrete formulation in both the $L^2(\Omega)$-and dG-norms. We further establish an a priori error estimate in a combined $L^2(\Omega)$-dG energy norm for sufficiently regular solutions for the semi-discrete formulation. Finally, numerical experiments confirm the predicted convergence behaviour and illustrate the robustness of the method under heterogeneous diffusion.

1 Introduction

Higher-order interactions can make pairwise population models inadequate, but cubic reactions complicate analysis and may permit finite-time blow-up. The paper proposes a hierarchical SWIP-PolyDG framework and establishes stability and error results under structural assumptions, with numerical tests of convergence and diffusion robustness.

1 Introduction

  • Higher-order interactions modify pairwise effects, and cubic reactions can affect ecological coexistence, stability, and persistence while also arising in epidemiological and biochemical models.The introduction identifies finite-time blow-up as a possibility for superlinear reaction–diffusion systems in general.Source: §1 · p. 1
  • In two dimensions, discrete Sobolev embeddings provide the higher-Lebesgue-norm bounds needed for cubic terms; extension of the analysis to three dimensions remains open.The cited obstacle is that the three-dimensional embedding only reaches p ≤ 6 and may not provide the required higher integrability.Source: §1 · p. 2
  • The method assumes each population block depends only on itself and preceding blocks, with dissipative diagonal cubic self-interactions to control reactions recursively.The analysis addresses the resulting semi-discrete system under these structural hypotheses.Source: §1 · p. 2
  • The SWIP-PolyDG discretization handles reactions up to cubic order on general polygonal meshes with heterogeneous diffusion tensors, and its diffusion form is coercive and continuous for sufficiently large harmonic-average-based penalties.Source: §1 · p. 2
  • The authors derive local-in-time L2- and dG-norm stability estimates recursively across blocks and an a priori combined L2(Ω)–dG error estimate for sufficiently regular solutions.For purely quadratic nonlinearities, the error result extends the Fisher–Kolmogorov analysis in [15, Theorem 2] to vector-valued systems while removing a previously required relation between diffusion coercivity and reaction terms.Source: §1 · p. 2
  • Numerical tests cover a three-species Lotka–Volterra system, a manufactured two-species cubic-HOI problem, and a coupled Fisher–Kolmogorov system with discontinuous diffusion.The last test illustrates robustness under heterogeneous diffusion.Source: §1 · p. 2

2 Mathematical model

The model is a multi-species reaction–diffusion system with linear, quadratic, and cubic reaction components and population-specific diffusion tensors. Its block hierarchy restricts reactions to current and preceding blocks, and the weak formulation expresses the system against test functions.

2 Mathematical model

  • The model describes n populations with diffusion tensor D, external forcing γ, and reaction components collected as G(w)w = (L + Qw + (Cw)w)w.L, Q, and C encode linear, quadratic, and cubic reaction coefficients, respectively.Source: §2 · p. 3

2.1 The block structure

  • The reaction structure is divided into population blocks, with each block depending only on its own and preceding blocks; each diagonal cubic block is diagonal with nonnegative coefficients.The block variables and forcing are partitioned consistently across the population system.Source: §2.1 · p. 3

2.2 Weak formulation

  • The weak formulation seeks w in the stated time-dependent function spaces so that the time derivative and diffusion terms balance reaction and forcing against every test function v.Source: §2.2 · p. 4

3 Polytopal discontinuous Galerkin semi-discrete formulation

The semi-discrete formulation uses discontinuous polynomial spaces on polytopal meshes, with diffusion-weighted face operators and a penalized diffusion bilinear form. Its analysis establishes coercivity and continuity and introduces inequalities used in the stability argument.

3.1 Discrete setting and preliminary estimates

  • The discrete setting uses a polytopal mesh with element diameters and face interfaces, alongside mesh regularity and submesh conditions.The paper distinguishes two-dimensional line-segment interfaces from three-dimensional interfaces decomposed into planar triangles.Source: §3.1 · p. 4
  • The discontinuous finite element space is built from elementwise polynomials of total degree ℓ ≥ 1, with trace, weighted-average, and jump operators defined on faces.The weighted averages are chosen according to the physical diffusion parameter.Source: §3.1 · pp. 4–5
  • Elementwise constant symmetric diffusion tensors may jump across interfaces, and the formulation uses diffusion-dependent weights and harmonic averages in its face terms.Source: §3.1 · p. 5, §3.2 · p. 5

3.2 Semi-discrete formulation

  • The semi-discrete PolyDG equation combines the time derivative, diffusion bilinear form, nonlinear reaction, and forcing for every discrete test function.Source: §3.2 · p. 5
  • The diffusion bilinear form is coercive when the penalty parameter is sufficiently large, and its continuity is also established.The stated coercivity condition is η̃ ≥ 4Cinv, where Cinv is the inverse-trace-inequality constant.Source: §3.2 · pp. 5–7
  • A two-dimensional Gagliardo–Nirenberg-type interpolation inequality and an extended Perov-type inequality are provided for the subsequent stability analysis.The Perov-type result accommodates a non-constant forcing term a.Source: §3.2 · p. 7

4 Stability analysis

The stability analysis proves local-in-time bounds recursively across population blocks in two dimensions. Later blocks inherit cumulative dependence on earlier blocks, and the estimates support existence of the discrete solution over the stated validity interval.

Stability analysis

  • In two dimensions, the first-block theorem establishes stability under Assumptions 1 and 2 when the penalty parameters are sufficiently large.Source: §4 · p. 8
  • The first-block result also gives L2- and DG-norm stability for semi-discrete systems with purely quadratic interactions, including several classical biological and chemical models.Examples include reaction–diffusion Lotka–Volterra systems, mass-action chemical kinetics, and epidemic models with spatial diffusion.Source: §4 · p. 9
  • The proof does not directly extend to three dimensions: H1 embeds only into Lq for q ≤ 6, not L8, which is needed for the higher-norm control.A three-dimensional extension may require adapting interpolation exponents and imposing corresponding admissibility conditions; the paper does not pursue it here.Source: §4 · pp. 9, 11
  • Each later block’s stability estimates depend on all preceding blocks, and its admissible time interval is generally no longer than those of earlier blocks.Instability or blow-up in an early block restricts the validity horizon of higher blocks that interact with it.Source: §4 · p. 12
  • The unique discrete solution is guaranteed to exist up to t∗_pM, as the stability estimates rule out blow-up on their validity intervals.Source: §4 · p. 12

5 Error analysis of the semi-discrete formulation

The section establishes an a priori error estimate for the semi-discrete problem under regularity and penalty-parameter conditions. Its proof uses an interpolant-based error decomposition and bounds the nonlinear terms using stability estimates and interpolation inequalities.

5 Error analysis of the semi-discrete formulation

  • The analysis extends the continuous bilinear form by assuming element-wise H2 regularity and continuity of the flow across interfaces.Continuity of the bilinear form is stated separately.Source: §5 · p. 12
  • For sufficiently regular solutions and a sufficiently large penalty parameter, Theorem 3 establishes an a priori error estimate for the semi-discrete problem.The hidden constant and α(τ) may depend on the solution, initial data, domain, and polynomial degree, but not on h.Source: §5 · p. 12
  • The proof splits the error into interpolation and discrete components, then controls cubic reaction differences using solution bounds and standard inequalities.The estimates use stability results, Hölder and Young inequalities, and a Sobolev–Poincaré inequality.Source: §5 · p. 13
  • For purely quadratic reactions, the estimate extends the two-dimensional Fisher–Kolmogorov result to vector-valued systems and removes its reaction–diffusion coercivity relation [15, Theorem 2].The stated improvement concerns the structural relation previously required between reaction terms and the diffusion coercivity constant.Source: §5 · p. 13

6 Time discretization: fully-discrete formulation

The semi-discrete problem is expressed as a matrix evolution system in basis coefficients, then advanced with an IMEX–RK scheme. The time discretization is used for numerical experiments, without a fully discrete stability or error analysis.

6 Time discretization: fully-discrete formulation

  • The semi-discrete formulation is written as MẆ(t) + KW(t) = G(W(t)) + Γ for the vector W(t) of basis-expansion coefficients.The basis is arranged in population-sized blocks, with each block corresponding to one population.Source: §6 · p. 14
  • The IMEX–RK scheme treats linear diffusion implicitly and the nonlinear reaction term G(W) + Γ explicitly.Its internal stages are defined by Butcher tableaux for a high-order optimized method.Source: §6 · p. 14
  • Fully discrete stability and error analysis is outside the paper’s scope, and convergence-study time steps are chosen so temporal error is negligible relative to spatial error.The IMEX–RK discretization is introduced for the numerical experiments.Source: §6 · p. 14

7 Numerical results

The numerical studies test spatial convergence on a three-species competition system and a manufactured cubic-interaction system, then examine qualitative robustness for heterogeneous diffusion. The reported convergence rates match predictions, while the heterogeneous case resolves asymmetric propagation and sharp layers without visible spurious oscillations.

7 Numerical results

  • The experiments cover predicted spatial convergence, p-refinement, and qualitative robustness under discontinuous diffusion using three representative reaction–diffusion cases.The cases include a Lotka–Volterra system, a manufactured two-species cubic-interaction problem, and a coupled Fisher–Kolmogorov model.Source: §7 · p. 14

7.1 Test case 1: Three-species competition-diffusion system

  • The three-species test matches the predicted O(h^ℓ) energy-norm rate and observes O(h^(ℓ+1)) L2 convergence, which is not covered by the analysis.The ℓ = 4 approximation is visually indistinguishable from the exact travelling wave, while ℓ = 1 loses accuracy near the transition region.Source: §7.1 · p. 15, Figure 1 · p. 15, Figure 2 · p. 16

7.2 Test case 2: Two-species system with higher-order interaction terms

  • The cubic-interaction test shows O(h^ℓ) convergence in dG and energy norms and O(h^(ℓ+1)) in L2; errors also decay rapidly with polynomial degree.The h-rates agree with predictions, while p-convergence is not covered by the theory.Source: §7.2 · p. 16, Figure 3 · p. 17

7.3 Test case 3: A heterogeneous double Fisher–Kolmogorov system

  • In the heterogeneous Fisher–Kolmogorov test, fronts advance faster where diffusion is 10 times larger, and sharp interface layers show no visible spurious oscillations.The snapshots depict two populations at t = 0, 5, 10, 15, and 20.Source: §7.3 · p. 18, Figure 4 · p. 18

8 Conclusions

The paper analyzes a weighted discontinuous Galerkin method for hierarchically coupled reaction–diffusion systems with up to cubic reactions, using dissipative diagonal structure to control intra-block cubic terms. It establishes semi-discrete stability and error results, while identifying three-dimensional analysis and fully discrete theory as open directions.

8 Conclusions

  • The method accommodates multiple populations and heterogeneous, possibly discontinuous or anisotropic diffusion tensors on general polygonal meshes.Its intra-block cubic interactions are assumed to have a dissipative diagonal structure.Source: §8 · p. 18
  • Recursive blockwise analysis yields local-in-time semi-discrete stability estimates and an a priori error estimate in a combined L2(Ω)–dG energy norm for sufficiently regular solutions.Numerical experiments support expected convergence for classical and higher-order systems and illustrate robustness with discontinuous diffusion in a coupled Fisher–Kolmogorov test.Source: §8 · p. 18
  • Extending the analysis to three dimensions remains open because the proof's Sobolev embeddings do not control the higher Lebesgue norms required by cubic reactions.The present theory is semi-discrete; fully discrete stability and error analysis for IMEX–Runge–Kutta time discretization is identified as a further direction.Source: §8 · p. 18
  • Further directions include adaptive hp-PolyDG methods for localized reaction fronts and heterogeneous media, and applications to more detailed ecological, epidemiological, and neurodegenerative-disease models.These are proposed developments rather than results established in this work.Source: §8 · p. 19

A Recursive coefficients in the block stability estimates

The appendix defines nonnegative coefficient functions used in the stability estimates for each population block. Their definitions expose dependence on preceding blocks and simplify the proof of Theorem 2.

A Recursive coefficients in the block stability estimates

A.1 Coefficients for the L2-energy estimate

The L2-energy estimate's coefficient definitions collect terms arising from Gagliardo–Nirenberg and Young inequalities. The listed definitions combine specific coefficient functions used in the estimate.

A.1 Coefficients for the L2-energy estimate

  • The coefficient functions collect terms produced by Gagliardo–Nirenberg and Young inequalities in the proof of Theorem 2's first estimate.Source: §A.1 · p. 19
  • One defined coefficient is the sum of A(p)2,1(t), A(p)3,1(t), A(p)4,1(t), A(p)3,2(t), and A(p)4,2(t).Source: §A.1 · p. 19
  • Another listed coefficient is A(p)4,3(t).Source: §A.1 · p. 19

A.2 Coefficients for the higher-energy estimate

The higher-energy estimate collects terms that depend only on previously controlled blocks. Those terms are nonnegative and known when the recursive argument reaches block p.

A.2 Coefficients for the higher-energy estimate

  • The higher-energy estimate collects terms depending only on preceding blocks.Source: §A.2 · p. 19
  • The functions in (23) and (24) are nonnegative and depend only on already controlled blocks 1 through p−1.Source: §A.2 · p. 19
  • Because those functions depend only on blocks already controlled, their values are known at the p-th recursive step.Source: §A.2 · p. 19
Loading 2609.34625v1…