A simple and fast finite difference method for the integral fractional Laplacian of variable order

The constant-order operator (1.1) may be less sufficient for the heterogeneous effect due to the spatial variability of complex medium. To account for heterogeneity, the variable-order operators depending on the spatial location variable have been alternatively proposed. In this work, we consider the following $\alpha(x)$-th order fractional Laplacian

The variable-order extension variable order enables us to truly capture the non-smooth effects such as fractures by prescribing variable degree of smoothness across the scales. For example, the authors in [4] advocated the variable-order operator to perform regularization in image processing. They used large value of the order $\alpha(x)$ in the flat region and the smaller value in the neighborhood of edges. In [16, 7], the authors pointed out that the usage of the variable order is partially motivated from the theoretical study of galaxy rotation curves and the dark matter. Their study requires the field equation for the gravitational potential to become a variable-order fractional Poisson’s equation. Another example is from the groundwater flow. The authors in [25] justified that the super-diffusion depends on spatial variables due to the local variation of aquifer properties. We also refer the interested readers to other applications including biology [11], geophysics [28], spatial statistics [32] and so on.

When restricted to the bounded domain, equations with the fractional Laplacian of constant order as the leading operator usually have non-physical weak boundary singularity (boundary layer) unless boundary conditions are special. The boundary layer is often undesirable for the real world applications. As a remedy, the authors in [45] suggested to use the Caputo-type variable-order operator and analyzed a one-dimensional (1D) two-point boundary value problem. In Section 3.2, we use a variable-order fractional Laplacian requiring variable order $\alpha(x)=2$ along the boundary and observe in Example 3.3 that a full convergence order can be achieved on uniform meshes, which may suggest full regularity of solutions to such equations. This observation may lead to more practical fractional calculus in real applications.

Changing directly the constant-order into the variable-order may increase not only the model capability but also the complexity in the computation. For example, one can use the variable-order function $\alpha(x,y)$ in place of $\alpha(x)$ in (1.2), which depends on both the spatial variable $x$ and the integration variable $y$ in the singular integral in the recent work [8, 29]. The wellposedness of corresponding fractional Poisson equations or other partial differential equations can be straightforwardly established from the variational formulation. However, the definition does not allow a simple Fourier transform of the fractional Laplacian as in the constant-order fractional Laplacian and Fourier-transformed based computational methods are not readily applicable.

The authors in [7] used the Fourier transform to extend the variable-order fractional Laplacian. However, the Fourier-type definition therein (see also Remark 2.2) is different from the one considered in this work. Here we prove that the singular integral type definition (1.2) can be rewritten via Fourier transforms but is different from the definition in [7]; see Theorem 2.1 for the details. One desirable feature of their definition is that the kernel is convolution-type but its explicit form is not available and have to be calculated through the forward and backward Fourier transforms.

Except for the above commonly used approaches to define the variable-order fractional Laplacian, one can also consider the spectral definition (see [44]) and harmonic-extension type definition (see [4]). For the spectral definition, the associated computation is straightforward when the eigen-values or eigen-functions are explicitly available. In this paper, we will not go through every definition but limit our attention into the singular-integral type (1.2).

The close or explicit form of the $\alpha(x)$-th order fractional Laplacian of functions seldom exists. In this work, we aim for an efficient numerical evaluation for the fractional Laplacian of variable order (1.2) with its applications in numerical solutions to nonlocal PDEs. In 1D case, we note that the authors in [46] showed the equivalence between the Fourier-type definition and two-sided Riemann-Liouville derivatives reformulation. Following upon their work, finite difference methods and fast solvers have been developed for the fractional Laplacian of coordinate-dependent type (see [30, 40]). However, to the best of our knowledge, we are not aware of any numerical results in 2D/3D in the literature due to the complexity of the operator discretization. Even if one can extend the finite difference scheme such as the popular Grunwald approximation to compute 1D variable-order fractional operators, the rigorous convergence analysis is not provided yet. The current research aims to fill this computational and theoretical gap.

For the constant-order fractional Laplacian, some progress has been made and extensive numerical methods have been developed among the computational mathematics community particularly during the past around five years. Finite or boundary element methods are developed in the context of the boundary integral equations, and they can be straightforwardly extended to solve the integral equations involving the volume integral (see [1]). However, it is extremely difficult to implement the finite element methods (FEM) particularly in 3D as it entails complicated numerical quadrature for the double integrals [12]. To avoid the expensive computation of the double integrals, meshless methods or collocation methods have been developed. For special domains such as a disk, the accurate and efficient spectral method using the eigen-functions was proposed in [43] and [19]. For a general domain, radial basis functions (RBF) methods using the standard Gaussian functions [6] and other generalized multi-quadratic functions have been developed. However, RBF methods suffer from the notorious ill-condition issue which makes it harder for solving the large-scale size problem.

For the large-scale computation and easy implementation in multi-dimensions for fractional Laplacian of constant order, finite difference methods or collocation methods have been proposed recently. Using the singularity subtraction technique, the authors in [27] proposed a simple solver via the Taylor expansion under the high-regularity assumption that the target function $u(x)\in C^{6}(\Omega)$ but the stability of the finite difference scheme is unknown. The authors in [9] converted the fractional Laplacian with the strong singularity into the one with the weak singularity through the integration by parts, and then constructed a quadrature-based finite difference approximation. However, the stability estimates were not provided and it is not clear if the associated finite difference scheme is stable when applied to the numerical solution of PDEs involving the fractional Laplacian operator. The authors in [2] presented the collocation method based on the sinc basis functions and rigorous analysis of the method was later provided in [3]. Nonetheless, the analysis therein relies on the properties of Fourier transform and cannot be directly extended to the problem considered in this work.

We remark that the Fourier transform of fractional Laplacian of constant-order allows designs of various and efficient numerical methods. However, this is not the case for the variable-order fractional Laplacian as the semi-group and symmetry properties cannot extend, and the explicit and elementary representation of the Fourier transform of variable order does not exist. Without these properties, straightforward extensions of previous working numerical methods are highly nontrivial including various discretization methods such as spectral methods in [38, 37], the fast finite element implementation [34], Dunford-Taylor reformulation based finite element methods and spectral methods and so on.

In [20] we establish an easy-to-implement and efficient finite difference method based on the generating functions approximation theory. The computation of the coefficients or weights in finite difference approximation can be efficiently calculated through the fast discrete Fourier transform and a fast solver is implemented for the structured resulting linear systems. Along this line of research, we aim to extend the idea of simple finite difference discretization in [20] to the variable-order context, and propose an efficient solver for the resulting unstructured matrix induced by the non-convolution kernel in the variable-order operator (1.2).

The significance and novelty of the paper is that it is the first work addressing a fast and accurate finite difference approximation for the numerical evaluation of the multi-dimensional fractional Laplacian of variable order (1.2). Our approximation is robust and stable in the sense that we can get a second-order approximation even for the nonsmooth piecewise variable-order function $\alpha(x)$ (see Tables 1, 2 and 3).

The approximation is derived from the discrete Fourier transform. Although the explicit form of the Fourier transform for the variable-order case does not exist, we find that it can be defined through the Fourier transform as in the constant-order case. Such equivalence allows us to discretize the singular operator in the frequency space which avoids the complicated quadrature.

Compared to the constant-order counterpart, the variable-order definition in multi-dimensions uses the non-convolution kernel and thus the corresponding discrete counterpart becomes non-Toeplitz, which brings up the extra computational difficulty regarding the efficient solver part. Nonetheless, by the Fourier transformation and the low-rank approximation of the Fourier symbol $|\xi|^{\alpha(x)}$ (see Lemma 2.4), we are able to overcome the difficulty and construct a fast solver based on the observation that the unstructured resulting matrix can be approximately decomposed into the linear combination of the structured matrices.

Several numerical results are demonstrated to show the accuracy and efficiency of the proposed scheme. Numerically, we find that the appropriately selected variable-order model can avoid the non-physical boundary singularity as in the constant-order model.

In a nutshell, the main contributions of this work are described as follows.

• •

  We establish the equivalence and provide its proof between Fourier-type definition and singular-integral form.

• •

  We derive a second-order finite difference approximation for the variable-order fractional Laplacian in multi-dimensions and present a rigorous convergence analysis.

• •

  We provide a fast solver with quasi-linear complexity to efficiently compute the variable-order operator and related PDEs.

• •

  We apply the developed finite difference discretization to solve various fractional PDEs with variable-order fractional Laplacian and carry out the extensive numerical experiments to show the accuracy and efficiency of our theory and algorithm.

We organize this paper as follows. In Section 2, we present a fractional finite difference approximation based on the equivalence between the Fourier-type definition and singular integral form. Moreover, we describe its implementation with a quasilinear fast solver, and show the accuracy and robustness of our approximation. In Section 3, we apply the discrete approximation to the numerical solution for the nonlocal elliptic PDEs. We carry out the extensive numerical experiments and demonstrate the accuracy and efficiency of our scheme. We put the proofs of the theoretical results in Section 4 for the interested readers, including the second-order approximation and the convergence analysis for the finite difference method in 1D case. Finally, we conclude this work in Section 5.

In our discussion, we assume that the variable-order function in the definition (1.2) is continuous and bounded, that is $\alpha(x)\in C(\mathbb{R}^{d})$, and

The Fourier transform and its inverse are defined respectively as follows

and

where $\mathbf{i}$ is the complex symbol $\mathbf{i}^{2}=-1$. Following the argument in the constant case as in [5], we show that the equivalence of two definitions still holds in the variable-order context.

The equivalence of the definitions inspires us to discretize the operator in the Fourier or frequency spaces rather than the direct discretization in the physical space using the complicated finite difference-quadrature.

Denote $h$ as the spatial stepsize of the grids. Let $h\mathbb{Z}^{d}$ be the set of infinite grids, i.e.,

Throughout the paper, the summation index $j$ is multi-index, so is the index $k$.

For a function $u_{h}:h\mathbb{Z}^{d}\longrightarrow\mathbb{R}$, define the semi-discrete Fourier transform as

Denote the parameterized domain $D_{h}=\big{[}-\frac{\pi}{h},\frac{\pi}{h}\big{]}^{d}$. Then we have the discrete Fourier inversion formula

For a continuous function $u(x)$, let $u_{h}(x)=u(x)$ for $x\in h\mathbb{Z}^{d}$. Inspired by the idea for 1D in [21] and our previous work [20] for multi-dimensional cases, we propose the discrete variable-order Laplacian operator as follows:

with

Introduce the following spectral Barron space (see [10, 26])

and the induced norm $\|\cdot\|_{\mathcal{B}^{s}(\mathbb{R}^{d})}$. It is clear that $\|u\|_{\mathcal{B}^{s}(\mathbb{R}^{d})}\leq c\|u\|_{\mathcal{B}^{t}(\mathbb{R}^{d})}$ for $s<t$.

For any grid function $v_{h}$ on $\Omega_{h}$, define the discrete maximum norm as

We state the main theoretical result in this work.

We postpone the proof of second-order approximation properties and put it in the section 4.

In the practical computation, we want to numerically evaluate the operator at the grid points, that is, $x=x_{j}=jh\in\Omega_{h}$. Thus, we have

where the weights

Making a change of variable $\eta=h\xi$ leads to

In the practical computation, for each fixed grid point $x_{j}$, we need to compute the coefficients $a_{n}^{(\alpha(x_{j}))}$ for $|j_{p}|\leq N_{p}$ for $p=1,\dots,d$ (Here $n=k-j$ and $N_{p}$ denotes the number of grid points in $x_{p}$ direction. For the sake of simplicity, let us assume they all equals to $N$).

To alleviate the burden of tedious notations, we drop the dependence on the variable $x_{j}$ and abbreviate it as $a_{n}^{(\alpha_{j})}$ without the confusion.

In 1D, the analytical expression of the above integrals can be derived. In fact, recall the formula [17, Page 399, 3.631.9]

Here $B(\cdot,\cdot)$ denotes the Beta function. Take $\nu=\alpha_{j}+1$, $a=2n$ and make a change of variable $x=\pi/2-t$. Then the formula (2.16) reduces to

Due to the symmetry of the integrand, the coefficients (2.15) can be reduced to

In multi-dimension cases, since the explicit formula is difficult to derive, we will use the standard numerical integration to compute the integrals. Specifically, take an integer number $M>\max\{N_{p},\,p=1,...d\}$ (here $N_{p}$ denote the number of grid points in $x_{p}$ direction) and step-size $\delta=2\pi/M$. Denote $\varphi(\eta)=[\sum_{p=1}^{d}4\sin^{2}(\frac{\eta_{p}}{2}))]^{\alpha_{j}/2}$. Applying the trapezoidal rule, we have

With the expression above we can use Matlab built-in function ‘ifftn’ to compute the coefficients $\tilde{a}^{(\alpha_{j})}_{n}$ for $0\leq n_{p}\leq M-1$ efficiently with the computational cost $\mathcal{O}(M^{d}\log^{d}M)$.

In this subsection we present a fast solver for the numerical evaluation of the fractional Laplacian. In general, we are interested in the numerical evaluation of the target function $u(x)$ restricted in the bounded computational domain $\Omega$. Suppose that $\Omega$ is a square (or a rectangular domain if we take a nonisotropic different mesh size in each spatial direction in the last subsection). If not, we can embed the general bounded domain $\Omega$ into a larger square (or rectangular) domain $\tilde{\Omega}$, consider the evaluation on the domain $\tilde{\Omega}$ instead and do the restriction in the domain $\Omega$. The following lemma is a key to develop a fast algorithm.

As direct consequence of the above lemma, we can use the superposition principle, the linear combination of the constant-order discrete fractional operators to approximate the variable-order fractional Laplacian to achieve a fast computation. That is

Based on the above discussion, we describe our algorithm for the numerical evaluation as follows:

Note that, in Step 2, we refer the readers to [20] for implementation details.

In the following example, we take the Gaussian function as our test function to illustrate the accuracy of the theoretical result.

In this experiment, the variable-order Laplacian of the Gaussian function is given as (see the Appendix for the derivation)

The error is measured as $E_{\infty}(h)=||(-\Delta_{h})^{\alpha(\cdot)/2}u-(-\Delta)^{\alpha(\cdot)/2}u||_{L_{h}^{\infty}}$, and the convergence order is estimated by $\log_{2}({E_{\infty}(2h)}/{E_{\infty}(h)})$.

Note that the function $u(x)$ is smooth and satisfies the regularity requirement in Theorem 2.3. Here we examine the effect of the variable order $\alpha(x)$ on the accuracy of the finite difference approximation. Tables 1-3 show the $\operatorname{\mathit{L_{h}^{\infty}}-\mathit{norm}}$ errors and convergence orders for different $\alpha(x)$ in 1D/2D/3D, respectively. Three different $\alpha(x)$ with different range and smoothness are tested. From the tables, we observe that our approximation is of second order, which agrees with the predicted convergence order in Theorem 2.3. In particular, the approximation is very robust in the sense that we are still able to get the second-order convergence even for the non-smooth piece-wise function $\alpha(x)$.

In this section, we will apply the proposed discrete operator into solving the elliptic equations and parabolic equations with the variable-order fractional Laplacian. For the parabolic equations, after the semi-discretization in temporal direction, they are reduced to the elliptic one at the each time step. Thus we only discuss the elliptic case as below:

where the reaction coefficient $b(x)$ is non-negative, the right-hand-side function $f(x)$ is given, $\Omega$ is the standard bounded domain with the sufficiently smooth boundary and $\Omega^{c}$ is the complement of $\Omega$ in $\mathbb{R}^{d}$.

Consider Equation (Eq.) (3.1) at the grid points $x_{j}=jh$ and we obtain

Replacing the fractional Laplacian by the difference operator defined in Eq. (2.9), Eq. (3.3) becomes

with the truncation error $T_{u}$ depending on the solution $u$. By Theorem 2.3, we assume $u\in\mathcal{B}^{s+\alpha_{\max}}(\mathbb{R}^{d})$ with a positive constant $s\leq 2$, such that there exists a constant $c_{u}$ satisfying

Omitting $T_{u}$ in Eq. (3.4) and denoting by $u_{j}$ the numerical approximation of $u(x_{j})$, and $\alpha(x_{j})=\alpha_{j}$ (the short notation $\alpha_{j}$ will be frequently used throughout the following without the confusion), $b_{j}=b(x_{j})$ and $f_{j}=f(x_{j})$, we get the finite difference scheme

In this subsection, several numerical examples are presented to show the efficiency and accuracy of the schemes. We first apply the proposed finite difference scheme to solve the elliptic equation in Example 3.3 and the time-dependent problem in Example 3.4, respectively. To show the wide application of finite difference method, the finite difference scheme is applied to solve Allen-Cahn equation in Example 3.5. In addition, we consider the coexistence of anomalous diffusion problem in Example 3.6. In the last Example 3.7, we turn to the time-dependent problem in 3D.

Throughout the examples, the low-rank approximation terms $r$ is taken as $7$ to compute the Chebyshev points and the quadrature number $M$ is taken as $2^{14}$ to compute the coefficients $a^{(\alpha_{j})}_{n}$ in Step 2 so that the accuracy of numerical solution will not be polluted. The tolerance of the BiCGSTAB (see [39]) method is set as $10^{-16}$ and the initial guess is fixed as zero in our simulations. All the simulations are performed in MATLAB 2021b on a 64-bit Ubuntu with Intel Xeon(R) Gold 6240 CPU processor at 2.60 GHz and 126 GB of RAM.

In Case 1, we take $\beta=4$. Since the right hand side $f(x)$ is unknown, we need to compute it with very fine step-size which is set as $f^{ref}\approx f_{h}=(-\Delta_{h})^{\alpha(\cdot)/2}u+u$ with small $h=2^{-9}$. In Table 4, the $\operatorname{\mathit{L_{h}^{\infty}}-\mathit{norm}}$ errors and convergence orders are listed with different $\alpha(x)$. Here, the $\operatorname{\mathit{L_{h}^{\infty}}-\mathit{norm}}$ error is measured as $E_{\infty}(h)=||u_{h}-u||_{L_{h}^{\infty}}$. Furthermore, from Table 4, we can observe that the convergence order of our scheme is of second order.

Next we consider Case 2, where $f(x)=1$ and the exact solution $u(x)$ is unknown. Table 5 shows the errors and convergence orders with different $\alpha(x)$. Here, the $\operatorname{\mathit{L_{h}^{\infty}}-\mathit{norm}}$ error is measured as $E_{\infty}(h)=||u_{h}-u_{h/2}||_{L_{h}^{\infty}}$. From the table, we observe that the convergence order is less than second order and is around $\tilde{\alpha}$, where $\min\limits_{x\in\Omega}\alpha(x)/2\leq\tilde{\alpha}\leq\max\limits_{x\in\Omega}\alpha(x)/2$.

In particular, if $\alpha(x)$ is a non-smooth piece-wise function, Table 4 shows that our scheme (3.6)-(3.7) still maintains second-order convergence when the exact solution $u(x)$ is smooth enough. But the convergence order reduces to $\min\limits_{x\in\Omega}\alpha(x)/2$ when the right-side data $f(x)=1$ and the exact solution $u(x)$ is unknown; see Table 5. The reason for the reduced convergence is that the solution may have the interior singularity inherited from the variable-order function.

To recover the second-order convergence rate, we take $\alpha(x)\in C(\mathbb{R}^{2})$ satisfying $\alpha(x)=2$ on $\partial\Omega$ and $\alpha(x)<2$ in $\Omega$. Figure 1 shows the profile of the function $\alpha(x)=0.8+1.2\max{(\left|x_{1}\right|,\left|x_{2}\right|)}$. From Table 6, we observe that the appropriately selected variable-order function $\alpha(x)$ enables us to recover the second-order convergence. We remark that for the case 1 corresponding to the first column results in Table 6, we further test the smaller case $h=1/128$ and $1/256$ and the convergence orders are 1.69 and 3.02, respectively. Thus the average order is around 2. For the case 2 in the second column, we also take smaller stepsize $h=1/128$ and the convergence order is 2.1. This implies that setting the variable-order $\alpha(x)$ can change the regularity of the solution and avoid the non-physical boundary singularity which reduces the convergence order as in Table 5.

In this example, we test the convergence order of spatial direction at the final time $T=N\Delta t=0.5$, where $\Delta t$ is time step-size. Since the solution decays to zero, we truncate the finite domain in space, $\Omega=[-4,4]^{2}$ as our computational domain. For time discretization we use Crank-Nicolson scheme. Since the exact solution is unknown, the error is measured as $\|u^{N}(\Delta t,h)-u^{N}({\Delta t/2,h/2})\|_{L_{h}^{\infty}}$.

Table 7 shows the errors and convergence orders in $\operatorname{\mathit{L_{h}^{\infty}}-\mathit{norm}}$ for different $\alpha(x)$, from which the second-order convergence can be observed. This is in agreement with the second-order approximation of our theoretic prediction.

The computational domain is $\Omega=[0,1]^{2}$, and $u$ is the phase field function. The initial condition is chosen as

where the constant $\kappa$ describes the diffuse interface width, the function $d_{i}(x)=\sqrt{(x_{1}-a_{i})^{2}+(x_{2}-b_{i})^{2}}$ and $a_{1}=b_{1}=0.42,a_{2}=b_{2}=0.58$. Here, we apply our method to study the benchmark problem coalescence of two ”kissing” bubbles in the phase field models. Letting $\overline{u}=u+1$, we can reformulate the problem (3.12) as an equation of $\overline{u}$ with the extended homogeneous boundary conditions. Three-level linearized finite-difference scheme (see [41]) is used in the simulation.