Source-linked AI summary
Numerical Methods for the Fractional Laplacian: a Finite Difference-quadrature Approach
Yanghong Huang, Adam Oberman
TL;DR
The paper addresses the limited accuracy and convergence theory of numerical methods for the fractional Laplacian, especially for the extended Dirichlet problem on bounded domains. It develops a finite difference/quadrature discretization with positive convolution weights, proves convergence, and reports O(h3−α) accuracy for quadratic interpolation, alongside improved far-field treatment and nonlinear obstacle-problem computations.
Problem
Accuracy and convergence were not firmly established for many numerical methods for the fractional Laplacian, whose bounded-domain formulations also have unsettled interpretations.
Method
The paper combines the singular-integral representation with finite differences and exact interpolation-based quadrature to construct a discrete convolution operator with positive weights.
Results
O(h3−α) accuracy is obtained with quadratic interpolation, while the scheme is proved convergent and far-field asymptotics improve truncation accuracy by an order of magnitude.
Takeaways & Limitations
The resulting consistent, stable positive-weight scheme supports computations for nonlinear equations involving the fractional Laplacian, including the obstacle problem.
Takeaways & Limitations
Higher accuracy for singular solutions remains an open problem, and the reported convergence-rate theory for weak solutions has an unknown exponent.
Abstract
from arXiv · showhide
The fractional Laplacian $(-Δ)^{α/2}$ is a non-local operator which depends on the parameter $α$ and recovers the usual Laplacian as $α\to 2$. A numerical method for the fractional Laplacian is proposed, based on the singular integral representation for the operator. The method combines finite difference with numerical quadrature, to obtain a discrete convolution operator with positive weights. The accuracy of the method is shown to be $O(h^{3-α})$. Convergence of the method is proven. The treatment of far field boundary conditions using an asymptotic approximation to the integral is used to obtain an accurate method. Numerical experiments on known exact solutions validate the predicted convergence rates. Computational examples include exponentially and algebraically decaying solution with varying regularity. The generalization to nonlinear equations involving the operator is discussed: the obstacle problem for the fractional Laplacian is computed.
1. Introduction.
The paper addresses numerical analysis challenges for the nonlocal fractional Laplacian on bounded domains by deriving and analyzing a finite difference/quadrature method for the extended Dirichlet problem.
- The fractional Laplacian models non-local diffusion with long-range interactions.
- Accuracy and convergence were previously unknown or insufficiently established for many numerical methods for the operator.
- The method targets a proved, accurate treatment of far-field truncation and supports later computation of nonlinear problems such as the obstacle problem.
- On bounded domains, different operator representations can produce different operators, while physical and probabilistic interpretations remain unsettled.
- The extended Dirichlet problem requires unknowns only inside D but prescribed values on the entire complement of D.
- The proposed one-dimensional discretization uses a discrete convolution of grid values with positive weights derived from the singular integral representation.
2. Quadrature and finite difference discretization of the fractional Laplacian.
The discretization splits the singular integral into a locally finite-difference part and a tail treated by exact weighted interpolation quadrature. The resulting schemes have error orders determined by the interpolation degree.
- Tail quadrature: The tail integral uses exact integration of the kernel multiplied by an interpolant of the regularized unknown function.Direct quadrature is avoided because its errors become too large as h approaches zero.
- Singular part: The singular part is approximated by a rescaled second derivative and then by centered finite differences on a uniform grid.
- Discrete weights: Linear and quadratic interpolation produce discrete weights wT and wQ, respectively, with symmetric, index-independent weights under the stated construction.
- Tail quadrature: The tail quadrature is built from piecewise Lagrange polynomials on subintervals outside the singular region.
- Error estimates: O(h2−α) accuracy is obtained with linear interpolation, while quadratic interpolation yields O(h3−α) accuracy for u ∈ C4.
- Error estimates: Higher-order interpolation can improve formal accuracy but may produce negative weights and sacrifice properties such as the maximum principle.
3. Explicit calculation of the weights.
The method derives explicit convolution weights by combining piecewise linear or quadratic interpolation with semi-exact quadrature against the fractional-Laplacian kernel. The resulting weights are positive, scale as h^-α, decay algebraically, and recover the central difference scheme as α→2−.
- Interpolation bases: Piecewise linear tent functions and piecewise quadratic Lagrange interpolants provide the interpolation bases used to derive explicit weights.The quadratic construction uses Qh and Rh, with local support extending over finite overlapping intervals.
- Weight construction: Primitives F and G of the kernel weight function convert the semi-exact quadrature formulas into explicit linear and quadratic weights.The formulas treat the singular interval and boundary-adjacent intervals using one-sided integrals.
- Weight properties: The weights from linear and quadratic interpolation are positive, supporting the maximum principle, stability, and convergence of the discrete extended Dirichlet problem.Higher-order interpolation is avoided because it can produce negative weights and lose these properties.
- Weight properties: The weights scale as h^-α and decay as j^-1-α, reflecting the fractional derivative order and the algebraic tail of the Lévy kernel.After removing h^-α, the weight sums remain uniformly bounded for α in (0, 2).
- Consistency: As α→2−, both interpolation-based discretizations recover the standard three-point central difference scheme for −∂xx.The construction therefore connects the fractional operator discretization to the classical second-derivative scheme.
- Time stepping: An explicit parabolic time integrator imposes a time-step restriction determined by the sum of the discrete weights.This links the weight normalization directly to the CFL condition.
4. Convergence of solutions to the Dirichlet boundary value problem.
The convergence proof combines positive discrete weights, a maximum principle, and a discrete super-solution to control the solution error by the truncation error. Consequently, the solution inherits the global accuracy O(h2−α) for wT and O(h3−α) for wQ under sufficient smoothness.
- Error control: A discontinuous quadratic super-solution is sampled on the grid to control errors inside the computational interval.Continuity is unnecessary for the stability results, although continuous barriers are used for existence proofs.
- Discrete problem: The discrete extended Dirichlet problem places the unknowns on Dh and prescribes exterior values gi on DC.The scheme is written as (−∆h)^α/2ui = fi in Dh with ui = gi in DC.
- Maximum principle: The discrete super-solution satisfies −δjvi ≥ min(2, 2(jh)^2) for i∈Dh and j>0, with separate cases handling interior and exterior grid points.These cases establish the bound needed for the discrete stability estimate.
- Maximum principle: Positive weights imply a maximum principle for the discrete operator, which supplies the stability mechanism for the convergence proof.The proof uses the maximum principle together with the sampled super-solution.
- Convergence theorem: The maximum solution error is bounded by a constant times the maximum truncation error when the exact solution is sufficiently smooth.This transfers the consistency estimate to the global solution estimate.
- Convergence theorem: The convergence rates are O(h2−α) for wT and O(h3−α) for wQ.The rate is inherited from the global accuracy of the corresponding discrete operators.
5. Far field boundary conditions.
Because the discrete operator contains an infinite sum, finite-domain computations require approximating the far field. Analytical and asymptotic treatments of the omitted tail preserve accuracy and are especially important for slowly decaying solutions.
- Domain truncation: Directly truncating the infinite operator can accumulate large errors because the solution and weights may decay slowly.The tail weights satisfy wj ≈ j^-1-α, so omitted contributions can remain significant.
- Truncation correction: Analytical integration of the truncated kernel terms reduces the domain-truncation error while retaining O(h3−α) or O(h2−α) accuracy for wQ and wT.The bound is independent of the truncation index M.
- Boundary data: For extended Dirichlet data, the exterior contribution is known when LW≥2L and vanishes identically when g≡0.Its evaluation may require numerical quadrature but is independent of the finite-difference scheme.
- Algebraic tails: For algebraically decaying whole-space solutions, asymptotic tail expansions approximate the dominant far-field contribution that cannot be neglected on practical domains.The method uses hypergeometric-function expressions for the leading tail terms.
- Numerical validation: Including the asymptotic boundary terms substantially reduces error on the small domain [−2, 2], where the interior terms alone are inaccurate.The experiment uses the crude algebraic-tail approximation with exact exponent β=1−α.
- Numerical validation: Before saturation, the L∞ error converges at rate 3−α; for fixed domain size, further grid refinement eventually stops reducing the error.The saturation highlights the importance of far-field boundary treatment.
- Extensions and cost: The approach extends to far-field conditions u(x)→c± with algebraic rates, and the convolution term can be evaluated in O(N log N) time.The computational cost is stated for fast convolution on the truncated domain.
6. Numerical experiments.
Numerical experiments test operator accuracy, equation-solution convergence, and the obstacle problem across decay rates, regularity levels, domain sizes, and discretization weights. Results generally confirm predicted rates, while limited regularity and finite domains reduce or eventually saturate accuracy.
- Numerical experiments: The experiments test operator accuracy and solution convergence, including exponentially and algebraically decaying functions, nonsmooth solutions, the extended Dirichlet problem, and the obstacle problem.Tests compare grid sizes, domain sizes, and the wT and wQ weights.
- Accuracy for smooth functions: O(h2−α) accuracy is obtained for wT and O(h3−α) for wQ on smooth exponentially decaying data.The wQ rule is more accurate than wT and the method in.
- Accuracy for smooth functions: O(h2−α) for wT and O(h3−α) for wQ are recovered for algebraically decaying data when the computational domain is large enough to control truncation error.For fixed domain sizes, the error eventually saturates.
- Accuracy for non-smooth functions: O(h1−α/2) accuracy is observed for a continuous solution with discontinuous derivatives, with the largest error near x = ±1.The reduced rate is consistent with the solution's lower regularity.
- Accuracy for non-smooth functions: For a C1 solution with discontinuous second derivatives, wT achieves O(h2−α), while wQ achieves O(h2−α/2), below the smooth-data rate.The wQ rate improves with smoothness but still requires more regularity to reach O(h3−α).
7. Conclusions.
The paper develops and validates a one-dimensional finite difference/quadrature scheme for the fractional Laplacian, proving convergence and achieving higher formal accuracy with quadratic tail weights. Far-field truncation is addressed asymptotically, while extensions to higher dimensions and singular solutions remain limited.
- Discretization: O(h3−α) formal accuracy is obtained with quadratic interpolation tail weights, compared with O(h2−α) for linear weights.The linear-weight scheme is retained mainly because its derivation is simpler.
- Discretization: Positive discrete weights preserve key properties of the singular-integral operator, and the scheme recovers centered finite differences as α→2−.The weight sum also provides a CFL stability condition for explicit time discretizations.
- Boundary treatment: An asymptotic approximation for far-field Dirichlet data improves accuracy by an order of magnitude when truncating the infinite computational domain.Without this treatment, truncation error can dominate, especially for algebraically decaying solutions requiring large domains.
- Numerical validation: Known smooth solutions validate the predicted convergence rates, while algebraically decaying solutions require larger domains than exponentially decaying solutions.The study also computes the singular Getoor solution on the line and for the extended Dirichlet problem on an interval.
- Nonlinear problems: The scheme enables computation of the fractional-Laplacian obstacle problem, with a computed convergence rate of O(h1+α/2).This uses a consistent, stable scheme with positive weights to compute a nonlinear viscosity solution.
- Limitations: For singular solutions, rigorous higher-order convergence techniques remain an open problem, although convergence for smooth solutions is bounded by truncation error.Weak-solution convergence rates can be fractional with an unknown exponent.
- Comparisons: Quadratic weights are the most accurate among the stable compared schemes, whereas spectral methods become unstable for α>1.The method was also compared with fractional-derivative, Fourier-based, and schemes.
- Limitations: Higher-dimensional extension remains challenging because the treatment of the operator’s tail does not extend, except for truncated tails or exponentially decaying extended boundary data.Singular-integral treatment and numerical quadrature extend more naturally than the tail treatment.