Rational spectral methods for PDEs involving fractional Laplacian in unbounded domains

TL;DR

This paper develops and analyzes rational spectral methods for solving PDEs with fractional Laplacians in unbounded domains, providing optimal error estimates and demonstrating superior numerical performance over existing Hermite-based approaches.

Contribution

The paper introduces rational basis spectral methods for fractional PDEs in unbounded domains, including explicit formulas, error analysis, and numerical validation.

Findings

Rational spectral methods achieve optimal error estimates.

The proposed methods outperform Hermite function approaches.

Numerical results confirm high accuracy and efficiency.

Abstract

Many PDEs involving fractional Laplacian are naturally set in unbounded domains with underlying solutions decay very slowly, subject to certain power laws. Their numerical solutions are under-explored. This paper aims at developing accurate spectral methods using rational basis (or modified mapped Gegenbauer functions) for such models in unbounded domains. The main building block of the spectral algorithms is the explicit representations for the Fourier transform and fractional Laplacian of the rational basis, derived from some useful integral identites related to modified Bessel functions. With these at our disposal, we can construct rational spectral-Galerkin and direct collocation schemes by pre-computing the associated fractional differentiation matrices. We obtain optimal error estimates of rational spectral approximation in the fractional Sobolev spaces, and analyze the optimal convergence of the proposed Galerkin scheme. We also provide ample numerical results to show that the rational method outperforms the Hermite function approach.

Tables

Table 1. Table 1: Rate of convergence using the generalized Hermite function, MMGFs with λ = 0 𝜆 0 \lambda=0 , and MMGFs with λ = 0.5 , 𝜆 0.5 \lambda=0.5, α = 1 𝛼 1 \alpha=1 and f ​ ( x ) = exp ⁡ ( − x 2 / 2 ) ​ ( 1 + x ) . 𝑓 𝑥 superscript 𝑥 2 2 1 𝑥 f(x)=\exp(-{x^{2}}/{2})(1+x).

	Hermite		MMGF $\lambda =0$		MMGF $\lambda =0.5$
$N$	$L^2$ error	Order	$L^2$ error	Order	$L^2$ error	Order
60	7.26e-4	–	2.32e-5	–	1.90e-5	–
80	5.72e-4	0.83	1.50e-5	-1.52	1.23e-5	1.50
100	4.77e-4	0.82	1.06e-5	-1.53	8.79e-6	1.51
120	4.09e-4	0.84	8.02e-6	-1.55	6.65e-6	1.53
140	3.58e-4	0.87	6.30e-6	-1.57	5.23e-6	1.56
160	3.17e-4	0.90	5.09e-6	-1.59	4.23e-6	1.59
180	2.84e-4	0.94	4.21e-6	-1.62	3.50e-6	1.62
200	2.57e-4	0.98	3.53e-6	-1.66	2.94e-6	1.66
220	2.32e-4	1.03	3.01e-6	-1.70	2.49e-6	1.71
240	2.12e-4	1.08	2.58e-6	-1.75	2.14e-6	1.77

Table 2. Table 2: Rate of convergence using the generalized Hermite function, MMGFs with λ = 0 𝜆 0 \lambda=0 , and MMGFs with λ = 0.5 . 𝜆 0.5 \lambda=0.5. α = 1 𝛼 1 \alpha=1 and f ​ ( x ) = 1 ( 1 + x 2 ) 2 𝑓 𝑥 1 superscript 1 superscript 𝑥 2 2 f\left(x\right)=\frac{1}{(1+x^{2})^{2}} .

	Hermite		MMGF $\lambda =0$		MMGF $\lambda =0.5$
$N$	$L^2$ error	Order	$L^2$ error	Order	$L^2$ error	Order
60	3.36e-3	–	2.36e-5	–	4.23e-5	–
80	2.61e-3	0.87	9.36e-6	-3.21	2.75e-5	1.49
100	2.17e-3	0.84	6.66e-5	-1.53	1.96e-5	1.51
120	1.86e-3	0.84	5.02e-6	-1.55	1.49e-5	1.53
140	1.63e-3	0.86	3.95e-6	-1.56	1.17e-5	1.56
160	1.45e-3	0.88	3.19e-6	-1.59	9.46e-6	1.59
180	1.30e-3	0.91	2.64e-6	-1.62	7.81e-6	1.62
200	1.18e-3	0.94	2.21e-6	-1.66	6.56e-6	1.66
220	1.08e-3	0.97	1.88e-6	-1.70	5.57e-6	1.71
240	9.84e-4	1.01	1.62e-6	-1.75	4.78e-6	1.77

Equations

$$\{\begin{array}{ll}(-\Delta {)}^{\alpha /2}u(x)+\rho u(x)=f(x),& x\in {\mathbb{R}}^d,\\ u(x)=0,& |x|\to \infty ,\end{array}$$

$$(-\Delta)^{\alpha/2}u(x):=C_{d,\alpha}{\rm p.v.}\int_{\mathbb{R}^{d}}\dfrac{u(x)-u(y)}{\lvert x-y\rvert^{d+\alpha}}dy\;\;\textmd{with}\;\;C_{d,\alpha}:=\dfrac{{\alpha}2^{\alpha-1}\Gamma\Big{(}\dfrac{\alpha+d}{2}\Big{)}}{\pi^{d/2}\Gamma\Big{(}\dfrac{2-\alpha}{2}\Big{)}}.$$

$$(-\Delta)^{\alpha/2}u(x):={\mathscr{F}}^{-1}\big{[}|\xi|^{\alpha}\mathscr{F}{\left[u\right]}({\xi})\big{]}(x).$$

$$J_{\mu }(x)=\sum _{m=0}^{\infty }\frac{(-1{)}^m}{m!\Gamma (m+\mu +1)}{\left(\frac{x}{2}\right)}^{2m+\mu }.$$

$$I_{\mu }(x)={\mathrm{i}}^{-\mu }{J}_{\mu }(\mathrm{i}x),K_{\mu }(x)=\frac{\pi }{2}\frac{I_{-\mu }(x)-I_{\mu }\left(x\right)}{\sin (\mu \pi )},$$

$${\int }_0^{\infty }{x}^{\lambda }{K}_{\mu }(ax)\cos (bx)dx$$

$$\displaystyle\;\;\;\;\times{}_{2}F_{1}\Big{(}\frac{\mu+\lambda+1}{2},\frac{1+\lambda-\mu}{2};\!\frac{1}{2};-\frac{b^{2}}{a^{2}}\Big{)},$$

$${\int }_0^{\infty }{x}^{\lambda }{K}_{\mu }(ax)\sin (bx)dx$$

$$\displaystyle\;\;\;\;\times{}_{2}F_{1}\Big{(}\frac{2+\mu+\lambda}{2},\frac{2+\lambda-\mu}{2};\frac{3}{2};-\frac{b^{2}}{a^{2}}\Big{)}.$$

$${}_2{F}_1\left(a,b;c;x\right)=\sum _{k=0}^{\infty }\frac{{\left(a\right)}_k{\left(b\right)}_k}{{\left(c\right)}_k}\frac{x^k}{k!},\mathrm{for}|x|<1,$$

$$(a{)}_0=1,(a{)}_k=a(a+1)...(a+k-1)=\frac{\Gamma (a+k)}{\Gamma (a)},k\in \mathbb{N}.$$

$$\underset{z\to 1^{-}}{\lim }\frac{{}_2{F}_1(a,b;a+b;z)}{-\ln (1-z)}=\frac{\Gamma (a+b)}{\Gamma (a)\Gamma (b)}.$$

$$\underset{z\to 1^{-}}{\lim }\frac{{}_2{F}_1(a,b;c;z)}{(1-z{)}^{c-a-b}}=\frac{\Gamma (c)\Gamma (a+b-c)}{\Gamma (a)\Gamma (b)}.$$

$$\frac{d^k}{d{x}^k}{}_2{F}_1(a,b;c;x)=\frac{(a{)}_k(b{)}_k}{(c{)}_k}{}_2{F}_1(a+k,b+k;c+k;x).$$

$$(2a-c-ax+bx){}_2{F}_1$$

$$+a(x-1){}_2{F}_1(a+1,b;c;x)=0.$$

$${}_{2}F_{1}(a,b;c;x)=(1-x)^{-a}\,{}_{2}F_{1}\Big{(}a,c-b;c;\frac{x}{x-1}\Big{)},$$

$${}_2{F}_1(-n,b;c;x)=\frac{(c-b{)}_n}{(c{)}_n}{}_2{F}_1(-n,b;b-c-n+1;1-x).$$

$$\displaystyle\int_{0}^{\infty}\cos(\mu x){}_{2}F_{1}\Big{(}a,b;\frac{1}{2};-c^{2}x^{2}\Big{)}dx=2^{-a-b+1}\pi c^{-a-b}\mu^{a+b-1}\frac{K_{a-b}({\mu}/{c})}{\Gamma(a)\Gamma(b)},$$

$$\displaystyle\int_{0}^{\infty}x\sin(\mu x){}_{2}F_{1}\Big{(}a,b;\frac{3}{2};-c^{2}x^{2}\Big{)}dx=2^{-a-b+1}\pi c^{-a-b}\mu^{a+b-2}\frac{K_{a-b}({\mu}/{c})}{\Gamma(a)\Gamma(b)}.$$

$$\begin{array}{rl}& n{C}_n^{\lambda }(t)=2t(n+\lambda -1)C_{n-1}^{\lambda }(t)-(n+2\lambda -2)C_{n-2}^{\lambda }(t),n\ge 2,\\ & C_0^{\lambda }(t)=1,C_n^{\lambda }(t)=2\lambda t.\end{array}$$

$${\int }_{-1}^1{C}_n^{\lambda }(t)C_m^{\lambda }(t){\omega }_{\lambda }(t)dt={\gamma }_n^{\lambda }{\delta }_{nm},{\gamma }_n^{\lambda }=\frac{\pi 2^{1-2\lambda }\Gamma (n+2\lambda )}{n!\left(n+\lambda \right){\Gamma }^2(\lambda )}.$$

$$\begin{split}&C_{2n}^{\lambda}(t)=\frac{(-1)^{n}}{(\lambda+n)B(\lambda,n+1)}\,{}_{2}F_{1}\Big{(}\!-n,n+\lambda;\frac{1}{2};t^{2}\Big{)};\\[4.0pt] &C_{2n+1}^{\lambda}(t)=\frac{(-1)^{n}2t}{B(\lambda,n+1)}\,{}_{2}F_{1}\Big{(}\!-n,n+\lambda+1;\frac{3}{2};t^{2}\Big{)},\end{split}$$

$$B(x,y)=\frac{\Gamma (x)\Gamma (y)}{\Gamma (x+y)}.$$

$$\begin{split}&C_{2n}^{\lambda}(t)=a_{n}^{\lambda}\,\,{}_{2}F_{1}\Big{(}\!-n,n+\lambda;\lambda+\frac{1}{2};1-t^{2}\Big{)},\\[4.0pt] &C_{2n+1}^{\lambda}(t)=b_{n}^{\lambda}\,t\,{}_{2}F_{1}\Big{(}\!-n,n+\lambda+1;\lambda+\frac{1}{2};1-t^{2}\Big{)},\end{split}$$

$$a_n^{\lambda }=\frac{{\left(\lambda \right)}_n}{{\left(1\right)}_n}\frac{{\left(\lambda +\frac{1}{2}\right)}_n}{{\left(\frac{1}{2}\right)}_n},b_n^{\lambda }=\frac{2\lambda {\left(\lambda +1\right)}_n}{{\left(1\right)}_n}\frac{{\left(\lambda +\frac{1}{2}\right)}_n}{{\left(\frac{3}{2}\right)}_n}.$$

$$T_n\left(t\right)=\frac{n}{2}\underset{\lambda \to 0}{\lim }\frac{C_n^{\lambda }(t)}{\lambda },n\ge 1.$$

$$\begin{split}&T_{2n}(t)={}_{2}F_{1}\Big{(}\!\!-n,n;\frac{1}{2};1-t^{2}\Big{)};\quad T_{2n+1}(t)=t\,{}_{2}F_{1}\Big{(}\!\!-n,n+1;\frac{1}{2};1-t^{2}\Big{)}.\end{split}$$

$$x=\frac{t}{\sqrt{1-t^2}}\mathrm{or}t=\frac{x}{\sqrt{1+x^2}},t\in I,x\in \mathbb{R}.$$

$$1-t^2=\frac{1}{1+x^2},\frac{dx}{dt}=\frac{1}{(1-t^2{)}^{3/2}}.$$

$$R_{n}^{\lambda}\left(x\right):=\big{(}1+x^{2}\big{)}^{-\frac{\lambda+1}{2}}C_{n}^{\lambda}\Big{(}\frac{x}{\sqrt{1+x^{2}}}\Big{)},\quad x\in{\mathbb{R}},$$

$$R_n^{\lambda }\left(x\right)=S(t)C_n^{\lambda }(t),S(t):=\sqrt{{\omega }_{\lambda }(t)\frac{dx}{dt}}=(1-t^2{)}^{\frac{\lambda +1}{2}},$$

Full text

\slugger

mmsxxxxxxxx–x

1 Introduction

Diffusion is a ubiquitous physical process, typically modeled by partial differential equations (PDEs) with usual Laplacian operators. Although they can describe the anisotropy of diffusion, many systems in science, economics and engineering exhibit anomalous diffusion, which can be more accurately and realistically modelled by PDEs with fractional Laplacian operators [4, 5, 13]. In the past decade, tremendous research attention has been paid to the analysis and numerical studies of fractional PDEs. The finite difference method and the finite element method are two widely studied methods in this direction (see, e.g., [18, 19, 20, 9, 40, 3, 43, 36, 44, 45, 41] and references therein). Most of efforts are devoted to dealing with the nonlocal nature or singularities of the fractional operators. Another powerful approach is the spectral method, which is more suitable for the non-local feature of the fractional operators (see, e.g., [47, 11, 10, 21, 25, 34, 35, 46, 37]). However, most of these works are for fractional problems in bounded domains. In particular, we refer to Bonito et al [6] for an up-to-date review of the various numerical methods for fractional diffusion based on different formulations of the fractional Laplacian.

It is known that many physically motivated fractional diffusive problems are naturally set in unbounded domains, but their investigation is still under-explored. For usual PDEs in unbounded domains, several approaches have been widely used in practice (see, e.g., [8, 32] and the original references cited therein). The first is direct domain truncation that works well for problems with rapidly decaying solutions, but is not feasible for fractional problems as the underlying solutions usually decay slowly, subject to certain power laws at infinity. On the other hand, the naive truncation introduces nonphysical singularities at the interface where the unbounded domain is terminated. The second is to design a suitable transparent boundary condition or artificial sponge layer, but this appears highly nontrivial for the fractional Laplacian. The third is the use of orthogonal functions in unbounded domains, which has been successfully applied to many usual PDEs (see, e.g., [38, 8, 14, 27, 32, 31]). Very recently, spectral methods for fractional PDEs on the half line are proposed by [22, 25] – using the generalized Laguerre functions as bases – extending the idea of [47]. A two-domain spectral approximations by Laguerre functions is developed in [12] for tempered fractional PDEs on the whole line. Mao and Shen [28] proposed both spectral-Galerkin and collocation methods using Hermite functions for fractional PDEs in unbounded domains. However, the collocation method therein relies on an equivalent formulation in frequency space by the Fourier transforms, and performs collocation methods to the equivalent formulation that involve forward/backward Hermite transforms. Tang, Yuan and Zhou [39] developed direct Hermite collocation methods with explicit formulations for the differentiation matrices, which is therefore more robust for nonlinear problems. Lastly, spectral approximation using non-classical orthogonal functions in unbounded domains – image of classical Jacobi polynomials through a suitable mapping, has proven to be more viable for usual PDEs with solutions decaying algebraically (see, e.g., [7, 8, 16, 17, 30]), compared with approximation by Hermite/Laguerre functions. As such, the rational basis (or mapped Jacobi functions) should be more desirable for PDEs with fractional Laplacian, due to the slow decaying solution with long tails subject to certain power law. However, to the best of our knowledge, there is essentially no work available along this line. Moreover, the extension of the mapping technique to the fractional setting is far from trivial, as we elaborate on below.

In this paper, we intend to fill in this gap, and design rational spectral methods for a class of PDEs with fractional Laplacian in $\mathbb{R}^{d}$. To fix the idea, we consider the model equation:

[TABLE]

for $\alpha\in(0,2),$ where the fractional Laplace operator is defined as in [23]:

[TABLE]

Here, p.v. stands for the Cauchy principal value, and $C_{d,\alpha}$ is a normalization constant. Equivalently, the fractional Laplacian can be defined as a pseudo-differential operator via the Fourier transform:

[TABLE]

For any expansion-based method, a critical issue is how to accurately evaluate the point-wise value of the fractional Laplacian performing upon the basis. For example, the key to the Hermite spectral method in [28] is the use of the attractive property that the Hermtie functions are the eigenfunctions of the Fourier transform, so the algorithm is largely implemented in the frequency $\xi$-space. In contrast, some analytically perspicuous formulas of $(-\Delta)^{\alpha/2}$ on the Hermite functions were derived in [39], which led themselves to the construction of efficient collocation algorithms in the physical $x$-space. In the spirit of [39], we search for the analytic formulas for computing the fractional Laplacian of the rational basis functions – the modified mapped Gegenbauer functions (MMGFs), orthogonal with respect to a uniform weight. Although the formulas (see Theorem 4) are not as compact as those for the Hermite functions, we can accurately compute the fractional Laplacian of the rational basis up to the degree $\thicksim 10^{3}$ by using e.g., Maple or Mathematica. Moreover, with these analytic tools at our disposal, we are able to study their asymptotic behaviors and dependence of the parameters so that the basis can be tailored to the decay rate of underlying solution. We propose and analyze a spectral-Galerkin scheme, and obtain optimal estimates (see Theorem 9). We also implement a direct collocation scheme based on the associated fractional differentiation matrices with the aid of the aforementioned explicit formulas. However, its error analysis appears very challenging and largely open. This is mostly for the reason that the fractional Laplacian takes the rational basis to a class of functions of completely different nature, as opposite to the usual Laplacian. In the multi-dimensional case, we implement the collocation schemes in the frequency space (cf. [28]), which relies on the approximability of spectral expansions to ${\mathscr{F}}[f](\xi)/(|\xi|^{\alpha}+\rho).$ We show that the rational approach outperforms the Hermite method in accuracy. In fact, it is common that the Fourier transform of a functions decays much slower than the function itself, so the rational basis is more desirable in this context.

The rest of the paper is organized as follows. In section 2, we collect some useful properties of the Bessel functions and Gegenbauer polynomials. In section 3, we present the main formulas for computing the fractional Laplacian of the modified rational functions, and study the asymptotic properties. In section 4, we derive optimal error estimates of the approximation by the modified rational functions in fractional Sobolev spaces. We propose and analyse spectral-Galerkin methods using modified rational basis functions in section 5. Then we implement the collocation methods in both one dimension and multiple dimensions in section 6. The final section is for some concluding remarks.

2 Preliminaries

In this section, we make necessary preparations for the algorithm development and analysis in the forthcoming sections. More precisely, we review some relevant properties of the hypergeometric functions, Gegenbauer polynomials, Bessel functions and their interwoven relations.

2.1 Bessel functions

Recall that the Bessel function of the first kind of real order $\mu$ has the series expansion (cf. [29]):

[TABLE]

The modified Bessel functions of the first and second kinds are defined by

[TABLE]

where ${\rm i}=\sqrt{-1}$ is the complex unit. For the modified Bessel functions of the second kind $K_{\mu}(x)$, we have the following important integral identities (see [15, P. 738]): for $-\lambda\pm\mu<1$ and $a,b>0,$

[TABLE]

and for $-\lambda\pm\mu<2$ and $a,b>0,$

[TABLE]

Here, $\Gamma(\cdot)$ is the usual Gamma function, and ${}_{2}F_{1}$ is the hypergeometric function defined in (8) below.

2.2 Hypergeometric functions

For any real $a,b,c$ with $c\not=0,-1,-2,\cdots,$ the hypergeometric function is a power series defined by

[TABLE]

and by analytic continuation elsewhere (cf. [15, P. 1014] or [2, Ch. 2]). Here $(a)_{k}$ is the rising Pochhammer symbol, i.e.,

[TABLE]

It is known that the series ${}_{2}F_{1}(a,b;c;x)$ is absolutely convergent for all $|x|<1.$ Moreover, (i) if $c-a-b>0,$ the series ${}_{2}F_{1}(a,b;c;x)$ is absolutely convergent at $x=\pm 1;$ (ii) if $-1<c-a-b\leq 0,$ the series ${}_{2}F_{1}(a,b;c;x)$ is conditionally convergent at $x=-1,$ but it is divergent at $x=1;$ (iii) if $c-a-b\leq-1$, it diverges at $x=\pm 1.$ Its divergent behaviour at $x=1$ can be characterised as follows (cf. [2, Ch. 2]).

• •

If $c=a+b,$ then

[TABLE]

• •

If $c<a+b,$ then

[TABLE]

From the definition (8), we can easily obtain

[TABLE]

According to [15, P. 1019], there holds

[TABLE]

We also recall the property of hypergeometric functions related to transformations of variable (cf. [29, P. 390]):

[TABLE]

and the Pfaff’s formula on the linear transformation (cf. [2, (2.3.14)]): for integer $n\geq 0,$

[TABLE]

Like (2.1)-(2.1), the following integral formulas (cf. [15, P. 825]) play a very important role in the algorithm development: for real $\mu>0$ and real $a,b,c>0,$

[TABLE]

and for $a,b>{1}/{2},$

[TABLE]

2.3 Gegenbauer polynomials

Gegenbauer polynomials, denoted by $C_{n}^{\lambda}(t),$   $t\in I:=(-1,1)$ and $\lambda>-1/2$, generalise Legendre and Chebyshev polynomials. They are defined by the three-term recurrence relation (cf. [15, P. 1000]):

[TABLE]

They are orthogonal with respect to the weight function $\omega_{\lambda}(t)=(1-t^{2})^{\lambda-1/2}$:

[TABLE]

The Gegenbauer polynomials can be defined by the hypergeometric functions ([15, P. 1000]):

[TABLE]

where $B(\cdot,\cdot)$ is the Beta function satisfying (cf. [15, P. 918]):

[TABLE]

Using the linear transformation (13) and (18)-(19), we have

[TABLE]

where

[TABLE]

Remark 2.1**.**

Note that when $\lambda=0$, we understand the classical Chebyshev polynomials in the sense of

[TABLE]

Correspondingly, it follows from (20) that

[TABLE]

Here, we still denote $T_{n}(t):=C_{n}^{0}(t)$.

3 Fractional Laplacian of the modified mapped Gegenbauer functions

In this section, we introduce the rational basis functions through the Gegenbauer polynomials with a singular mapping. For convenience, we term the resulting mapped basis as modified mapped Gegenbauer functions (MMGFs), which are different from the usual mapped Gegenbauer functions by absorbing the weight function in the basis. We also present the explicit formulas for the evaluation of their fractional Laplacian, which plays an essential role in the spectral algorithms.

3.1 The mapping and MMGFs

Consider the one-to-one mapping between $t\in I=(-1,1)$ and $x\in\mathbb{R}=(-\infty,\infty)$ of the form:

[TABLE]

It is clear that

[TABLE]

Definition 1**.**

For $\lambda>-1/2,$ let $C_{n}^{\lambda}(t),t\in I=(-1,1),$ be the Gegenbauer polynomial of degree $n$ as defined in (16). We define the modified mapped Gegenbauer functions (MMGFs) as

[TABLE]

or equivalently,

[TABLE]

where $x,t$ are associated with the mapping (1).

One verifies readily from (17) and (3)-(4) that

[TABLE]

Thanks to (16) and (3), the MMGFs satisfy the three-term recurrence relation:

[TABLE]

Moreover, we can show that

[TABLE]

It is clear that by (8), (20) and (2), we have

[TABLE]

and

[TABLE]

It is seen that the MMGFs are expressed in terms of

[TABLE]

3.2 Formulas for computing fractional Laplacian of MMGFs

In view of (3.1)-(9), we first compute the fractional Laplacian of the simple functions in (9).

Theorem 2**.**

For real $s>0,$ we have that for any $\gamma>0,$

[TABLE]

and for any $\gamma>1/2,$

[TABLE]

where the factor

[TABLE]

Proof.

Recall the formula (cf. [29, 15.4.6]):

[TABLE]

Note that (13) also holds for $z<1,$ with the analytic extension by the transformation formula (12) (see [15, 9.130]). Thus, we have

[TABLE]

Then using (14) with $\mu=\xi,\alpha=\gamma$ and $\beta=1/2,$ we obtain that for $\xi>0,$

[TABLE]

Note that for $\xi<0,$ we have $\hat{v}(\xi)=\hat{v}(-\xi)$. Thus, from the definition (3) and (15), we obtain

[TABLE]

Then using the formula (2.1) with $\lambda=2s+\gamma-1/2,\mu=\gamma-1/2$ and $b=x,$ we find

[TABLE]

Hence, we derive (10).

The formula (11) can be derived in a similar fashion. Like (15), we obtain from (15) with $\mu=\xi,\alpha=\gamma$ and $\beta=3/2$ that for $\gamma>1/2$ and $\xi>0,$

[TABLE]

Note that in this case, ${\mathscr{F}}[xv](-\xi)=-{\mathscr{F}}[xv](\xi).$ Similar to (16), we find

[TABLE]

Thus, we derive from (2.1) with $\lambda=2s+\gamma-1/2,\mu=\gamma-3/2$ and $b=x$ that

[TABLE]

Finally, the formula (11) follows from the property: $\Gamma(z+1)=z\Gamma(z).$ ∎

Using the transformation formula (12), we can represent the formulas in Theorem 2 in the terms of the hypergeometric function defined through the series in (8). Note that the former is more convenient for computation, while the latter is more suitable for analysis.

Corollary 3**.**

For real $s>0,$ we have that for any $\gamma>0,$

[TABLE]

and for any $\gamma>1/2,$

[TABLE]

*where $A_{s}^{\gamma}$ is defined as in (12). *

Remark 3.1**.**

It is seen that if $s$ is a positive integer, then the hypergeometric functions in (21) and (22) become finite series. We can directly verify by using (8) that

[TABLE]

for $s=1,2,\cdots,$ and $|x|\to\infty.$ However, for non-integer $s>0,$ the hypergeometric functions may diverge as $|x|\to\infty.$ Indeed, we find from (9) and (10) that

• (i)

if $\gamma=1/2,$ then

[TABLE]

• (ii)

if $\gamma>1/2,$ then

[TABLE]

• (iii)

if $0<\gamma<1/2,$ it has the same behaviour as in (23).

Similarly, we can analyse the behaviour at infinity for (22) for three cases: (i) $\gamma=3/2;$ (ii) $\gamma>3/2$ and (iii) $1/2<\gamma<3/2.$

Remark 3.2**.**

In a distinctive difference with the integer case, we see that the decay rate in the fractional case in (24) is independent of $\gamma,$ if $\gamma>1/2.$

With the above preparations, we can now derive the explicit representation of the fractional Laplacian of $\{R_{n}^{\lambda}\}$.

Theorem 4**.**

For real $s>0$ and $\lambda>-1/2,$ the fractional Laplacian of the MMGFs can be represented by

[TABLE]

and

[TABLE]

where the constants $a_{n}^{\lambda},b_{n}^{\lambda}$ and $A_{s}^{\gamma}$ are the same as in (21) and (12).

Proof.

By (3.1), we have

[TABLE]

so substituting (10) with $\gamma=k+\frac{\lambda+1}{2}$ into the above leads to (26).

Similarly, we derive from (3.1) that

[TABLE]

so substituting (11) with $\gamma=k+\frac{\lambda}{2}+1$ into the above leads to (26). ∎

In view of the asymptotic results in Remark 3.1, we can analyse the decay rate of fractional Laplacian of the basis. Indeed, by virtual of (23)-(25), we obtain from (28)-(29) that (i) if $-1/2<\lambda<0,$ then

[TABLE]

(ii) if $\lambda=0,$ we have

[TABLE]

(iii) if $\lambda>0,$ we have

[TABLE]

Similar results are available for $(-\Delta)^{s}R_{2n+1}^{\lambda}(x).$ It is noteworthy from (3) and the above that the fractional Laplacian on the basis does not always lead to the gain in decay rate of $1/(1+x^{2})^{s}.$

Remark 3.3**.**

It is important to point out that the involved hypergeometric functions in (26) and (27) can be evaluated recursively by using (2.2). Denote

[TABLE]

Then by (2.2), we have

[TABLE]

for $k\geq 1.$ Similarly, we can efficiently compute the hyergeometric functions in (27).

Remark 3.4**.**

To enhance the resolution of the basis, one can also introduce a scaling parameter $\mu>0$ (cf. [31]). More precisely, the algebraic mapping in (1) turns to

[TABLE]

The corresponding modified rational function can be defined as

[TABLE]

*In fact, it is straightforward to extend the previous properties and formulas to the scaled basis. For simplicity, we omit the details. *

4 Estimates of MMGF approximation in fractional Sobolev spaces

In this section, we analyse the approximation property by the modified rational basis functions in fractional Sobolev spaces. We remark that there exist very limited results on the Legendre or Chebyshev rational approximations (see [17, 42, 33]). However, most of them are suboptimal. Here, we derive the optimal estimates in more general settings.

4.1 Fractional Sobolev spaces

For real $r\geq 0,$ we define the fractional Sobolev space (as in [24, P. 30] and [1, Ch. 1]):

[TABLE]

equipped with the norm

[TABLE]

We have the following space interpolation property (cf. [1, Ch. 1]).

Lemma 5**.**

For real $r_{0},r_{1}\geq 0,$ let $r=(1-\theta)r_{0}+\theta r_{1}$ with $\theta\in[0,1]$. Then for any $u\in H^{r_{0}}(\mathbb{R})\cap H^{r_{1}}(\mathbb{R}),$ we have

[TABLE]

Proof.

For the readers’ reference, we sketch the derivation of this interpolation property. It is clear that by (35),

[TABLE]

Using the Hölder’s inequality with $p=1/(1-\theta)$ and $q=1/\theta$, we obtain

[TABLE]

This completes the proof. ∎

4.2 Error estimate of orthogonal projections

Define the approximation space

[TABLE]

Consider the $L^{2}$-orthogonal projection $\pi_{N}^{\lambda}:L^{2}(\mathbb{R})\to V_{N}^{\lambda},$ i.e.,

[TABLE]

For notational convenience, we introduce the pairs of functions associated with the mapping (1):

[TABLE]

In what follows, the notation with or without “  $\breve{}$  ” has the same meaning.

In order to describe the approximation errors, we introduce new differential operators as follows

[TABLE]

where $a(x)={dx}/{dt}=(1+x^{2})^{\frac{3}{2}}.$ Correspondingly, we define the Sobolev space

[TABLE]

equipped with the norm and semi-norm

[TABLE]

Theorem 6**.**

For any $u\in H^{s}(\mathbb{R})\cap{\mathbb{B}}^{m}_{\lambda}(\mathbb{R})$ with integer $1\leq m\leq N+1,s\in(0,1),$ and $\lambda>-1/2,$ we have

[TABLE]

where $c$ is a positive constant independent of $N$ and $u.$

Proof.

We take two steps to carry out the proof.

Step 1: We first prove that

[TABLE]

For this purpose, we study the close relation between $\pi_{N}^{\lambda}$ and the orthogonal projection $\varPi_{N}^{\lambda}:L^{2}_{\omega_{\lambda}}(I)\to{\mathcal{P}}_{N},$ such that for any $\varPhi\in L^{2}_{\omega_{\lambda}}(I),$

[TABLE]

Recall the Gegenbauer polynomial approximation result (cf. [31, Thm 3.55]): if $\varPhi^{(l)}(t)\in L^{2}_{\omega_{\lambda+l}}(I)$ for $0\leq l\leq m,$ we have

[TABLE]

where the weight function $\omega_{a}(t)=(1-t^{2})^{a-1/2}.$

From (38) and (39), we find

[TABLE]

Therefore, we have

[TABLE]

As a result, there holds

[TABLE]

Thus, using (46) with $l=0$, we derive from (40)-(42) that

[TABLE]

Like (48), we can show

[TABLE]

Similar to (49), we derive from (46) with $l=0,1$ that

[TABLE]

Then the estimate (44) is a direct consequence of (50) and the above.

Step 2: It is evident that the result (44) implies

[TABLE]

and

[TABLE]

Then using the interpolation inequality in Lemma 5 with $r_{0}=0,r_{1}=1$ and $\theta=s,$ we obtain from (52)-(53) that

[TABLE]

This completes the proof. ∎

In the error analysis, it is necessary to consider the $H^{s}$-orthogonal projection. Define the bilinear form on $H^{s}({\mathbb{R}}):$

[TABLE]

Consider the orthogonal projection $\pi_{N,\lambda}^{s}:H^{s}(\mathbb{R})\to V_{N}^{\lambda}$ such that

[TABLE]

Then by the projection theorem, we have

[TABLE]

Taking $\phi=\pi_{N}^{\lambda}u,$ we immediately derive the following estimate.

Theorem 7**.**

For any $u\in H^{s}(\mathbb{R})\cap{\mathbb{B}}^{m}_{\lambda}(\mathbb{R})$ with integer $1\leq m\leq N+1,s\in(0,1),$ and $\lambda>-1/2,$ we have

[TABLE]

where $c$ is a positive constant independent of $N$ and $u.$

4.3 Error estimate of interpolation

Let $\{t_{j}^{\lambda},\rho_{j}^{\lambda}\}_{j=0}^{N}$ be the Gegenbauer-Gauss quadrature nodes and weights, where $\{t_{j}^{\lambda}\}$ are zeros of the Gegenbauer polynomial $C_{N+1}^{\lambda}(t).$ Define the mapped nodes and weights:

[TABLE]

Then by the exactness of the Gagenbauer-Gauss quadrature (cf. [31, Ch 3]), we have

[TABLE]

which, together with (37), implies the exactness of quadrature

[TABLE]

We now introduce the interpolation operator $I_{N}^{\lambda}u\,:\,C(\mathbb{R})\to V_{N}^{\lambda}$ such that

[TABLE]

As a consequence of (5) and (60), we have

[TABLE]

We have the following interpolation approximation result.

Theorem 8**.**

For any $u\in H^{s}(\mathbb{R})\cap{\mathbb{B}}^{m}_{\lambda}(\mathbb{R})$ with integer $1\leq m\leq N+1,s\in(0,1),$ and $\lambda>-1/2,$ we have

[TABLE]

where $c$ is a positive constant independent of $N$ and $u.$

Proof.

Recall the Gegenbauer-Gauss interpolation $I_{N}^{G}:C(-1,1)\to{\mathcal{P}}_{N},$ such that

[TABLE]

Then we have the expansion

[TABLE]

One verifies from (3), (39), (59) and (61)-(63) that

[TABLE]

Thus,

[TABLE]

where with a little abuse of notation, we still use the same notation as in (48). Following the lines as in (49)-(4.2), we can show that

[TABLE]

and

[TABLE]

According to [31, Thm. 3.41] on the Gegenbauer-Gauss interpolation error estimate, we have

[TABLE]

Then by the interpolation inequality in Lemma 5, we obtain from the above that

[TABLE]

This completes the proof. ∎

5 Modified rational spectral-Galerkin methods

In this section, we consider the spectral-Galerkin approximation to a model equation, and conduct the error analysis. We also present some numerical results to show our proposed method outperforms the Hermite approximations in [28, 39].

5.1 The scheme and its convergence

Consider the model equation

[TABLE]

for $\alpha\in(0,2),$ where $f\in L^{2}(\mathbb{R})$ and the constant $\rho>0.$

For notational convenience, let $s=\alpha/2.$ A weak form of (68) is to find $u\in H^{s}(\mathbb{R})$ such that

[TABLE]

The spectral-Galerkin scheme is to find $u_{N}\in V_{N}^{\lambda}$ (defined in (37)) such that

[TABLE]

Denote $e_{N}=u_{N}-\pi_{N,\lambda}^{s}u$ and $\tilde{e}_{N}=u-\pi_{N,\lambda}^{s}u.$ By a standard analysis, we find that for any $v_{N}\in V_{N}^{\lambda},$

[TABLE]

Taking $v_{N}=e_{N}$ and using the Cauchy-Schwarz inequality, we obtain

[TABLE]

Thus, by the triangle inequality, we derive

[TABLE]

In summary, we have the following convergence result.

Theorem 9**.**

For any $u\in H^{s}(\mathbb{R})\cap{\mathbb{B}}^{m}_{\lambda}(\mathbb{R})$ and $f\in{\mathbb{B}}^{k}_{\lambda}(\mathbb{R})$ with integer $1\leq m,k\leq N+1,s=\alpha/2\in(0,1),$ and $\lambda>-1/2,$ we have

[TABLE]

where $c$ is a positive constant independent of $N$ and $u,f.$

5.2 Numerical examples

We now present several examples to show the convergence behaviour of the above spectral Galerkin method. In all tests, we report the numerical errors in the $L^{2}$-norm, and set $\rho=1$. Here, we only consider the cases with $\lambda=0$ and $\lambda=0.5$, which correspond to the modified mapped Chebyshev rational functions and modified mapped Legendre functions, respectively.

Example 1: Exponential decay $f(x).$  We first consider equation (68) with $f(x)=\exp(-{x^{2}}/{2})(1+x).$ Since the closed-form exact solution is not available, we take the numerical solution with $N=600$ as the reference solution. The convergence results with MMGFs for $\alpha=0.4,1,1.6$ are presented in Figure 1 (middle and right). In the left plot, we have also presented the convergence results for the Hermite function approach in [28]. It is clearly seen that the MMGFs approach outperforms the Hermite approximations for all cases, namely, the MMGFs approach admits much higher convergence rates. This can also be seen form Table 1, where we have presented the order of convergence for both approaches.

Example 2: Algebraic decay $f(x).$  We next consider equation (68) with an algebraic decay source term: $f\left(x\right)=\frac{1}{(1+x^{2})^{2}}.$ The plots of the error decay for both Hermite functions and MMGFs are in Figure 2. Indeed, we observe the convergence behaviour similar to the previous example – the MMGF approach has a much better performance. The comparison in Table 2 also shows that the proposed approach converges much faster than the Hermite method.

To better understand the solution behaviours, we present in Figure 3 the asymptotic behavior of the “exact” solutions as $|x|\gg 1$ for the above two examples. We see that, for both examples with very different decay of $f(x)$, the solution $u(x)$ decays at the same rate: $\lvert x\rvert^{-\alpha-1}$. This testifies the solution decays at a rate of a power law, as opposite to the usual Laplacian. This also explains the reason why MMGFs have a better performance than the Hermite functions.

6 Modified rational spectral-collocation methods

With the formulas in Theorem 4 at our disposal, we can directly generate the spectral fractional differentiation matrices and develop the direct collocation methods, similar to the Hermite collocation methods in [39]. However, it seems nontrivial and largely open to analyse its convergence. In fact, we can also implement the collocation method in the Fourier transformed domain which turns to be more a natural way to extend the method to multiple dimensions.

6.1 Fractional differentiation matrices

Let $\{x_{j}^{\lambda},\omega_{j}^{\lambda}\}_{j=0}^{N}$ be the mapped Gegenbauer-Gauss collocation points and weights as given in (59). For any $u_{N}\in V_{N}^{\lambda},$ we write

[TABLE]

where $u_{j}=u_{N}(x_{j}^{\lambda}).$ Note that the corresponding Lagrange basis function $\{l_{j}(x)\}_{j=1}^{N}$ can be expressed as

[TABLE]

Consequently, we can easily derive the associated differential matrix $\mathcal{D}^{L,\alpha,\lambda}$ with Lagrange type bases

[TABLE]

where $(-\Delta)^{\alpha/2}R_{j}^{\lambda}(x_{i}^{\lambda})$ can be computed via (26) and (27).

6.2 Numerical examples

We now present several numerical examples to show the performance of the spectral collocation method based on MMGFs. Notice that the collocation method is more practical for problems with variable coefficients and nonlinear problems. Also, we shall carry out comparisons with the Hermite collocation method in [39].

6.2.1 A multi-term fractional model

We first consider the following multi-term fractional Laplacian equation:

[TABLE]

Here we set $J=4$ and

[TABLE]

Numerical results with two different souce terms are presented in Figure 4. It can be seen that, similar to the Galerkin methods, the MMGF approach has a much better performance than the Hermite function approach in all cases.

6.2.2 Fractional model with variable coefficients

We next consider the following problem

[TABLE]

where $r(x)=1+\exp(-x^{2})$ and $f(x)=\frac{1}{(1+x^{2})^{1.2}}$. The convergence results for $\alpha=0.4,1,1.6$ are provided in Figure 5 for both approaches. Again, the MMGFs spectral collocation method outperforms the Hermite collocation method.

6.2.3 An eigenvalue problem

Finally, we consider the following eigenvalue problem as in [39]:

[TABLE]

Notice that exact eigenvalues for the case of $\alpha=1$ are available in [26]. For this example, we shall compute the first three eigenvalues by the MMGFs spectral collocation method and the Hermite collocation methods for comparison. The numerical results are given in Figure 6, which shows that the MMGF collocation method is more accurate than the Hermite collocation method.

6.3 Spectral-collocation methods in multiple dimensions

To this end, we propose the modified rational collocation methods based on a formulation in the Fourier transformed domain in multiple dimensions and show that it is more accurate than the Hermite spectral collocation methods in [28].

To fix the idea, we consider the $d$-dimensional model problem:

[TABLE]

where we denote $x=(x_{1},\cdots,x_{d})$ and $|x|=\sqrt{x^{t}x}.$ It is known that in the Fourier transformed domain, it can be expressed as

[TABLE]

where $\hat{u},\hat{f}$ are the Fourier transform of $u,f,$ respectively. Thus, we have

[TABLE]

That is, the Fourier transform of solution $u$ can be expressed explicitly as above. This motivates the construction of the collocation method in the frequency space. To describe the algorithm, we denote

[TABLE]

and define the tensorial grids and tensorial MMGFs as

[TABLE]

As the first step, we approximate $f(x)$ by the multidimensional interpolation:

[TABLE]

where the coefficients $\{\tilde{f}_{n}\}_{n\in\Upsilon_{\!N}}$ can be computed from the samples $\{f(x_{j}^{\lambda})\}_{j\in\Upsilon_{\!N}}$ by the tensorial version of the quadrature (60). Then we have the approximation:

[TABLE]

where ${\mathscr{F}}[R_{n_{i}}^{\lambda}](\xi_{i})$ can be computed by the formulas in Theorem 10 below.

Then the direct collocation approximation of $u(x)$ in the frequency space is given by

[TABLE]

where $\{\xi_{j}^{\lambda}\}$ are the tensorial grids as in (80). With these samples, we can write the final approximation $u_{N}(x)$ by taking Fourier inverse transform of $\hat{u}_{N}(\xi)$ as follows:

[TABLE]

where the coefficients $\{\tilde{u}_{n}\}_{n\in\Upsilon_{\!N}}$ can be computed from $\{\widehat{u}_{N}(\xi_{j}^{\lambda})\}_{j\in\Upsilon_{\!N}}$ in (83) by the quadrature formula (cf. (60)) as before. Here, the inverse Fourier transforms can be computed by the formulas in Theorem 10 and Remark 6.1 below.

Like Theorem 4, we have the following formulas for computing the Fourier transform of the MMGFs.

Theorem 10**.**

For real $\lambda>-1/2,$ the Fourier transform of the MMGFs can be computed by

[TABLE]

and

[TABLE]

where the constants $a_{n}^{\lambda},b_{n}^{\lambda}$ are defined in (21).

Proof.

By (15), we have that for $\gamma>0,$

[TABLE]

Similarly, we derive from (18) that for $\gamma>1/2,$

[TABLE]

Consequently, the formulas (85)-(86) follow from (3.1)-(3.1) directly. ∎

Remark 6.1**.**

In (84), we need the inverse transform of $R_{n}^{\lambda}(\xi),$ which can be computed by the same formulas (3.1)-(3.1). Indeed, by definition, we have

[TABLE]

We now consider a two dimensional example with $f(x,y)=\exp(-\sqrt{x^{2}+y^{2}})$. Notice that the Fourier transform of this source term can be computed as

[TABLE]

The corresponding numerical results are presented in Figure 7. Once again, the MMGF collocation method is more accurate and converges faster than the Hermite collocation method.

7 Summary and concluding remarks

In this paper, we have developed accurate spectral methods using rational basis (or modified mapped Gegenbauer functions) for PDEs with fractional Laplacian in unbounded domains. The main building block of the spectral algorithms is some explicit formulas for the Fourier transforms and fractional Laplacian of the rational basis. With these, we can construct rational spectral-Galerkin and collocation schemes by pre-computing the associated fractional differentiation matrices. We obtain optimal error estimates of rational spectral approximation in the fractional Sobolev spaces, and analyze the optimal convergence of the proposed Galerkin scheme. Numerical results show that the rational method outperforms the Hermite function approach. Future studies along this line include the error estimates of the rational collocation methods in section 6, fast pre-conditioner/solvers for high dimensional problems, and applications of the MMGFs approach to tempered fractional PDEs.