Source-linked AI summary

A class of high-order discontinuous-Galerkin methods satisfying infinitely many entropy conditions with provable error estimates and strong convergence for general nonlinear conservation laws

Yuanzhe Wei, Chi-Wang Shu

arXiv:2609.04687v1math.NAphysics.comp-ph

TL;DR

High-order DG methods need to combine accuracy with entropy consistency and convergence for nonlinear conservation laws, especially when solutions develop discontinuities. The paper derives DG schemes through operator semigroups and nonlinear projections, extending them to systems and multiple dimensions. The resulting methods satisfy infinitely many entropy inequalities, achieve optimal smooth-solution error estimates, and strongly converge for discontinuous strictly convex scalar problems.

  • Problem

    High-order DG methods must retain accuracy while satisfying enough entropy conditions to select the physically relevant solution and support convergence for nonlinear conservation laws.

  • Method

    The paper uses operator semigroups and a nonlinear projection with cellwise damping to construct entropy-preserving DG schemes, extending the framework to systems and arbitrary space dimensions.

  • Results

    The schemes satisfy infinitely many local entropy inequalities, attain optimal L2 error estimates for smooth scalar solutions, and strongly converge for discontinuous strictly convex scalar solutions.

  • Takeaways & Limitations

    The framework provides a unified high-order DG construction combining entropy control with accuracy and convergence results across scalar, system, and multidimensional settings.

  • Takeaways & Limitations

    The analysis is semi-discrete, assumes a uniform L∞ bound for compensated-compactness convergence, and does not yet establish strong convergence or optimal-error results for systems and multidimensional problems.

Abstract

from arXiv · show

We propose a novel framework for deriving semi-discrete discontinuous-Galerkin (DG) methods using operator semigroups for scalar conservation laws, then apply it to construct a class of high-order OFDG-type schemes [13] satisfying infinitely many local entropy inequalities with general E-fluxes on non-uniform meshes. Such schemes are further generalized to systems of conservation laws in any number of space dimensions by using entropy stable numerical fluxes in the sense of [1]. Finally, we prove optimal error estimates for smooth solutions to nonlinear scalar conservation laws, and prove strong convergence for discontinuous solutions to strictly convex conservation laws via compensated compactness.

1. Introduction

The paper develops an operator-based high-order DG framework that enforces infinitely many entropy inequalities while retaining conservation, extending to systems and multiple dimensions. It proves optimal smooth-solution error estimates and strong convergence for discontinuous strictly convex scalar problems.

  • Hyperbolic conservation laws can admit infinitely many weak solutions, with entropy inequalities selecting the physically relevant entropy solution.
  • The operator framework interprets classical DG through exact entropy-solution evolution followed by projection onto the DG space.The exact evolution carries the entropy structure, while projection determines which properties survive discretization.
  • A nonlinear projection adds cellwise damping on nonconstant modes, preserving conservation while enforcing local inequalities for prescribed and, with uniform bounds, infinitely many convex entropies.
  • The construction extends to systems in one dimension and arbitrary space dimensions using strictly convex entropy pairs and entropy-stable numerical fluxes.
  • The schemes achieve optimal L2 error estimates for smooth scalar solutions and strong convergence for discontinuous solutions of strictly convex scalar conservation laws.The convergence proof uses compensated compactness and assumes a uniform L∞ bound on numerical solutions.
  • The paper derives multiple scalar, system, and multidimensional scheme variants, including simplified OFDG-related damping mechanisms, and reports numerical tests of accuracy and robustness.

2. An operator framework for the classical DG method

This section derives classical semi-discrete DG from an operator semigroup rather than directly from the weak form. The exact entropy-solution operator followed by cellwise projection recovers the Godunov-flux method and explains its stability and entropy properties.

  • The semi-discrete DG approximation evolves in a finite-dimensional polynomial space as an ODE with a locally Lipschitz right-hand side and flow map Φ_t.
  • The operator construction applies the exact entropy-solution operator S_t and then the cellwise L2 projection π onto the DG polynomial space.
  • As Δt→0, the composite operator πS_Δt recovers classical DG with the Godunov flux.
  • The derivation isolates interface effects by tracking short-time Riemann evolutions whose waves remain in localized regions near cell boundaries.
  • Because the exact evolution satisfies entropy inequalities and projection decreases the cellwise L2 norm, the resulting scheme admits a local square-entropy inequality.
  • Under periodic or compactly supported boundary conditions, the Godunov-flux classical DG flow is L2-stable: ∥Φ_t2u∥L2(I)≤∥Φ_t1u∥L2(I) for 0≤t1≤t2.

3. The Operator DG method for 1D scalar conservation laws

The paper derives a nonlinear projection-based Operator DG method whose damping preserves local entropy inequalities while retaining conservation and high-order accuracy. Suitable damping bounds yield infinitely many entropy inequalities, including for general E-fluxes and related schemes.

  • Operator construction: The nonlinear projection contracts square entropy and a prescribed strictly convex entropy, producing DG damping on nonconstant modes.The projection preserves cell averages and decreases high-frequency modes.
  • Operator construction: The resulting method is the classical DG spatial operator plus a cellwise damping remainder, with computable upper bounds for its coefficient.Replacing the exact remainder by an upper bound can only decrease the entropy production rate.
  • Accuracy: The damping bounds are chosen tightly enough that the resulting scheme retains optimal-order accuracy.The construction estimates the remainder using polynomial norm bounds and Jensen-type inequalities.
  • Entropy stability: For a consistent E-flux, Scheme B satisfies a local entropy inequality for the prescribed entropy, including monotone fluxes such as Godunov and local Lax-Friedrichs fluxes.The cited E-flux class also includes HLL fluxes with suitable wave-speed bounds.
  • Entropy stability: With the square entropy, Scheme B satisfies at least two local entropy inequalities.The square-entropy inequality follows as a separate corollary of the projection and damping construction.
  • Infinitely many entropy inequalities: For the analytic entropy choice in Corollary 3.13, Scheme B satisfies infinitely many entropy inequalities, and the same conclusion extends to Scheme C.Scheme D also has infinitely many entropy inequalities when the numerical solution remains bounded.

4. The Operator DG method for 1D systems of conservation laws

The Operator DG construction extends to one-dimensional systems by combining strictly convex entropy pairs with entropy-stable numerical fluxes. The resulting schemes satisfy local entropy inequalities, including simplified variants under a uniform boundedness condition.

  • System extension: For one-dimensional systems, the method uses a general strictly convex entropy pair and entropy-stable numerical fluxes.The entropy flux and entropy-stability condition are formulated for vector-valued conservation laws.
  • System extension: The system scheme adds damping to the deviation from the element average while retaining the DG formulation.The vector-valued construction uses estimates analogous to the scalar case.
  • Entropy stability: HLL fluxes, including local Lax-Friedrichs fluxes, belong to the entropy-stable class used by the system construction.This supplies concrete flux choices for the entropy-stable formulation.
  • Entropy stability: If the numerical flux is entropy stable, Scheme B′ satisfies a local semi-discrete entropy inequality.The theorem applies to any entropy-stable flux satisfying the stated condition.
  • Simplified schemes: The entropy result also holds for Scheme C′ and for Scheme D′ with sufficiently large damping parameters when the numerical solution is uniformly bounded.The extension follows from monotonicity of the remainder for systems.

5. The Operator DG method for higher dimensional systems of conservation laws

The framework generalizes to systems in arbitrary space dimensions on convex, shape-regular meshes. Directional entropy-stable fluxes supply interface entropy stability, while elementwise damping yields the corresponding local entropy inequality.

  • Multidimensional formulation: The multidimensional formulation treats systems of conservation laws on convex, shape-regular meshes using piecewise polynomial DG spaces.The mesh consists of convex elements with geometric regularity measured through element diameters and inscribed spheres.
  • Multidimensional formulation: The method adds elementwise damping based on the deviation of the solution from its element average.The damping coefficient depends on entropy derivatives and the minimum Hessian eigenvalue over the element’s state convex hull.
  • Flux conditions: Directional numerical fluxes are required to be consistent, conservative, and entropy stable for the selected entropy.The multidimensional entropy condition is defined relative to the outward normal of each element boundary.
  • Entropy stability: For any entropy-stable directional flux, Scheme B″ satisfies a local semi-discrete entropy inequality on each element.The inequality is expressed through the corresponding numerical entropy flux integrated over the element boundary.
  • Simplified schemes: Scheme D″ satisfies the same entropy theorem for sufficiently large damping parameters when the numerical solution is uniformly bounded.The multidimensional construction retains the boundary integration of the damping-control quantity.

6. Optimal error estimates to smooth solutions of general nonlinear scalar conservation laws

The paper proves optimal error estimates for the proposed schemes on quasi-uniform meshes when the exact scalar solution remains smooth, covering both upwind and general monotone fluxes. Numerical tests in one and two dimensions agree with the predicted optimal convergence orders.

  • Theoretical estimates: Theorem 6.1 establishes optimal error estimates for Schemes B, C, and D applied to smooth solutions of general nonlinear scalar conservation laws.The result assumes a quasi-uniform mesh, sufficiently smooth flux, smooth exact solution on a fixed time interval, and projected initial data.
  • Theoretical estimates: O(h^(k+1)) is obtained with an upwind numerical flux for polynomial degree k.The estimate is stated for sufficiently small mesh size h under the theorem’s smooth-solution assumptions.
  • Proof strategy: For the upwind-flux analysis, the modified Gauss–Radau projection selects interface conditions according to the sign of the initial characteristic speed.The projection definition is unchanged when the initial sign condition is replaced by the time-dependent characteristic speed while the exact solution remains smooth.
  • Numerical tests: One-dimensional Burgers tests for Scheme D on smooth periodic solutions match Theorem 6.1’s predicted orders across polynomial degrees k=1,2,3.The experiments use uniform meshes and compare cases r=k+1 and r=k+2; the reported L1, L2, and L∞ results are shown in Tables 3 and 4.
  • Numerical tests: Two-dimensional scalar Burgers tests also show optimal convergence orders for Scheme D” with polynomial degrees k=1,2,3.The tests use smooth periodic solutions on uniform rectangular meshes and report L1, L2, and L∞ errors for r=k+1 and r=k+2.

7. Strong convergence through compensated compactness for general strictly convex conservation laws

Under strict convexity, uniform L∞ bounds, suitable initial convergence, and the stated mesh, flux, and damping assumptions, the proposed scheme converges strongly to the unique entropy solution. The proof establishes entropy dissipation and compact entropy productions, then applies compensated compactness.

  • Strong-convergence theorem: Theorem 7.4 proves strong convergence of Scheme D solutions to the unique entropy solution as h→0.The result assumes a strictly convex flux, ε1 > 0, periodic or compactly supported boundary conditions, and strongly convergent initial projections.
  • Assumptions: The analysis assumes quasi-uniform meshes, fixed polynomial degree k≥1, and uniform L∞ boundedness of the numerical solutions.The initial projections must converge strongly in L2, and the boundary conditions are periodic or compactly supported.
  • Entropy estimates: The damping coefficient requires ε1 > 0 because the first term supplies estimates needed for the entropy-dissipation and consistency bounds.The same term is identified as necessary for the subsequent compactness argument.
  • Compensated compactness: The proof establishes local and global square-entropy inequalities together with compact entropy productions in negative Sobolev spaces.The production sequences are shown relatively compact in W^-1,p and then H^-1 using the stated lemmas and estimates.
  • Limit identification: Compensated compactness yields a distributional limit satisfying the limiting square-entropy inequality, while strict convexity guarantees uniqueness and convergence of the entire family.The argument identifies the Young-measure limit with a Dirac mass and uses Vitali’s theorem to obtain strong convergence.

8. Numerical results

The numerical experiments examine damping choices, non-convex scalar problems, and demanding Euler tests. They show that the flux-derivative contribution can be essential for correct entropy solutions in non-convex cases, while simplified damping robustly handles several strong-shock computations.

  • Damping selection: The first damping term is indispensable for the proved entropy inequalities and strong convergence, although its necessity varies across numerical tests for non-convex problems.The section explicitly distinguishes the theoretical requirement from cases where ε1 = 0 performs adequately computationally.
  • Non-convex scalar tests: For the Buckley–Leverett problem, ε1 = 0 with ε2 = 10 fails to produce the correct solution for k=1,2,3,4.The experiment uses Scheme D, N=275, and a global Lax–Friedrichs flux.
  • Non-convex scalar tests: Increasing ε1 restores convergence for the Buckley–Leverett problem with r=k+1 and r=k+2.The corresponding solutions are reported to converge correctly to the entropy solution.
  • Damping selection: For k=4, the parameters ε1 = 5.0 × 10^-6 and ε2 = 1.0 × 10^-6 are reported as sufficient in one reduced-damping test.The observed parameter ratios are consistent in order of magnitude with the theoretical relation discussed in the section.
  • System and multidimensional tests: In the Sod, Shu–Osher, double Mach reflection, and shock-diffraction computations, the simplified system schemes resolve strong shocks and fine-scale structures.The multidimensional Euler tests require no positivity-preserving limiter for the reported parameters.

9. Concluding remarks

The paper’s operator framework produces entropy-controlled high-order DG schemes that preserve conservation, extend to systems and multidimensional meshes, and support both accuracy estimates and strong convergence results. Numerical experiments further assess damping and shock-resolution behavior, while several theoretical and computational extensions remain open.

  • Framework and schemes: The framework derives classical semi-discrete DG through an evolution operator followed by L2 projection, then replaces linear projection with entropy-contracting nonlinear projection.The resulting correction damps only nonconstant polynomial modes and preserves cell averages and conservation.
  • Framework and schemes: The constructed schemes cover one-dimensional scalar and system equations and multidimensional systems on general convex, shape-regular meshes.The schemes are labeled B–D, B′–D′, and B″–D″ across these settings.
  • Analytical properties: Scalar schemes with E-fluxes satisfy infinitely many local entropy inequalities, while system schemes accommodate prescribed strictly convex entropies and entropy-stable fluxes.The added dissipation is compatible with optimal smooth-solution accuracy and strong convergence for discontinuous strictly convex scalar problems.
  • Numerical evidence: Numerical experiments indicate high-order accuracy in smooth tests and robust shock resolution, while damping components can determine correctness for non-convex scalar problems.The interface contribution alone sufficed in the reported multidimensional Euler tests.
  • Open problems: Open directions include fully discrete theory, removal of the uniform L∞ assumption, broader strong-convergence results, general non-convex entropy selection, alternative projections, and more efficient damping evaluation.Adaptive selection of ε1 and ε2 is also identified as a computational direction.

Disclosure on the use of artificial intelligence

The disclosure states that ChatGPT assisted with language editing, part of the compensated-compactness proof, Julia test codes, and coefficient computations, while the human authors independently verified the mathematical work.

  • Uses: ChatGPT was used to improve readability and language quality.
  • Uses: ChatGPT assisted with the compensated-compactness proof, including the proof of Lemma 7.3.
  • Uses: ChatGPT assisted with programming Julia codes for the numerical tests and computing optimal coefficients recorded in the tables.
  • Verification: The human authors independently verified and validated all mathematical derivations, results, and proofs.

Appendix A Operator for generating the classical DG method with HLL flux

The appendix extends the operator-semigroup interpretation of DG from the Godunov flux to HLL fluxes, including local Lax–Friedrichs, while preserving the entropy conclusions needed by the scheme.

  • Operator construction: The operator construction replaces the exact solution operator with a modified S_HLL^Δt that yields the HLL flux.The construction modifies the entropy solution only near interfaces over regions of length O(Δt).
  • Operator construction: For sufficiently small Δt, interface regions remain disjoint, allowing S_HLL^Δt to be defined locally from Riemann-solution wave-speed bounds.The relevant bounds are a− = min{0, λ−} and a+ = max{0, λ+}.
  • Entropy properties: Averaging the modified solution preserves the required entropy property because Jensen’s inequality does not increase total entropy.The construction also retains the properties used in the analogue of Theorem 2.1.
  • Consequences: The resulting DG method inherits the conclusions of Corollary 2.2 with the Godunov flux replaced by the HLL flux.The proof is unchanged once the modified operator satisfies the required properties.
  • Consequences: Because the scheme uses only each cell and its two neighbors, it also satisfies a local square-entropy inequality for the HLL flux.The appendix further states that Corollary 2.3 extends to this case.

Appendix B Improved estimate of 𝜷𝒌from Lemma 3.4

The appendix replaces the earlier β_k bound with an explicit computable estimate derived from a generalized eigenvalue problem, using symmetry to cover both endpoints.

  • Motivation: The previous bound β_k ≤ (k + 1)^4 combines two sharp inequalities whose optimizers do not coincide.Using a single inequality yields a sharper estimate.
  • Derivation: The improved constant is obtained from an endpoint evaluation problem in the finite-dimensional polynomial space with a zero-average constraint.The supremum is attained because the admissible space is finite dimensional.
  • Derivation: The same constant applies at the left and right endpoints by the transformation p(ξ) ↦ p(−ξ), which preserves the constraint and L2 norm.Thus the endpoint estimates are related by symmetry.
  • Result: β_k ≤ 2A_k, where A_k is computed from a k × k symmetric generalized eigenvalue problem.The eigenvalue problem uses matrices M(k) and G(k).
  • Result: Table B.1 compares the improved computable bounds β_k ≤ 2A_k with the earlier bound (k + 1)^4.The table reports values obtained by solving the generalized eigenvalue problem.
Loading 2609.04687v1…