A Positivity-Preserving Relaxation Algorithm

Abstract: We combine Patankar-type methods with suitable relaxation procedures that are capable of ensuring correct dissipation or conservation of functionals such as entropy or energy while producing unconditionally positive and conservative approximations. To that end, we adapt the relaxation algorithm to enforce positivity by using either ideas from the dense output framework when a linear invariant must be preserved, or simply a geometric mean if the only constraint is positivity preservation. The latter merely requires the solution of a scalar nonlinear equation while former results in a coupled linear-nonlinear system of equations. We present sufficient conditions for the solvability of the respective equations. Several applications in the context of ordinary and partial differential equations are presented, and the theoretical findings are validated numerically.

• The paper presents a unified relaxation algorithm that combines Patankar-type updates with correction steps to enforce unconditional positivity, conservation, and entropy stability.
• It employs a bootstrapping technique using NSARK schemes to maintain high-order accuracy while correcting for nonlinear invariants such as entropy and energy.
• Extensive numerical experiments demonstrate improved error control and preservation of key physical properties in simulations of stiff reaction networks, fluid flows, and other systems.

Positivity-Preserving Relaxation for High-Order Time Integration

Problem Context and Motivation

The integration of positivity preservation and conservation/dissipation of nonlinear functionals, especially entropy, is central to robust numerical simulation of ODEs and PDEs arising in fields such as chemical kinetics, population dynamics, and fluid mechanics. Standard time-stepping methods—including implicit/explicit Runge–Kutta (RK), multistep, and projection techniques—face significant trade-offs between positivity, conservation, stability, and order. Patankar-type methods and their modern extensions (e.g., MPRK, MPSSPRK, gBBKS, GeCo) offer unconditional positivity for production-destruction-rest systems (PDRS), but their incorporation with entropy conservation/dissipation, particularly at high order and for general nonlinear invariants, is nontrivial.

Overview of Methodology

This work presents a formal relaxation algorithm that couples Patankar-type methods with relaxation procedures to unconditionally enforce positivity and conservation or correct dissipation of targeted functionals—including (but not limited to) entropy and energy. It builds upon the relaxation framework for ODEs/PDEs introduced by Ranocha et al. [ranocha2020general], extending it to non-standard additive RK (NSARK) schemes whose coefficients may depend on the solution and time step, specifically encompassing MPRK and related schemes as prevalent in stiff PDRS integration.

Relaxation Framework and Algorithmic Structure

Given a high-order method advancing the state $u^n \mapsto u^{n+1}$, the relaxation step corrects $u^{n+1}$ to $\tilde{u}^{n+1}$ such that a chosen functional $\eta$ is exactly dissipated or conserved:

• For dissipative settings, $\eta(\tilde{u}^{n+1}) \leq \eta(u^n)$.
• For conservative settings, $\eta(\tilde{u}^{n+1}) = \eta(u^n)$. The update is constructed as an affine or, for stronger positivity, geometric or nonlinear interpolation between $u^n$ and $u^{n+1}$, i.e., $$\tilde{u}^{n+1} = (u^{n+1})^{\gamma} (u^n)^{1-\gamma}$$ or similarly structured as a Patankar-weighted output. The value of $\gamma$ is determined by solving a scalar or coupled linear-nonlinear system encoding the desired property for $u^{n+1}$0.

A key technical insight is that, within this relaxation, choices for the interpolation (or "dense output") must themselves be designed to preserve positivity unconditionally, even when $u^{n+1}$1—a limitation of prior affine-based approaches. The geometric mean is posited for pure positivity, while bootstrapped, linearly implicit update formulae with positivity-preserving weights are crafted when preservation of linear invariants is necessary.

Order Preservation and Bootstrapping

The authors prove that, through careful tracking of local expansion, the overall relaxation method retains the design order $u^{n+1}$2 of the underlying scheme, provided certain technical conditions on the functional $u^{n+1}$3 and the baseline NSARK scheme are met. To generalize to arbitrarily high order, a bootstrapping technique constructs the dense output for the relaxation step by nesting lower-order positive update rules, thus extending unconditional positivity and conservation/dissipation to high-order integrators (including third-order MPRK, MPSSPRK, etc.).

Sufficient Conditions for Solvability

The paper provides sufficient conditions for solvability of the scalar or coupled nonlinear relaxation equation. For convex $u^{n+1}$4 (e.g., entropy-like functionals), existence and uniqueness of a suitable $u^{n+1}$5 near unity are shown, ensuring both algorithmic practicality and stability.

Numerical Experiments

A series of tests on ODE and PDE systems underscore the practical implications:

• Lotka–Volterra System: The relaxation algorithm achieves linear error growth for invariant preservation, contrasting with the quadratic growth from positivity-only MPRK methods. Entropy is conserved up to solver tolerance.
• Stratospheric Reaction Problem: The relaxation-enhanced MPRK is the only considered scheme to simultaneously maintain both nonlinear positivity and all underlying linear invariants, with significant improvement in long-time accuracy.
• Linear Advection Equation: For several convex entropy fluxes (logarithmic, geometric, harmonic mean), the NSARK-relaxation pairing successfully maintains both positivity and the corresponding entropy conservation, unlike classical explicit/implicit central difference methods.
• Isothermal Euler and Porous Medium Equations: The algorithm robustly preserves positivity of density (Euler) and supports physically meaningful, nonnegative solutions for nonlinear diffusion (porous medium), even in challenging stiff regimes where standard IMEX or explicit methods fail.
• Higher-Order Methods: The bootstrapping approach is shown to yield third-order positivity-preserving, entropy-conservative updates for stiff PDRS and hyperbolic PDE examples.

Strong results in all cases include unconditionally positive, invariant-respecting solutions without observable loss of order, and with automatic adaptation of the interpolation parameter (Newton, secant, or Regula Falsi solvers for the nonlinear equation are compared).

Implications and Integration with Modern Structure-Preserving Discretization

This framework advances the program of structure-preserving time integration, and is compatible with spatial discretizations based on entropy-conservative/entropy-stable fluxes and summation-by-parts (SBP) finite difference, finite volume, or DG schemes. Critically, it complements recent work on high-order, fully discrete entropy-conservative methods (e.g., ADER-DG [GABURRO2023127644], SBP-DG [FWDGP2018]) by providing a general strategy for enforcing positivity without conditional time-step restrictions or order reduction.

Additionally, the bootstrapping and "flux-balanced" versions of MPRK/NSARK emerging from this work may be transferred to other settings where nonlinear invariants, such as multiple species mass or energy, must be exactly or approximately preserved—opening possibilities in plasma physics, reactive flows, and uncertain quantification under random perturbations.

Outlook and Directions for Future Research

This methodology establishes a flexible component for black-box structure preservation in stiff large-scale systems. Open directions include:

• Expansion to implicit or semi-implicit NSARK methods beyond the currently linearly implicit Patankar-type class.
• Systematic analytic and numerical stability analysis of the full suite of relaxation-updated NSARK schemes, with particular attention to long-time nonlinear dynamics, oscillatory Hamiltonian PDEs, and "stiff accuracy."
• Integration with modern adaptive time stepping and error control mechanisms leveraging, e.g., Bayesian optimization or embedded estimators [IR2023].
• Extension to more general classes of invariants, including multiple nonlinear constraints, via multi-dimensional relaxation systems.

Conclusion

The positivity-preserving relaxation algorithm for NSARK/multistage integrators described in the paper (2604.02308) rigorously addresses the joint requirement of unconditional positivity and conservation or dissipation of entropy-like functionals for high-order time integration. Its modular construction, theoretical guarantees, and robust numerical demonstration mark a significant advance in the field of structure-preserving numerical integration for ODE/PDE systems, with clear pathways for further generalization and integration into large-scale, multiphysics simulation codes.