Numerical Methods for Fractional Diffusion

Andrea Bonito; Juan Pablo Borthagaray; Ricardo H. Nochetto; Enrique Otarola; Abner J. Salgado

Source-linked AI summary

Numerical Methods for Fractional Diffusion

Andrea Bonito, Juan Pablo Borthagaray, Ricardo H. Nochetto, Enrique Otarola, Abner J. Salgado

arXiv:1707.01566v1Jul 5, 2017math.NA

TL;DR

Fractional diffusion requires numerical methods for nonlocal operators with distinct spectral and integral definitions. The paper presents three finite element strategies based on extension PDEs, singular integral formulations, and Dunford–Taylor representations, and analyzes their convergence and computational behavior. The methods achieve reported convergence rates, while their applicability is bounded by multidimensional singular-kernel difficulties and suboptimal or unresolved cases.

Problem
Fractional diffusion involves nonlocal operators, singular kernels, and limited boundary regularity, creating numerical difficulties especially in multidimensional settings.
Method
The paper presents three finite element methods: an extension-based PDE method for the spectral operator, a direct integral formulation, and Dunford–Taylor discretizations for both definitions.
Results
The analyzed schemes provide convergence estimates, including (#T_Y)^-1/(d+1) decay for the basic extension FEM and second-order convergence up to a logarithmic term for time discretization.
Takeaways & Limitations
The methods offer complementary ways to approximate spectral and integral fractional Laplacians, with the Dunford–Taylor algorithm reported as advantageous for multiprocessor scaling.
Takeaways & Limitations
The a posteriori error analysis remains problematic in dimensions d > 1 because evaluating the fractional Laplacian is singular near mesh skeletons.

Abstract

from arXiv · show

We present three schemes for the numerical approximation of fractional diffusion, which build on different definitions of such a non-local process. The first method is a PDE approach that applies to the spectral definition and exploits the extension to one higher dimension. The second method is the integral formulation and deals with singular non-integrable kernels. The third method is a discretization of the Dunford-Taylor formula. We discuss pros and cons of each method, error estimates, and document their performance with a few numerical experiments.

1. Introduction

The paper introduces fractional diffusion as a nonlocal model motivated by anomalous processes and distinguishes integral and spectral fractional Laplacians. It presents three finite element strategies: an extension-based PDE method, an integral formulation, and Dunford–Taylor approaches.

Classical Brownian-motion diffusion models can be inadequate for anomalous diffusion observed across physical and social environments.
Fractional diffusion is a nonlocal process requiring information from all spatial points.
Boundary regularity is limited, producing reduced convergence rates.
The integral and spectral definitions of the fractional Laplacian do not coincide and produce different boundary behavior.
The extension formulation converts the spectral operator into a local, higher-dimensional PDE, but introduces singular behavior near the extended boundary.
The paper presents three FEMs: an extension-based method for the spectral definition, an integral method for singular kernels, and Dunford–Taylor methods for both operators.

2. The Spectral Fractional Laplacian

The spectral fractional Laplacian is studied through an extension-based PDE formulation on a truncated cylinder over the domain. The framework also extends to more general symmetric elliptic operators.

The spectral fractional Laplacian is discretized using PDE techniques based on an extension formulation.
The analysis uses a semi-infinite cylinder and its finite-height truncation, with the extended variable distinguished from spatial variables.

2.1. Extension Property.

The extension property represents the spectral fractional Laplacian through a weighted, local divergence-form problem in one additional dimension. Its trace and conormal derivative recover the fractional operator, while the weight creates nonuniform ellipticity near the boundary.

The extension problem is a weighted divergence-form PDE with zero lateral boundary data and prescribed conormal derivative on the domain boundary.
The weight y^α belongs to the Muckenhoupt A2 class and requires weighted Lebesgue and Sobolev spaces.
The extension is locally variationally tractable but nonuniformly elliptic because y^α degenerates or blows up as y ↓ 0.
The weighted Sobolev setting provides trace control from the extension space to Hs(Ω).
The fractional Laplacian is recovered from the extension through the trace of the solution and its Dirichlet-to-Neumann relation.

2.2. Regularity.

The extended solution admits a spectral representation whose regularity is governed by one-dimensional mode functions. Regularity is substantially weaker in the added dimension, which shapes the finite element analysis.

The extension solution is represented through eigenfunction coefficients and mode functions solving a one-dimensional boundary-value problem.
For s = 1/2, the mode functions are exponential, while other fractional powers involve modified Bessel functions.
The regularity results support finite element error analysis, although additional estimates may be unnecessary for the subsequent development.
The analysis derives global regularity estimates for the α-harmonic extension when f ∈ H1−s(Ω).
The extension has substantially worse regularity in the added dimension than in the original spatial variables.

2.3. Truncation.

The unbounded extension domain is truncated to a finite cylinder because the extension solution decays exponentially in the extended variable. Homogeneous Dirichlet data at the truncation boundary therefore introduces an exponentially small error.

Exponential decay of the extension solution in y justifies truncating the unbounded cylinder to C_Y = Ω × (0, Y).The truncated problem imposes homogeneous Dirichlet data on y = Y.
The truncation error is exponentially small in the truncation parameter Y.The estimate is stated for the solutions of the original and truncated problems and depends on the first Dirichlet Laplace eigenvalue.

2.4. FEM: A Priori Error Analysis.

The a priori FEM analysis discretizes the truncated extension problem on tensor-product meshes and derives anisotropic interpolation and error estimates. Graded meshes in the extended direction recover near-optimal complexity, while quasi-uniform meshes are suboptimal in the original domain.

The extension problem is approximated on a tensor-product mesh T_Y = T × I_Y using continuous piecewise linear finite elements with homogeneous Dirichlet data.The mesh has #T_Y = M#T degrees of freedom and scales as M^(d+1) when #T ≈ M^d.
The Galerkin approximation is characterized by a best-approximation property in the weighted energy norm.Existence and uniqueness follow from the finite element space and the Lax–Milgram framework.
Anisotropic interpolation estimates separate errors in the spatial and extended directions and motivate grading the mesh in y.Uniform meshes exhibit a suboptimal rate because of the behavior of the second y-derivative near y = 0.
The graded mesh uses nodes y_m = m^γM^−γY with γ > 3/(1 − α) = 3/(2s), yielding a priori estimates in terms of degrees of freedom.The truncation parameter is selected as Y ≈ log(#T_Y).
Except for a logarithmic factor, quasi-uniform complexity decays as (#T_Y)^−1/(d+1), whereas geometric grading improves the decay to (#T_Y)^−1/d.The a priori estimates require f ∈ H^(1−s)(Ω) and a domain satisfying the stated regularity condition.

2.5. Numerical Experiments.

Numerical experiments on square and circular domains use graded meshes to test the extension FEM. The observed convergence rates agree with the a priori theory across several values of s, although the approximation of the original solution can remain suboptimal.

The experiments use graded meshes, Gaussian quadrature, and conjugate gradients with an ILU preconditioner and residual tolerance 10^−12.The computations are performed for d = 2 within the deal.II library.
Square Domain: The graded-mesh approximation is near-optimal for the extension U but suboptimal for the original solution u.
Square Domain: For the square-domain manufactured solution, graded meshes produce the convergence rate predicted by Theorem 2.4 for s = 0.2 and s = 0.8.The reported experimental rate is (#T_Y)^−1/3.
Circular Domain: For the circular-domain experiment, errors in the weighted extension norm agree with Theorem 2.4 for s = 0.3 and s = 0.7.The test uses eigenfunction-based data and the same graded-mesh construction.

2.6. FEM: A Posteriori Error Analysis.

Residual error estimators fail for the weighted extension problem because the boundary integral at y = 0 is not meaningful. An anisotropic estimator based on local problems on cylindrical stars is therefore developed and tested with adaptive refinement.

The a priori theory may fail when data or domain regularity assumptions are violated, producing x-variable singularities that make quasi-uniform refinement inefficient.Adaptive refinement driven by a posteriori estimation is needed in such cases.
Residual error estimators cannot be used because their derivation relies on a strong local residual and integration by parts.The boundary integral is meaningless at y = 0 when α ∈ (−1, 1), so the derivation fails at its first step.
The proposed estimator handles mesh anisotropy and the coefficient y^α through local problems on cylindrical stars.Under certain assumptions, it is equivalent to the error up to data-oscillation terms.
The ideal local spaces are infinite dimensional, so the practical estimator replaces them with continuous piecewise polynomial spaces of degree at least two.Quadratic tensor-product spaces are used on quadrilaterals, while simplices use quadratic polynomials augmented by a cubic bubble.
For Ω = (0, 1)^2 and f ≡ 1, adaptive refinement restores the optimal decay rate for s < 1/2 despite incompatible data for the a priori theory.The incompatibility arises because f lacks a vanishing trace and does not belong to H^(1−s)(Ω).

2.7. Extensions and Applications.

The paper extends its discussion to time-dependent, nonlinear, and optimal-control problems involving spectral fractional diffusion, while noting an unresolved regularity issue for incompatible data.

Extensions and Applications: The extension framework supports fractional parabolic problems with Caputo derivatives and yields fully discrete convergence rates consistent with solution regularity.For nonzero initial data, the solution is singular as t decreases to zero, and this affects regularity-based convergence analysis.
Extensions and Applications: Extending the analysis to space-time fractional wave equations with γ ∈ (1, 2] remains under investigation.
Extensions and Applications: Finite element methods have been developed for elliptic and parabolic obstacle problems with the spectral fractional Laplacian.These analyses require combining Sobolev and Hölder regularity with positivity-preserving interpolation on anisotropic meshes.
Extensions and Applications: Spectral fractional Laplacians have also been used in stationary and parabolic PDE-constrained optimization and optimal-control problems with convergence results for fully discrete schemes.
Extensions and Applications: Handling incompatible right-hand sides f /∈ H1−s(Ω) remains open because it produces boundary singularities governed by the stated regularity behavior.

3. The Integral Fractional Laplacian

The integral fractional Laplacian is formulated through a nonlocal energy bilinear form and approximated directly with finite elements, while singular kernels and limited boundary regularity drive the analysis and mesh design.

The Integral Fractional Laplacian: The integral formulation uses the Hs(Rd)-seminorm and a symmetric bilinear form over Q = (Ω × Rd) ∪ (Rd × Ω), yielding a well-posed weak problem.Well-posedness follows directly from the Lax–Milgram lemma.
Regularity: Smooth right-hand sides do not generally produce solutions smoother than Hs+1/2−ε(Ω), reflecting boundary-driven limited regularity and reduced convergence rates.The reduced regularity is also linked to the lower regularity of integral-Laplacian eigenfunctions.
Regularity: For Lipschitz domains, solution regularity depends on s and data assumptions, with Theorem 3.3 providing Hs+1/2−ε regularity under corresponding conditions.For s > 1/2, the theorem guarantees at least H1(Ω) regularity.
FEM: A Priori Error Analysis: A direct continuous finite element method applies in arbitrary dimensions, with energy- and L2-error estimates established for quasi-uniform and graded meshes.For s > 1/2 in two dimensions, graded meshes exploit weighted regularity near the boundary.
Numerical Experiments: Numerical experiments agree with predicted convergence rates, while graded meshes show an additional increased order for s ≤ 1/2 not covered by Theorem 3.10.The L2 experiments for s = 0.25 and s = 0.75 also agree with the predicted rates.
FEM: A Posteriori Error Analysis: The a posteriori error analysis has restrictions for d > 2 and encounters problematic evaluation of (−∆)sφv near the mesh skeleton when d > 1.The latter computation is tractable in one dimension, but improvement remains open.

4. Dunford-Taylor Approach for Spectral and Integral Laplacians

The Dunford-Taylor approach represents spectral and integral fractional Laplacians through equivalent operator formulas, then combines sinc quadrature with finite element approximations on general Lipschitz domains.

Dunford-Taylor Approach: The method uses the Dunford-Taylor representation for the spectral fractional Laplacian and an equivalent representation for the integral fractional Laplacian.
Dunford-Taylor Approach: Both variants use sinc quadratures and finite element approximations of the resulting integrands at each quadrature point.
Dunford-Taylor Approach: The approach is proved efficient on general Lipschitz domains, but the integral representation produces a non-conforming method based on approximating stiffness-matrix actions.The spectral representation permits direct approximation of the solution.
Finite Element Discretization: The approximation is expected to deliver the same rate of convergence as the underlying finite element discretization, subject to the stated condition f ∈ L2(Ω).
Finite Element Discretization: Finite element approximation rates depend on elliptic regularity of index α and on the regularity parameters governing the fractional solution and data.Theorem 4.2 introduces γ := max{r + 2α*−2s, 0} in its rate characterization.

4.2. Sinc Quadrature.

The Dunford-Taylor method approximates the fractional operator through sinc quadrature, whose exponential convergence is balanced against finite element error. The approach extends to broader operators, geometries, and space-time fractional problems, with stated convergence and implementation trade-offs.

4.2. Sinc Quadrature.: Sinc quadrature uses uniformly spaced points after a logarithmic change of variables, with N+ and N− chosen to balance quadrature and finite element errors.The choice of N+ and N− is designed to improve robustness with respect to s.
4.2. Sinc Quadrature.: Exponential convergence of the sinc quadrature follows from the decay and holomorphic properties of the Dunford-Taylor integrand.This property underpins the quadrature approximation used by the method.
4.2. Sinc Quadrature.: For quasi-uniform meshes, the method achieves quasi-optimal linear order for f ∈H1−s(Ω) and rate 2α∗ = 2−s > 1 for smoother data.The optimal smoother-data rate applies with σ = 2(1−s) when s ≤ 2/3 and σ = s when s > 2/3.
4.2. Sinc Quadrature.: The Dunford-Taylor method appears advantageous for multiprocessor computation because it exhibits good strong and weak scaling properties.Strong scaling fixes the degrees of freedom while increasing processors; weak scaling keeps degrees of freedom per processor constant as problem size grows.
4.2. Sinc Quadrature.: The algorithm requires N+ + N− + 1 independent standard Laplacian finite element solves per quadrature point, which are then aggregated into U^k.Its implementation is described as straightforward when starting from a Poisson finite element solver.
4.2.1. Extensions.: The framework extends to symmetric elliptic and regularly accretive operators, surface problems, non-symmetric advection-diffusion, and space-time fractional diffusion.For space-time problems, the contour is parametrized hyperbolically because the Mittag-Leffler function prevents deformation onto the negative real axis; the method does not require time stepping for the initial-value problem.
4.2.1. Extensions.: For u0 = 0, midpoint quadrature on a graded time mesh yields second-order convergence up to a logarithmic term in the number of time intervals.The time mesh is graded toward t = 0 and uses N log(N) time steps.

4.3. Integral Laplacian.

The integral Laplacian method reformulates the nonlocal bilinear form through auxiliary whole-space problems, then approximates quadratures, truncated domains, and finite elements. Its error analysis combines consistency estimates with Strang’s lemma, while the approach faces mesh-generation and convergence-rate limitations.

The spectral-Laplacian strategy does not apply to the integral Laplacian because formulas like (1.10) are not well defined.
The method represents the integral bilinear form using auxiliary solutions v(ψ, µ) of whole-space problems and a sinc quadrature.The quadrature uses spacing k and truncation indices N+ and N−.
Quadrature consistency is established under regularity conditions s < δ ≤ min(2 − s, σ), with σ strictly smaller than 3/2.
The whole-space auxiliary problems are truncated to bounded domains whose size depends on the quadrature point, then discretized with standard finite elements.The resulting meshes must match the original mesh inside Ω, and constructing these partitions is identified as the method’s bottleneck.
Truncation consistency decays exponentially with the truncation parameter M, while the fully discrete scheme combines quadrature, truncation, and finite-element approximations.The truncation estimate uses a constant independent of M and k.
The final error estimate assumes u ∈Hβ(Ω) for β ∈(s, 3/2) and yields an optimal convergence rate up to a logarithmic factor, but graded-mesh applicability remains open.For quasi-uniform meshes, the rate is about h1/2; balancing errors uses M ≈ log(1/h) and N+ ≈ N− ≈ |log h|2.