PaperlayerView exact revision on arXiv

Source-linked AI summary

An Optimal IPDG Scheme for the Biharmonic Equation

Bohua Zhang, Xia Ji, Shuo Zhang

arXiv:2609.12675v1math.NAcs.NA

TL;DR

Large penalties in standard IPDG methods can cause locking by overconstraining low-degree approximations, while optimal conforming spaces are unavailable for k = 3, 4 on general grids. The paper uses vertex continuity and degree-matched projected penalties, whose infinite-penalty limit is an optimal nonconforming space. It proves optimal error estimates and reports locking-free behavior for very large penalties.

  • Problem

    Standard IPDG methods can lock when large penalties force solutions into overconstrained spaces, and optimal C1-conforming spaces for k = 3, 4 are not generally available on unstructured grids.

    Source: §1 · p. 2

  • Method

    For k = 3, 4, the method imposes vertex continuity and penalizes projected function and normal-derivative jumps in P^{k-3} and P^{k-2}, respectively.

    Source: §1 · p. 3, §2.1 · p. 4

  • Results

    O(h^{k-1}) energy-norm and O(h^{k+1}) L2 error estimates are proved, and experiments show locking-free behavior on convex and non-convex domains for penalties up to 10^5.

    Source: §1 · p. 3, §5 · p. 22

  • Takeaways & Limitations

    As the penalties grow, the discrete solution converges to a nonconforming method on the optimal space B_k, supporting locking-free accuracy uniformly in penalty magnitude.

    Source: §5 · p. 22

    The work treats the planar problem; extending the framework to three-dimensional polyhedral meshes remains nontrivial, and a full 3D stability analysis remains open.

    Source: §5 · p. 22

Abstract

from arXiv · show

This paper presents an optimal interior penalty discontinuous Galerkin (IPDG) scheme for the planar biharmonic equation using piecewise polynomials of degree $k=3$ or $4$. In standard IPDG methods, large penalty parameters force the discrete solution into an overconstrained space, severely degrading accuracy---phenomenon known as numerical locking. To overcome this, our method enforces only vertex continuity and projects the jumps of the function and its normal derivative onto $\mathcal{P}^{k-3}$ and $\mathcal{P}^{k-2}$, respectively. Consequently, as the penalty parameters tend to infinity, the discrete solution is forced to lie in a constrained subspace $V_{h,\infty}^k$, which we identify as the optimal nonconforming finite element space $B_h^k$. This intrinsic connection fundamentally eliminates numerical locking. We prove optimal error estimates of $\mathcal{O}(h^{k-1})$ in a mesh-dependent energy norm and $\mathcal{O}(h^{k+1})$ in the $L^2$ norm. Numerical experiments on both convex and L-shaped domains confirm that the proposed scheme is robust and entirely locking-free, even for extremely large penalty parameters.

1 Introduction

Standard IPDG methods can lock when large penalties force low-degree approximations into overconstrained spaces, while conforming alternatives for cubic and quartic elements are not optimal on general grids. This paper introduces a reduced-penalty IPDG method whose limit is an optimal nonconforming space and proves optimal error estimates.

1 Introduction

  • Degree-k conforming discretizations of fourth-order problems have an energy-norm rate bounded by O(h^{k-1}), making lower-degree elements practically useful for simpler element structures and reduced computational overhead.Source: §1 · p. 2
  • Large penalties can cause numerical locking by forcing standard IPDG solutions into overconstrained spaces, and optimal C1-conforming spaces are rarely available for cubic or quartic polynomials on unstructured grids.The paper focuses on k = 3, 4; for k ≥ 5, the Argyris element already provides an optimal C1-conforming space.Source: §1 · p. 2
  • The method uses vertex continuity and penalizes projected function and normal-derivative jumps in P^{k-3} and P^{k-2}, respectively.As the penalties tend to infinity, the discrete solution converges to a nonconforming method on the optimal space B_k.Source: §1 · p. 3
  • The constrained limit space has the optimal O(h^{k-1}) energy-norm approximation property, supporting accuracy uniformly in the penalty magnitude for cubic and quartic elements.The paper identifies this property as eliminating numerical locking.Source: §1 · p. 3
  • O(h^{k-1}) energy-norm and O(h^{k+1}) L2 error estimates are proved, with locking-free behavior on convex and non-convex domains for penalties up to 10^5.The estimates apply for k = 3, 4 and solutions in H^{k+1}(Ω).Source: §1 · p. 3

2 IPDG schemes for biharmonic equation

The method addresses the clamped planar biharmonic equation on polygonal domains with a partially continuous discontinuous Galerkin space and reduced-order projected penalties. Its stability analysis establishes control equivalent to the full-jump DG norm and coercivity for sufficiently large penalties.

2 IPDG schemes for biharmonic equation

  • The model problem is the planar biharmonic equation with clamped boundary conditions on a polygonal domain and an L2 source term.Source: §2 · p. 3

2.1 The Unified Discrete IPDG Formulation

  • Trial and test functions are piecewise degree-k polynomials, k ∈ {3, 4}, constrained to be continuous at mesh vertices.At boundary vertices, this constraint enforces zero function values.Source: §2.1 · p. 4
  • The bilinear form uses reduced-order projected penalties to avoid artificial stiffness, with penalty parameters controlling the strength of continuity constraints.Its unified form is assembled from four components.Source: §2.1 · p. 4

2.2 Discrete Stability Analysis

  • The projected mesh-dependent energy norm controls discrete functions equivalently to the standard full-jump DG norm, with equivalence constants independent of mesh size.The proof uses control of unpenalized higher-frequency jump components by tangential derivatives.Source: §2.2 · pp. 5–6
  • For sufficiently large penalty parameters, the bilinear form is coercive, and the penalty thresholds remain uniformly bounded as the mesh is refined.As both penalties tend to infinity, the coercivity constant tends to a fixed value independent of the penalties.Source: §2.2 · pp. 6–7

3 Uniformly robust convergence of the scheme

This section develops the discrete Stokes complex and quasi-interpolation results used to establish convergence and penalty-robust error estimates. Under the stated regularity and coercivity conditions, it gives optimal O(h^{k+1}) L² and O(h^{k−1}) energy convergence, while describing the method’s approach to its locking-free limit as penalties increase and the instability possible with insufficient penalties.

3 Uniformly robust convergence of the scheme

Characterization of the Constrained Subspace and the Discrete Stokes Complex

Interpolation Operator

  • For k = 3 and 4, the discrete Stokes complex is exact without supplementary macro-element assumptions; quartic exactness uses the Crouzeix-Falk pair and its unconditional stability proved by Carstensen and Sauter.The cubic case uses the Fortin-Soulie pair; the quartic case extends the interpolation mechanism to k = 4.Source: Interpolation Operator · pp. 8–9

Optimal Quasi-Interpolation Error Estimates

3.2 Penalty-robust A Priori Error Estimates

  • For k = 3, 4 and sufficiently large penalties, the energy-norm error estimate holds under u ∈ H^(k+1)(Ω); its continuity bound is independent of γ1 and γ2.The interpolation error has zero projected jumps, making its penalty contribution vanish.Source: §3.2 · pp. 11–12
  • Under full H4 elliptic regularity for the dual biharmonic problem and sufficiently large penalties, the method has an optimal L2 error estimate.The stated regularity holds on smooth domains but generally fails on polygons.Source: §3.2 · p. 12
  • O(h^k+1) L2 convergence and O(h^k−1) energy convergence are optimal, with the energy rate uniform for fixed penalties above the coercivity thresholds.For large penalties, the L2 error is dominated by best approximation and the penalty contribution becomes negligible.Source: §3.2 · p. 14
  • On convex polygonal domains, dual regularity limits the rigorous L2 rate to O(h^k+0.739), although square-domain experiments show O(h^k+1) superconvergence.For a convex square, the leading corner singularity behaves as r^λ with Re(λ) ≈ 2.739, restricting the dual solution to H^3.739−ε(Ω); the reported superconvergence on structured grids does not contradict locking-free behavior.Source: §3.2 · p. 15

3.3 Evolution of the Penalty Parameter and Locking-Free Mechanism

  • For penalties below the coercivity thresholds, the discrete system may become unstable; experiments show oscillations for λ < 10^2.Source: §3.3 · p. 15
  • Just above the coercivity thresholds, the error-bound constant can grow large, but the O(h^k−1) mesh-refinement rate remains unaffected.The degenerate quasi-optimal constant is independent of mesh size h.Source: §3.3 · p. 15
  • As penalties grow, coercivity approaches c_eq/2 while continuity remains penalty-independent, making the error penalty contribution negligible and the method locking-free.The dimensionless stabilization parameters need not scale with mesh size; the O(γ^−1/2) bound describes convergence toward the limit solution on a fixed mesh as γ1, γ2 increase.Source: §3.3 · p. 15

4 Numerical Experiments

The experiments test convergence on regular meshes, behavior on a non-convex domain, and sensitivity to penalty parameters. They use a convex square and a non-convex L-shaped domain, with uniform refinements for convergence tests.

4 Numerical Experiments

  • The numerical experiments assess optimal convergence on regular meshes, robustness on non-convex domains, and sensitivity to penalty parameters.Source: §4 · p. 15
  • Tests use a convex square and a non-convex L-shaped domain, with successive uniform refinements for convergence analysis.Source: §4 · p. 15

4.1 Experimental Setup and Error Measures

The experiments use MATLAB and a direct solver for the sparse systems. Accuracy is measured in discrete L2, broken H1, and broken H2 norms; the penalty parameters are unified and the meshes are uniformly refined.

4.1 Experimental Setup and Error Measures

  • The implementation uses MATLAB and solves the discretization’s large sparse linear systems with a direct solver.Source: §4.1 · p. 15
  • Accuracy is measured using the discrete L2 norm and broken H1 and H2 semi-norms.Source: §4.1 · p. 16
  • The convergence order is computed by the standard logarithmic ratio of errors on adjacent meshes.Source: §4.1 · p. 16
  • The implementation sets γ1 = γ2 = λ and fixes λ = 1000 unless otherwise specified, exceeding the stated coercivity threshold.Source: §4.1 · p. 16
  • The coarsest mesh has h = 1 for both domains, and each refinement subdivides every triangle into four congruent subtriangles.Source: §4.1 · p. 16

4.2 Accuracy and Optimal Convergence on the Square Domain

Square-domain tests examine polynomial and oscillatory exact solutions. The P3 and P4 schemes achieve the reported optimal convergence rates for both tests, including the high-frequency trigonometric case.

4.2 Accuracy and Optimal Convergence on the Square Domain

Test 1: Smooth Polynomial Solution

Test 2: Smooth Trigonometric Solution

4.3 Performance on Non-Convex L-Shaped Meshes

On the non-convex L-shaped domain, the method retains near-optimal convergence for a smooth polynomial solution, while a singular test approaches the rates allowed by limited regularity.

4.3 Performance on Non-Convex L-Shaped Meshes

Test 1: Smooth Polynomial Solution

Test 2: Singular Solution

  • The singular solution’s reentrant-corner regularity limits theoretical rates to O(h^α) in H2, O(h^(1+α)) in H1, and roughly O(h^(2α)) in L2.Here α = 0.544483736782464.Source: Test 2: Singular Solution · p. 20
  • At the finest mesh levels for the singular test, H2 error convergence stabilizes near 0.55 and L2 convergence approaches 1.09–1.23.These rates approach limits attributed to the solution’s low regularity and the domain’s compromised dual regularity.Source: Test 2: Singular Solution · p. 20

4.4 Sensitivity Analysis of the Penalty Parameters

The penalty sensitivity analysis finds a stable regime at medium-to-large penalties, while penalties below the coercivity threshold produce fluctuations and spurious modes.

4.4 Sensitivity Analysis of the Penalty Parameters

  • The P4 sensitivity study fixes h = 0.0625 and scans λ from 1 to 100000; convergence rates are computed between successive uniformly refined meshes.The P3 sensitivity analysis is provided in the appendix.Source: §4.4 · p. 21
  • For λ below the coercivity threshold, weak interfacial constraints cause error and rate fluctuations and can produce non-physical spurious modes.The passage identifies strict positive-definiteness only when stabilization parameters exceed a critical lower bound.Source: §4.4 · p. 21
  • For λ ≥ 10^3, errors and convergence rates remain stable as the penalty increases, including on the L-shaped domain.The analysis compares the square and L-shaped domains and reports stability up to λ = 10^5.Source: §4.4 · pp. 21–22
  • As λ grows, penalization drives projected jumps toward zero and constrains the solution to an O(λ^-1/2) neighborhood of the constrained subspace.The approximation error remains dominated by the best approximation error in the optimal nonconforming space.Source: §4.4 · p. 21

5 Conclusion

The conclusion links the reduced-penalty IPDG framework to an optimal nonconforming space, reporting locking-free behavior for cubic and quartic elements and optimal convergence rates.

5 Conclusion

  • As the penalty grows, the discrete solution converges to the nonconforming method on B_h^k, separating approximation capability from penalty magnitude and eliminating locking for k = 3, 4.The framework combines a partially continuous finite element space with degree-matched reduced-order projections.Source: §5 · p. 22
  • For k = 3, 4, the analysis proves O(h^(k−1)) energy-norm and O(h^(k+1)) L2 convergence; experiments confirm stability up to λ = 10^5 on convex and non-convex geometries.The conclusion contrasts this pure-IPDG approach with weakly over-penalized methods and mixed formulations.Source: §5 · p. 22
  • Extending the framework to three-dimensional polyhedral meshes remains open, with face-based jumps, discrete Stokes complex construction, and full stability analysis posing challenges.The construction is especially delicate for quartic elements on general tetrahedral partitions.Source: §5 · p. 22

A Supplementary Test: Smooth Trigonometric Solution on the L-Shaped Domain

On the non-convex L-shaped domain, a supplementary trigonometric test at fixed penalty shows the expected energy-norm convergence rates for P3 and P4.

A Supplementary Test: Smooth Trigonometric Solution on the L-Shaped Domain

  • For the trigonometric exact-solution test on the non-convex L-shaped domain, with fixed penalty parameter λ = 1000, P3 and P4 energy-norm errors decrease as h is halved and converge approximately at O(h2) and O(h3), respectively.Source: A Supplementary Test: Smooth… · p. 26
  • With λ = 1000, the trigonometric test’s energy errors converge at approximately O(h^2) for P3 and O(h^3) for P4.The reported results cover both polynomial degrees on the L-shaped domain.Source: A Supplementary Test: Smooth… · p. 26, Table 9 · p. 26, Table 10 · p. 26
  • Tables 9 and 10 report the P3 and P4 computational results, respectively, at fixed λ = 1000.The tables present numerical errors and convergence rates.Source: Table 9 · p. 26, Table 10 · p. 26

B Sensitivity Analysis for the P3 Element

The P3 sensitivity analysis examines numerical errors and convergence rates as the penalty parameter varies on square and L-shaped domains. The L-shaped domain shows the same stable behavior for large penalty parameters as the square domain.

B Sensitivity Analysis for the P3 Element

Loading 2609.12675v1…