# Numerical Approximation of the Fractional Laplacian on R Using Orthogonal Families

In this paper, using well-known complex variable techniques, we compute explicitly, in terms of the _2F_1 Gaussian hypergeometric function, the one-dimensional fractional Laplacian of the Higgins functions, the Christov functions, and their sine-like and cosine-like versions. After discussing the numerical difficulties in the implementation of the proposed formulas, we develop a method using variable precision arithmetic that gives accurate results.

• 2 publications
• 2 publications
• 4 publications
08/24/2019

### A Pseudospectral Method for the One-Dimensional Fractional Laplacian on R

In this paper, we propose a novel pseudospectral method to approximate a...
04/19/2022

### Efficient Monte Carlo Method for Integral Fractional Laplacian in Multiple Dimensions

In this paper, we develop a Monte Carlo method for solving PDEs involvin...
03/05/2021

### Highly accurate operator factorization methods for the integral fractional Laplacian and its generalization

In this paper, we propose a new class of operator factorization methods ...
06/04/2020

### Metastable Speeds in the Fractional Allen-Cahn Equation

We study numerically the one-dimensional Allen-Cahn equation with the sp...
01/14/2022

### Explicit formulas for concatenations of arithmetic progressions

The sequence (Sm(n))_n⩾ 0: 1, 12, 123, … formed by concatenating the fir...
02/13/2020

### Generalised Hermite functions and their applications in Spectral Approximations

In 1939, G. Szegö first introduced a family of generalised Hermite polyn...
08/06/2021

### Computing solution landscape of nonlinear space-fractional problems via fast approximation algorithm

The nonlinear space-fractional problems often allow multiple stationary ...

## 1 Introduction

In this paper, we obtain explicit expressions for the one dimensional fractional Laplacian of the Higgins and Christov functions. We also explain how to implement these results efficiently in Matlab [16], and give numerical examples as an application. More precisely, we work with the following definition [12] of the fractional Laplacian:

 (−Δ)α/2u(x)=cα∫∞−∞u(x)−u(x+y)|y|1+αdy, with cα=α2α−1Γ(1/2+α/2)√πΓ(1−α/2), α∈(0,2). (1)

We recall that the Fourier transform of

is , and that the inverse Fourier transform is ; hence, another definition of the fractional Laplacian consistent with (1) is given by associating the operator with the Fourier symbol:

 [(−Δ)α/2u]∧(ξ)=|ξ|α^u(ξ). (2)

We note that, when , according to (2), . We also observe, however, that, when the Fourier transform of does not exist in the classical sense, the limit is singular. For instance, take , whose Fourier transform is the Dirac delta distribution, i.e., ; then

 (−Δ)α1={1,α=0,0,α>0. (3)

On the other hand, when , ; and, when , , where denotes the Hilbert transform [9]:

 H(u)(x)=1π∫∞−∞u(y)x−ydy, (4)

which has the Fourier symbol . Thus, formally, .

The fractional Laplacian appears in a number of applications (see, for instance, [13, Table 1] and the references therein). In what regards its numerical computation, most references in the literature involve truncation of the domain (see for instance the works [10], [15], and [1], where, although a truncation of the domain is not explicitly given, the method relays on the spectral definition of the fractional Laplacian on a bounded domain). However, in a recent paper [6], we proposed a pseudospectral method to compute the fractional Laplacian of a bounded function on without truncation. The idea was to map into the finite interval by means of the change of variable , and then obtain the trigonometric Fourier series expansion of . Therefore, the central point of [6] was the efficient and accurate numerical computation of the fractional Laplacian of , for ; observe that the resulting expressions for are quite different, depending on whether

is even or odd. At this point, a crucial observation is that the equivalent on

of the functions (i.e., the even case) are precisely the complex Higgins functions defined in [4] (see also [8, 14]):

 λn(x)=(ix−1ix+1)n,n∈Z, (5)

because . Indeed, in this paper, thanks to the structure of , and after rewriting (1) as

 (−Δ)α/2u(x)=cαα∫∞0ux(x−y)−ux(x+y)yαdy, (6)

we obtain an explicit expression of their fractional Laplacian by means of contour integration:

 (−Δ)α/2λn(x)=⎧⎪ ⎪⎨⎪ ⎪⎩0,n=0,−2|n|Γ(1+α)(isgn(n)x+1)1+α2F1(1−|n|,1+α;2;2isgn(n)x+1),n∈Z∖{0}, (7)

where is the Gaussian hypergeometric function (see for instance [2, Ch. 15]). From this result, we also derive several other related ones outlined below.

The structure of this paper is as follows. In Section 2, we prove (7), which constitutes the main result. We observe that is a complete orthogonal system in with weight , because is a complete orthonormal system in , normalized by . Therefore, the related family of functions

 μn(x)=(ix−1)n(ix+1)n+1,n∈Z, (8)

known as the complex Christov functions [4, 7, 14, 18], form a complete orthogonal system in (normalized by the factor ).

Throughout this paper, we use the definition and notation from [4] for the following families of functions:

• Cosine-like Higgins functions:

 CH2n(x)=λn(x)+λ−n(x)2,n=0,1,2,… (9)
• Sine-like Higgins functions:

 SH2n+1(x)=λn+1(x)−λ−n−1(x)2i,n=0,1,2,… (10)
• Cosine-like Christov functions:

 CC2n(x)=μn(x)−μ−n−1(x)2,n=0,1,2,… (11)
• Sine-like Christov functions:

 SC2n+1(x)=−μn(x)+μ−n−1(x)2i,n=0,1,2,… (12)

In Section 3, starting from (7), we calculate the fractional Laplacian of (8)-(12).

Even if the fractional Laplacian of all the families considered here can be computed accurately with the technique explained in [6], expressions like (7) have the advantage of being very compact and, hence, it is effortless to use them in numerical applications, provided that fast accurate implementations of the Gaussian hypergeometric function are available. Therefore, in Section 4, using Matlab, we test their adequacy from a numerical point of view, comparing the numerical results with those in [6]. On the one hand, for moderately large values of , the use of variable precision arithmetic seems unavoidable; on the other hand, our implementation of largely outperforms that of Matlab. Finally, even if the method in [6] is faster, the method in this paper is much easier to implement and still competitive for not too large values of .

## 2 Fractional Laplacian of the complex Higgins functions

Before we proceed, let us recall some well-known definitions. Given , the generalized binomial coefficient is defined by

 (zn)=⎧⎨⎩z(z−1)…(z−n+1)n!,n∈N,1,n=0.

where . Therefore, if is a nonnegative integer, and , then . Furthermore, it is immediate to check that, for all and ,

 (zn)=(−1)n(n−1−zn).

We will also need the Pochhammer symbol, which represents the rising factorial, and is defined by

 (z)n={z(z+1)…(z+n−1),n∈N,1,n=0.

Observe that, when is not zero or a negative integer, an equivalent definition is

 (z)n=Γ(z+n)Γ(z);

in particular, when ,

 (z)n=(z+n−1)!(z−1)!.

Remark that, if is a negative integer or zero, and , then . Moreover, the following identities will be useful, too:

 (−z)n =(−1)n(z−n+1)n, (zn) =(z−n+1)nn!=(−1)n(−z)nn!.

The Pochhammer symbol also appears in the definition of the Gaussian hypergeometric function (see for instance [2, Ch. 15]). Let ; then, is defined by

 2F1(a,b;c;z)=∞∑k=0(a)k(b)k(c)kzkk!. (13)

In general, the infinite series converges for . However, in our case, we take to be a negative integer, so the sum is finite, because of the properties of the Pochhammer symbol. More precisely, in this paper, we are interested in the following two particular cases:

 2F1(−m,1+α;1;z)=m∑k=0(mk)(−1−αk)zk, (14) 2F1(−m,1+α;2;z)=−1αm∑k=0(mk)(−αk+1)zk=1m+1m∑k=0(m+1k+1)(−1−αk)zk, (15)

for . Observe that the identities also hold after replacing by in the sums, a fact that we will also use below. On the other hand, if , .

Bearing in mind the previous arguments, let us prove (7), which is the main result of this paper. Note that we work with the binomial coefficient rather than with the Pochhammer symbol, because we think that the former is more intuitive.

###### Theorem 2.1.

Let be defined as in (5). Then, is given by (7).

###### Proof.

The case with is trivial, because . Assume . The derivative of (5) is

 λ′n(x)=−2ni(x−i)2(x+ix−i)n−1.

Introducing this expression in (6), we get

 (−Δ)α/2λn(x)=2ni2α−1Γ(1/2+α/2)√πΓ(1−α/2)∫∞0gn(y;x)dy, (16)

with integrand

 gn(z;x)=(z+x+i)n−1(z−x+i)n+1−(z−x−i)n−1(z+x−i)n+1zα(z+x−i)n+1(z−x+i)n+1,

where, in this notation, we regard as a parameter, rather than as an independent variable.

We next compute (16) for every , by integrating along certain integration contours in , and using Cauchy’s integral theorem. Since , has a branch cut. In what follows, we consider the principal branch of the logarithm, which corresponds to ; in particular, , unless the branch cut is crossed. The branch choice determines also how we choose the contours. In Figure 1, we have depicted one such contour for , which consists of four parts, avoids the branch cut, but encloses the poles and , for every . Then, by the residue theorem, the integral along it is equal to the sum of the residues. The pieces of the contour that run parallel to the branch cut will give the approximation of the integral from to ; the other pieces will give integrals that tend to zero, when tends to one contour that encloses , except for the brach cut. More precisely, for every , we take , such that , and , such that , for instance, . We also take , such that, with and fixed, . Then, we define

 C=C1∪CR∪C2∪Cr,

with

 C1 ={−y−iδ:y∈(−rcosθ1,−Rcosθ1)}, CR ={Reθi: θ∈(θ1,θ2)}, C2 ={y+iδ:y∈(Rcosθ2,rcosθ2)}, Cr ={re−θi: θ∈(−θ2,−θ1)}.

Here, we have omitted the dependecy on for simplicity of notation. Now, by Cauchy’s residue theorem, we have:

 ∫Cgn(z;x)dz =∫C1gn(z;x)dz+∫Crgn(z;x)dz+∫C2gn(z;x)dz+∫CRgn(z;x)dz (17) =2πi[Res(gn(z;x),i−x)+Res(gn(z;x),x−i)]. (18)

In order to compute the residues of at and , we use the general Leibniz rule:

 dndzn(p(z)q(z))=n∑k=0(nk)p(k)(z)q(n−k)(z).

In this way, we obtain

 Res(gn(z;x),i−x) (19) =1n!dndznz−α(z+x+i)n−1∣∣∣z=i−x (20) =1n!n∑k=1(nk)[k!(−αk)z−α−k][(n−1)!(n−1−(n−k))!(z+x+i)n−1−(n−k)]∣∣ ∣∣z=i−x (21) =n∑k=1(n−1k−1)(−αk)(i−x)−α−k(2i)k−1. (22)

Likewise, we have

 Res(gn(z;x),x−i)=−n∑k=1(n−1k−1)(−αk)(x−i)−α−k(−2i)k−1. (23)

In order to compute (17), we observe that the first and third integrals tend to zero, as tends to infinity (this can be easily seen by changing to polar coordinates and recalling that ). In what regards the first and third integrals, we have that tends to , parameterized by , with , and tends to , parameterized by , with . We observe that , where the argument of is determined by the curve before taking the limit . Thus, it is on , and on . Then,

 limR→∞∫C1gn(z;x)dz=−∫∞0gn(−y;x)dy=eiπα∫∞0gn(y;x)dy,

and

 limR→∞∫C2gn(z;x)dz =∫0−∞gn(y;x)dy=∫∞0gn(−y;x)dy=−e−iπα∫∞0gn(y;x)dy.

Therefore,

 ∫Cgn(z;x)dz=(eiπα−e−iπα)∫∞0gn(y;x)dy=2πi[Res(gn(z;x),i−x)+Res(gn(z;x),x−i)].

Finally, after some manipulations, we have

 ∫∞0gn(y;x)dy =2πi2isin(πα)[Res(gn(z;x),i−x)+Res(gn(z;x),x−i)] =πsin(πα)n∑k=1(n−1k−1)(−αk)[(i−x)−α−k(2i)k−1−(x−i)−α−k(−2i)k−1] =−π((−1)−α+1)sin(πα)n∑k=1(n−1k−1)(−αk)(x−i)−α−k(−2i)k−1 =παe−iπα/2(i)1+αsin(πα/2)(ix+1)1+α2F1(1−n,1+α,2,2ix+1),

where we have used (15). Substituting this last expression into (16), and using Euler’s reflection formula,

 Γ(w)Γ(1−w)=πsin(πw),

at , we obtain

 (−Δ)α/2λn(x)=2ni2α−1Γ(1/2+α/2)√πΓ(1−α/2)iαΓ(α/2)Γ(1−α/2)(ix+1)1+α2F1(1−n,1+α,2,2ix+1). (24)

Then, (7) for follows from (24), using Legendre’s duplication formula , at . In order to finish the proof, we notice that the case where is a negative integer follows from the symmetry

 λn(x)=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯λ−n(x),n∈N. (25)

## 3 Fractional Laplacian of other families of functions

In this section, we compute the fractional Laplacian of the families of functions defined by (8)-(12). Let us obtain first the Fractional Laplacian of the complex Christov Functions:

###### Proposition 3.1.

Let be defined as in (8). Then,

 (−Δ)α/2μn(x)=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩Γ(1+α)(ix+1)1+α2F1(−n,1+α;1;2ix+1),n=0,1,2,…,−Γ(1+α)(−ix+1)1+α2F1(1+n,1+α;1;2−ix+1),n=−1,−2,−3,… (26)

Moreover, the expression for reduces to

 (−Δ)α/2μ0(x)=Γ(1+α)(ix+1)1+α.
###### Proof.

First, we note that can be expressed in terms of as follows:

 μn(x)=λn(x)1(ix+1)=12λn(x)(ix+1)−(ix−1)(ix+1)=λn(x)−λn+1(x)2;

therefore,

 (−Δ)α/2μn(x)=(−Δ)α/2λn(x)−(−Δ)α/2λn+1(x)2, (27)

so we just have to use (7) and simplify the resulting expression. First we assume that . Substituting (7) into this last equation, and using that , we get

 (−Δ)α/2μn(x) =Γ(1+α)(ix+1)1+α[n∑k=0(n+1k+1)(−1−αk)(2ix+1)k−n−1∑k=0(nk+1)(−1−αk)(2ix+1)k] =Γ(1+α)(ix+1)1+αn∑k=0(nk)(−1−αk)(2ix+1)k,

which is (26), for . On the other hand, when , again from (7),

 (−Δ)α/2μ0(x)=−(−Δ)α/2λ1(x)2=Γ(1+α)(ix+1)1+α2F1(0,1+α;2;2ix+1)=Γ(1+α)(ix+1)1+α,

and (26) also holds. Finally, when is a negative integer, we observe that

 μn(x)=−¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯μ−1−n(x), (28)

so, in that case, we use , which concludes the proof.

Before we compute the fractional Laplacian of the cosine-like and sine-like Higgins functions (9) and (10), we first express (7) as a polynomial on times a negative power of .

###### Lemma 3.1.

Let be defined as in (5), and . Then, (7) can be expanded as

 (−Δ)α/2λn(x)=−2|n|Γ(1+α)(isgn(n)x+1)|n|+α|n|−1∑l=0(|n|−1l)(isgn(n)x)|n|−1−l2F1(−l,1+α;2;2). (29)
###### Proof.

Assume first that . We express in (24) as a sum; then, replacing the index by in the sum, and expanding by Newton’s binomial formula, we get

 (−Δ)α/2λn(x) =−2Γ(1+α)(ix+1)1+αn−1∑k=0(nk+1)(−1−αk)(2ix+1)k (30) =−2nΓ(1+α)(ix+1)n+αn−1∑k=0k∑l=012k(nk)(kl)(−1−αn−1−k)(ix)l. (31)

We now interchange the order of the sums and replace the indices by and by ; then, after some rewriting, we get,

 (−Δ)α/2λn(x) =−2nΓ(1+α)(ix+1)n+αn−1∑l=0(ix)ln−1∑k=l12k(nk)(kl)(−1−αn−1−k) (32) =−2nΓ(1+α)(ix+1)n+αn−1∑l=0(ix)n−1−ll∑k=012n−1−k(nn−1−k)(n−1−kn−1−l)(−1−αk) (33) =−2nΓ(1+α)(ix+1)n+αn−1∑l=0(n−1l)(ix)n−1−l1l+1l∑k=0(l+1k+1)(−1−αk)2k, (34)

which yields (29). The case with a negative integer follows from (25).

###### Remark.

If we only consider in (29), we get the asymptotic behavior of :

 (−Δ)α/2λn(x)∼−2|n|Γ(1+α)(isgn(n)x)|n|−1(isgn(n)x+1)|n|+α∼−2|n|Γ(1+α)(isgn(n)x)1+α,x→±∞,

i.e., decays at infinity as for all , whereas

 λn(x)−1=(1−2ix+1)n−1∼−2nix+1∼−2nix,x→±∞,

i.e., decays at infinity as , for all . Therefore, the operator introduces an extra decay of .

As a consequence of Lemma 3.1, we also have the following result.

###### Proposition 3.2.

The fractional Laplacian of (9) and (10) is given respectively by

 (−Δ)α/2CH2n(x)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩0,n=0,−2Γ(1+α)(1+x2)(1+α)/2sin((1+α)arccot(x)−πα2),n=1,−2nΓ(1+α)(1+x2)(n+α)/2sin((n+α)arccot(x)−πα2)×⌊n−12⌋∑l=0(−1)l(n−12l)xn−1−2l2F1(−2l,1+α;2;2)+2nΓ(1+α)(1+x2)(n+α)/2cos((n+α)arccot(x)−πα2)×⌊n−22⌋∑l=0(−1)l(n−12l+1)xn−2−2l2F1(−2l−1,1+α;2;2),n=2,3,4,…, (35)

and

 (36)
###### Proof.

From the definition of (9) and (10), together with (25) and Euler’s formula , we have immediately, for , , ,

 CH2n(x)=R(λn(x)),SH2n(x)=I(λn+1(x)),

 CH2n(cot(s))=cos(2ns),SH2n+1(cot(s))=sin(2(n+1)s).

Since and are real, their fractional Laplacian is also real. We thus want to obtain expressions of and that avoid the use of complex numbers. To do this, we express all the complex numbers in (32) in their polar form: e.g., ,

 1(x−i)n+α=(x+i)n+α(1+x2)n+α=ei(n+α)arccot(x)(1+x2)(n+α)/2,

etc. Note that, since , with , we have chosen the definitions of the arctangent and the arccotangent such that in the last expression. Then, (32) becomes

 (−Δ)α/2λn(x)=−2nΓ(1+α)(1+x2)(n+α)/2n−1∑l=0(n−1l)xn−1−lei((n+α)arccot(x)−π(l+1+α)/2)2F1(−l,1+α;2;2). (37)

With respect to (35), the case is trivial, because ; otherwise, taking the real part of (37),

 (−Δ)α/2CH2n(x) =−2nΓ(1+α)(1+x2)(n+α)/2n−1∑l=0(n−1l)xn−1−