    # Well-Conditioned Galerkin Spectral Method for Two-Sided Fractional Diffusion Equation with Drift

In this paper, we focus on designing a well-conditioned Glarkin spectral methods for solving a two-sided fractional diffusion equations with drift, in which the fractional operators are defined neither in Riemann-Liouville nor Caputo sense, and its physical meaning is clear. Based on the image spaces of Riemann-Liouville fractional integral operators on L_p([a,b]) space discussed in our previous work, after a step by step deduction, three kinds of Galerkin spectral formulations are proposed, the final obtained corresponding scheme of which shows to be well-conditioned—the condition number of the stiff matrix can be reduced from O(N^2α) to O(N^α), where N is the degree of the polynomials used in the approximation. Another point is that the obtained schemes can also be applied successfully to approximate fractional Laplacian with generalized homogeneous boundary conditions, whose fractional order α∈(0,2), not only having to be limited to α∈(1,2). Several numerical experiments demonstrate the effectiveness of the derived schemes. Besides, based on the numerical results, we can observe the behavior of mean first exit time, an interesting quantity that can provide us with a further understanding about the mechanism of abnormal diffusion.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

In this paper, we target on investigating a well-condtioned Galerkin spectral methods for the following two-sided fractional diffusion equation with drift

 (1.1) {−(¯p⋅Dαu(x)+¯q⋅Dα∗u(x))+¯d⋅Du(x)=h(x),  x∈(a,b),u(a)=u(b)=0,

where , satisfying , and , are neither the Riemann-Liouville operators nor the Caputo ones ; rather, in general, for ,

 (1.2) Dγu(x):=DaIn−γxDn−1u(x),
 (1.3) Dγ∗u(x):=(−1)nDxIn−γbDn−1u(x),

with and , , denote separately the left Riemann-Liouville fractional integral

 aIβxu(x)=1Γ(β)∫xa(x−s)β−1u(s)ds,

and the right Riemann-Liouville fractional integral operator

 xIβbu(x)=1Γ(β)∫bx(s−x)β−1u(s)ds.

Here presents the Euler gamma function.

The fractional Dirichlet problem and variants thereof appear in many applications, in particular in physical settings where anomalous dynamics occur and where the spread of mass grows faster than linearly in time. Examples include turbulent fluids, contaminant transport in fractured rocks, chaotic dynamics and disordered quantum ensembles; see [15, 16, 26]. The authors in  believe that problem (1.1), which can be interpreted as the steady-state equation for a time dependent advection and anomalous diffusion problem, is a more physical model than the corresponding Riemann-Liouville or Caputo fractional equation. During the derivation of Eq. (1.1), the authors in  point out that besides obeying the conservation of energy principle, the physical interpretation of the flux at a given cross section , is that “there is a nonlocal effect from a flux originating at a cross section , proportional to distance for that point”. In other words, the contribution to the flux at cross section , from points to its left and right is given by

 −k∫xa(x−s)1−α∂u(s,t)∂sds,

and

 −k∫bx(s−x)1−α∂u(s,t)∂sds,

respectively, where is a dispersion coefficient. In this way, when considering the case of Dirichlet boundary conditions , after by changing the unknown , the simulation of the model equation would require the same energy source as for the case , which physically makes sense.

Besides can be viewed as the steady-state equation for a time dependent advection and anomalous diffusion problem, we shall see that when , under the framework of the image spaces of Riemann-Liouville fractional integral operators on space , Problem (1.1) itself can also be used to describe the mean first exit time of a stochastic process never leaving a fixed region in the state space —an interesting deterministic quantity that can provide us with a further understanding about the mechanism of the anomalous diffusion.

Another topic we want to note is that from a mathematical view, under suitable assumptions on , fractional Laplacian operator 

 (1.4) (−△)su(x)=22s−1Γ(s+12)π12Γ(1−s)∫Ru(x+y)+u(x−y)−2u(x)|y|1+2sdy.

is equivalent to

 (−△)α2u(x)=12cos(πα2)(RL−∞Dαx+RLxDα∞)u(x)=12cos(πα2)(D2−∞Iαx+D2xIα∞)u(x),

where . Actually, the proof in  also ensures that

 (−△)α2u(x)=12cos(πα2)(D−∞IαxD+DxIα∞D)u(x).

Therefore, mathematically, the following one-dimensional Poisson problem with generalized Dirichlet boundary condition:

 (1.5) {(−△)α/2u(x)=h(x), x∈Ω=(a,b),u(x)=0, x∈R∖Ω,

can be changed as

 (1.6) {−12(Dαu(x)+D% α∗u(x))=−cos(πα2)h(x),  x∈(a,b),u(a)=u(b)=0,

which is a special case of (1.1), where (not only limited to (1,2) as in (1.1)).

From the view of stochastic processes, the physical meaning of the fractional Laplacian defined in above way with Dirichlet boundary conditions is the negative infinitesimal generator of stopped subordinated Brownian motion (i.e., stopped -stable Lévy motion), which represents particles that are stopped upon exiting the domain via a jump over the boundary [7, 8, 21]. Here, we do not concern the detailed conditions under which (1.5) and (1.6) are equivalent. Instead, we mainly focus on the spectral methods that are effective for them, and leave the theoretical part in our future work.

Comparing with the classical differential equations, one of the big challenges we have to face is the expensiveness of its computation cost besides its complexity, since fractional operators are pseudodifferential operators which are non-local. Finite difference methods and finite elements methods are not easy to apply when solving especially a two-sided fractional problems, because the information on the whole domain is needed which results in a huge computational cost. In this case, spectral method, as a global method, appears to be a natural choice. There are existing spectral work, used to solve one-sided or two-sided fractional differential equations with Riemann-Liouville or variable order fractional operators [5, 17, 20, 30, 31]

. Early spectral collocation methods for fractional problems using classical interpolation basis functions with Legendre-Gauss-Lobatto or Chebyshev-Gauss-Lobatto collocated points are proposed in

 and 

. Eigenfunctions of a fractional Sturm-Liouville operator are derived in

. Spectral approximation results in weighted Sobolev spaces involving fractional derivatives are derived in , including also rigorous convergence analysis. The authors in  introduce fractional Birkhoff interpolation basis functions into collocation methods to reduce the condition numbers when solving the one-sided Caputo fractional equations.

As for the problem (1.1), a variational formulation is studied in , together with a finite element error analysis. The regularity of (1.1) is studied, also a finite elements method and a spectral type approximation method are proposed in it.

As far as we know, there is little literature to discuss the weak formulation of the two-sided fractional diffusion problems with drift, in which the fractional operators are physically well-defined. Also, there has been no relevant work to talk about the corresponding well-conditioned scheme.

This paper mainly proposes three kinds of Galerkin spectral schemes for solving Eq. (1.1). These three Galerkin spectral schemes are based on different weak variational formulations and have different regularity requirements, all of which shows to be effective to this kind of two-sided fractional diffusion equation with drift, even when the solution has a low regularity. In special, based on the former two formulations, the third one, named as mixed Galerkin spectral formulation, is designed by splitting the Eq. (1.1) into three subequations. In this way, the trial and test functions are more flexible to choose, so that the coefficient matrices can be expressed in a simpler way. Besides, compared with the condition numbers of the stiff matrices in the other two schemes, the condition number in mixed Galerkin spectral scheme can be reduced from to about , where is the degree of the polynomials used in the approximation.

The rest of this paper is organized as follows. Section 2 reviews some important definitions and results about the image spaces of Riemann-Liouville fractional integrals on Sobolev space space, which are the framework of the weak formulation in this paper. Three different weak formulations and Galerkin spectral methods are presented step by step in Sections 3, where the differences among them are discussed. Section 4 provides the numerical results for solving problems (1.1) and (1.5), in which one can observed that the condition numbers are substantially decreased in the mixed Galerkin spectral method. Finally, the main results are summarized in Section 5.

## 2. Preliminaries

In this section, we outline the definition and some results about the image spaces of -order Riemann-Liouville fractional integral operators on or , which is called “spaces of fractional integrals” for short , where is a given classical Sobolev space and .

As we all know that the concept of fractional calculus is almost as old as their more familiar integer order counterparts, and many mathematical results about fractional operators are also discussed in the early days [23, 22, 24, 25]. Until recently, fractional derivatives have been widely and successfully explored as a tool for developing more sophisticated mathematical models. Here, we borrow (not simply copy but sometimes have to flip through pages) the space, which we call as the image space of Riemann-Liouville fractional integral operators on space, introduced in , and some results given in  and , to begin our discussion. The reason we choose this kind of space, not only because it comes from a “non state of the art” references, but also because the key difficulty of the fractional operators that are widely used, such as Riemann-Liouville derivative or Caputo deriavative, are actually come from the pseudo-differential or Riemann-Liouville fractional integral operator in them. Since the space of fractional integrals of functions can catch this characteristic very well, it is a natural way to begin our discussion from it.

Denote as space on . The set of -th order left and right Riemann-Liouville fractional integrals of functions, , are firstly given in Definition 2.3 of . We rearrange them as follows:

###### Definition 2.1.
 (2.1) Iα[Lp(Ω)]:={f:f(x)=aIαxφ(x),φ(x)∈Lp(Ω),x∈Ω}, α>0,

and

 (2.2) Iα∗[Lp(Ω)]:={f:f(x)=xIαbφ(x),φ(x)∈Lp(Ω),x∈Ω}, α>0.

Now we only list some results about ; similar results can be derived for .

In , Corollary 2.10 shows that actually if , then . Therefore, the following lemmas hold .

###### Lemma 2.2.

If , , then

 (2.3) aIαxDαu(x)=u(x).
###### Lemma 2.3.

Let . If , , then

 (2.4) (Dαu(x),v(x))=(u(x),Dα∗v(x)).
###### Lemma 2.4.

Let , , . If , then

 (2.5) Dαu(x)=Dα1Dα2u(x),

and

 (2.6) Dα2u(x)∈Iα1[Lp(Ω)].

If , then there exists a unique [25, 34], such that . Using Lemma 2.2, we have

 (2.7) ∫bau(x)⋅Dα∗ϕ(x)dx=∫baaIαxφ(x)⋅Dα∗ϕ(x)dx=∫baφ(x)⋅xIαbDα∗ϕ(x)dx=∫baφ(x)ϕ(x)dx∀ϕ(x)∈C∞c(Ω),

where the integrals make sense because of the Hölder inequality .

Because , , so, Eqs. (2.3)-(2.7) still hold for , . Therefore, we can say that is a Sobolev space.

Since for , (Theorem 2.6 in ), i.e., (similar to Poincaré inequality )

 ∥aIαxφ(x)∥p≤(b−a)αΓ(α+1)∥φ(x)∥p   ∀φ(x)∈Lp(Ω).

We can introduce the norm in by

###### Remark 2.5.

In the later sections, we can see that actually, for , , , functions and belong to and , respectively, where denote the Jacobi polynomials, which are defined by Rodrigues’ formula

 (1−x)σ(1+x)ηJσ,ηn(x)=(−1)n2nn!dndxn[(1−x)n+σ(1+x)n+η],

and they are orthogonal on with respect to when , .

Next, the Sobolev space with higher regularity can be defined :

###### Definition 2.6.

The image space of -th order left Riemann-Liouville fractional integrals on is defined as

 (2.9) Iα[Wm,p(Ω)]:={f:f(x)=aIαxφ(x),φ(x)∈Wm,p(Ω),x∈Ω},

and with norm

 ∥f(x)∥Iα[Wm,p(Ω)]:=∥Dαxf(x)∥Wm,p(Ω),

where is a given classical integer Sobolev space.

The relationships between the image spaces and are briefly listed in the following lemmas, in which besides the case , the most interested cases is when and .

###### Lemma 2.7.

 When , and , then

 (2.10) Hα,p(Ω)=^Iα[Lp(Ω)]:=Iα[Lp(Ω)]=Iα∗[Lp(Ω)].

When , then

 (2.11) Hα,p0(Ω)=Iα[Lp(Ω)]∩Iα∗[Lp(Ω)],

where

 Hα,p0(Ω)={f:f(x)∈Hα,p(Ω),~{}and~% {}f(a)=f(b)=0},
 Hα,p(Ω)={f:∃ g(x)∈Hα,p(R), s.t. g(x)∣∣Ω=f(x)},
 Hα,p(R)={f(x)∈Lp(R):F−1[(1+|ξ|2)α2F[f]]∈Lp(R)}.
###### Lemma 2.8.

 If , then

 (2.12) Iα[Hm(Ω)]∩Iα∗[Hm(Ω)] = {f:f(x)∈Wm,q(Ω),f(x)=o((x−a)m+α−12), as x→a, f(x)=o((b−x)m+α−12), as x→b},  q=21−2α.

If , then

 (2.13) Iα[Hm(Ω)]∩Iα∗[Hm(Ω)] = {f:f(x)∈Wm+1,q(Ω),f(x)=o((x−a)m+α−12), as x→a, f(x)=o((b−x)m+α−12), as x→b},  q=23−2α.

Denote as the polynomials spaces of degree less than or equal to on . Then is a subspace of .

Denote as the orthogonal projection operator from onto . Then the following approximation property holds:

###### Lemma 2.9.

 If , and , then there exists a constant , such that

where .

## 3. Variational formulations and spectral methods

We use the spaces of fractional integrals introduced above to design Galerkin spectral methods for solving problem (1.1). Without loss of generality, we now restrict our attention to the interval .

### 3.1. Variational formulations

In order to derive a variational form of (1.1

), we firstly assume for the moment that

is a sufficiently smooth solution. By multiplying an arbitrary , it can be obtained that

 (3.1) ∫1−1−(¯p⋅Dαu(x)+¯q⋅Dα∗u(x))⋅v(x)dx+¯dDu(x)⋅v(x)dx=∫1−1h(x)v(x)dx.

#### 3.1.1. Variational formulation I

Taking integration by parts for the left hand of (3.1), and noting that for smooth with by the definition of in(1.2), we can obtain

 (3.2) −¯p∫1−1Dα2u(x)⋅Dα2∗v(x)dx−¯q∫1−1Dα2∗u(x)⋅Dα2v(x)dx+¯d∫1−1D12u(x)⋅D12∗v(x)dx=∫1−1h(x)v(x)dx.

Denote

 Φα21(Ω):=Iα2[L2(Ω)]∩Iα2∗[L2(Ω)].

Now we define the associated bilinear form for (1.1) as

 (3.3) B1(u,v):=−¯p(Dα2u,Dα2∗v)−¯q(Dα2∗u,Dα2v)+¯d⋅(D12u,D12∗v).

For a given function , which belongs to the dual space of , and be denoted as , where , , we define the associated linear functional as

 (3.4) F1(v):=⟨h,v⟩,

where is the duality pair of and .

By Lemma 2.4 and formula (2.13) in Lemma 2.8, we can check that both (3.3) and (3.4) make sense.

Thus, the corresponding variational formulation of (1.1) can be defined as follows.

###### Definition 3.1 (Variational Formulation I).

A function is a variational solution of problem (1.1) provided

 (3.5) B1(u,v)=F1(v)∀v(x)∈Φα21(Ω).

Denote

 Φα21,N(Ω)=Iα2[PN(Ω)]∩Iα2∗[PN(Ω)].

Then the Galerkin approximation of (3.5) is: find , such that

 (3.6) B1(u1,N,v1,N)=F1(v1,N)∀v1,N(x)∈Φα21,N(Ω).

#### 3.1.2. Variational formulation II

Actually, for smooth solution with , and an arbitrary given , instead of Eq. (3.2), we can get another formula by taking integration by part for the left side of Eq. (3.1), as follows:

 (3.7) −¯p∫1−1Dα−12u(x)⋅Dα+12∗v(x)dx−¯q∫1−1D% α−12∗u(x)⋅Dα+12v(x)dx            −¯d∫1−1u(x)⋅Dv(x)dx=∫1−1h(x)v(x)dx.

Denote

 Φα−122(Ω):={f:f∈^Iα−12[L2(Ω)],~{}and~{}f(−1)=f(1)=0}.

We now define another type of bilinear form for (1.1) as

 (3.8) B2(u,v):=−¯p(Dα−12u,Dα+12∗v)−¯q(Dα−12∗u,Dα+12v)−¯d⋅(u,Dv).

For a given source term , which belongs to the dual space of , and be denoted as , where , , we define the associated linear functional as

 (3.9) F2(v):=⟨h,v⟩,

where is the duality pair of and .

By Lemma 2.4, formula (2.11) in Lemma 2.7, and formula (2.12) in Lemma 2.8, we can check that both (3.8) and (3.9) make sense.

Thus, the corresponding variational formulation of (1.1) can be defined as follows.

###### Definition 3.2 (Variational Formulation II).

A function is a variational solution of problem (1.1) provided

 (3.10) B2(u,v)=F2(v)∀v(x)∈Φα+121(Ω).
###### Remark 3.3.

It is not difficulty to see that the weak solution as well as the linear functional in (3.10) lie in weaker spaces than the weak solution and the linear functional of (3.5) do; the classical solution can be recovered from both (3.5) and (3.10) if is smooth enough.

Denote

 Φα−122,N(Ω):=^Iα−12[PN(Ω)].

We can see that if , then .

The Galerkin approximation of (3.10) is: find , such that

 (3.11) B2(u2,N,v2,N)=F2(v2,N)∀v2,N(x)∈Φα+121,N(Ω).

#### 3.1.3. Variational formulation III

Since by Lemma 2.7, when , , it is not simple to manipulate during the numerical realization. One way to get rid of using it during the computation, is based on the following splitting formula, which is equivalent to problem (1.1):

 (3.12) ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩l(x)=¯p −1I2−αxDu(x),r(x)=¯q xI2−α1Du(x),−D[l(x)+r(x)]+¯d⋅Du(x)=h(x),u(−1)=u(1)=0.

Similarly to the above discussions, by assuming for the moment that is a sufficiently smooth solution, then multiplying the first three equalities of (3.12) separately by , , , and taking integration by parts, we can get

 (3.13) ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩∫1−1l(x)ψ1(x)dx=¯p ∫1−1Dα−12u(x)⋅% Dα−12∗ψ1(x)dx,∫1−1r(x)ψ2(x)dx=¯q ∫1−1Dα−12∗u(x)⋅Dα−12ψ2(x)dx,∫1−1[l(x)+r(x)]Dψ3(x)dx−¯d⋅∫1−1u(x)Dψ3(x)dx=∫1−1h(x)ψ3(x)dx.

If , and belong to , then and belong to , , and , where , . In this case, (3.13) is actually the same as Variational Formulation II (3.2).

For the convenience of computation and implementation, we partially yield to the requirements on the regularity in (3.7). Specifically, we define the third type of mixed variational formulation of (1.1) in the following way:

###### Definition 3.4 (Variational Formulation III).

Find , , and , such that

 (3.14) ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩(l(x),ψ1(x))−¯p (Dα−12u(x),Dα−12∗ψ1(x))=0∀ψ1(x)∈^Iα−12[L2(−1,1)],(r(x),ψ2(x))−¯q (Dα−12∗u(x),Dα−12ψ2(x))=0∀ψ2(x)∈^Iα−12[L2(−1,1)],(l(x)+r(x),Dψ3(x))−¯d⋅(u(x),Dψ3(x))=F3(ψ3)∀ψ3(x)∈H10(−1,1),

where for a given , is a linear functional defined as

 (3.15) F3(v):=⟨h,v⟩,

and is the duality pair of and .

We can see that the main difference between the weak formulae (3.10) and (3.14) is that the linear functional of later one lies in a bit smaller space than that of former one , and as a sequence, the weak solution of (3.14) lies in a smaller space.

Denote

 Ψ3,N(Ω)={f:f∈PN(Ω),~{}and~{}f(−1)=f(1)=0}.

Then the Galerkin approximation of (3.14) is: find , , and , such that

 (3.16) ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩(lN,ψ1,N)−¯p (Dα−12u3,N,Dα−12∗ψ1,N)=0∀ψ1,N∈^Iα−12[PN(Ω)],(rN,ψ2,N)−¯q (Dα−12∗u3,N,Dα−12ψ2,N)=0∀ψ2,N∈^Iα−12[PN(Ω)],(lN+rN,Dψ3,N)−¯d⋅(u3,N,Dψ3,N)=F3(ψ3,N)∀ψ3,N∈Ψ