PaperlayerView exact revision on arXiv

Source-linked AI summary

Normal Equations and Discrete Energy Structures for High-Order Streamline Diffusion

Erik Burman, Peter Hansbo, Mats G. Larson

arXiv:2609.12506v1math.NAcs.NA

TL;DR

Transient SUPG methods for skew-symmetric first-order systems need stability and error analysis beyond second-order linear multistep schemes. This paper develops a common discrete-energy framework and dual-residual estimate, establishing stability across several time integrators and reporting improved localisation in acoustic tests.

  • Problem

    For multidimensional first-order hyperbolic systems, stability and a priori error analysis for strongly consistent Streamline Diffusion finite elements with linear multistep schemes above order two was not previously available.

    Source: §1 · p. 2

  • Method

    The paper combines a common operator-and-energy framework for fully discrete SUPG with method-dependent discrete energies and a direct dual-residual error estimate.

    Source: §3 · p. 8, §8 · p. 36

  • Results

    Away from moving jumps, normal-equation SUPG has d’Alembert local P1 convergence rate 2.0006 versus 0.4995 for the unstabilised solution; stability and error estimates are also established for the stated time schemes.

    Source: §8 · p. 36

  • Takeaways & Limitations

    The experiments provide evidence of improved localisation with stabilisation, but do not establish a general theorem for nonsmooth solutions.

    Source: §8 · p. 36

    The analysis assumes a time-independent operator; space- or time-dependent coefficients, which require additional commutator bounds, are not covered.

    Source: §2 · p. 4

Abstract

from arXiv · show

We analyse fully discrete Streamline Diffusion finite element methods, also known as SUPG methods, for time-dependent first-order systems with skew-symmetric spatial operators. To the best of our knowledge, this is the first stability and a priori error analysis of strongly consistent Streamline Diffusion finite elements combined with linear multistep schemes of order greater than two for multidimensional first-order hyperbolic systems. Local-in-time $L^2$-residual minimisation singles out $\delta=b_0\tau$, where $\delta$ is the stabilisation parameter, $\tau$ the time step, and $b_0$ the current-time coefficient in the method average. For this choice, each implicit linear system is a symmetric positive definite normal equation in a graph norm. A direct dual-residual estimate yields errors of order $O(h^{k+1/2}+\tau^\nu)$ when $\delta$ is proportional to $h$ and $\tau=O(h)$, where $h$ is the mesh size, $k$ the polynomial degree, and $\nu$ the temporal order. Stability is established for the $\theta$-method and the Adams-Moulton methods AM3-AM5. The general-degree AM3 and AM4 estimates use strengthened CFL conditions; for continuous piecewise affine elements, an elementwise cancellation recovers the standard hyperbolic regime. AM5 and the explicit Adams-Bashforth methods AB3 and AB4 are stable under a standard hyperbolic CFL condition, though the explicit schemes do not have normal-equation structure. Order-matched acoustic tests recover orders two through five and material-residual rates one half-order lower. In a d'Alembert test with discontinuous initial data, normal-equation SUPG recovers second-order local $\mathbb P_1$ convergence away from the wave fronts, whereas the unstabilised error is close to half order. A two-dimensional compact-wave test shows improved localisation with modest dissipation.

1 Introduction

The section frames transient SUPG analysis as a coupled energy and algebra problem for multidimensional first-order systems. It introduces a common framework for high-order linear multistep methods, with method-dependent stability results, a dual-residual error bound, and supporting acoustic experiments.

Introduction

  • Skew-symmetry supports L2 energy stability, but unstabilised error estimates are typically suboptimal and may allow dispersive errors to pollute spatially separated regions.Stabilised methods commonly recover half an order and better preserve local propagation of perturbations.Source: §1 · p. 2
  • Transient consistency tests the full space–time residual against perturbed tests, producing a cross term that complicates energy analysis and general SUPG steps are nonsymmetric.Prior rigorous results either stop at second-order linear multistep schemes or use Runge–Kutta, discontinuous-in-time, or coupled space–time formulations; this work addresses higher-order linear multistep schemes for multidimensional first-order systems.Source: §1 · pp. 2–3

Main Contributions

  • Choosing δ = b0τ from local residual minimisation makes each implicit step a symmetric positive definite normal equation in a graph norm.For general δ the system is not symmetric in general; explicit Adams–Bashforth methods have b0 = 0 and cannot use this structure with positive δ.Source: §1 · p. 2
  • The dual-residual estimate gives O(h^k+1/2 + τ^ν) error when δ ∼ h and τ ≲ h, under the stated regularity assumptions.The proof uses the ordinary L2-projection, skew-symmetry, and discrete summation by parts.Source: §1 · p. 2
  • AM3–AM4 require strengthened CFL conditions for general degree, while affine elements recover the standard regime and AM5, AB3, and AB4 satisfy standard hyperbolic CFL conditions.The framework establishes stability for the θ-method and AM3–AM5; separate energy identities establish it for AB3 and AB4.Source: §1 · pp. 2–3
  • Order-matched acoustic tests recover temporal orders two through five, while the material-residual errors follow a half-order-lower sequence.The d’Alembert test finds second-order local P1 convergence away from wave fronts for normal-equation SUPG, versus order close to one-half unstabilised; the compact-wave test shows improved localisation.Source: §1 · p. 3

2 Abstract Setting and Formulation of the Method

The formulation sets up a time-independent skew-symmetric operator on a Hilbert space and builds discrete time averages, residuals, and test operators around it. The distinguished residual-minimising choice yields a graph-norm normal equation for implicit steps; the analysis also specifies approximation assumptions and acoustic-wave structure.

2 Abstract Setting and Formulation of the Method

  • The abstract setting uses a real Hilbert space H and dense subspace V, with a time-independent skew-symmetric spatial operator; space- or time-dependent coefficients are outside the stated analysis.Such coefficients require additional commutator bounds.Source: §2 · p. 4
  • For acoustic waves, scaling velocity as u = cv gives the system’s product inner product; periodic or reflecting-wall conditions make the boundary pairing vanish and yield skew-symmetry.The text notes that the powers of c, operator, and energy are consistent.Source: §2 · p. 4

2.1 Time Differences, Averages, and Their Adjoints

  • The time grid uses t_n = nτ, and the method average has stencil length ℓ+1 with b0 as the current-value coefficient; earlier stencil values are prescribed start-up data.The finite-interval average is not self-adjoint, so its discrete adjoint is defined using test sequences extended by zero outside computed levels.Source: §2.1 · p. 5
  • Discrete summation by parts retains both endpoint terms, and the time-shifted sequence cannot be replaced by the method average without additional symmetry and endpoint terms.The average and its adjoint are related by an exact change-of-index identity.Source: §2.1 · p. 5
  • The approximation setup uses a shape-regular mesh, conforming degree-k finite elements, and the L2-orthogonal projection to define the local truncation residual.The spatial approximation index is r = min(s, k + 1).Source: §2.1 · p. 5
  • Temporal consistency is expressed through local truncation estimates on stencil intervals, whose uniformly bounded overlap controls the corresponding global time-derivative norms.The dual-residual estimate uses differentiated consistency bounds requiring one additional time derivative and endpoint regularity.Source: §2.1 · pp. 5–6

2.2 Shift and Method Operators

  • The shift and method operators are defined for δ > 0, and skew-symmetry supplies a graph-inner-product identity and associated graph norm.Source: §2.2 · p. 6

2.3 SUPG and Residual-Minimisation Formulations

  • In the residual-minimisation formulation, the Euler–Lagrange equation gives the update, whose left-hand side is symmetric positive definite in the graph norm; general δ need not give symmetry.The formulation allows δ to be chosen independently of b0.Source: §2.3 · pp. 6–7

3 Abstract Stability and Error Analysis

The section builds a method-independent energy and residual framework for stability and error analysis, then uses a dual-residual estimate to obtain a half-order spatial gain. Its parameter balance makes the spatial rate explicit, but accuracy balance alone does not establish stability.

3.1 Sequence Norms and Abstract Stability

  • The abstract energy hypothesis yields a stability bound uniform in the terminal horizon, with start-up data, lifted-dual forcing, and any strong forcing defect treated separately.The energy and dissipation norms are method-dependent; the θ-method has no strong forcing defect.Source: §3.1 · p. 9

3.2 Error Equation and Residual Bounds

  • The error analysis splits approximation and discrete errors, then estimates the residual functional in a lifted dual norm to gain half a spatial order.Projection orthogonality and energy-controlled quantities support the dual estimate.Source: §3.2 · pp. 9–11

3.3 Parameter Balance

  • For the standard stabilised choice δ ∼h and τ ≲h, the spatial part is order h^(r−1/2); with r = k +1, smooth degree-k approximation gives order h^(k+1/2) under stated regularity assumptions.Source: §3.3 · p. 12
  • For the normal-equation choice δ = b0τ, balancing the dual contributions is an accuracy choice, so stability must still be checked for the time-stepping method.The Adams–Moulton theorem includes this choice under method-dependent CFL conditions.Source: §3.3 · p. 12

4 An A-Stable Scheme: The Theta Method

The θ-family obtains an energy balance from a one-step identity without a CFL condition. Its temporal order and resulting error rate depend on θ.

4 An A-Stable Scheme: The Theta Method

  • For 1/2 ≤ θ ≤ 1, the θ-method satisfies the abstract energy hypothesis, yielding stability without a strong forcing defect.Its energy–dissipation identity converts increments into telescoping energy differences and nonnegative dissipation.Source: §4 · p. 13

4.1 Consistency and Error Estimate

  • The method has temporal order one for fixed θ > 1/2 and order two for θ = 1/2 or θ = 1/2 + O(τ).These orders carry through to the error estimate under the stated regularity and projected-start-up assumptions.Source: §4.1 · pp. 13–14
  • Under δ ∼ h and τ ≲ h, the graph-norm error is O(h^(k+1/2) + τ^ν), with ν set by the θ-dependent temporal order.The θ-family includes Crank–Nicolson at θ = 1/2 and backward Euler at θ = 1.Source: §4 · p. 12, §4.1 · p. 14

5 Adams–Moulton Methods of Orders Three to Five

The Adams–Moulton analysis establishes consistency and method-specific energy estimates for AM3–AM5. Stability and error rates depend on the method, polynomial degree, and CFL regime.

5.1 Consistency of the Adams–Moulton Methods

  • AMν approximates the exact time derivative by interpolating it with a polynomial of degree at most ν−1, giving consistency orders ν = 3, 4, 5.The residual and its first forward difference also receive the uniform and square-summable controls required by the dual-residual estimate.Source: §5.1 · pp. 14–15

5.2 Stability Assumptions and Increment Bounds

  • For general-degree AM3 and AM4, the proof uses strengthened CFL conditions; AM5 requires no strengthened condition.For continuous piecewise affine elements with constant-coefficient G, elementwise constant gradients allow cancellation that avoids the recurrence producing stronger CFL powers.Source: §5.2 · pp. 15, 17

5.3 Exact AM3–AM5 Energy Decompositions

  • The energy decompositions give AM3 and AM4 negative bulk terms in their highest relevant increments, whereas AM5 has a positive third-increment term.For AM5, that positive term can absorb the negative fourth-increment term under the standard CFL restriction.Source: §5.3 · pp. 19–20

5.4 Lower Bounds for the Energy Decompositions

  • The lower-bound analysis controls the bulk terms under the relevant conditions and bounds the endpoint terms when c0 is sufficiently small.These bounds are used in the subsequent stability estimates.Source: §5.4 · pp. 20–21

5.5 Stability and Error Estimates

  • AM3 and AM4 require strengthened CFL conditions for general degree, while AM5 needs only the standard hyperbolic condition; AM3–AM4 also reach that regime for continuous piecewise affine elements.The piecewise affine relaxation uses an elementwise cancellation; the error estimate has temporal order three, four, or five for AM3, AM4, or AM5, respectively.Source: §5.5 · pp. 21–23
  • Under δ ∼ h and τ ≲ h, the error rate is O(h^(k+1/2) + τ^ν) under the applicable general-degree stability condition or piecewise affine relaxation.The relaxation depends on affine elements and constant-coefficient first-order operators; higher degrees retain strengthened conditions in this proof.Source: §5.5 · p. 23

6 Explicit Schemes: Adams–Bashforth Methods

AB3 and AB4 are explicit, lack the normal-equation structure, and are analyzed through consistency estimates, energy decompositions, and increment control. Under the stated scaling and stability conditions, the analysis yields graph-norm error rates for both methods.

6.4 Stability and Error Estimates

  • AB3 and AB4 are stable under the standard hyperbolic restriction τ ≲ h and 0 < δ ≲ h; absorption controls graph-energy maxima and the material-residual sum without discrete Gronwall.The theorem assumes the inverse inequality and the stated condition with sufficiently small c_AB.Source: §6 · p. 23, §6.4 · p. 27

6.1 Consistency of the Adams–Bashforth Methods

  • AB3 and AB4 have consistency orders three and four, respectively, with integral-remainder bounds providing the residual control used in the dual-residual estimate.The estimates control both the residual and its first forward difference.Source: §6.1 · p. 24

6.2 Exact AB3 and AB4 Energy Decompositions

  • The exact AB3 and AB4 energy decompositions pair positive lower-increment bulk terms with negative highest-increment terms that the increment estimates must absorb.For AB4, a relatively large negative fourth-increment coefficient is balanced by a positive third-increment term.Source: §6.2 · p. 24, §6 · p. 24

6.3 Increment and Boundary Control

  • Because explicit updates have no current spatial term, differencing controls the highest increment from the preceding increment with a small factor τ/h.The boundary analysis bounds start-up contributions and absorbs final increment terms into the leading final energy.Source: §6.3 · pp. 25–27

6.4 Stability and Error Estimates

  • O(h^(k+1/2) + τ^ν) under δ ∼ h and τ ≲ h, with ν = 3 for AB3 and ν = 4 for AB4.The error estimate follows from the stability result and the stated regularity hypotheses with projected start-up data.Source: §6.4 · p. 28
  • For AB2, arbitrary-degree stability requires τ ≲ h^(4/3), while continuous piecewise affine elements recover τ ≲ h and balanced second-order space-time errors.The piecewise affine result uses the elementwise cancellation described in the paper.Source: §6.4 · p. 28

7 Numerical Experiments

The experiments test formal temporal accuracy, order-matched convergence, and error localisation for smooth and discontinuous acoustic waves. Their refinement rates are empirical and do not establish universal stability thresholds or a general localisation theorem.

7 Numerical Experiments

  • The experiments ask whether integrators attain formal order, whether matched temporal and spatial orders produce predicted rates, and whether stabilisation confines errors from discontinuities and compact waves.Observed refinement rates are empirical and do not establish universal stability thresholds or a general localisation theorem.Source: §7 · p. 28

7.1 Implementation and Test Problems

  • The experiments use periodic acoustic problems, consistent mass matrices, and sufficient-order Gaussian quadrature; no mass lumping, filtering, limiter, added viscosity, or refinement-path retuning is used.When δ = b0τ, assembly reproduces the normal-equation cancellation to floating-point precision.Source: §7.1 · pp. 28–29

7.2 Temporal Accuracy

  • Every tested time integrator attains its formal temporal order; the AM5 final refinement reaches a floating-point error floor and is excluded from the mean rate.Table 1 reports mean rates over the three finest refinements unaffected by roundoff.Source: §7.2 · p. 29, Table 1 · p. 29

7.3 Order-Matched One-Dimensional Convergence

  • Final L2 and maximum graph errors recover orders 2, 3, 4, and 5 with order-matched one-dimensional methods; material-residual rates are 3/2, 5/2, 7/2, and 9/2.The method–space pairs are Crank–Nicolson/P1, AM3/P2, AM4/P3, and AM5/P4.Source: §7.3 · p. 30, Table 2 · p. 30

7.4 Two-Dimensional Convergence

  • P1 two-dimensional tests yield graph rates near 2 and material-residual rates near 3/2; multistep final-L2 rates near 2.87 are structured-mesh superconvergence, not general P1 rates.The rates are reported for the oblique plane wave on slash-diagonal production meshes.Source: §7.4 · p. 31, Table 3 · p. 31
  • Across tested mesh orientations, CN graph rates are 2.015–2.095 and material-residual rates 1.502–1.506; AM5 rates are 2.056–2.168 and 1.493–1.498.Error levels depend on orientation, but the tested graph and material-residual rate trends remain robust.Source: §7.4 · p. 31

7.5 One-Dimensional Localisation of Discontinuities

  • Local error rates are 0.4995 without stabilisation and 2.0006 with normal-equation SUPG; finest-mesh local errors are 5.6971 × 10^-3 and 3.9288 × 10^-7, respectively.The smooth-region result is numerical evidence for this linear one-dimensional configuration, not a general theorem for nonsmooth solutions.Source: §7.5 · p. 32, Table 4 · p. 32

7.6 Two-Dimensional Localisation of a Compact Wave

  • The stabilised compact-wave run has exterior-energy fraction 1.2126 × 10^-7 versus 4.7772 × 10^-6 unstabilised, with final energy ratio 0.993225 versus 1.0000.This configuration shows improved localisation with modest dissipation; no general multidimensional result is claimed.Source: §7.6 · p. 34, Figure 5 · p. 35

8 Conclusions

The conclusions establish a common energy framework for high-order time integrators, with method-dependent stability conditions and a dual-residual error estimate. Experiments show order recovery and improved localisation in the tested settings, while several theoretical extensions remain open.

8 Conclusions

  • For implicit methods, δ = b0τ yields a symmetric positive definite normal equation associated with graph-norm residual minimisation; explicit Adams–Bashforth methods do not have this structure.The common formulation separates the method average, backward time difference, material residual, and perturbed test operator.Source: §8 · p. 36
  • Stability is established for the θ-method, AM3–AM5, AB3, and AB4, with strengthened CFL restrictions for general-degree AM3 and AM4 and the standard hyperbolic regime recovered for continuous piecewise affine elements.AM5, AB3, and AB4 use standard hyperbolic CFL restrictions.Source: §8 · p. 36
  • The dual-residual argument gives O(h^k+1/2 + τ^ν) error for δ ∼ h and τ ≲ h, subject to stability and regularity assumptions.The proof uses the ordinary L2-projection, skew-symmetry, and discrete summation by parts.Source: §8 · p. 36
  • Order-matched acoustic tests recover orders two through five in final-time and graph errors, while material-residual rates follow the half-order sequence 3/2, 5/2, 7/2, and 9/2.The two-dimensional P1 tests give graph and material-residual rates close to two and 3/2, respectively.Source: §8 · p. 36
  • The discontinuous d’Alembert test gives local P1 rates 0.4995 without stabilisation and 2.0006 with normal-equation SUPG away from moving jumps; the stabilised finest-mesh local error is about 1.45×10^4 times smaller.In the smooth two-dimensional compact-wave test, the stabilised exterior-energy fraction is 39.4 times smaller, with modest dissipation; these experiments are not a general theorem for nonsmooth solutions.Source: §8 · p. 36
  • The stated extensions include variable coefficients, nonperiodic boundary conditions, and quantitative local error estimates for nonsmooth solutions in multidimensional and nonlinear systems.The conclusions also identify robust multistep start-up procedures and efficient solvers and preconditioners as computational priorities.Source: §8 · p. 36

A Proof of the Dual Residual Estimate

The appendix proves the dual-residual estimate using the ordinary L2-projection, skew-symmetry, and discrete summation by parts.

A Proof of the Dual Residual Estimate

A.1 Projection Contribution

The projection contribution uses the time-commutation of the spatial L2-projection and skew-symmetry, retaining both endpoint terms in the subsequent bound.

A.1 Projection Contribution

  • Because Πh is the spatial L2-projection, it commutes with the time operators.Source: §A.1 · p. 37
  • Skew-symmetry of G is then used in the projection-contribution argument.Source: §A.1 · p. 37
  • Both endpoint terms are retained before applying Cauchy–Schwarz and the stated bounds.Source: §A.1 · p. 38

A.2 Truncation Contribution

The truncation contribution is bounded by expanding Pδ, applying discrete summation by parts, and combining the resulting estimates with consistency bounds.

A.2 Truncation Contribution

  • The truncation contribution is treated separately by setting its relevant quantity.Source: §A.2 · p. 38
  • Expanding Pδ and applying the same summation-by-parts formula yields the truncation estimate.Source: §A.2 · p. 38
  • The consistency estimates and 0 < δ ≤ 1 bound the resulting contributions; combining the projection and truncation residual bounds and taking the supremum proves (3.32)–(3.33).Source: §A.2 · p. 38
Loading 2609.12506v1…