Numerical Approximation of the Fractional Laplacian on R Using Orthogonal Families

by   Jorge Cayama, et al.

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.


page 1

page 2

page 3

page 4


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

In this paper, we propose a novel pseudospectral method to approximate a...

Efficient Monte Carlo Method for Integral Fractional Laplacian in Multiple Dimensions

In this paper, we develop a Monte Carlo method for solving PDEs involvin...

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 ...

Metastable Speeds in the Fractional Allen-Cahn Equation

We study numerically the one-dimensional Allen-Cahn equation with the sp...

Explicit formulas for concatenations of arithmetic progressions

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

Generalised Hermite functions and their applications in Spectral Approximations

In 1939, G. Szegö first introduced a family of generalised Hermite polyn...

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:


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:


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


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


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]):


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


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


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


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:

  • Sine-like Higgins functions:

  • Cosine-like Christov functions:

  • Sine-like Christov functions:


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

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

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

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

in particular, when ,

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

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


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:


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).


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

Introducing this expression in (6), we get


with integrand

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


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

Figure 1: An example of integration contour for .

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

In this way, we obtain


Likewise, we have


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,



Finally, after some manipulations, we have

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

at , we obtain


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


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,


Moreover, the expression for reduces to


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



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

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

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


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


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


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


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


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

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

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




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

and this implies (see also [4]) that

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., ,

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


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