Computer Modeling in Engineering & Sciences

DOI: 10.32604/cmes.2021.015310

ARTICLE

A Pseudo-Spectral Scheme for Systems of Two-Point Boundary Value Problems with Left and Right Sided Fractional Derivatives and Related Integral Equations

I. G. Ameen1, N. A. Elkot2, M. A. Zaky3,*, A. S. Hendy4,5 and E. H. Doha2

1Department of Mathematics, Faculty of Science, Al-Azhar University, Cairo, Egypt
2Department of Mathematics, Faculty of Science, Cairo University, Giza, 12613, Egypt
3Department of Applied Mathematics, Physics Division, National Research Centre, Dokki, Cairo, 12622, Egypt
4Department of Computational Mathematics and Computer Science, Institute of Natural Sciences and Mathematics, Ural Federal University, Yekaterinburg, 620002, Russia
5Department of Mathematics, Faculty of Science, Benha University, Benha, 13511, Egypt
*Corresponding Author: M. A. Zaky. Email: ma.zaky@yahoo.com; zaky.nrc@gmail.com
Received: 08 December 2020; Accepted: 04 March 2021

1  Introduction

Fractional-order differential operators have recently risen to prominence in the modelling of several processes. The mathematical models involving these operators have also attracted much attention, a survey of recent activity is given in [1–3]. The issue we address in this paper is to construct and analyse a spectral collocation method to solve the following nonlinear system of Caputo fractional two-point boundary value problems:

{−1CDzμu(z)=fL(z,u(z),v(z)),zCD1μv(z)=fR(z,u(z),v(z)),u(−1)=u0,u′(−1)=u1,v(1)=v0,v′(1)=v1,μ∈(1,2), (1)

where fR and fL:[−1,1]×R→R are continuous functions and satisfy the Lipschitz condition (64), and −1CDzμ and zCD1μ are the left- and right-sided Caputo fractional derivatives, respectively (see definition 2). In case of μ=2, then −1CDzμ and zCD1μ coincide with the usual second order derivative u″(z) and v″(z), and the system (1) recovers the integer-order system of two point boundary value problems.

Because the fractional-order differential operators are nonlocal with weakly singular kernels, the numerical discretization of the fractional models is more change than the classical schemes. There are several analytical schemes to solve fractional differential equations, such as the Green’s function method, the Mellin transform method, the Laplace transform method, the Fourier transform method, and so on [4–7]. However, analytical methods are rare for most of fractional differential equations, e.g., with non linearities or linear equations with time-dependent coefficients. Hence, constructing efficient numerical approaches is of great importance in practical applications.

Many numerical schemes have been developed to solve the fractional differential equations, mostly with the finite element methods (e.g., [8–11] and the references therein) and the finite difference methods (e.g., [12–18] and the references therein). Since spectral methods are capable of providing high-order accurate numerical approximations with less degrees of freedoms [19–23], they have been widely used for numerical approximations of fractional differential equations [24–29] or its related integral equations [30–37]. In particular, well designed spectral methods appear to be particularly attractive to tackle the difficulties associated with the weakly singular kernels of the fractional differential operators and the integral equations [38,39].

The system of fractional two-point boundary value problems (1) can be converted to an equivalent weakly singular nonlinear system of Volterra integral equations (15). The key idea of the presented approach is to solve (15) using the Legendre spectral collocation scheme. The aim of that convert to (15) is to approximate the related integral terms by the Gauss quadrature formula. The presented method has spectral convergence. This theoretical estimate is confirmed by two numerical test examples. Specifically, our strategies and contributions are highlighted as follows:

     i)  The system of fractional two-point boundary value problems is recast into an equivalent weakly singular nonlinear system of Volterra integral equations.

    ii)  The Legendre spectral collocation method is applied to the transformed equation.

   iii)  The convergence analysis of the Legendre collocation method under the L2-norms is derived.

The structure of the paper is as follows. In Section 2, we introduce some necessary definitions, notations and lemmas. The Legendre spectral collocation scheme is presented in Section 3. The convergence analysis is provided in Section 4. In Section 5, numerical examples are performed to confirm the efficiency of the numerical method. A brief conclusion is highlighted in Section 6.

2  Mathematical Preliminaries

In this section, we provide some notations, definitions, and some useful lemmas about the fractional differential and integral operators [4] and the Jacobi polynomials.

Definition 1. Let t∈[−1,1], for α>0, the left and right Riemann-Liouville fractional integrals of order μ are defined, respectively, as:

−1Izμu(z)=1Γ(μ)∫−1z(z−τ)μ−1u(τ)dτ,zI1μυ(z)=1Γ(μ)∫z1(τ−z)μ−1υ(τ)dτ, (2)

where Γ(.) is the usual Gamma function.

Definition 2. The left- and right-sided Caputo fractional derivatives are defined as:

(−1CDzμu)(z)=1Γ(m−μ)∫−1z(z−τ)m−μ−1u(m)(τ)dτ,(zCD1μu)(z)=(−1)mΓ(m−μ)∫z1(τ−z)m−μ−1u(m)(τ)dτ, (3)

where , m∈N.

Definition 3. The left- and right-sided Riemann-Liouville fractional derivatives are defined as:

(−1RLDzμu)(z)=1Γ(m−μ)dmdzm∫−1z(z−τ)m−μ−1u(τ)dτ,(zRLD1μu)(z)=(−1)mΓ(m−μ)dmdzm∫z1(τ−z)m−μ−1u(τ)dτ. (4)

It is worthy to mention here that the left- and right-sided Caputo fractional derivatives satisfy the following fundamental properties

Theorem 1. There hold [4]

aIzα(aCDzμu(z))=aCDzμ−αu(z)−∑j=⌈μ−α⌉⌈μ⌉−1u(j)(a)Γ(j+α−μ+1)(z−a)j+α−μ, (5)

zITα(zCDTμu(z))=zCDTμ−αu(z)−∑j=⌈μ−α⌉⌈μ⌉−1u(j)(T)Γ(j+α−μ+1)(T−z)j+α−μ, (6)

where μ≥α.

The following formulas introduce the relationship between the Riemann-Liouville and the Caputo fractional derivatives [5].

aDzμu(z)=aCDzμu(z)+∑j=0m−1u(j)(a)Γ(j−μ+1)(z−a)j−μ, (7)

zDTμu(z)=zCDTμu(z)+∑j=0m−1u(j)(T)Γ(j−μ+1)(T−z)j−μ. (8)

Let θ>0 and ϑ>−1, then

and for θ∈(m−1,m) with m∈N and ϑ≠0,

Now, we give some basic properties of the Jacobi polynomials and related Jacobi–Gauss interpolation. For ν, υ>−1, the Jacobi polynomials Jiν,υ(ζ), ζ∈Λ=[−1,1] of degree i form a complete Lων,υ2(Λ) orthogonal system with the weight function ων,υ=(1−ζ)ν(1+ζ)υ, i.e.,

∫ΛJiν, υ(ζ)Jjν, υ(ζ)ων, υ(ζ)dζ=γiν, υδi, j (9)

where, δi,j is the Kronecker function, and

γiν, υ=2ν+υ+1Γ(i+ν+1)Γ(i+υ+1)(2i+ν+υ+1)i!Γ(i+ν+υ+1). (10)

Denote PN(Λ) the space of all polynomials of degree less than or equal to N and {ω¯iν,υ,xiν,υ}i=0N are the set of weights and nodes of the Gauss-Jacobi interpolation. The associated Gauss-Jacobi integration formula can be written as:

∫ΛQ(x)ων,υ(x)dx≈∑i=0NQ(xi)ω¯iν,υ. (11)

The formula (11) is exact for any Q(x)∈P2N+1(Λ). Accordingly,

∑k=0NJiν,υ(xkν,υ)Jjν,υ(xkν,υ)ω¯kν,υ=γiν,υδi,j,∀0≤i+j≤2N+1. (12)

For any u∈C[−1,1], the Gauss-Jacobi interpolating operator Jx,Nν,υ:C[−1,1]→PN(Λ) is determined uniquely by

Jx,Nν,υu(xjν,υ)=u(xjν,υ),0≤j≤N. (13)

The above condition indicates that Jx,Nν,υu=u, ∀u∈PN. Consequently, since Jx,Nν,υu∈PN, then we can write

Jx,Nν,υu(x)=∑i=0Nυ^iν,υJiν,υ(x),where υ^iν,υ=1γiν,υ∑j=0Nu(xj)Jjν,υ(xj)ω¯jν,υ. (14)

The Legendere polynomials can be obtained directly from the properties of the Jacobi polynomials by setting ν=υ=0 as Li(x)=Ji0,0(x). In the following sections, we drop the parameters ν and υ whenever ν=υ=0.

3  The Pseudo-Spectral Method

We consider the system of two-point fractional boundary value problems (1) with homogeneous boundary conditions. There is no loss of generality since this can always be accomplished by a simple change of variables.

{−1CDzμu(z)=fL(z,u(z),v(z)),zCD1μv(z)=fR(z,u(z),v(z)),u(−1)=u′(−1)=v(1)=v′(1)=0,μ∈(1,2). (15)

The above system is equivalent to the following system of weakly singular integral equations:

{u(x)=1Γ(μ)∫−1x(x−σ)μ−1fL(σ,u(σ),υ(σ))dσ,υ(x)=1Γ(μ)∫x1(σ−x)μ−1fR(σ,u(σ),υ(σ))dσ. (16)

The variable transformations σ1(x,ζ)=x+12ζ−1−x2 and σ2(x,ζ)=x+12−x−12ζ are used in the first and second equations of the system (16) to convert the intervals (−1,x) and (x,1) to the unit interval Λ as follows

u(x)=1Γ(μ)(1+x2)μ∫−11(1−ζ)μ−1fL(σ1(x,ζ),u(σ1(x,ζ)),v(σ1(x,ζ))dζ, (17)

v(x)=1Γ(μ)(1−x2)μ∫−11(1+ζ)μ−1fR(σ2(x,ζ),u(σ2(x,ζ)),v(σ2(x,ζ)))dζ. (18)

The Legendre spectral collocation scheme for (17) and (18) is to seek uN(x) and vN(x)∈PN(Λ) with N≥1, such that

uN(x)=1Γ(μ)Jx,N[(x+12)μ∫−11(1−ζ)μ−1Jζ,Nμ−1,0fL(σ1(x,ζ),uN(σ1(x,ζ)),vN(σ1(x,ζ))dζ] (19)

vN(x)=1Γ(μ)Jx,N[(1−x2)μ∫−11(1+ζ)μ−1Jζ,N0,μ−1fR(σ2(x,ζ),uN(σ2(x,ζ)),vN(σ2(x,ζ)))dζ]. (20)

We now describe the implementation procedure of (19), (20) in detail. We consider the following Legendre approximations

uN(x)=∑r=0Nu^rLr(x),Jx,NJζ,Nμ−1,0((1+x2)μfL(σ1(x,ζ),uN(σ1(x,ζ)),vN(σ1(x,ζ)))=∑q=0N∑q′=0Nρ^q,q′Lq(x)Jq′μ−1,0(ζ), (21)

and

vN(x)=∑s=0Nv^sLs(x),Jx,NJs,N0,μ−1((1−x2)μfR(σ2(x,ζ),uN(σ2(x,ζ)),vN(σ2(x,ζ))))=∑r=0N∑r′=0Nξ^r,r′Lr(x)Jr′0,μ−1(ζ). (22)

Then, by (21) and (9) direct computations lead to

1Γ(μ)∫−11(1−ζ)μ−1Jx,NJζ,Nμ−1,0((1+x2)μfL(σ1(x,ζ),uN(σ1(x,ζ)),vN(σ1(x,ζ))))dζ=1Γ(μ)∑q=0N∑q′=0Nρ^q,q′Lq(x)∫−11(1−ζ)μ−1Jq′μ−1,0(ζ)dζ=2μΓ(μ+1)∑q=0Nρ^q,0Lq(x). (23)

Applying (12) to (21), one can verify readily that

ρ^q,0=μ(2q+1)2μ+1∑i=0N∑j=0N(xi+12)μ×fL(σ1(xi,ζjμ−1,0),uN(σ1(xi,ζjμ−1,0)),vN(σ1(xi,ζjμ−1,0)))Lq(xi)ω¯iω¯jμ−1,0. (24)

Similarly

1Γ(μ)∫−11(1+s)μ−1Jx,NJζ,N0,μ−1((1−x2)μfR(σ2(x,ζ),uN(σ2(x,ζ)),vN(σ2(x,ζ)))dζ=1Γ(μ)∑r=0N∑r′=0Nξ^r,r′Lr(x)∫−11(1+ζ)μ−1Jr′0,μ−1(s)dζ=2μΓ(μ+1)∑r=0Nξ^r,0Lr(x). (25)

Using (12)–(22) yields

ξ^r,0=μ(2r+1)2μ+1∑i=0N∑j=0N(1−xi2)μ×fR(σ2(xi,ζj0,μ−1),uN(σ2(xi,ζj0,μ−1)),vN(σ2(xi,ζj0,μ−1))Lr(xi)ω¯iω¯j0,μ−1. (26)

Hence, by using (17)–(26) we deduce that

∑i=0Nu^iLi(x)=2μΓ(μ+1)∑i=0Nρ^i,0Li(x), (27)

∑i=0Nv^iLi(x)=2μΓ(μ+1)∑i=0Nξ^i,0Li(x). (28)

Finally, using (9) we obtain

u^i=2μΓ(μ+1)ρ^i,0,0≤i≤N, (29)

v^i=2μΓ(μ+1)ξ^i,0,0≤i≤N. (30)

This system of equations can be solved for u^i and v^i. Then by using (21) and (22), we obtain an approximate solution for the problem (1).

4  Convergence Analysis

4.1 Auxiliary Lemmas

In this section, we introduce some functional spaces. We denote ∂xmg(x) to be the mth derivative of g i.e., ∂xmg(x):=dmgdxm(x). We also denote the Lων,μ2(Λ) inner product and norm by

(g,h)ων,μ:=∫Λg(x)h(x)ων,μdx, (31)

‖g‖ων,μ:=(∫Λ|g(x)|2ων,μdx)12. (32)

Definition 4. Let s≥1 be an integer. The Sobolev space Hων,μs is defined as

Hων,μs(Λ):={g∈Lων,μ2(Λ):∂xmg∈Lων,μ2(Λ),0≤m≤s}, (33)

with the inner product and norm

(g,v)Hων,μs=∑m=0s(∂xmg,∂xmv)ων,μ, (34)

‖g‖Hων,μs=(g,v)Hων,μs12. (35)

Definition 5. For a non-negative integer s. The weighted Jacobi non-uniformly Sobolev space Bων,μs is defined as

Bων,μs(Λ):={g:∂xmg∈Lων+m,μ+m2(Λ),0≤m≤s}. (36)

with the inner product, norm, and semi-norm

(g,v)Bων,μs=∑m=0s(∂xmg,∂xmv)ων+m,μ+m, (37)

‖g‖Bων,μs=(g,g)Bων+m,μ+ms12. (38)

In particular, L2(Λ)=Bω0,00(Λ) and ‖.‖=‖.‖ω0,0. It is obvious that Hων,μs(Λ) is a subspace of Bων,μs(Λ), that is

The space L∞(Λ) is defined with the norm

∥u∥∞=esssupx∈Λ|u(x)|. (39)

Lemma 1. (cf. [40]) Let ν,μ>−1, then for any u∈Bων,μs(Λ) with s≥1 and 0≤k≤s≤N+1,

‖∂xk(u−Jx,Nν,μu)‖ων+k,μ+k≤cNk−s‖∂xsu‖ων+k,μ+k, (40)

where Jx,Nν,μ is the Gauss-Jacobi interpolation operator.

Lemma 2. (cf. [40]) For any u∈Hs(Λ) with 0≤s≤N+1,

‖u−Jx,Nu‖∞≤cN34−s‖∂xsu‖∞. (41)

Lemma 3. (cf. [41]) Let {Fi(x)}i=0N be the Nth interpolation Lagrange polynomials related to the N + 1 Gauss points of the Jacobi polynomial. Then,

‖JNν,μ‖∞:=maxx∈Λ∑i=0N|Fi(x)|={O(log⁡N),−1<ν,μ≤−12,O(Nγ+12),γ=max(ν,μ),otherwise. (42)

Let ζiμ−1,0 be the Gauss-Jacobi nodes in Λ and σ1,iμ−1,0=σ1(x,ζiμ−1,0). The mapped Gauss-Jacobi interpolating operator xJ~σ,Nμ−1,0:C(−1,x)→PN(−1,x) is given as

xJ~σ,Nμ−1,0u(σiμ−1,0)=u(σiμ−1,0),0≤i≤N. (43)

Hence

xJ~σ,Nμ−1,0u(σiμ−1,0)=u(σiμ−1,0)=u(σ1(x,ζiμ−1,0))=xJσ,Nμ−1,0u(σ1(x,ζiμ−1,0)), (44)

and

xJ~σ,Nμ−1,0u(σ)=Jζ,Nμ−1,0u(σ1(x,ζ))|ζ=2σ11+x+1−x1+x. (45)

It is not difficult to obtain the following results

∫−1x(x−σ)μ−1xJ~σ,Nμ−1,0fL(σ,u(σ),v(σ))dσ=(1+x2)μ∫−11(1−ζ)μ−1Jζ,Nμ−1,0fL(σ1(x,ζ),u(σ1(x,ζ)),v(σ1(x,ζ)))dζ=(1+x2)μ∑j=0NfL(σ1(x,ζjμ−1,0),u(σ1(x,ζjμ−1,0)),v(σ1(x,ζjμ−1,0)))ω¯jμ−1,0=(1+x2)μ∑j=0NfL(σ1,jμ−1,0,u(σ1,jμ−1,0),v(σ1,jμ−1,0))ω¯jμ−1,0. (46)

Similarly

∫−1x(x−σ)μ−1(xJ~σ,Nμ−1,0fL(σ,u(σ),v(σ)))2dσ=(1+x2)μ∑j=0NfL2(σ1,jμ−1,0,u(σ1,jμ−1,0),v(σ1,jμ−1,0))ω¯jμ−1,0. (47)

We denote J to be the identity operator. Then, for any 1≤s≤N+1 we have

∫−1x(x−σ)μ−1|J−xJ~σ,Nμ−1,0u(σ)|2dσ=(1+x2)μ∫−11(1−ζ)μ−1|(J−Jζ,Nμ−1,0)u(σ1(x,ζ))|2dζ≤cN−2s(1+x2)μ∫−11(1−ζ)μ+s−1(1+ζ)s|∂ζsu(σ1(x,ζ))|2dζ=cN−2s∫−1x(x−σ)μ+s−1(1+σ)s|∂ζsu(σ)|2dσ (48)

Let ζi0,μ−1 be the Jacobi-Gauss nodes in Λ and ζ2,i0,μ−1=σ2(x,ζi0,μ−1). The mapped Jacobi-Gauss interpolation operator xJ~ζ,N0,μ−1:C(x,1)→PN(x,1) is defined by

xJ~σ,N0,μ−1v(ζi0,μ−1)=v(ζi0,μ−1),0≤i≤N. (49)

Hence,

xJ~σ,N0,μ−1v(ζi0,μ−1)=v(ζi0,μ−1)=v(σ2(x,ζi0,μ−1))=xJζ,N0,μ−1v(σ2(x,ζi0,μ−1)) (50)

and

xJ~ζ,N0,μ−1v(ζ)=Jζ,N0,μ−1v(σ2(x,ζ))|ζ=2σ21−x−1+x1−x (51)

We can also derive the following results

∫x1(σ−x)μ−1xJ~σ,N0,μ−1fR(σ,u(σ),v(σ))dσ=(1−x2)μ∫−11(1+ζ)μ−1Jζ,N0,μ−1fR(σ2(x,ζ),u(σ2(x,ζ)),v(σ2(x,ζ)))dζ=(1−x2)μ∑j=0NfR(σ2(x,ζj0,μ−1),u(σ2(x,ζj0,μ−1),v(σ2(x,ζj0,μ−1)))ω¯j0,μ−1=(1−x2)μ∑j=0NfR(σ2,j0,μ−1,u(σ2,j0,μ−1),v(σ2,j0,μ−1))ω¯j0,μ−1. (52)

Similarly

∫x1(σ−x)μ−1(xJ~σ,N0,μ−1fR(σ,u(σ),v(σ)))2dζ=(1−x2)μ∑j=0NfR2(σ2,j0,μ−1,u(σ2,j0,μ−1),v(σ2,j0,μ−1))ω¯j0,μ−1, (53)

and

∫x1(σ−x)μ−1|J−xJ~σ,N0,μ−1v(σ)|2dσ=(1−x2)μ∫−11(1+ζ)μ−1|(J−Jζ,N0,μ−1)v(σ2(x,ζ))|2dζ≤cN−2s(1−x2)μ∫−11(1+ζ)μ+s−1(1−ζ)s|∂ζsv(σ2(x,ζ))|2dζ=cN−2s∫x1(σ−x)μ+s−1(1−σ)s|∂ζsv(σ)|2dσ. (54)

4.2 Error Analysis in L2-Norm

In this section, we analyze the numerical errors of the systems (19) and (20). Let E = |Eu|+|Ev|, where Eu = u − uN, Ev = v − vN and denote by J to be the identity operator.

Lemma 4. The following inequality holds

‖E‖≤‖Eu‖+‖Ev‖≤∑i=16‖Ei‖ (55)

where