Numerical simulation method of lines together with a pseudospectral method for solving space-time partial differential equations with space left- and right-sided fractional derivative

Abstract

Purpose

This study aims to use new formula derived based on the shifted Jacobi functions have been defined and some theorems of the left- and right-sided fractional derivative for them have been presented.

Design/methodology/approach

In this article, the authors apply the method of lines (MOL) together with the pseudospectral method for solving space-time partial differential equations with space left- and right-sided fractional derivative (SFPDEs). Then, using the collocation nodes to reduce the SFPDEs to the system of ordinary differential equations, which can be solved by the ode45 MATLAB toolbox.

Findings

Applying the MOL method together with the pseudospectral discretization method converts the space-dependent on fractional partial differential equations to the system of ordinary differential equations.

Originality/value

This paper contributes to gain choosing the shifted Jacobi functions basis with special parameters a, b and give the authors this opportunity to obtain the left- and right-sided fractional differentiation matrices for this basis exactly. The results of the examples are presented in this article. The authors found that the method is efficient and provides accurate results, and the authors found significant implications for success in the science, technology, engineering and mathematics domain.

Citation

Ali, M., Almoaeet, M. and Albuohimad, B.K. (2022), "Numerical simulation method of lines together with a pseudospectral method for solving space-time partial differential equations with space left- and right-sided fractional derivative", Arab Journal of Mathematical Sciences, Vol. 28 No. 2, pp. 168-183. https://doi.org/10.1108/AJMS-02-2021-0052

1. Introduction

The recent development in the last few decades has shown that most of the complex system in engineering and other several phenomena can be accurately modeled using partial differential equations with fractional order. This contributed to a great development in several areas such as biotechnology, chemistry, signal and image processing, finance and many others [1–5]. The main aim of this paper is to introduce an efficient numerical method to approximate the fractional partial differential equation (FPDE) of the form

(1)∂u(x,t)∂t=s(x,t)+c+(x,t)0RDxαu(x,t)+c−(x,t)xRDℓαu(x,t),

we also assume the initial and Dirichlet boundary conditions:

(2)u(x,0)=F(x),

(3)u(0,t)=u(ℓ,t)=0,

on the domain of space 0 < x < l and time 0 ≤ t ≤ T, and we consider that the parameter α is the fractional order where 1 < α ≤ 2. And a source or sink term is the function s(x, t). The functions c₊(x, t) ≥ 0 and c−(x,t)≥0 represented the interpreted as transport-related coefficients or the advection, and the diffusion coefficients.

Many approximation methods in the numerical analysis have been a survey to solve space-dependent on fractional partial differential equations (SFPDEs), and the target of the main subject of these methods in terms of convergence to real solutions, the stability of methods, and order of accuracy and value of error.

One of the best types of methods for solving SFPDEs numerically is by discretization of the space variable without the time variable. These kinds of method are referred to as method of lines (MOL). In this method, the spatial dimensions of the space variable can be discretized by using diverse techniques such as mesh methods, or meshless methods [6–12]. In general, these methods convert the SFPDEs to a system of ordinary differential equations, or differential-algebraic equations based on the type of boundary conditions [13]. There are many types of partial differential problems that have been solved by the MOL, we refer to [7, 10, 14–19], and therein. Also, many different methods have been discussed for SFPDEs, we refer to [20–24], and therein.

Our target of this work is to use the advantage of the pseudospectral method based on the shifted Jacobi functions together with MOL and employed the collocation method to approximate the solution to the SFPDE (1)–(3).

The structure of this paper was arranged in the following way: In Section 2, preliminaries, the definitions of fractional derivatives and some notations of Jacobi–Gauss (JG) nodes. In Section 3, the new numerical technique for solving SFPDE (1)–(3) is presented. In Section 4, the illustrative examples were included to demonstrate the validity and applicability of the proposed method. In Section 5, a brief conclusion and some remarks.

2. Preliminaries and definitions

We introduce several important basic definitions and properties of fractional Riemann–Liouville integrals and derivatives and JG nodes.

Definition 2.1.

(Left and right Riemann–Liouville fractional integral).

Let α is real number where 0 ≤ α ≤ 1 and g: I → R is a continuous function and where the bounded interval I = [a, b], then the left and right Riemann–Liouville fractional integrals of order α are defined as:

IxαaRg(x)=1Γ(α)∫ax(x−t)α−1g(t)dt,x∈I,

IbαxRg(x)=1Γ(α)∫xb(t−x)α−1g(t)dt,x∈I,

respectively.

Definition 2.2.

(Left and right Riemann-Liouville fractional derivative).

If n − 1 ≤ α < n such that n is a positive integer number and the continuous function g: I → R then the left and right Riemann-Liouville fractional derivatives of order α are defined by:

Dxαag(x)=dndxnIx(n−α)aRg(x)=1Γ(n−α)dndxn∫ax(x−t)n−α−1g(t)dt,

Dbαxg(x)=(−1)ndndxnIx(n−α)bRg(x)=(−1)nΓ(n−α)dndxn∫xb(t−x)n−α−1g(t)dt.

We can verify that the following properties hold true.

Dxαa(x−a)r=Γ(r+1)Γ(r+1−α)(x−a)r−α,

and

Daαx(b−x)r=Γ(r+1)Γ(r+1−α)(b−x)r−α,

where r is any real number. For more information see [25–28].

Definition 2.3.

The fractional Riesz operator of the order α is defined as [29]:

∂α∂xαu(x)=−12 cosπα2Dxα0u(x)+Dlαxu(x),

where x belongs to the finite interval 0 ≤ x ≤ l and n − 1 < α ≤ n such that α ≠ 1.

2.1 Jacobi–Gauss nodes

The Jacobi polynomials Pj(a,b)(x) have been applied in a wide range of engineering disciplines and used in system analysis, optimal control, numerical analysis, signal analysis for representation solution of problems [26]. The basis {Pj(a,b)(x)}j=0m,a,b>−1,x∈[−1,1] are orthogonal with respect to weight function w^a,b(x) = (1 − x)^a(1 + x)^b as follows:

(4)∫−11Pj(a,b)(ξ)Pi(a,b)(ξ)w(a,b)(ξ)dξ=γja,bδij,

where

(5)γja,b=2a+b+1Γ(j+a+1)Γ(j+b+1)(2j+a+b+1).j!.Γ(j+a+b+1).

The JG quadrature formula

(6)∫−11g(τ)wa,b(τ)dτ=∑i=0n−1wig(τi),

is exact for any polynomials g ∈ P_2n−1 = span{1,τ,…τ2n−1}, where

wi(a,b)=2a+b+1Γ(n+a+1)Γ(n+b+1)n!Γ(n+a+b+1)(1−τi2)[(Pn(a,b))′(τi)],

and τ₀, τ₁, …, τ_n−1 are the JG nodes to zeros of the classical Jacobi polynomials Pn(a,b)(τ), see [26].

According to the definition of Jacobi polynomials, various types of functions with fractional order have been constructed based on them. For more information see [26].

We insert some properties of the classical Jacobi polynomials in the following:

(7)Pk(a,b−1)(2u−1)=k+a+b2k+a+bPk(a,b)(2u−1)+k+a2k+a+bPk−1(a,b)(2u−1),

(8)Pk(a−1,b)(2u−1)=k+a+b2k+a+bPk(a,b)(2u−1)−k+b2k+a+bPk−1(a,b)(2u−1).

And the following theorem plays a cornerstone to establish our proposed method.

Theorem 1.

Let α, a, b are real numbers with conditions − 1 < b − α < b and − 1 > a, b > 0. Then for all x belong to the interval [0, 1] we have:

(9) Dxα0RxbPk(a,b)2x−1=Γ(k+b+1)Γ(k+b−α+1)xb−αPk(a+α,b−α)2x−1.

It is clear from combined properties and theorems that have been inserted as it hold the proof of theorem in References [30–32]. And there are many applications of the above theorem that can be found in References [33–40].

3. Th presented pseudospectral discretization method

The important step of our method includes to discretize the space variable of the unknown function u(x, t) that appears in SFPDE (1)–(3). To get started, we use the pseudospectral method based on the shifted Jacobi function of the parameters a = 1, b = 1, together with the collocation JG nodes. Let the domain space of the function u(x, t) belongs to the interval [0, 1], then the JG nodes will be correspondent to the interval [0, 1] as:

(10)τ^i≔12τi(1, 1)+1,i=1,…,n.

Let us begin to approximate the unknown function u(x, t) by ũn(x,t) as follows:

(11)ũn(x,t)=∑j=0naj(t)ϕj(x),

where

(12)ϕj(x)≔x(1−x)τ^j1−τ^jℓj(x),x∈(0,1),

such that ℓ_j(t), j = 1, 2, 3, …, n, are the Lagrange basis polynomials based on the JG nodes τ^ii=1n, then:

(13)ℓj(x)≔∏k=1k≠jnx−τ^kτ^j−τ^k.

From the Kronecker properties in the JG nodes, the ℓ_j(x) and ϕ_j(t) are satisfied as follows:

(14)ϕj(τ^i)=ℓj(τ^i)=0ifi≠j1ifi=j,j=1,2,3…,n.

Let us denote the shifted Jacobi polynomials P^k(a,b)(t), k = 0, 1, …, n, which can be written as:

P^k(a,b)(t)=Pk(a,b)(2t−1),k=0,1,2,3,…,n.

Since the degree of polynomial ℓ_j(x) is n − 1. Then, the expansion of ℓ_j(x) in terms of the shifted Jacobi polynomials P^k−1(1,1)(x), k = 1, …, n, can be written as follows:

(15)ℓj(x)=∑k=1nλ^j,kP^k−1(1,1)(x).

Now, by multiplying both sides of Eqn (15) by x(1−x)P^k−1(1,1)(x), and using the orthogonality property (4) to the shifted Jacobi polynomials in the interval (0, 1), get that

(16)λ^jk≔(2k+3)(k+2)k+1∫01x(1−x)ℓj(x)P^k−1(1,1)(x)dt.

By using the JG quadrature rule (6), and then using the Kronecker property (14), we get

λ^jk≔(2k+3)(k+2)k+1∑i=1nωi(1, 1)ℓj(τ^i)P^k−1(1,1)(τ^i).

(17)λ^jk≔(2k+3)(k+2)k+1ωj(1, 1)P^k−1(1,1)(τ^i).

By using above Eqn (17), then we can define ϕ_j(t) as:

(18)ϕj(x)=∑k=1nλjk(1−x)xP^k−1(1,1)(x),

where

(19)λjk≔(2k+3)(k+2)k+1ωj(1,1)τ^j1−τ^jP^k−1(1,1)(τ^i), wherej=1,2,3,…,n.

Now dependent on the fact of Eqns (19), (18), (15) and (11) the approximation of the function u(x, t) is complete. Also we can write ∂∂tũn(x,t) as follows:

(20)∂∂tũn(x,t)≃∑k=0nȧk(t)ϕk(x).

Similarly, we approximate the initial and boundary conditions as follows:

(21)u(x,0)≃ũn(x,0)=∑k=0nak(0)ϕk(x)=F(x),

(22)u(0,t)≃ũn(0,t)=0,u(ℓ,t)≃ũn(ℓ,t)=0.

By substituting the approximations functions (20), (21) and (22) into the problems (1)–(3), we get:

(23)∑j=0nȧj(t)ϕj(x)=s(x,t)+c+(x,t)∑j=0naj(t) 0Dxαϕj(x)+c−(x,t)∑j=0naj(t)xD 1αϕj(x).

Let τi(1,1),i=0,2,3,…,n be the zeros of P^n+1(1,1)(x). By using collocating nodes at x=τi(1,1),i=0,1,2,3,…,n the above Eqn (23), will change to the system of algebra equations dependent on the time variable:

∑j=0nȧj(t)ϕj(τi(1,1))=s(τi(1,1),t)+c+(τi(1,1),t)∑j=0nDxα0ϕj(τ^i)aj(t)+

(24)c−(τi(1,1),t)∑j=0nD1αxϕj(τ^i)aj(t)

In general, Eqn (24) can be rewritten in the following matrix form:

(25)Pȧ(t)=s(t)+C+(t)D+a(t)+C−(t)D−a(t),

where

P=ϕ0(τ0(1,1))ϕ0(τ1(1,1))…ϕ0(τn(1,1))ϕ1(τ0(1,1))ϕ1(τ1(1,1))…ϕ1(τn(1,1))⋮⋮⋱⋮ϕn(τ0(1,1))ϕn(τ1(1,1))…ϕn(τn(1,1)),

and

a(t)=a0(t)a1(t)⋮an(t),s(t)=s(τ0(1,1),t)s(τ1(1,1),t)⋮s(τn(1,1),t),D+=d00+d01+…d0n+d10+d11+…d1n+⋮⋮⋱⋮dn0+dn1+…dnn+,

D−=d00−d01−…d0n−d10−d11−…d1n−⋮⋮⋱⋮dn0−dn1−…dnn−,

C+(t)=diagnolc+(τ0(1,1),t),c+(τ1(1,1),t),…,c+(τn(1,1),t)and

C−(t)=diagnolc−(τ0(1,1),t),c−(τ1(1,1),t),…,c−(τn(1,1),t).

Collocating the initial condition (21) at x=τi(1,1),i=0,1,2,3,…,n, then we get:

∑j=0naj(0)ϕj(τi(1,1))=f(τi(1,1)),i=0,1,2,3…,n,

or

Pa(0)=F,

where the vector F=[F(τ0(1,1)),F(τ1(1,1)),…,F(τn(1,1))]T and a(0)=[a0(0),a1(0),a2(0),…,an(0)]T. In the final result, the main problems (1)–(3) will be reduced to the following system of ODEs with initial conditions:

(26a)Pȧ(t)=s(t)+C+(t)D+a(t)+C−(t)D−a(t),

(26b)Pa(0)=F.

3.1 On the derivation of the left and right fractional differentiation matrices

The coefficients dij+ and dij−, appear in the system of Eqn (25), we need to compute them accurately and efficiently. In the following, we present an efficient method to compute the coefficients dij+ and dij−. It is interesting to point out that, we choose ϕ_k(x) as the basis functions in the pseudospectral method. The best feature of these bases is that, we can obtain a closed form of the left and right fractional derivatives, by the use of the next two theorems 2 and 3, which may simplify the discretization stage.

Theorem 2.

The left fractional differentiation matrix D+α=[dij+] can be obtained as the element of it in the explicit form as follows:

dij+≔Dxα0ϕj(τ^i)=∑k=1nλjkζk(τ^i),

and

ζk(x)≔Γ(1+k)Γ(k+1−α)x1−αP^k−1(1+α,1−α)(x)−Γ(k+2)Γ(k+2−α)k+22k+1x2−αP^k−1(1+α,2−α)(x)−Γ(k+1)Γ(k+1−α)k2k+1x2−αP^k−2(1+α,2−α)(x),

where τ^i are the shifted JG nodes defined in Eqn. (10).

Proof: Depending on the fact we get in Eqn (19), we have

(27)Dxα0ℓ^j(x)=∑k=1nλjkDxα0((1−x)xP^k−1(1,1)(x)),

where P^k−1(1,1)(x)=Pk−1(1,1)(2x−1) is the shifted Jacobi polynomial, and by using the property (7) with special parameters a = 1 and b = 2 gives

(1−x)xP^k−1(1,1)(x)=xP^k−1(1,1)(x)−x2P^k−1(1,1)(x)=xP^k−1(1,1)(x)−k+22k+1x2P^k−1(1,2)(x)−k2k+1x2P^k−2(1,2)(x).

Using the above equation and Theorem 1, we get

Dxα0((1−x)xP^k−1(1,1)(x))=Γ(k+1)Γ(k+1−α)x1−αP^k−1(1+α,1−α)(x)−Γ(k+2)Γ(k+2−α)k+22k+1x2−αP^k−1(1+α,2−α)(x)−Γ(k+1)Γ(k+1−α)k2k+1x2−αP^k−2(1+α,2−α)(x).

By substituting the above equation in (27), with using the nodes τ^i, i = 1, 2, 3, …, n , the proof of theorem is complete.♠

Theorem 2.

The right fractional differentiation matrix D−α=[dij−] can be obtained as the element of it in the explicit form as:

(28) dij−≔D1αxϕj(τ^i)=∑k=1nλjkζk(1−τ^i),

where

ζk(x)≔Γ(k+1)Γ(k+1−α)x1−αP^k−1(1+α,1−α)(x)−Γ(k+2)Γ(k+2−α)k+22k+1x2−αP^k−1(1+α,2−α)(x)−Γ(k+1)Γ(k+1−α)k2k+1x2−αP^k−2(1+α,2−α)(x).

Proof: Depending on the fact we get it in Eqn (19), and by using property (8) for the special parameter a = 2 and b = 1, we obtain:

x(1−x)P^k−1(1,1)(x)=(1−x)P^k−1(1,1)(x)−(1−x)2P^k−1(1,1)(x)=(1−x)P^k−1(1,1)(x)−k+22k+1(1−x)2P^k−1(2,1)(x)−k2k+1(1−x)2P^k−2(2,1)(x).

Now, combining Theorem 1 with the above equation, we obtain:

D1αx((1−x)xP^k−1(1,1)(x))=Γ(1+k)Γ(k+1−α)(1−x)1−αP^k−1(1−α,1+α)(x)−Γ(k+2)Γ(k+2−α)k+22k+1(1−x)2−αP^k−1(2−α,1+α)(x)−Γ(k+1)Γ(k+1−α)k2k+1(1−x)2−αP^k−2(2−α,1+α)(x)=ζk(1−x).

By substituting the nodes τ^i, i = 1, 2, 3, …, n in the above equation, the proof of theorem is complete.♠

The inverse of matrix P in the explicit form can be found by using the next theorem.

Theorem 3.

The element of the matrix P can be written in the explicit form as follows:

(29) P−1s,j=ws(1,1)τs(1,1)(1−τs(1,1))P^j(1,1)(τs(1,1)),s,j=0,1,…,n,

where ws(1,1),s=0,1,2,…,n are the JG weights with respect to the weight function w^(1,1)(x) = x(1 − x) on [0, 1].

Proof: Start the proof by noting the fact that:

(30)∫01t(1−t)ρkPk(1,1)(2t−1)Pj(1,1)(2t−1)dx=δkj,j,k=0,1,…,n.

On the other hand, we can approximate the above integral by using the JG quadrature rule with respect to the weight function t(1 − t) on [0, 1] as follows:

(31)∫01t(1−x)ρkPk(1,1)(2t−1)Pj(1,1)(2t−1)dt=∑s=0nws(1,1)ρkP^k(1,1)(τs(1,1))P^j(1,1)(τs(1,1)),

j, k = 0, 1, …, n to reach the mass matrix P, we can rewrite the above equation as:

(32)∫01t(1−t)ρkPk(1,1)(2t−1)Pj(1,1)(2t−1)dt=∑s=0nℓ°k(τs(1,1))ws(1,1)τs(1,1)(1−τs(1,1))P^j(1,1)(τs(1,1)),

j, k = 0, 1, …, n. Comparing the left sides of (3.3) and (32), yields:

(33)∑s=0nℓ°k(τs(1,1))ws(1,1)τs(1,1)(1−τs(1,1))P^j(1,1)(τs(1,1))=δkj,j,k=0,1,…,n.

Thus by rewriting the above equation in the matrix form, we have:

(34)PP−1=I,

which concludes the proof.♠

The above theorem shows that the matrix P is nonsingular, and we can rewrite the implicit initial value problem (26) to the following explicit initial value problem:

(35a)ȧ(t)=P−1s(t)+P−1C+(t)D++C−(t)D−a(t),

(35b)a(0)=P−1F.

It should be noted that the IVP (35) can be solved by various well-known software. We solve this IVP by the ode45 MATLAB toolbox. In the next section, we provide some numerical examples to check the efficiency of the proposed method.

4. Numerical results and comparisons

Our goal to solve problems (1)–(3) and find the unknown function u(x, t), by approximate the space variables via pseudospectral method based on JG nodes, and then solving the new IVP by the ode45 MATLAB toolbox. We will begin in the first step by transforming interval (0, ℓ) to interval [0, 1], where the function u(x, t) defined on the interval x ∈ (0, ℓ), t ∈ (0, T), and the collocation points {τ_i, i = 0, 1, 2, …, m} belong to the interval [−1, 1], then x=12[(ℓ)τi+ℓ] are the corresponding collocation points on [0, ℓ]. Also if we transform this interval to [0, 1], then,

x=4τ−3ℓℓ,t∈[0,1].

In order to confirm the utility of the presented method, we apply the method to solve some IVPs. The proposed method has been implemented with MathWorks MATLAB 2017a in a personal computer 3.5 GHz Core i7 PC with 8 GB of RAM.

Example 1.

The following two-sided FPDE has been solved in Reference [41]:

(36) ∂u(x,t)∂t=s(x,t)+c+(x,t)0RDx1.8u(x,t)+c−(x,t)xRDℓ1.8u(x,t),0<x<2,

with initial and boundary conditions

(37) u(x,0)=4x2(2−x)2,u(0,t)=u(2,t)=0,

has the exact solution u(x, t) = 4e^−tx²(2 − x)², where the coefficient functions,

c+(x,t)=Γ(1.2)x1.8,c−(x,t)=Γ(1.2)(2−x)1.8,

and the forcing function

s(x,t)=−32e−tx2+(2−x)2−2.5(x3+(2−x)3)+2522(x4+(2−x)4)−u(x,t).

By applying the presented method with n = 40, the obtained function ũ(x,t) is plotted in Figure 1. Moreover, the error function Ei,j(u)≔u(τ^i,tj)−ũ(τ^i,tj) is plotted in this figure too. From that, this figure can be seen that our numerical solution is in good agreement with the analytic solution. Consumed CPU time for an accurate solution is obtained in just 21.357 seconds.

Example 2.

Consider the following Riesz space fractional diffusion equation [23]:

(38) ∂u(x,t)∂t=∂1+αu(x,t)∂x1+α+s(x,t),

subject to

(39) u(x,0)=x2(1−x)2,0≤x≤1,u(0,t)=u(1,t)=0,0<t≤T,

where 0 < α < 1

s(x,t)=−x2(1−x)2e−t+e−t2 cos(1+α)π224Γ(4−α)[x3−α+(1−x)3−α]−e−t2 cos(1+α)π212Γ(3−α)[x2−α+(1−x)2−α]+2Γ(2−α)[x1−α+(1−x)1−α].

and the exact solution is u(x, t) = x²(1 − x)²e^−t.

By applying the presented method with n = 20 and α = 0.2, the obtained function ũ(x,t) is plotted in Figure 2. Moreover, the error of the obtained function ũ(x,t) is plotted in this figure too. Consumed CPU time for an accurate solution is obtained in just 23.427 seconds. In Table 1, En,m2(u) is the two-norms of the error for the obtained function ũ(x,t), which is defined as,

En,m2(u)≔∑i=1n∑j=1mũ(τ^i,tj)−u(τ^i,tj)212.

From Table 1, it can be seen that norms of errors for various values of n and α appear that the small numbers of n, with that an accurate solution.

Example 3.

We consider the following SFPDEs [29]:

(40) ∂u(x,t)∂t=Kα∂α∂xαu(x,t),0<x≤π,1<α≤2,

with initial and boundary conditions:

(41)u(x,0)=sin(4x),u(0,t)=u(π,t)=0,

Figure 3 shows the approximation function of ũ(x,t), for order of fractional derivative α = 1.25 and ℵ = 1.5 and the time for 0 < t < 1. Observed that if value of α is decreased from 2 to 1 the amplitude of the sinusoidal solution behavior is increased. In Table 2, the values of ũ(x,1) in various values of x and n with the value of α = 1.5

Example 4.

We consider the following FPDEs [42, 43]

(42) ∂u(x,t)∂t+u(x,t)=∂α∂xαu(x,t)+s(x,t),1<α≤2,

with initial and boundary conditions:

(43) u(x,0)=x2(L−x)2,0<x≤L,u(0,t)=u(L,t)=0,0≤t≤T,

s(x,t)=4e−tcosπα23x4−α+(L−x)4−αΓ(5−α)−x3−α+(L−x)3−αΓ(4−α)+x2−α+(L−x)2−αΓ(3−α).

The exact solution is of the form u(x, t) = x²(L − x)²e^−t The problem (42) has been solved in [43] by proposed a class of difference scheme based on the parameter spline function and improved matrix transform method. Also, the authors in [42] have solved this problem (42) by proposing a space-time spectral algorithm based on the shifted Jacobi tau technique. To demonstrate the accuracy of our proposed method, in Table 3, compare the absolute errors |u(x_i, 0.1) − u(x_i, 0.1)|, i = 1, …, 9 with the numerical method proposed in [42, 43] and our results for various choices of α.

5. Conclusion

In the present paper, we developed an efficient and accurate method for solving SFPDE. Applying the MOL together with the pseudospectral discretization method converts the SFPDE to the system of ordinary differential equations. Choosing the shifted Jacobi functions as a test basis with special parameters a, b gives us this opportunity to obtain the left- and right-sided fractional differentiation matrices for this basis exactly.

Four examples have been solved, and the results are reported. These results show that our method is efficient and provides accurate results, whereas a small number of JG nodes are used based on the collocation method. Obtaining some theoretical estimates for the approximation errors would be desirable.

Figures

Figure 1

The approximation function ũ(x,t) (left) and the error function (right) with n = 40 for Example 1

Figure 2

The approximation function ũ(x,t) (left) and the error function (right) with n = 20 and α = 0.2 for Example 2

Figure 3

The approximation function ũ(x,t) with n = 40 and α = 1.25, 1.5 for Example 3

Table 1

Example 2: Norms of errors E_n,100(u) for various values of n and α

Α	E4,1002(u)	E8,1002(u)	E12,1002(u)	E16,1002(u)
α = 1.2	3.4816e − 15	4.1010e − 15	5.7264e − 15	9.4558e − 15
α = 1.5	3.4191e − 15	4.1058e − 15	6.4447e − 15	9.5015e − 15
α = 1.8	2.4772e − 15	4.3510e − 15	7.4810e − 15	1.1811e − 14

Table 2

Example 3: The values of ũ(x,1) in various values of x and n with the value of α = 1.5

Values of n	n = 11	n = 16	n = 21	n = 26	n = 31
ũ(π/5,1)	0.099371	0.099576	0.099791	0.099962	0.100035
ũ(2π/5,1)	−0.124557	−0.124638	−0.124584	−0.124602	−0.124594
ũ(3π/5,1)	0.124557	0.124638	0.124584	0.124602	0.124594
ũ(4π/5,1)	−0.099371	−0.099576	−0.099791	−0.099962	−0.100035

Table 3

Example 4: Comparison of absolute errors of our scheme with scheme in [42, 43]

Method of [42]				Method of [43]			The presented method
x	α = 1.2	α = 1.5	α = 1.8	α = 1.2	α = 1.5	α = 1.8	α = 1.2	α = 1.5	α = 1.8
0.2	1.54 e-12	3.78 e-11	4.34 e-11	1.61 e-3	1.42 e-3	7.24 e-4	0.1055e-12	0.9843e-12	1.9413e-12
0.4	1.44 e-10	8.64 e-11	1.37 e-10	1.42 e-3	1.02 e-3	5.74 e-4	1.47905e-12	3.9820e-12	2.3070e-12
0.6	1.63 e-10	1.54 e-10	1.48 e-10	8.75 e-4	7.49 e-4	4.12 e-4	3.3255e-12	5.0410e-12	9.1279e-13
0.8	1.30 e-10	2.10 e-10	1.55 e-10	7.50 e-4	6.16 e-4	3.22 e-4	1.1055e-12	3.0410e-12	2.3986e-12
1.0	1.89 e-10	2.21 e-10	2.19 e-10	7.13 e-4	5.76 e-4	2.95 e-4	2.1152e-12	1.0010e-12	8.3765e-12
1.2	2.48 e-10	1.94 e-10	2.53 e-10	7.50 e-4	6.16 e-4	3.22 e-4	7.0055e-12	1.0410e-13	5.1819e-12
1.4	1.73 e-10	1.51 e-10	1.74 e-10	8.75 e-4	7.49 e-4	4.12 e-4	1.3498e-12	1.2430e-12	2.2781e-12
1.6	5.11 e-11	9.89 e-11	6.21 e-11	1.42 e-3	1.02 e-3	5.74 e-4	6.245e-12	3.0410e-12	5.2410e-12
1.8	3.43 e-11	3.72 e-11	3.08 e-11	1.61 e-3	1.42 e-3	7.24 e-4	1.1055e-13	3.0410e-12	1.3452e-12