Finite element approximation for the Bessel $(p,s)$-Laplacian

Juan Pablo Borthagaray; José Camilo Rueda Niño

Source-linked AI summary

Finite element approximation for the Bessel $(p,s)$-Laplacian

Juan Pablo Borthagaray, José Camilo Rueda Niño

arXiv:2609.37747v1Sep 29, 2026math.NAcs.NA

TL;DR

The paper addresses the computational cost of direct finite element discretization for the Bessel (p, s)-Laplacian Dirichlet problem. It develops an augmented Lagrangian formulation using a projected discrete gradient and analyzes its consistency and convergence. The convergence analysis depends on a discrete inf-sup condition whose uniform validity remains under investigation.

Problem
Direct finite element implementation is expensive because the stiffness matrix is dense and must be reassembled at every nonlinear iteration.
Source: §1 · p. 2
Method
The paper develops an augmented Lagrangian formulation based on a projected discrete problem that replaces Riesz fractional gradients with elementwise averages.
Source: §1 · pp. 2–3
Results
The analysis estimates the consistency residual, derives a Strang-type estimate, and establishes convergence rates; numerical experiments illustrate performance across nonlinear regimes and geometries.
Source: §1 · p. 3
Takeaways & Limitations
For f in the dual of the zero-exterior Lions–Calderón space, the Dirichlet problem has a unique minimizer and weak solution.
Source: §2.2 · p. 7
Uniform validity of the discrete inf-sup condition remains under investigation; explicit validity is shown for p = 2 on uniform one-dimensional meshes.
Source: §4.1 · p. 13

Abstract

from arXiv · show

We study the finite element approximation of the Dirichlet problem for the Bessel $(p,s)$-Laplacian, $\operatorname{div}^{s} \left(|\nabla^{s} u|^{p-2}\nabla^{s} u\right)$, built from the Riesz fractional gradient $\nabla^s$ and posed on the Bessel potential space $X^{s,p}$, obtained by complex interpolation between $L^p$ and $W^{1,p}$. For $p = 2$ this operator reduces to the fractional Laplacian, while for $p \neq 2$ it differs from the fractional $p$-Laplacian arising from real interpolation. Since a direct Galerkin discretization requires reassembling a dense stiffness matrix at every nonlinear iteration, we propose an augmented Lagrangian formulation based on the projected fractional gradient $B_h = \mathbf{\Pi}_h \nabla^s$, in which the only dense matrix is assembled once and the nonlinearity decouples elementwise. We analyze the resulting consistency error, establish a priori convergence rates in $X^{s,p}$, conditional on a discrete inf-sup condition for the augmented Lagrangian scheme, and present numerical experiments that illustrate the convergence of the method and its ability to deal with degenerate and strongly nonlinear problems.

1. Introduction

The introduction frames the Bessel (p, s)-Laplacian Dirichlet problem within the broad study of nonlocal operators. It motivates an augmented Lagrangian finite element approach by the cost of direct discretization and outlines its analysis and numerical evaluation.

1. Introduction

Fractional-order operators model nonlocal phenomena in applications including phase transitions, neural networks, biology, elasticity, transport, and image processing.The introduction focuses on the Riesz fractional gradient, which is connected with the fractional Laplacian.Source: §1 · p. 1
Direct Galerkin implementation is costly because its stiffness matrix is dense and must be reassembled at every nonlinear iteration.Each matrix entry also involves an integral over Rd.Source: §1 · pp. 2–3
The proposed augmented Lagrangian approach uses a projected discrete problem that replaces Riesz fractional gradients with elementwise averages.The paper also develops an iterative solver and extends the approach to variable diffusivity.Source: §1 · p. 3
The analysis estimates the projected formulation’s consistency residual, derives a Strang-type estimate, and establishes convergence rates.Numerical experiments examine different nonlinear regimes and geometrical configurations.Source: §1 · p. 3

2. The Bessel (p, s)-Laplacian: Variational Setting

This section establishes the variational setting for the Bessel (p, s)-Laplacian and formulates its Dirichlet problem as a convex minimization problem. It distinguishes the Lions–Calderón framework from fractional Sobolev spaces and identifies regularity estimates used in the later error analysis.

2.1. Fractional calculus operators and Lions–Calderón spaces.

The Riesz fractional gradient and its associated divergence provide the operators for the variational framework, with the zero-exterior space norm given by the Lp norm of the gradient.The section also relates the Bessel operator to complex interpolation and contrasts it with the fractional p-Laplacian arising from real interpolation.Source: §2.1 · pp. 3–4, §2.2 · p. 7
The Lions–Calderón spaces are identified with Bessel potential spaces, while fractional Sobolev spaces arise through real interpolation and are strictly distinct except when p = 2.The paper uses the Lions–Calderón space as the variational setting associated with the Riesz fractional gradient.Source: §2.1 · pp. 4–5

2.2. Dirichlet problem.

For data f in the dual of the zero-exterior Lions–Calderón space, strict convexity and weak lower semicontinuity yield a unique minimizer and weak solution.The energy functional combines the fractional-gradient p-energy with the duality pairing against f.Source: §2.2 · p. 7
The section develops continuity and ellipticity-type estimates for the nonlinear operator to obtain quantitative stability estimates for the solution.The solution regularity estimates used for a priori error analysis are expressed in the Besov scale and are drawn from [13, Theorems 4.1 and 4.2].Source: §2.2 · pp. 7–8

3. Finite element approximation

The section develops a finite-element approximation of the variational problem and analyzes its error in the nonlocal setting. It establishes discrete well-posedness and convergence estimates, with the rates depending on solution regularity and deteriorating near the endpoints of the range of p.

3.1. Setting.

The Galerkin problem uses linear Lagrange elements on a shape-regular simplicial mesh, with its finite-dimensional space contained in the variational space.The section assumes Ω is a polytope and defines the mesh with maximum element diameter h.Source: §3.1 · pp. 8–9
The discrete problem has a unique solution because it minimizes the restriction of the convex energy, and it satisfies a discrete stability estimate.The section sets aside practical implementation when analyzing approximation quality.Source: §3.1 · p. 9

3.2. Interpolation and error bounds.

Because the Riesz fractional gradient is nonlocal, elementwise error contributions cannot be treated independently; the analysis instead uses localization and quasi-interpolation estimates.The interpolation framework draws on quantitative contiguity and localized fractional Sobolev seminorms.Source: §3.2 · pp. 9–10
A priori convergence rates are established for the finite-element solution, with constants independent of h under the stated assumptions.The corollary’s hidden constants may depend on f, d, s, p, Ω, and mesh regularity; suitable regularity of the right-hand side yields rates.Source: §3.2 · p. 10
The error exponents degenerate as p approaches 1 or infinity; p-robust rates in quasi-norms are left for future work.The cited remark contrasts this open analysis with robust quasi-norm rates in the local setting.Source: §3.2 · p. 11

4. Augmented Lagrangian formulation

The section presents a projected augmented Lagrangian scheme that replaces direct nonlinear discretization, with an exact discrete constraint and a reusable nonlocal matrix, and discusses its conditional well-posedness, iterative convergence, and extension to variable diffusivity. It also reports numerical experiments for the p = 2 case with A = χΩ, outside the stated hypotheses, indicating boundary behavior different from the fractional Laplacian.

4. Augmented Lagrangian formulation

Direct assembly would require dense stiffness matrices to be recomputed at each nonlinear iteration, motivating a constrained augmented Lagrangian reformulation.Introducing w = ∇^s u makes the whole-space constraint explicit; saddle points recover the original solution and identify the multiplier as |∇^s u|^(p−2)∇^s u.Source: §4 · pp. 11–12

4.1. An approximate energy

Replacing the continuous energy with the projected discrete energy introduces a nonlocal consistency error; for s = 1, the projection agrees with the gradient on V_h.The discrete constraint uses B_h = Π_h∇^s, so w_h = B_hu_h.Source: §4.1 · p. 12
For p > 1 and f ∈ X−s,p′(Ω), the discrete functional Jh: Vh → R defined in (4.3) admits a unique minimizer uh ∈ Vh.The operator Bh is injective; whether the condition is satisfied uniformly remains under research, although explicit calculations show validity of (4.6) for p = 2 over uniform one-dimensional meshes.Source: §4.1 · pp. 12–13

4.2. Discrete augmented Lagrangian formulation

Critical points of the discrete augmented Lagrangian satisfy the projected problem with w_h = B_hu_h and λ_h = |B_hu_h|^(p−2)B_hu_h.The unprojected alternative is an exact reformulation of the Galerkin problem only when p = 2 and r = 1; convergence with its consistency error is unclear.Source: §4.2 · pp. 14–15
Only B requires nonlocal-integral assembly, and the nonlinearity reduces to independent elementwise algebraic equations without matrix recomputation during iteration.The associated bilinear-form matrix can be expressed using B.Source: §4.2 · p. 16

4.3. Iterative solution of the discrete system

Algorithm 4.9 is an Uzawa-type iteration whose iterates converge to the unique discrete saddle point.Injectivity of B_h ensures coercivity of the primal subproblem; the other two iteration steps are computed elementwise.Source: §4.3 · pp. 16–17

4.4. Extension to variable diffusivity

The scheme accommodates space-dependent diffusivity with little modification because the Bessel energy has a pointwise coefficient, unlike the two-point coefficient required for the fractional p-Laplacian.The variable coefficient enters through a diagonal matrix, leaving B independent of the coefficient and reusable across diffusivities.Source: §4.4 · p. 17
For degenerate diffusivity, a different weighted fractional Sobolev framework is needed; general coefficients also require projection and incur an unpursued consistency term.Explicit error bounds are unavailable for piecewise constant diffusivity because the cited regularity argument is not expected to carry over.Source: §4.4 · pp. 17–18

4.4. Extension to variable diffusivity.

For p = 2, Section 6 numerically explores the variable-diffusivity configuration A = χΩ.Source: §4.4 · p. 18
In this configuration, a0 = 0, so it falls outside the hypotheses above.Source: §4.4 · p. 18
The configuration may resemble the Dirichlet problem for the censored fractional Laplacian.Source: §4.4 · p. 18
For the censored fractional Laplacian, the energy stems from the Sobolev seminorm |·|Hs(Ω).Source: §4.4 · p. 18
The operator considered here is different because its energy involves the L2(Ω)-norm of the s-gradient.Source: §4.4 · p. 18
Numerical experiments indicate that its solution boundary behavior differs from that of the fractional Laplacian.Source: §4.4 · p. 18

5. Convergence of the discrete augmented Lagrangian formulation

The analysis bounds the consistency error from projecting the fractional gradient and truncating the computational domain, then combines it with approximation error to derive convergence under a uniform stability condition. The projected approximation preserves the Galerkin rate for p = 2, while its bounds are less sharp for p ≠ 2.

5. Convergence of the discrete augmented Lagrangian formulation

The projection defect satisfies δ_h(v) ≤ C h^(σ-s) ||v||_{B^σ_{p,∞}(Ω)} on shape-regular, quasi-uniform meshes.For functions smoother than the variational space, the defect tends to zero at a rate determined by their regularity.Source: §5 · p. 19

5.1. Consistency of the projected operator.

The consistency residual separates contributions from projection onto L_h and truncation outside the computational domain.The truncation contribution must be balanced against discretization error by choosing the growth of R_ext with h.Source: §5.1 · pp. 20–21, §5.2 · p. 23
At the interpolant π_hu, the approximation error is O(h^β |log h|) and the solution’s projection defect is O(h^β).These estimates support the consistency bound for sufficiently small h under Proposition 5.4’s domain and mesh assumptions.Source: §5.1 · p. 21

5.2. A Strang-type estimate and convergence rates.

A Strang-type argument decomposes the solution error into approximation and consistency terms, whose bounds yield error rates.The analysis applies the consistency and approximation estimates at π_hu.Source: §5.2 · pp. 21, 23
The convergence-rate corollary assumes uniform stability of the discrete scheme, in addition to the stated solution, mesh, and domain hypotheses.The stability assumption is required for the derived rate result.Source: §5.2 · pp. 22–23
For p = 2, the projected approximation preserves the Galerkin convergence rate; for p ≠ 2, its error bounds are less sharp.The authors attribute the reduced rate to their continuity-and-monotonicity-based Strang argument; the Galerkin rate exploits minimization of the same energy at continuous and discrete levels.Source: §5.2 · p. 23
For 1 < p ≤ 2, an alternative monotonicity-and-continuity argument gives a sharper Galerkin bound by using Galerkin orthogonality.The paper contrasts that argument with the projected scheme’s Strang estimate, which precludes the same absorption step.Source: §5.2 · pp. 23–24

6. Numerical Experiments

The experiments validate the linear case, examine how nonlinearity and fractional order affect solutions, and compare a localized variable-diffusivity operator with the fractional Laplacian.

6.1. Validation in the linear case.

Observed orders of 0.505, 0.518, and 0.560 for s = 0.25, 0.5, and 0.75 match the expected 1/2 rate.The figure reports energy-norm errors for p = 2.Source: §6.1 · p. 24, Figure 6.1 · p. 25
For s = 0.5, projected orders were unchanged across r = 1, 2, and 5, while unprojected orders were 0.526, 0.423, and 0.281, respectively.The projected scheme's observed order was independent of r; the unprojected scheme showed reduced rates for r ≠ 1.Source: §6.1 · p. 24, Figure 6.1 · p. 25

6.2. Effect of Nonlinearity. (

At p = 10, solutions were markedly asymmetric at s = 0.1 and became progressively more symmetric as s increased.The authors attribute the asymmetry to the nonlocal exterior weighting, which pulls the solution toward zero less strongly near the reentrant corner.Source: §6.2 · p. 25, Figure 6.2 · p. 26

6.3. Variable Diffusivity. (

Suppressing the fractional-gradient contribution outside Ω changes the solution particularly near the boundary.The comparison uses A = χΩ versus A ≡ 1 on the same mesh configuration.Source: §6.3 · pp. 25–26
For s = 0.5, fitted boundary exponents were γ ≈ 0.48 (A ≡ 1) and γ ≈ 0.34 (A = χΩ), with localized-operator behavior differing from the censored fractional Laplacian.The censored fractional Laplacian is described as defined only for s > 1/2, with boundary behavior u(x) ∼ dist(x, ∂Ω)^(2s−1).Source: §6.3 · p. 26

Appendix A. On the uniform bound for Bh

Appendix A studies whether the constant Ch can be bounded uniformly in h and identifies a tractable special case.

The appendix's general question is whether the constant Ch from Remark 4.4 remains bounded independently of h.Source: Appendix A · p. 29
Uniform boundedness of Ch reduces to controlling the relative Lp size of the projection defect (I − Πh)∇s vh.The defect separates into a projection part and an exterior tail, with the projection part identified as the difficult component.Source: Appendix A · p. 29
For uniform meshes in d = 1 with p = 2, the authors prove Ch is uniformly bounded, with the bound tending to β*^-1 = 1.0805 as Ωext grows.Source: Appendix A · p. 29

A.1. The approximation question.

The uniform-bound question becomes a strict relative-approximation condition on the projected fractional gradient, while scale invariance limits the usefulness of ordinary h-order estimates.

The exterior tail can be controlled uniformly in h by taking the computational domain sufficiently large.Source: §A.1 · p. 29
Uniform control requires the Lp projection defect to be bounded by λ times the Xs,p norm, with λ < 1 independent of h.This is the condition stated in (A.4).Source: §A.1 · p. 29
No estimate with a positive power of h establishes the required bound: the gain of 1 − s derivatives is canceled by mesh-scale roughness, yielding only O(1).Source: §A.1 · p. 30
For self-similar mesh families, dilation invariance makes uniformity in h automatic, leaving a dimensionless quantity determined by mesh pattern, s, and p.Source: §A.1 · p. 30

A.2. The case p = 2 and uniform meshes in d = 1.

For p = 2 on uniform one-dimensional meshes, Fourier analysis characterizes the projection defect condition underlying a uniform bound for Ch.

For p = 2, orthogonality and Fourier analysis turn the sufficient projection-defect condition into a characterization and sharpen the constant.Source: §A.2 · p. 30
A uniform bound for Ch holds if and only if the projection defect satisfies (A.4) with λ < 1 independent of h.Source: §A.2 · p. 30
R(s,t) extends to a continuous, strictly positive function, so βs = inf R(s,t) is positive for every s in [0,1].Source: §A.2 · p. 32
The Fourier-based estimate transfers to bounded domains, but does not account for truncation of the mesh at Ωext, which leaves an exterior-tail contribution.Source: §A.2 · p. 32
The proof specializes to d = 1, p = 2, and uniform meshes, using translation invariance and the representation of discrete functions by hat functions.Source: §A.2 · pp. 30, 32