PaperlayerView exact revision on arXiv

Source-linked AI summary

A Matrix-free Augmented High Order Compact Solver for Variable-Coefficient Biharmonic Problems

Jin Li, Kejia Pan, Xu Qian, Li-Lian Wang

arXiv:2609.16478v1math.NA

TL;DR

The paper addresses the global coupling caused by unavailable auxiliary-variable boundary values in mixed discretizations of large-scale variable-coefficient biharmonic problems. It augments those boundary values as unknowns, producing decoupled second-order subproblems linked through a lower-dimensional Schur complement solved matrix-free. The resulting framework retains high-order accuracy and is reported to solve problems with more than one billion unknowns within practical runtimes, while extending to related high-order elliptic applications.

  • Problem

    Unavailable auxiliary-variable boundary values leave mixed biharmonic discretizations coupled, limiting high-order, compact, and efficient solution methods.

  • Method

    The method augments unknown auxiliary-variable boundary values, decouples second-order subproblems, solves them with FFT-based methods, and handles the Schur complement using matrix-free GMRES.

  • Results

    Numerical experiments confirm accuracy, robustness, and efficiency, with more than one billion unknowns solved within a practical runtime and applications demonstrated for related high-order elliptic problems.

  • Takeaways & Limitations

    The augmented framework supports efficient biharmonic computation and applies to triharmonic equations, incompressible flows, and plate bending.

  • Takeaways & Limitations

    The reported complexity is a conditional cost estimate rather than a uniform complexity bound under grid refinement.

Abstract

from arXiv · show

We propose an augmented high-order compact finite difference method for biharmonic equations with clamped boundary conditions and variable coefficients. Standard mixed-type formulations introduce an auxiliary variable, but its boundary values are unavailable, leaving the resulting discrete systems globally coupled and difficult to solve at large scales. Our key contribution is the development of a new augmented formulation that treats these unavailable boundary values as additional unknowns, reduces the global coupling to a lower-dimensional Schur complement system, and yields decoupled second-order subproblems. The Schur complement is solved by matrix-free GMRES, while the subproblems are handled by FFT-based fast solvers. The method achieves fourth-order accuracy using compact stencils, and has $O(n\log n)$ computational complexity, enabling the solution of the biharmonic equation with $1024^3$ degrees of freedom within several minutes. To the best of our knowledge, this level of computational efficiency has not previously been achieved in either the literature or practice. Using energy estimates and Fourier analysis, we derive a new $L^2$-estimate for Poisson equations with inexact Dirichlet boundary and then prove the convergence of the proposed scheme. We provide ample numerical experiments to confirm the accuracy, efficiency, and further apply the fast and accurate solver to triharmonic equations, high-wavenumber problems, Stokes flow, and plate bending problems.

1 Introduction

Biharmonic problems are difficult to solve with methods that are simultaneously high-order, compact, and efficient, especially under clamped boundary conditions and at large scales. The paper addresses this through an augmented mixed formulation that enables decoupled second-order solves and matrix-free computation.

  • 1 Introduction: Mixed reformulations reduce biharmonic equations to second-order systems, but unavailable auxiliary-variable boundary values typically leave the discrete system globally coupled.This coupling limits the simultaneous achievement of high-order accuracy, compactness, and computational efficiency.
  • 1 Introduction: The paper introduces an augmented HOC solver for variable-coefficient biharmonic equations with clamped boundary conditions.
  • 1 Introduction: Treating auxiliary-variable boundary values as additional unknowns confines the global coupling to a lower-dimensional Schur complement system.The augmented variables reside on the boundary, which has one fewer spatial dimension than the primary unknown.
  • 1 Introduction: After the augmented variables are computed, the original mixed system becomes a set of decoupled second-order subproblems.This exposes Poisson- or Helmholtz-like subproblems to fast solver technology.
  • 1 Introduction: The matrix-free implementation evaluates Schur-complement matrix-vector products through FFT-based subproblem solves, yielding overall O(n log n) complexity.The method avoids explicitly assembling the Schur complement and uses fast iterative or direct solvers for the subproblems.
  • 1 Introduction: The convergence analysis derives an L2 estimate for discrete Poisson problems with inexact Dirichlet data and uses it to support the scheme’s convergence proof.The estimate is established using discrete sine expansions.

2 The augmented HOC solver

The augmented HOC solver introduces unavailable auxiliary boundary values as lower-dimensional unknowns, reducing global coupling to a Schur complement while retaining compact fourth-order discretization. Matrix-free GMRES and fast subproblem solvers yield conditional O(n log n) complexity.

  • 2.1 Augmented method: The biharmonic equation is split into Poisson and Helmholtz-type second-order equations for u and v, coupled through the augmented boundary constraint.The auxiliary variable is defined by v = α∆u, and the resulting equations are discretized separately.
  • 2.1 Augmented method: The method augments unavailable boundary values of the auxiliary variable, producing an equivalent formulation with three coupled subsystems and a lower-dimensional boundary unknown.The augmented variable represents the otherwise unavailable trace of v on the clamped boundary, rather than an additional modeling quantity.
  • 2.2 HOC finite difference discretization: The scheme uses fourth-order compact finite-difference operators and high-order one-sided boundary approximations for the normal derivatives entering the augmented constraint.Tangential derivatives are obtained from the Dirichlet data, while the normal second derivative is approximated with a one-sided formula using the Neumann condition.
  • 2.3 Efficient implementation of the augmented HOC discretization: Matrix-free Schur-complement products are evaluated through decoupled subproblem solves, making GMRES suitable without explicitly assembling the Schur complement.Explicit assembly becomes expensive under refinement; the matrix-free action instead invokes the discrete subproblem solves.
  • 2.3 Efficient implementation of the augmented HOC discretization: O(n log n) computational cost is obtained when inner and outer iteration counts remain bounded under refinement, but this complexity estimate is conditional.For N^3 grids, the augmented boundary has O(N^2) unknowns, while constant-coefficient subproblems use discrete sine transforms.

1 Main procedure

The main procedure first solves the augmented Schur-complement system by matrix-free GMRES, then recovers v and u through FFT-based solves.

  • 1 Main procedure: Matrix-free GMRES solves the Schur-complement system for the converged augmented variable Qm.Each matrix-vector product is evaluated through the SCHURMATVEC procedure.
  • 1 Main procedure: The method computes V from Qm and then computes U from V using FFT-based solvers for the two discrete second-order equations.The output is the fourth-order numerical solution pair V and U.

6 End

The Schur right-hand side is constructed by solving the v- and u-equations with zero augmented data and evaluating the resulting boundary residual.

  • 6 End: The Schur right-hand side b is formed by computing V(0), then U(0), and evaluating the residual R(0).The construction uses FFT-based solves for both zero-augmented-data subproblems.

14 End

The Schur matrix-vector product is obtained by solving both second-order subproblems for a given Q and combining the residual with the Schur right-hand side.

  • 14 End: SCHURMATVEC computes V(Q), then U(Q), evaluates R(Q), and returns R(Q) + b as the product ShQ.Because b = −R(0), this residual combination yields the Schur-complement action.

3 Convergence analysis

The convergence analysis establishes conditional error bounds for the augmented HOC scheme in two and three dimensions using discrete stability, elliptic regularity, Sobolev estimates, and sine expansions. The rigorous analysis gives O(h7/2) convergence, while numerical results indicate fourth-order behavior and motivate sharper treatment of boundary contributions.

  • Stability estimate: The analysis derives stability through discrete elliptic regularity, positive definiteness, and Sobolev embedding for grid functions vanishing on the boundary.These ingredients apply on uniform Cartesian grids in both two and three dimensions.
  • Stability estimate: The central technical estimate treats fourth-order compact problems with boundary data on one edge or face and provides a mesh-independent stability constant.The boundary forcing is prescribed on Γ_h and vanishes on the remaining boundary grid points.
  • Error analysis: Two-dimensional error analysis decomposes the auxiliary-variable error into boundary-driven and truncation-error components before applying the stability estimate to the primary-variable error.The same framework extends to three dimensions, where boundary forcing lies on a face and a tensor-product sine expansion is used in both transverse directions.
  • Convergence theorem: The convergence theorem assumes the augmented boundary-variable error satisfies an O(h7/2) discrete L2 bound, which numerical results support but the analysis does not yet prove.Under this assumption, the remaining convergence analysis is rigorous, although the resulting estimate is weaker than the observed fourth-order rate.
  • Three-dimensional extension: O(h7/2) convergence is rigorously established in three dimensions, but the bound does not recover the fourth-order rate observed numerically.The loss of one-half order is attributed to estimating the boundary-driven component without exploiting possible cancellation or additional structure.

4 Numerical experiments

Numerical experiments assess convergence, efficiency, robustness, and broader applicability of the augmented HOC scheme across variable-coefficient, oscillatory, singularly perturbed, triharmonic, flow, and plate problems.

  • Convergence tests: Fourth-order convergence is recovered on sufficiently fine grids for both ε = 10−4 and ε = 10−5 despite thinner boundary layers at smaller ε.The ε = 10−5 case has slightly lower coarse-grid rates because its boundary layer near x = 0 is thinner.
  • Convergence tests: Variable-coefficient Example 2 maintains near-fourth-order convergence, with ∥Eu∥∞ decreasing monotonically and GMRES counts between 4 and 10.The more strongly varying (q, p) = (10, 0.2) solution produces larger same-mesh errors, while asymptotic convergence remains essentially unchanged.
  • Convergence tests: Example 3 approaches fourth-order convergence under both frequency settings, although (k1, k2, k3) = (25, 5, 25) has larger errors and lower coarse-grid orders.The 256 × 256 × 256 solution captures anisotropic oscillations, stronger in the x- and z-directions than in y.
  • Computational efficiency: 362.76 seconds and 161.88 seconds suffice for 1024 × 1024 × 1024 meshes in Examples 2 and 3, respectively.GMRES iteration counts remain very small across meshes, indicating mild overhead from the augmented solve.
  • Triharmonic equation: The triharmonic extension approaches fourth-order convergence on fine meshes, including the more oscillatory setting where same-mesh errors are larger.Grid-refinement tests report approximately fourth-order maximum-norm error decay for sufficiently fine meshes across tested wave numbers.
  • Applications: The method agrees closely with literature benchmarks for lid-driven cavity flow and with classical central-deflection and bending-moment values for a clamped plate.The plate quantities converge to 0.005612 and −0.12577, compared with reference values 0.00560 and −0.1257.

5 Conclusions

The proposed augmented HOC scheme decouples biharmonic subproblems through boundary augmentation and combines FFT-based solvers with matrix-free GMRES. Experiments support its accuracy, efficiency, robustness, and applicability to related high-order elliptic problems.

  • Boundary augmentation reduces the biharmonic problem to decoupled second-order subproblems coupled through a lower-dimensional Schur complement, solved with FFT-based methods and matrix-free GMRES.The scheme applies in two and three dimensions with clamped boundary conditions, and its convergence is proved using discrete stability estimates and discrete sine expansions.
  • Numerical experiments confirm accuracy, robustness, and efficiency for constant- and variable-coefficient problems, with more than one billion unknowns solved within a practical runtime.Only a small number of GMRES iterations are required.
  • Applications to triharmonic equations, incompressible flows, and plate bending demonstrate the framework’s applicability to related high-order elliptic problems.

A Proofs of Auxiliary Lemmas and Corollaries

The auxiliary proofs analyze discrete sine expansions, orthogonality, and recurrence relations for the discrete operator. These tools establish explicit modal representations and bounds used in the convergence analysis.

  • Discrete sine expansion and orthogonality reduce the homogeneous discrete problem to independent modal recurrences and provide a Parseval identity.The reduction follows from substituting the expansion into the discrete operator and using linear independence of the sine basis.
  • Each modal coefficient satisfies a second-order recurrence with parameter relation cosh α_k = 5 − 2 cos θ_k.The recurrence is obtained after rewriting the modal equation in hyperbolic form.
  • Boundary conditions determine the recurrence solution and yield the stated modal representation.The representation follows by solving the recurrence subject to the boundary conditions.
  • The auxiliary estimates combine bounds for α_k regimes with discrete orthogonality to complete the lemma proofs.The argument treats α_k ≤ 1 and α_k > 1 separately before concluding the estimate.
Loading 2609.16478v1…