Hermite spectral collocation methods for fractional PDEs in unbounded domains
Tao Tang Huifang Yuan Tao Zhou
math.NA


























 Tang et. al.]Tao Tang1,
       Huifang Yuan2, and Tao Zhou3
 1 Department of Mathematics, Southern University of Sciences and Technology, Shenzhen, China. 

           2 Department of Mathematics, Hong Kong Baptist University, Hong Kong, China.
           3LSEC, Institute of Computational Mathematics, Academy of Mathematics and Systems Science,
Chinese Academy of Sciences, Beijing, China.
 tangt@sustc.edu.cn (T. Tang), 13480510@life.hkbu.edu.hk (H. Yuan),
          tzhou@lsec.cc.ac.cn (T. Zhou)







This work is concerned with spectral collocation methods for fractional PDEs in unbounded domains. The method consists of expanding the solution with proper global basis functions and imposing collocation conditions on the Gauss-Hermite points. In this work, two Hermite-type functions are employed to serve as basis functions. Our main task is to find corresponding differentiation matrices which are computed recursively. Two important issues relevant to condition numbers and scaling factors will be discussed. Applications of the spectral collocation methods to multi-term fractional PDEs are also presented. Several numerical examples are carried out to demonstrate the effectiveness of the proposed methods.





33C45,34K37,35R11,65M70



Hermite spectral collocation methods for fractional PDEs in unbounded domains
    [
    December 30, 2023
=============================================================================








§ INTRODUCTION


Many systems in science and engineering can be more accurately described by  using fractional partial differential equations (PDEs) rather than the traditional approaches <cit.>. This leads to an intensive investigation over the past two decades on efficient numerical methods for fractional PDEs. Among others, the finite difference method and the finite element method are two widely investigated methods in this direction, see, e.g., <cit.> and references therein.

Another powerful approach for fractional models is the spectral methods. In this approach, the key is to construct suitable basis functions to handle to the solution singularities. Along this direction, a recent advance is brought by Karniadakis and co-authors who proposed the so-called Jacobi poly-fractonomials based spectral methods <cit.>. These bases are eigenfunctions of the corresponding fractional and tempered fractional Sturm-Liouville problems. Another approach that employs the generalized Jacobi functions (GJFs) is proposed by Shen et al. <cit.>. Those bases are adapted to the fractional operator, as a fractional derivative of poly-fractonomials/GJFs is simply another poly-fractonomials/GJFs with different parameter. Consequently, fractional derivatives become a local operator in the physical space spanned by poly-fractonomials/GJFs, and this property leads to very efficient spectral methods for fractional PDEs in bounded domains. The poly-fractonomials/GJFs have been successfully applied to various fractional models <cit.>. However, compared to fractional PDEs in bounded domain, little works have been done for fractional PDEs defined on unbounded domains. Very recently, a spectral method for fraction differential equations in the half line is proposed in <cit.> – using the generalized Largurre functions as bases – extending the idea of <cit.>.

 When this paper is prepared, we noticed very recent work of Mao and Shen <cit.> who proposed both the spectral Galerkin and collocation method for fractional PDEs in unbounded domains. However, the collocation method therein relies on an equivalent formulation in frequency space by the Fourier transform, and performs collocation methods to the equivalent formulation that involve forward/backward Hermite transform. In contrast, our collocation methods are direct methods that based on the derivation of explicit DMs. Moreover, our approach can be easily applied to nonlinear problems as the differentiation matrices are constructed explicitly.

In this work, we aim at designing spectral collocation methods for fractional PDEs in unbounded domains (the whole space ℝ^d ).
To better demonstrate our idea, we consider the following model equation:

    (-Δ)^α/2u(𝐱) + γ f(u) = g(𝐱),      𝐱∈ℝ^d,
    
      u(𝐱)=0,     |𝐱| →∞,

where f(u) is a linear/nonlinear function of u, and the fractional Laplace operator is defined as <cit.>

    (-Δ)^α/2u(𝐱)=C_n,α∫_ℝ^du(𝐱)-u(𝐲)|𝐱-𝐲|^n+αd𝐲,     
    C_n,α=α2^α-1Γ(α+n2)π^n/2Γ(2-α2).

Notice that the fractional Laplace operator (-Δ)^α/2, where 0<α<2, recovers the standard Laplace operator as α→ 2.

For such problems that are defined on the whole space ℝ^n, there are alternative (equivalent) ways to define the fractional Laplace operator. For example, it can be defined as a pseudo-differential operator via the Fourier transform:

    ℱ[(-Δ)^α/2u](ξ)=|ξ|^αℱ[u](ξ).


Our spectral collocation methods handling the above equation consist of two parts: expanding the solution with the basis functions and imposing collocation conditions on the Guass-Hermite points. In particular, we shall consider two types of expansion bases, that is, the normalized Hermite functions {e^-x^2/2H_n(x)}_n and the over-scaled bases {e^-x^2H_n(x)}_n. For both approaches, we shall derive explicit formulas for the associated differential matrix (DM), for which the components can be computed efficiently by using a recurrence formula. To deal with solutions with different decay rate, a scaling factor will be included in the expansion. Application to multi-term fractional Laplace equations will also be discussed. It is noticed that both methods admit spectral convergence for solutions with exponential decay in infinity.

The rest of the this paper is organized as follows. The next section provides some preliminaries for some special functions. Our spectral collocations methods are presented and discussed in Sections 3 and 4, for the over-scaled bases {e^-x^2H_n(x)}_n  and the normalized Hermite functions, respectively. In Section 5 we shall discuss an equivalent collocation scheme by using the Lagrange type bases. Numerical examples are presented in Section 6 to demonstrate the effectiveness of the proposed spectral methods. We finally give some concluding remarks in Section 7.




§ PRELIMINARIES



This section will provide some preliminaries useful for designing our spectral collocation methods. For ease of notations, in this section we shall focus our attention to the one-dimensional case.



 §.§ Confluent hypergeometric functions

We first introduce the definition of confluent hypergeometric function of the first kind which is defined by the following power series <cit.>

    _1F_1(a,b;x)=∑_k=0^∞(a)_k(b)_kx^kk!,

where (a)_k is the Pochhammer symbol defined as

    (a)_0=1,   (a)_k=a(a+1)(a+2)...(a+k-1).

By the above definitions, we can easily get

    d^kdx^k_1F_1(a,b;x)=(a)_k(b)_k_1F_1(a+k,b+k;x).

We also have the following integral representation <cit.>

    _1F_1(a,b;x)=Γ(b)Γ(b-a)Γ(a)∫_0^1e^txt^a-1(1-t)^b-a-1dt,    a,   b >0.

It is noticed that the confluent hypergeometric function satisfies the Kummer's transformation formula

    _1F_1(a,b;-x)=e^-x_1F_1(b-a,b;x).

It it easy to check that the following recurrence formula holds

    (2a-b+x)_1F_1(a,b;x)=a_1F_1(a+1,b;x)-(b-a)_1F_1(a-1,b;x).




 §.§ Hermite polynomials/functions

The Hermite polynomials, denoted by H_n(x), n ⩾ 0, x∈ℝ, are defined by the following three-term recurrence relation (see e.g., <cit.>):

    H_0(x)=1,    H_1(x)=2x,   
    H_n+1(x)=2xH_n(x)-2nH_n-1(x),    n≥ 1.

The Hermite polynomials are orthogonal with respect to the weight function ω(x)=e^-x^2, namely,

    ∫_ℝH_m(x)H_n(x)e^-x^2dx=γ_nδ_mn,    γ_n=√(π)2^nn!.

It is well known that Hermite polynomials and the confluent hypergeometric function satisfy the following formulas:

    H_2n(x)=(-1)^n(2n)!n!_1F_1(-n, 1/2;x^2);
       H_2n+1(x)=(-1)^n(2n+1)!n!2x_1F_1(-n, 3/2;x^2).

The corresponding normalized Hermite functions are defined as

    H_n(x)=1√(2^nn!)e^-x^2/2H_n(x).

and they are orthogonal with respect to the weight function ω(x)=1, i.e.,

    ∫_ℝH_m(x)H_n(x)dx=√(π)δ_mn.

We shall discuss the spectral collocation methods based on the above Hermite functions in Section 4.

We shall also discuss the spectral collocation methods based on the over-scaled bases that defined as follows

    H_n(x)= e^-x^2/2H_n(x).

It is easy to see that the generalized Hermite defined here is orthogonal with respect to the weight function ω(x)=e^x^2, i.e.,

    ∫_ℝH_m(x)H_n(x)e^x^2dx=√(π)δ_mn.

Such bases were first proposed by Brinkman in <cit.> and have been well studied
in physics, see e.g., <cit.>.
The above over-scaled basis was first proposed by Brinkman when studying the so-called Fokker-Planck equations,
where the velocity part of the probability distribution function was expanded in Hermite functions
(<ref>).
His approach has become one of the most popular methods used for solving the Fokker-Planck equation, see, e.g., <cit.>.




 §.§ Bessel functions

We shall also use properties of generalized Bessel functions.  Recall that the Bessel function of order μ is defined as

    J_μ(x)=∑_m=0^∞ (-1)^m/m!  Γ(m+μ+1)(x/2)^2m+μ.

In particular we have

    J_-1/2(x)=√(2/π x)cosx.

For the Bessel functions it holds in <cit.> that

    ∫_ℝ^+J_μ(bt)exp(-p^2 t^2)t^ν-1dt = (b/2p)^μΓ(μ+ν/2)/2p^νΓ(μ+1)_1F_1(μ+ν/2, μ+1, -b^2/4p^2).

It is easy to check that for Bessel function with integer parameter it holds

    d/dx{x^-nJ_n(x)}=-x^-nJ_n+1(x).

Let m, n be two integers, and let ν≥ -n-1,   μ≥ -m-1 be two real numbers,  we denote

    μ + ν + m + n := δ;   μ + ν - m - n := ζ.

Then for a> 0 it holds (see e.g, <cit.>, P.216)

    ∫_0^ax^v+2n+1(a^2-x^2)^m+μ/2J_μ(b√(a^2-x^2))J_v(cx)dx 
       =a^ζ+1b^μc^v(∂b∂ b)^m(∂c∂ c)^n[(b^2+c^2)^-δ+12J_δ+1(a√(b^2+c^2))].







§ SPECTRAL COLLOCATION METHODS BASED ON {H_N(X)}_N


In this section, we first consider the spectral collocation method based on the {H_n(x)}_n.
We assume that the solution admits an exponential decay in infinity, and we shall approximate u(x) with a finite sum of the the basis {H_n(x)}_n, i.e.,

    u(x)≈ u_N(x)=∑_n=0^N-1c_nH_n(x).

By inserting the above expansion into the fractional PDE (<ref>), we obtain

    (-Δ)^α/2u_N(x) + γ f(u_N(x)) = g(x).

Let {x_i}_i=0^N-1 be the roots of the N-th order Hermite polynomials, we then impose the collocation conditions on these collocation points, which yields

    ∑_n=0^N-1c_n(-Δ)^α/2H_n(x_i) + γ f(u_N(x_i)) =g(x_i),     i=0, 1,..., N-1.

Then we can write the above equations into the following system

    𝒟^α𝐜 + γ F(𝐜)=𝐠,

where 𝐜=(c_0, ..., c_N-1)^T is the unknown coefficient vector, and 𝒟^α∈ℝ^N× N is the differential matrix with components

    𝒟_i,j^α=(-Δ)^α/2H_j(x_i),    i , j = 0, 1,..., N-1.

By solving the above linear system one gets an approximated solution u_N(x). Next, we shall derive explicit formulas for the components of the differential matrix.



 §.§ The one dimensional case

We first consider the even terms H_2n. Precisely, we have the following theorem

For 0< α <2 , we have

    (-Δ)^α/2H_2n(x)=2^α(-1)^n√((2n)!)2^nn!Γ(n+α/2+1/2)Γ(n+1/2)_1F_1(n+α2+12,12,-x^2).



Consider the forward Fourier transform, we have

    ℱ[H_2n](ξ)
       =1/√(2π)√(2^2n(2n)!)∫_ℝexp(-x^2)H_2n(x)e^-xξdx
       =2/√(2π)√(2^2n(2n)!)∫_ℝ^+exp(-x^2)H_2n(x)cos(xξ)dx
       =(-1)^n/√(2)√(2^2n(2n)!)ξ^2ne^-ξ^2/4.

Then by considering the inverse Fourier transform we have

    (-Δ)^α/2H_2n(x)
       =1/√(2π)∫_ℝ|ξ|^αℱ[H_2n](ξ)e^xξdξ
       =1/√(2π)(-1)^n/√(2)√(2^2n(2n)!)∫_ℝ|ξ|^αξ^2ne^-ξ^2/4e^xξdξ
       =2/2√(π)(-1)^n/√(2^2n(2n)!)∫_ℝ^+ξ^2n+αe^-ξ^2/4cos(xξ)dξ
       =(-1)^n/√(π)√(2^2n(2n)!)2^2n+αΓ(2n+α+1/2)_1F_1(2n+α+1/2,1/2;-x^2)
       =2^α(-1)^n√((2n)!)2^nn!Γ(n+α/2+1/2)Γ(n+1/2)_1F_1(n+α2+12,12,-x^2).

This completes the proof.

.3cm

In the above derivation, we have adopted some integral formulas in <cit.>. Using similar arguments, for the odd terms {H_2n+1(x)}, we have the following theorem

For 0<α<2, it holds

    (-Δ)^α/2H_2n+1(x) 
       =    2^α+1(-1)^n√((2n+1)!)2^n+1/2n!Γ(n+α/2+3/2)Γ(n+3/2)x_1F_1(n+α2+32,32,-x^2).



By Theorems 3.1-3.2, we get the following explicit formula for the components of the differential matrix (DM)

    𝒟_ij^α=(-Δ)^α/2H_j(x_i) 
       =   
    2^α(-1)^n√((2n)!)2^nn!Γ(n+α/2+1/2)Γ(n+1/2)_1F_1(n+α2+12,12,-x_i^2),    j=2n; 
    
    2^α+1(-1)^n√((2n+1)!)2^n+1/2n!Γ(n+α/2+3/2)Γ(n+3/2)x_i_1F_1(n+α2+32,32,-x_i^2),   j=2n+1.

We now summarize the procedure for computing the components of the DM:


  * For each x_i, compute the quantities 𝒟^α_ij with j=0,1,2,3, by the above formula. Notice that one has to deal with confluent hypergeometric function _1F_1 and this may be non-trivial, as by definition (<ref>), this is an infinite expansion.  Nevertheless, one can find a fast & accurate algorithm for example in <cit.>.


  * Compute the quantities 𝒟^α_ij for 4<j ≤ N-1 in a recurrence way using the recurrence formula (<ref>).

In general, the above DM is easy to construct, and in fact, the matrix can be stored in priori (offline). We provide in Fig. <ref> the condition number of the differential matrix with respect to the number of collocation points N. It is noticed that the condition number grows very fast with respect to N. Thus, efficient pre-conditioners should be designed for practice applications.

Notice that we have used the definition (<ref>) to derive the differential matrix, however, one can also adopt the definition (<ref>) to derive a similar differential matrix.



 §.§ The two dimensional case


For multi-dimensional cases, the bases used will be the tensorized 1D bases. Here we take the two dimensional case as an example. To this end, the two dimensional bases take the following  form {H_n(x)H_m(y)}_n,m.


By (<ref>), we know that

    ℱ[H_nH_m](ξ,η)
       =1/2π∫_ℝ∫_ℝH_n(x)H_m(y)e^-xξe^-yηdxdy 
       =1/2π∫_ℝH_n(x)e^-xξdx∫_ℝH_m(y)e^-yηdy=ℱ[H_n](ξ)ℱ[H_m](η).

Then, using the inverse Fourier transform gives

    (-Δ)^α/2[H_n(x)H_m(y)] 
       =1/2π∫_ℝ∫_ℝ(ξ^2+η^2)^α/2ℱ[H_n](ξ)ℱ[H_m](η)e^xξe^yηdξ dη.

We take the double-even term {H_2n(x)H_2m(y)} as an example. We denote

    L_nm=(-1)^n/√(2)√(2^2n(2n)!)(-1)^m/√(2)√(2^2m(2m)!).

Then we have

    (-Δ)^α/2[H_2n(x)H_2m(y)] 
       =L_nm/2π∫_ℝ∫_ℝ(ξ^2+η^2)^α/2ξ^2ne^-ξ^2/4η^2me^-η^2/4e^xξe^yηdξ dη
       =4L_nm/2π∫_ℝ^+∫_ℝ^+(ξ^2+η^2)^α/2ξ^2ne^-ξ^2/4η^2me^-η^2/4cos(xξ)cos(yη)dξ dη
       =4L_nm/2π ∫_ℝ^+ℐ(x,y,ρ)  ρ^2n+2m+α+1e^-ρ^2/4dρ,

where

    ℐ(x,y,ρ)=∫_0^π/2(cosθ)^2n(sinθ)^2mcos(xρcosθ)cos(yρsinθ)dθ.

By the definition of the Bessel function (<ref>) and the property (<ref>)  we have

    ∫_0^π/2(cosθ)^2n(sinθ)^2mcos(pcosθ)cos(qsinθ)dθ
       =π/2(pq)^1/2∫_0^1(1-x^2)^n-1/4x^2m+1/2J_-1/2(p√(1-x^2))J_-1/2(qx)dx
       =π/2(∂p∂ p)^n(∂q∂ q)^m[(√(p^2+q^2))^-(n+m)J_n+m(√(p^2+q^2))].

Notice that by property (<ref>) we have

    (∂p∂ p)^n(∂q∂ q)^m[(√(p^2+q^2))^-(n+m)J_n+m(√(p^2+q^2))]
       =(-1)^n+m(√(p^2+q^2))^-(2n+2m)J_2n+2m(√(p^2+q^2)).

Now, by replacing p and q with ρ x and ρ y, we obtain

    ℐ(x,y,ρ)=π/2(-1)^n+mρ^-(2n+2m)(√(x^2+y^2))^-(2n+2m)J_2n+2m(ρ√(x^2+y^2)).

We then do inverse Fourier transform to obtain (for t=m+n)


    (-Δ)^α/2[H_2n(x)H_2m(y)]
    
       =(-1)^tL_nm(x^2+y^2)^-t∫_ℝ^+ρ^α+1e^-ρ^2/4J_2t(ρ√(x^2+y^2))dρ
       =2^αΓ(t+α/2+1)/√(2^2n(2n)!)√(2^2m(2m)!)Γ(2t+1)_1F_1(t+α/2+1;2t+1;-(x^2+y^2)).


In the above derivation, we have used property (<ref>). In a similar way, we can derive similar results for the bases of forms  H_2n+1(x)H_2m(y), H_2n(x)H_2m+1(y) and H_2n+1(x)H_2m+1(y). To summaries these results in an unified way, we introduce the parameter δ_1 and δ_2, with δ_1,δ_2 ∈{0,1}. And let

    a=n+m+α/2+δ_1+δ_2+1;     b=2n+2m+δ_1+δ_2+1.

Using similar arguments as above, we can derive that

    (-Δ)^α/2[H_2n+δ_1(x)H_2m+δ_2(y)]= C_x,y(α, m,n,δ_1,δ_2) Γ(a)/Γ(b)_1F_1(a;b;-(x^2+y^2)),

with

    C_x,y(α,m,n,δ_1,δ_2)=2^α+δ_1+δ_2x^δ_1y^δ_2/√(2^2n+δ_1(2n+δ_1)!)√(2^2m+δ_2(2m+δ_2)!).






We finally provide the explicit formulas for the components of the DM as following

    𝒟   _(i-1)*N+j,(p-1)*N+q^α=(-Δ)^α/2{H_p(x_i)H_q(y_j)}
       =2^αΓ(n+m+α/2+1)_1F_1(n+m+α/2+1;2n+2m+1;-(x_i^2+y_j^2))√(2^2n(2n)!)√(2^2m(2m)!)Γ(2n+2m+1), 
                                                    p=2n, q=2m;  
    2^α+1Γ(n+m+α/2+2)y_j_1F_1(n+m+α/2+2;2n+2m+2;-(x_i^2+y_j^2))√(2^2n(2n)!)√(2^2m+1(2m+1)!)Γ(2n+2m+2), 
                                                   p=2n, q=2m+1;  
    2^α+1Γ(n+m+α/2+2)x_i_1F_1(n+m+α/2+2;2n+2m+2;-(x_i^2+y_j^2))√(2^2n+1(2n+1)!)√(2^2m(2m)!)Γ(2n+2m+2), 
                                                   p=2n+1, q=2m;  
    2^α+2Γ(n+m+α/2+3)x_iy_j_1F_1(n+m+α/2+3;2n+2m+3;-(x_i^2+y_j^2))√(2^2n+1(2n+1)!)√(2^2m+1(2m+1)!)Γ(2n+2m+3), 
                                                    p=2n+1, q=2m+1.

Again, the above components can be computed in a similar procedure as in the one-dimensional case. We remark that one may also derive such formulas for the three dimensional case, and we omit it here as it may involves complex notations.



 §.§ The use of scaling factors

It is well known that for Hermite-type spectral methods,  the convergence rate deteriorates if the decay rates between the solution and the bases function have a relatively large gap.
A remedy to fix this problem is to use the so-called
scaling factor  <cit.>. We now introduce the basic idea of the scaling factor and choose the one-dimensional case as an example for illustration.
To this end, let
u(x) be a function that decay exponentially, namely,

    |u(x)|∼ 0,    ∀   |x|>M,

where M>0 is some constant.
The idea of using the scaling factor is to expand u as

    u(x)=∑_n=0^N-1 c_n H_n(r x)   ⇔  
    u(x/r)=∑_n=0^N-1 c_n H_n(x),

where r >0 is a scaling factor. The key point of using
r is to scale the Hermite-Gauss nodes {x_k}_k=0^N-1 so that the collocation points {x_k/r}_k=0^N-1
are well within the effective support of u. This suggests the following choice

    max_0 ≤ k ≤ N-1{|x_k|}/r ≤ M   ⇒   r = max_0≤ k ≤ N-1{|x_k|}/M.

We remark however, in practice, finding the quantity M may be non-trivial and thus the optimal scaling factor is hard to obtain in general.

We now show how to include the scaling factor in our spectral collocation methods for the fractional PDEs. Let us illustrate the idea in 1D. We now seek to the following expansion

    u_N(x)=∑_n=0^N-1c_nH_n(rx).

By inserting the expansion to the fractional PDE (<ref>) and imposing the collocation condition we obtain

    𝒟^α,r𝐜 + γ F^r(𝐜)=𝐠,

where the 𝒟^α,r is the differential matrix with components

    𝒟_ij^α,r =(-Δ)^α/2H_j( rx_i),    i, j =0, ... , N-1.

We need to deal with the fractional Laplacian of H_j(rx). To this end, suppose that

    (-Δ)^α/2v(x)=ϕ(x),

Then the fractional Laplacian of v_r(x)=v(rx) is

    (-Δ)^α/2v_r(x)=r^αϕ(rx).

This can be simply proved by using the definition (<ref>). By using together the above argument and the DM formula (<ref>), we have

    𝒟_ij^α,r   =(-Δ)^α/2H_j( rx_i) 
       =
    (2r)^α(-1)^n√((2n)!)2^nn!Γ(n+α/2+1/2)Γ(n+1/2)_1F_1(n+α/2+1/2,1/2,-z_i^2),    j =2n ; 
    
    2^α+1r^α(-1)^n√((2n+1)!)2^n+12n!Γ(n+α2+32)Γ(n+32)z_i_1F_1(n+α2+3/2,32,-z_i^2),    j =2n+1.

where z_i = r x_i for i=0,1, ..., N-1.



 §.§ Applications to multi-term fractional PDEs

In this section, we claim that our spectral collocation method can be easily applied to the multi-term fractional PDEs.  We shall still take the one-dimensional case as an example. Consider the following multi-term fractional PDEs

    ∑_j=1^J (-Δ)^α_j/2u(x) + γ f(u) = g(x),    x∈ℝ.

The above problem is motivated by the approximation of distributed order fractional models
using a quadrature rule, see e.g. <cit.>.

For the above multi-term models, by inserting the Hermite expansion (<ref>) into the equation and
imposing the collocation condition, one gets the following system:

    𝒟^𝒥𝐜 + γ F(𝐜)=𝐠        𝒟^𝒥= ∑_j=1^J 𝒟^α_j.

Notice that the components of each differential matrix 𝒟^α_j can be computed by the explicit formula (<ref>).




§ SPECTRAL COLLOCATION METHODS BASED ON THE NORMALIZED HERMITE FUNCTIONS {H_N}_N.


In the last section, we have proposed the spectral methods with the over-scaled bases {H_n}_n. While the associated DM is easy to compute, its condition number grows fast with respect to the number of collocation points, and this is due to the poor property of the bases. In this section, we shall discuss the spectral collocation methods based on the normalized Hermite functions {H_n}_n. The main task is still to derive the explicit formula for the differential matrix.




 §.§ The one dimensional case

The aim is the present the explicit formula for 𝒟_mn^α=(-Δ)^α/2H_n(x_m). This also relies on the forward/inverse Fourier transform as done in the last section.

By the forward Fourier transform, we Notice that

    ℱ{H_n(x)}=1/√(2π)∫_ℝH_n(x)e^-ξ xdx=(-)^nH_n(ξ).

That is, the normalized Hermite functions {H_n}_n are eigenfunctions of the Fourier transform with eigenvalues {(-)^n}_n with =√(-1).


A more difficult part is to compute the inverse Fourier transform. By the above fact and (<ref>), we have to deal with the inverse fourier transform of  {(-)^n |ξ|^αH_n(ξ)}_n. Notice that we can write

    H_n(ξ)=∑_k=0^nâ_n,kexp(-ξ^2/2)ξ^k,      n=0, 1, ..., N-1

with a_n,k=1/√(2^nn!)a_n,k for k≤ n, where a_n,k can be computed in a recursion way

    a_0,0=1,    a_1,0=0,    a_1,1=2;
        a_n+1,k=-a_n,k+1,    k=0;
        a_n+1,k=2a_n,k-1-(k+1)a_n,k+1,    k>0.

Let us first consider the even terms with k=2m, the inverse Fourier transform yields
 
    ℱ^-1[exp(-ξ^2/2)ξ^2m|ξ|^α](x)
       =1/√(2π)∫_ℝexp(-ξ^2/2)ξ^2m|ξ|^αe^ξ xdξ
       =2/√(2π)∫_ℝ^+exp(-ξ^2/2)ξ^2m+αcos(ξ x)dξ
       =2^2m+α/2/√(π)Γ(2m+1+α/2)_1F_1(2m+1+α/2,1/2,-x^2/2).

Then for the odd terms with k=2m+1, it holds
 
    ℱ^-1[exp(-ξ^2/2)ξ^2m+1|ξ|^α](x) 
       =1/√(2π)∫_ℝexp(-ξ^2/2)ξ^2m+1|ξ|^αe^ξ xdξ
       =2/√(2π)∫_ℝ^+exp(-ξ^2/2)ξ^2m+1+αsin(ξ x)dξ
       =2^2m+2+α/2/√(π)Γ(2m+3+α/2)x_1F_1(2m+3+α/2,3/2,-x^2/2).

For ease of notations, we denote

    F_k(x)=ℱ^-1[exp(-ξ^2/2)ξ^k|ξ|^α](x),    k=0, 1, ..., N-1.

Then by (<ref>) the components of the differentiation matrix yield

    𝒟_mn^α=(-Δ)^α/2H_n(x_m)=(-)^n
    ∑_k=0^N-1â_n,kF_k(x_m),    0 ≤ n, m ≤ N-1.

Notice that when k>n, we set a_n,k=0. In Fig. <ref> we present the condition number of this differential matrix with respect to N.

It is noticed from Fig. <ref> that the condition number grows algebraically with N – much well behaved than the previous case (where the over scaled bases are used). In fact the computation complexity of the DM above is almost the same as the formula (<ref>) since each components here is simply a evaluation of a linear combination of the confluent hypergeometric functions.




 §.§ The two dimensional case

We now consider the two dimensional case. By the forward Fourier transform we have

    ℱ[H_nH_m](ξ,η)
       =1/2π∫_ℝ∫_ℝH_n(x)H_m(y)e^-xξe^-yηdxdy 
       =1/2π∫_ℝH_n(x)e^-xξdx
    ∫_ℝH_m(y)e^-yηdy 
       =(-)^n+mH_n(ξ)H_m(η).

Then, by the inverse Fourier transform we obtain

    (-Δ)^α/2   [H_n(x)H_m(y)]=(-)^n+m/(2π)∫_ℝ∫_ℝ(ξ^2+η^2)^α/2H_n(ξ)H_m(η)e^xξe^yηdξ dη.

Similar as in the one dimensional case, we expand H_n(ξ)H_m(η) as a combination of ξ^kη^l. Inspired by (<ref>) we have

    H_n(ξ)H_m(η)=∑_k=0^N-1∑_l=0^N-1â_n,kâ_m,lexp(-(ξ^2+η^2)/2)ξ^kη^l.

Next, we should deal with the inverse Fourier transform of terms like  
    (ξ^2+η^2)^α/2exp(-(ξ^2+η^2)/2)ξ^kη^l,    0≤ k,l ≤ N-1.

The derivation is similar as in the one dimensional case, and thus we omit the details here. To summarize, we introduce the matrix F(ξ,η) as

    F_k,l=ℱ^-1{(ξ^2+η^2)^α/2exp(-(ξ^2+η^2)/2)ξ^kη^l}
    
    =   2^α/2(-1)^p+qΓ(p+q+α/2+1)_1F_1(p+q+α/2+1;2p+2q+1;-x^2+y^2/2)2^p+qΓ(2p+2q+1), 
                                                    k=2p,   l=2q;  
    2^(α+1)/2(-1)^p+qiΓ(p+q+α/2+2)y_1F_1(p+q+α/2+2;2p+2q+2;-x^2+y^2/2)2^p+qΓ(2p+2q+2), 
                                                   k=2p,  l=2q+1;  
    2^(α+1)/2(-1)^p+qiΓ(p+q+α/2+2)x_1F_1(p+q+α/2+2;2p+2q+2;-x^2+y^2/2)2^p+qΓ(2p+2q+2), 
                                                   k=2p+1,   l=2q;  
    -2^α/2+1(-1)^p+qΓ(p+q+α/2+3)xy_1F_1(p+q+α/2+3;2p+2q+3;-x^2+y^2/2)2^p+qΓ(2p+2q+3), 
                                                    k=2p+1,   l=2q+1.

Finally, we provide the explicit formula for the components of the differentiation matrix as follows

    𝒟_(p-1)*N+q,(n-1)*N+m^α   =(-Δ)^α/2{H_n(x_p)H_m(y_q)}
       =(-i)^n+m∑_k=0^N-1∑_l=0^N-1â_n,kâ_m,lF_k,l(x_p,y_q).


Similar as in Section 3, for the spectral collocation methods with the normalized Hermite functions, one can easily include the scaling factors with slight modifications to the differential matrix. And furthermore, the application to multi-term fractional PDEs is also straightforward.




§ SPECTRAL COLLOCATION METHODS WITH LAGRANGE TYPE BASES


In the above sections, we have derived explicit formulas for the differential matrices by using the Hermite-type bases. In this section, we shall discuss the equivalent spectral collocation methods with Lagrange type bases. For notation simplicity, we shall only present the one dimensional case. Notice that for this type of bases, we still search an (equivalent) approximated solution in the same finite space as in the above sections, however, such an approach indeed results in different differential matrices. By using the Lagrange type bases, we approximate the solution in the following way

    u_N(x)=∑_j=0^N-1u_jh_j(x),      h_j(x_k) = δ_jk,    0≤ j,k≤ N-1.

Here the Lagrange type bases {h_j(x)}_j=1^N are defined as

    h_j(x)=e^-x^2/2/e^-x_j^2/2∏_i=0,i≠ j^N-1x-x_j/x_i-x_j,    0≤ j≤ N-1.

The associated points {x_j}_j=0^N-1 are the Gauss-Hermite points. It is clear that we can express each Lagrange-type basis with the normalized Hermite functions, i.e.,

    h_j(x)=∑_k=0^N-1  b_k^jH_k(x),       b_k^j=1/√(π)H_k(x_j)ω̂_j,     0≤ j,k≤ N-1,

where {ω̂_j}_j=0^N-1 is the weights of the Gauss quadrature rule associated with the Hermite functions, which are defined as:

    ω̂_j=√(π)/NH_N-1^2(x_j),   j=0,1,...,N-1.

Consequently,  we can easily derive the associated differential matrix D^L,α with Lagrange type bases

    D_i,j^L,α=(-Δ)^α/2h_j(x_i)=∑_k=0^N-1b_k^j(-Δ)^α/2H_k(x_i).

The quantities (-Δ)^α/2H_k(x_i) can be obtained via equation (<ref>). The condition number of the above differential matrix with respect to the order N is presented in Fig. <ref>. Again, we can see that the condition number grows algebraically respect to the order N.





§ NUMERICAL EXAMPLES


In this section, we shall present several constructive examples to show the convergence property of our spectral collocation methods.  In all our computations, we shall report the numerical error both in the weighted norm e_w and in the maximum norm e_m, which are defined respectively as

    e_w=u(x)-u_N(x)_ω^2,    e_m=max_j|u(x_j)-u_N(x_j)|.

Here ω(x)=e^x^2 for the bases {H_n(x)}_n and ω(x)=1 for the normalized Hermite functions.



 §.§ The fractional Laplace equation

Our first example is the fractional Laplace equation

    (-Δ)^α/2u(x)  = g(x),      x∈ℝ,
    
      u(x)=0,     |x| →∞.

The right hand side is chosen such that the exact solution is u(x)=exp(-x^2)x. Notice that for a given exact solution u(x), the right hand side does not necessary have explicit formula, and in such cases, we shall compute it by expanding u(x) with a large enough number of the basis functions. This is also true for all other examples.


We perform computations with different fractional order, i.e., α=0.4, 1 and 1.6. The numerical errors against the numerical of points N with the over scaled bases and the normalized Hermite functions are presented in Fig. <ref> and Fig. <ref>, respectively. For the over-scaled bases, no scaling factor is needed, i.e, r=1, while for the normalized Hermite functions, we choose r=√(2). In Fig. <ref> and Fig. <ref>, the numerical error in weighted norm and in maximum norm are reported in the left plot and right plot, respectively. It is clear that spectral convergence is obtained for all cases of α. While no scaling factor is needed for the over scaled bases, it is noticed that the convergence is polluted when a larger number of collocation points are used, due to the fast grow of the condition number.





 §.§ A linear fractional PDE


Next, we consider the following fractional PDE

    (-Δ)^α/2u(x)+2u(x)=f(x),     x∈ℝ
    
    u(x)=0,     x →∞

We first test the performance of the over-scaled bases. We set u(x)=exp(-x^2/2)x^2cos(x), and the right hand side can be computed accordingly.  Again, we consider α=0.4,1 and 1.6. The numerical results with a scaling factor r=1/√(2) and without the scaling factor (r=1) are reported in Fig. <ref> and Fig. <ref>, respectively. It is clear seen that using a proper scaling factor results in faster convergence  rate for both the weighted error and the maximum error.









Now we test the normalized Hermite functions. For the same equation we choose the right hand side such that the solution yields u(x)=exp(-2x^2)x^2cos(x). It is clear that the optimal scaling factor is r=2. The numerical results with a scaling factor r=2 and without a scaling factor (i.e., r=1) are reported Fig. <ref> and Fig. <ref>, respectively. Again, it is shown that a proper scaling factor can be useful to speed up the convergence.





 §.§ A two-dimensional example

We now consider a two dimensional example, and the equation considered is

    (-Δ)^α/2u(x,y)+2 u(x,y) = g(x,y).

The exact solution is chosen as u(x,y)=exp(-(x^2+y^2))sin(x+y). We also perform the computations with α=0.4,  1 and 1.6. To avoid too much pictures, here we only test the performance of the over-scaled bases. The numerical errors in weighted and maximum norm against the numerical of number of collocation points N are presented in Fig. <ref>. Spectral convergence is again observed.


















 §.§ A multi-term fractional model

Our next example the multi-term Laplacian equation:

    ∑_j=1^J (-Δ)^α_j/2u(x) = g(x),    x∈ℝ.

Here we set J=4 and {α_j}_j=1^J are chosen as the transformed Legendre-Gauss points:

    α_1= 0.139,    α_2= 0.660,   α_3= 1.340,   α_4= 1.861.

We set the exact solution to be u(x)=(- 3x^2/2)(sinx+x^6+x^2cosx) and the right hand side can be computed accordingly. Again, we test the performance of the over-scaled bases. In this example, we consider scaling factors r=√(1.5), √(1.3), and the approach without a scaling factor, i.e., r=1. The corresponding numerical results are reported in Fig. <ref>.

We can see that both the weighted error and maximum error decay fast for the case of r=√(1.5). And this indicates the effectiveness of using a scaling in improving convergence rate.



























 §.§ A nonlinear example

Our next example is a nonlinear fractional PDE

    (-Δ)^α/2u(x)+ u^2(x)=g(x)

For the over scaled bases H_n(x), we set the exact solution as u(x)=exp(-x^2)(sin(x)+x^2).
For the normalized Hermite functions H_n(x), the exact solution is chose to be u(x)=exp(-x^2/2)(sin(x)+x^2). In our computations, for each expansion number N, we use the Newton iteration method with a tolerance 10^-16 to deal with the nonlinear term. The performance of collocation method with the normalized Hermite functions are presented Fig. <ref>, and it is shown that the method yields a spectral convergence rate.










 §.§ An eigenvalue problem


Finally we consider the following eigenvalues problem

    ((-Δ)^α/2+x^2)u(x)=λ u(x)

The above eigenvalue problem with α=1 has been analyzed in <cit.>. In particular, the eigenvalues of this problems is given by

    λ_2k-1=-a_k',   λ_2k=-a_k,    k=1,2,...,

where a_k and a_k' are the roots of the following Airy function and its derivative (in the decreasing order)

    A(x)=1/π∫_0^∞cos(t^3/3+xt)dt.


In this example, we shall compute the first three eigenvalues by the spectral collocation method. The exact eigenvalues are

    λ_1≈1.01879297164747,   λ_2≈2.33810741045976,   λ_3≈3.24819758217983.

Numerical result are presented in Fig. <ref> with loglog scale. An algebraic decay is observed and this is due to the algebraic decay (non-exponential decay) of eigenvalues.







§ SUMMARY AND CONCLUSION


In this work, we have proposed a spectral collocation method based on Hermite functions for fractional PDEs in unbounded domain. One basis function is associated with the over-scaled weight exp(-x^2)H_n(x) and another is with the standard normalized-weight exp(-x^2/2)H_n(x), each has some advantages or popularity in practice. For both approaches, explicit differentiation matrices are derived. To deal with solutions with fast or  slow decay rate, a scaling factor in spectral implementation is discussed.

Although the numerical experiments indicate the spectral rate of convergence, there are still several issues requiring future investigations:



  * Rigorous convergence analysis is not covered in this work, which is still an ongoing work.


  * For over-scaled Hermite functions (which was used extensively in physics, see <cit.>) the condition numbers of the associated matrices grow fast with respect to N. Thus, it is useful to investigate some efficient pre-conditioners in this case.


  * We only provide some ad-hoc discussions on the scaling factors. It will be more meaningful to provide some more practical guidance on the optimal choice of the scaling factors, see, e.g., <cit.>.







§ ACKNOWLEDGMENTS

This work is partially supported by the National Natural Science Foundations of China under grant numbers 91630312, 91630203, 11571351, and 11731006. The second author is supported by a Hong Kong PhD Fellowship. The last author is supported by the science challenge project (No. TZ2016001), NCMIS, and the youth innovation promotion association (CAS).




99
physics_II E. Berkai, R. Metzler, and J. Klafter, From continuous time random walks to the fractional Fokker-Planck equation, Phys. Rev. E(3) 61(2000), no. 1, 132-138.

BFW98  P. Biler, T. Funaki, and W. A. Woyczynski. Fractal Burgers equations. J. Differential Equations, 148(1):9-46, 1998.




Brinkman H. C. Brinkman, Brownian motion in a field of force and the diffusion theory of chemical
reactions, Physica, 22 (1956), pp. 29-34.




distributed_I M. Caputo. Distributed order differential equations modelling dielectric induction and diffusion. Fract. Calc. Appl. Anal., 4:421-442, 2001.

Chen_Mao_Shen  S. Chen, Z. Mao and J. Shen,  Efficient and accurate spectral methods using general Jacobi functions for solving Riesz fractional differential equations. Appl. Numer. Math., 106:165-181, 2016.


spectralshen S. Chen, J. Shen and L. Wang, Generalized Jacobi functions and their applications to fractional differential equations, Math. Comp. 85(2016), no. 300, 1603-1638.




distributeDiethelm Kai Diethelm, and Neville J. Ford, Numerical analysis for distributed-order differential equations. J. Comput. Appl. Math. 225(2009), no. 1, 96-104.

Fok
J. C. M. Fok, B.-Y. Guo, and T. Tang,
Combined Hermite spectral-finite difference method for the Fokker-Planck equations, Math. Comp. 71 (2002), 1497-1528.

tableofintegrals I. S. Gradshteyn, I. M. Ryzhik, Table of Integrals, Series, and Products, Translated from the Russian. Sixth Edition. Translation edited and with a preface by Alan Jeffrey and Daniel Zwillinger. Academic Press, Inc., San Diego, CA, 2000.


integraltransform B, Harry, Tables of Integral Transforms, Vol.I, McGraw-Hill Book Company , New York, 1954, ISBN 07-019549-8.

differenceHuang Y. Huang and A. M. Oberman, Numerican methods for the fractional Laplacian: a finite difference-quadrature approach, SIAM J. Numer. Anal. (2014), 52(6), 3016-3084.

recurrenceof1F1 M. Humi, W. Miller, Second Course in Ordinary Differential Equations for Scientists and Engineers, Universitext. Springer-Verlag, New York, 1988. x+441 pp. ISBN: 0-387-96676-5.

differenceJi C. Ji and Z. Sun, A high-order compact finite difference scheme for the
fractional sub-diffusion equation, J. Sci. Comput. 64(2015), no.3, 959-985.

elementJin B. Jin, R. Lazarov and Z. Zhou, Error estimates for a semi discrete finite element method for fractional order parabolic equations, SIAM J. Numer. Anal. 51(2013)445-466.

eigenvalue Lörinczi, József, and J. Ma  -0.37em lecki. Spectral properties of the massless relativistic harmonic oscillator, J. Differ. Equ., 253.10(2012):2846-2871.

KZK16 E. Kharazmi, M zayernouri, and G. E. Karniadakis. Petrov-Galerkin and spectral collocation methods for distributed order differential equations, SIAM J. Sci. Comput., 39(3), pp. A1003-A1037, 2017.

spectralArab H. Khosravian-Arab, M. Dehghan, M.R.Eslahchi, Fractional Sturm-Liouville boundary value problems in unbounded domains: theory and applications, J. Comput. Phys. 299(2015), 526-560.

Jiao16 Y.-J. Jiao, T.-J. Wang and Q. Zhang,
    A fully discrete spectral method for fisher's equation on the whole line,
    East Asian J. Appl. Math., 6 (2016), pp. 400-415.

Lan72  N. S. Landkof. Foundations of modern potential theory. Springer-Verlag, New York, 1972. Die
Grundlehren der mathematischen Wissenschaften, Band 180.

NMTMA G. Li, C. Sun, X. Jia, and D. Du, Numerical solution to the multi-term time fractional diffusion equation in a finite domain, Numer. Math. Theor. Meth. Appl., 9(3), pp.337-357, 2016.

LZK16 A. Lischke, M. Zayernouri, and G. E. Karniadakis. A tunably-accurate Laguerre Petrov-Galerkin spectral method for multi-term fractional differential equations on the half line,  arXiv: 1607.08571, 2016.

Sun H. Ma, W. Sun and T. Tang,
Hermite spectral methods with a time-dependent scaling for parabolic equations in unbounded domains, SIAM J Numer Anal., 43(2005), 58-75.

Mao_Shen Z. Mao and J. Shen, Hermite spectral methods for fractional PDEs in unbounded domains, to appear in SIAM J. Sci. Comput., 2017.

dPQRV11  A. de Pablo, F. Quirós, A. Rodríguez, and J. L. Vázquez. A fractional porous medium equation. Adv. Math., 226(2):1378-1409, 2011.

computent1F1 J. Pearson, Computation of hypergeometric function, master's thesis, University of Oxford, UK, Sep. 2009.

integralsandseriesA. P. Prudnikov, Y. A. Brychkov, O. I. Marichev, Integrals and series: special functions, Gordon and Breach Science Publishers, 1992, ISBN 2-88124-090-9(v.2).

Ren15
J. Ren and Z. Sun, Efficient numerical solution of the multi-term time fractional diffusion-wave equation, East Asian J. Appl. Math., 5 (2015), pp. 1-28

Risken H. Risken, The Fokker-Planck Equation: Methods of solution and Applications, 2nd ed., Springer-Verlag, Berlin, 1989.

SWT J. Shen, L.-L. Wang and T. Tang, Spectral Methods: Algorithms, Analysis and Applications, Springer, 2011.

SS C. Sheng and  J. Shen, A hybrid spectral element method for fractional
two-point boundary value problems, Numer. Math. Theor. Meth. Appl., 10(2), pp.437-464, 2017.

Scaling T. Tang, The Hermite spectral method for Gaussian type functions, SIAM J. Sci. Comput., 14 (1993), 594-606.

TD15 W. Tian, H. Zhou, and W. Deng, A class of second order difference approximations for solving space fractional diffusion equations,  Math. Comput., Vol. 84, pp. 1703-1727, 2015.

differenceWang H. Wang, T.S. Basu, A fast finite difference method for two dimensional space fractional diffusion equations, SIAM J. Sci. Comput. 34 (2012), no. 5, A2444-A2458.

Wang16
T. Wang, C. Zhang, and Q. Zhang, Mixed spectral method for heat transfer using generalised Hermite functions and Legendre polynomials. East Asian J. Appl. Math., 6 (2016),  pp. 448-465

WZ16 S. Wu and T. Zhou, Fast parareal iterations for fractional diffution equations, J. Comput. Phys., 329:210-226, 2017.





spectralgeorge M. Zayernouri and G.E.Karniadakis, Fractional Sturm-Liouville eigen-problems: Theory and numerical approximations, J. Comput. Phys. 252 (2013), 495-517.



