    # Optimal Petrov-Galerkin spectral approximation method for the fractional diffusion, advection, reaction equation on a bounded interval

In this paper we investigate the numerical approximation of the fractional diffusion, advection, reaction equation on a bounded interval. Recently the explicit form of the solution to this equation was obtained. Using the explicit form of the boundary behavior of the solution and Jacobi polynomials, a Petrov-Galerkin approximation scheme is proposed and analyzed. Numerical experiments are presented which support the theoretical results, and demonstrate the accuracy and optimal convergence of the approximation method.

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

Of interest in this paper is the approximation of the solution to the fractional diffusion, advection, reaction equation

 Lαru(x) + b(x)Du(x) + c(x)u(x) = f(x),  x∈I, (1.1) subject to u(0)=u(1) =0, (1.2) where  Lαru(x) := −D(rD−(2−α) + (1−r)D−(2−α)∗)Du(x), (1.3)

and , , , , denotes the usual derivative operator, the -order left fractional derivative operator, and the -order right fractional derivative operator, defined by:

 Dαu(x) := D0D−(2−α)xDu(x) = D1Γ(2−α)∫x01(x−s)α−1Du(s)ds, (1.4) Dα∗u(x) := DxD−(2−α)1Du(x) = D1Γ(2−α)∫1x1(s−x)α−1Du(s)ds. (1.5)

In recent years fractional differential equations have received increased attention as they have been used in modeling a number of physical phenomena such as contaminant transport in ground water flow , viscoelasticity , image processing [7, 15], turbulent flow [28, 34], and chaotic dynamics .

The are two important properties that distinguish a fractional order differential equations from its integer order counterpart. Firstly, as can be noted from (1.3), fractional differential equations are nonlocal in nature. Secondly, the solution of fractional differential equations (typically) have a lack of regularity at the boundary of the domain. Finite difference methods [10, 26, 33, 36, 37], finite element methods [14, 22, 27, 38], discontinuous Galerkin methods , and mixed methods [8, 25], have all been developed for fractional differential equations. These methods typically exhibit slow convergence due to the lack of regularity of the solution at the boundary. In [21, 23] an enriched subspace was given for one sided fractional differential equations, where the boundary behavior of the solution was included in the finite element trial space. Mao and Shen in  extended the work of Gui and Babuška in  to establish that for an assumed boundary behavior of the solution a geometrically spaced mesh with increasing polynomial degree trial function on the subintervals resulted in an exponential rate of converge for the approximation. For a special class of self-adjoint fractional differential equations a spectral approximation scheme was presented in  using a special class of functions, polyfractonomials. Spectral methods, exploiting a special property satisfied by fractional diffusion operator applied to Jacobi polynomials (see (2.16)) has been particularly effective for the approximation of the solution to fractional diffusion equations [9, 13, 24, 29, 31, 30, 42, 43].

Two recent papers have established the explicit form of solutions to fractional diffusion, advection, reaction equations on a bounded domain in . In , Hao and Zhang studied the case for , for which is a symmetric operator. Their work was extended in  to the general case of . The solution was shown to have the form , where is contained in the weighted Sobolev space (defined in Section 2), where and are explicit functions of , and the regularity of the right hand side function, (see Theorems 2.2 and 2.3 below). Of particular note is that for the fractional diffusion, reaction problem, and the fractional diffusion, advection, reaction problem, the regularity of the solution is bounded, regardless of the regularity of . This boundedness in the regularity of is not the case for the fractional diffusion, advection, reaction equation on , as was recently established by Ginting and Li in .

The numerical approximation scheme presented below is accurate as, using , the precise boundary behavior of the solution is incorporated into the approximate solution. Additionally, using the special property of the fractional diffusion operator applied to Jacobi polynomials (see (2.16))

 Lαrω(x)ˆG(α−β,β)k(x) = λkˆG(β,α−β)k(x),

and that is a basis for , the approximation scheme using Jacobi polynomial is efficient in that if the solution is (very rarely the case) the approximation converges exponentially. If the solution has bounded regularity (typically the case) the approximation converges optimally at an algebraic rate of convergence.

This paper is organized as follows. In the following section definitions, notation, and several known results are summarized. Section 3 contains the Petrov-Galerkin weak formulation for (1.1),(1.2), and establishes the existence and uniqueness of its solution. The analysis follows the work of Jin, Lazarov and Zhou in , wherein the lower order terms are handled using the Petree-Tartar Lemma. The approximation scheme is given in Section 4

, and associated error estimates derived. Numerical experiments are presented in Section

5.

## 2 Notation and Properties

Jacobi polynomials have an important connection with fractional order diffusion equations [2, 13, 30, 29]. We briefly review their definition and some of their important properties [1, 35].

Usual Jacobi Polynomials, , on .
Definition: , where

 pn,m := 12n(n+am)(n+bn−m). (2.1)

Orthogonality:

 ∫1−1(1−t)a(1+t)bP(a,b)j(t)P(a,b)k(t)dt = {0,k≠j|∥P(a,b)j|∥2,k=j, where   |∥P(a,b)j|∥ = (2(a+b+1)(2j+a+b+1)Γ(j+a+1)Γ(j+b+1)Γ(j+1)Γ(j+a+b+1))1/2. (2.2)

In order to transform the domain of the family of Jacobi polynomials to , let and introduce . From (2.2),

 ∫1−1(1−t)a(1+t)bP(a,b)j(t)P(a,b)k(t)dt = ∫102a(1−x)a2bxbP(a,b)j(2x−1)P(a,b)k(2x−1)2dx = 2a+b+1∫10(1−x)axbG(a,b)j(x)G(a,b)k(x)dx = {0,k≠j,2a+b+1|∥G(a,b)j|∥2,k=j. where   |∥G(a,b)j|∥ = (1(2j+a+b+1)Γ(j+a+1)Γ(j+b+1)Γ(j+1)Γ(j+a+b+1))1/2. (2.3)

From [29, equation (2.19)] we have that

 dkdtkP(a,b)n(t) = Γ(n+k+a+b+1)2kΓ(n+a+b+1)P(a+k,b+k)n−k(t). (2.4)

Hence,

 dkdxkG(a,b)n(x) = Γ(n+k+a+b+1)Γ(n+a+b+1)G(a+k,b+k)n−k(x). (2.5)

Note that, from Stirling’s formula, we have that

 limn→∞Γ(n+σ)Γ(n)nσ = 1, for σ∈R. (2.6)

For compactness of notation, let

 ω(a,b)=ω(a,b)(x):=(1−x)axb. (2.7)

We let and use to denote that there exists constants and such that, as , . Additionally, we use to denote that there exists a constant such that .

For , is used to denote the largest integer that is less than or equal to , and is used to denote the smallest integer that is greater than or equal to .

Function space .
For , let

 L2σ(I):={f(x):∫10σ(x)f(x)2dx < ∞}. (2.8)

Associated with is the inner product, , and norm, , defined by

 (f,g)σ:=∫10σ(x)f(x)g(x)dx,and% ∥f∥σ:=(⟨f,f⟩σ)1/2.

The set of orthogonal polynomials form an orthogonal basis for , and for , form an orthonormal basis for .

Without a subscript, denotes the usual inner product.

Function space .
The weighted Sobolev spaces differ from the usual spaces in that the associated norms apply a polynomial weight at each endpont of , namely, and . These weights increase with the order of the derivative. We give two equivalent definitions for the spaces. In the first definition the spaces , for , are defined by the

- method of interpolation. The second definition is based on the decay rate of the Jacobi coefficients of a function expanded in terms of the Jacobi polynomials

. Both definitions are useful, and used in the analysis below. The equivalence of the spaces is discussed in .

Definition: Using Interpolation
Following Babuška and Guo , and Guo and Wang , we introduce the weighted Sobolev spaces .

###### Definition 2.1

Let , , . Then

 Hsω(a,b)(I):={v:∥v∥2s,ω(a,b):=s∑j=0∥∥Djv∥∥2ω(a+j,b+j)<∞}. (2.9)

Definition (2.9) is extended to using the - method of interpolation. For the spaces are defined by (weighted) duality.

Definition: Using the decay rate of Jacobi coefficients
Next we define function spaces in terms of the decay property of the Jacobi coefficients of their member functions.

Given , let

 vj = ∫10ω(a,b)(x)v(x)ˆG(a,b)j(x)dx. (2.10)

Note that for ,

 v(x) = ∞∑j=0vjˆG(a,b)j(x). (2.11)
###### Definition 2.2

Let , , , and be given by (2.10). Then, define

 Hs(a,b)(I):={v:∞∑j=0(1+j2)sv2j<∞} (2.12)

as the -weighted Sobolev space of order .

###### Theorem 2.1

[12, Theorem 4.1] The spaces and coincide, and their corresponding norms are equivalent.

With the structure of the spaces, and properties (2.5) and (2.3), it is straight forward to show that is a bounded mapping from onto .

###### Lemma 2.1

[12, Lemma 4.5] For , , the differential operator is a bounded mapping from onto .

For convenience, from hereon we use to represent the spaces and .

Definition: Condition A
The parameters , , and and constant satisfy: , ,

 c∗∗ = sin(πα)sin(π(α−β))+sin(πβ), (2.13)

where is determined by

 r = sin(πβ)sin(π(α−β))+sin(πβ). (2.14)

For compactness of notation, for and defined in (1.1) and defined in (2.14) we introduce

 ω(x):=ω(α−β,β)(x)=(1−x)α−βxβ,  and  ω∗(x):=ω(β,α−β)(x)=(1−x)βxα−β. (2.15)

Additionally, we use to denote the weighted duality pairing between functions if and .

From [13, 20],

 Lαrω(x)ˆG(α−β,β)k(x) = λkˆG(β,α−β)k(x),  % where   λk = −c∗∗Γ(k+1+α)Γ(k+1), k=0,1,2,…, (2.16)

and given by (2.13). Also, using (2.6), .

Let denote the space of polynomials of degree less than or equal to . We define the weighted orthogonal projection by the condition

 (v−PNv , ϕN)ω = 0,  ∀ϕN∈SN. (2.17)

Note that , where .

###### Lemma 2.2

[18, Theorem 2.1] For and , with , there exists a constant , independent of and such that

 ∥∥v−PNv∥Hμω(I) ≤ CNμ−t∥v∥Htω(I). (2.18)

Remark: In  (2.18) is stated for . The result extends to using interpolation.

The regularity of the solution to (1.1) can be influenced by the regularity of the coefficients and . The following lemma enables us to insulate the influence of these terms.

Introduce the space and its associated norm, defined for , as

 Wk,∞w(I) := {f: (1−x)j/2xj/2Djf(x)∈L∞(I), j=0,1,…,k}, (2.19) ∥f∥Wk,∞w := max0≤j≤k∥(1−x)j/2xj/2Djf(x)∥L∞(I). (2.20)

The subscript denotes the fact that is a weaker space than in that the derivative of functions in may be unbounded at the endpoints of the interval.

###### Lemma 2.3

[12, Lemma 7.1] Let , , , and . For

 (i) g ∈Hs(α,β)(I) then  fg∈Hs(α,β)(I),  and for (2.21) (ii) g ∈H−s(α,β)(I) then  fg∈H−s(α,β)(I). (2.22)
###### Theorem 2.2

[12, Theorem 7.1] Let , be determined by Condition A, satisfying and

 f∈H−α/2(I)∩Hs(β,α−β)(I). (2.23)

Then there exists a unique solution , with , to

 Lαru(x) + c(x)u(x) = f(x), x∈I,  subject to u(0)=u(1)=0. (2.24)

The inclusion of an advection term can significantly reduced the regularity of the solution.

###### Theorem 2.3

[12, Theorem 7.2] Let , be determined by Condition A, satisfying , and

 f∈H−α/2(I)∩Hs(β,α−β)(I). (2.25)

Then there exists a unique solution , with , to

 Lαru(x) + b(x)Du(x) + c(x)u(x) = f(x), x∈I,  subject to u(0)=u(1)=0. (2.26)

Introduce defined by

 ˜s:={min{s,α+(α−β)+1,α+β+1}, if b=0  (see Theorem ???)min{s,α+(α−β)−1,α+β−1}, if b≠0  (see Theorem ???). (2.27)

## 3 Weak Formulation

In place of (1.1), (1.2), we consider the following problem.

Given , and and satisfying the hypothesis of Theorem 2.3, determine such that satisfies

 ⟨Lαru + bDu + cu,ψ⟩ω∗ = ⟨f,ψ⟩ω∗,  ∀ψ∈Hα/2ω∗(I). (3.1)

Note that the formulation (3.1) has different test and trial spaces. With this in mind we recall the Banach-Nečas-Babuška theorem.

###### Theorem 3.1

[11, Pg. 85, Theorem 2.6] Let and denote two real Hilbert spaces, a bilinear form, and a bounded linear functional on . Suppose there are constants and such that

 (i) |B(w,v)|≤C1∥w∥H1∥v∥H2,  for all w∈H1, v∈H2, (3.2) (ii) sup0≠v∈H2|B(w,v)|∥v∥H2 ≥ C2∥w∥H1, for all w∈H1, (3.3) (iii) supw∈H1|B(w,v)| > 0, for all v∈H2, v≠0. (3.4)

Then there exists a unique solution satisfying for all . Further, .

For , and and satisfying the hypothesis of Theorem 2.3, let , and be defined by

 B(ϕ,ψ) := ⟨Lαru + bDu + cu,ψ⟩ω∗, (3.5) F(ψ) := ⟨f,ψ⟩ω∗. (3.6)

### 3.1 Continuity of B(⋅,⋅)

In order to establish that is well defined and continuous we need to determine which space lies in.

The space a function lies in is determined by its behavior at: (i) the left endpoint (), (ii) the right endpoint (), and (iii) away from the endpoints. In order to separate the consideration of the endpoint behaviors, following , we introduce the following function space . Let , and

 Λ∗ := {(x,y):23x

Introduce the semi-norm and norm

 |f|2Hs(γ)(J) := ∬Λxγ+s|D⌊s⌋f(x)−D⌊s⌋f(y)|2|x−y|1+2(s−⌊s⌋)dydx + ∬Λ1xγ+s|D⌊s⌋f(x)−D⌊s⌋f(y)|2|x−y|1+2(s−⌊s⌋)dydx :=|f|2Hs(γ)(Λ) + |f|2Hs(γ)(Λ1), and  ∥f∥2Hs(γ)(J) := ⎧⎪ ⎪⎨⎪ ⎪⎩∑sj=0∥Djf∥2L2(γ+j)(J), for s∈N0∑⌊s⌋j=0∥Djf∥2L2(γ+j)(J) + |f|2Hs(γ)(J), for s∈R+∖N0, where  ∥g∥2L2(γ)(J) := ∫Jxγg2(x)dx.

Then, .

Note: A function is in if and only if and .

From  we have the following theorem.

###### Theorem 3.2

[12, Theorem 6.4] Let , , , , and . Then provided

 0≤t≤s,  σ+2p≥μ,  σ+2p−t>−1,   and   σ+2p+t≥μ+s. (3.7)

Additionally, when (3.7) is satisfied, there exists (independent of ) such that .

###### Lemma 3.1

The terms , and are well defined. Additionally, there exists such that for and

 |B(ϕ,ψ)| = |⟨Lαrωϕ+bDωϕ+cωϕ,ψ⟩ω∗| ≤ C∥ϕ∥Hα/2ω(I)∥ψ∥Hα/2ω∗(I). (3.8)

Proof: We begin by considering the term.

From Theorem 3.2, with , , , and choosing we have that . Hence for , , with .

Again, using Theorem 3.2, with , , , and choosing we have that . Hence for , with .

Combining the above two applications of Theorem 3.2 we have that for , with

 ∥ωϕ∥Hα/2(β−1,α−β−1)(I) ≲ ∥ϕ∥Hα/2ω(I). (3.9)

A similar application of Theorem 3.2 establishes that for , with

 ∥ωϕ∥Hα/2ω∗(I) ≲ ∥ϕ∥Hα/2ω(I). (3.10)

From (3.9) and Lemma 2.1 we have that with