Numerical study of multidimensional fractional time and space coupled Burgers’ equations

Numerical study of multidimensional fractional time and space coupled Burgers’ equations

HODA F AHMED, M S M BAHGAT^∗and MOFIDA ZAKI

Mathematics Department, Faculty of Science, Minia University, El-Minia, Egypt

∗Corresponding author. E-mail: msmbahgat66@hotmail.com

MS received 1 August 2019; revised 17 October 2019; accepted 4 November 2019

Abstract. This paper declares a new spectral collocation technique to provide accurate approximate solutions of the one- and two-dimensional time and space fractional coupled Burgers’ equations (TSFCBEs) in which the fractional derivatives are defined according to Caputo’s definition. The suggested method is based on the shifted Gegenbauer polynomials (SGPs) for approximating the solution of TSFCBEs. The suggested technique reduces the considered problems to the solution of nonlinear algebraic equations (NLAEs). Moreover, the accuracy and reliability of the proposed method are confirmed through numerical examples. Finally, the obtained numerical results are compared with those previously reported in the literature.

Keywords. Collocation method; Caputo derivative; fractional-order Burgers’ equations; Gegenbauer polynomials.

PACS Nos 02.60.−x; 02.30.Jr; 02.30.Mv

1. Introduction

Spectral methods are useful approaches of discretisa- tion for many kinds of differential equations (DEs). In the last four decades, spectral methods [1,2] have been widely used in many fields. The speed of convergence and the high level of accuracy are their main advantages.

There are three popular techniques for spectral methods;

they are the collocation, Tau and Galerkin methods [3–

5]. The spectral collocation method (SCM) is very easy to apply and adjust to various problems because it pro- vides highly accurate solutions. The estimated solution as a finite sum of basis functions, which may be orthog- onal polynomials or a mixture of them, is the major step in spectral methods and then the coefficients are deter- mined in a way that decays the difference between the exact and approximate solutions as potential. In collo- cation technique, the coefficients of this expansion are determined by making the approximate solution to sat- isfy the DE at some applicably selected points from the domain. These points are identified as collocation points. The selection of these points is vital for the effi- ciency of the SCM [3,6–8], and consequently it has become more and more widespread for solving frac- tional differential equations (FDEs) [7,8]. Orthogonal polynomials such as Chebyshev, Legendre, Gegenbauer, Jacobi are employed for solving various kinds of DEs and FDEs [3–5,7,9,10]. The Gegenbauer polynomials

(GPs) [10–12] are an appropriate basis for estimations because they get speedy rates of convergence for a small range of the spectral expansion terms. So they have received much attention for their vital properties.

In 1995, Esipov [13] presented the coupled viscous Burgers’ equations as a simple model of sedimenta- tion or evolution of scale volume concentrations of two kinds of particles in fluid suspensions or colloids, under the effect of gravity. Many scientists have studied one- dimensional (1D) coupled viscous Burgers’ equations analytically and numerically [3,6,14–16]. The analytic solution of two-dimensional (2D) coupled viscous Burg- ers’ equations were proposed by Fletcher using the Hope–Cole transformation [17]. The two-dimensional coupled viscous Burgers’ equations have been con- sidered in many articles analytically and numerically [16,18–21]. In recent years, the study of the 1D and 2D TSFCBEs has attracted more attention (see [22–25]), because the reality of nature is translated superbly into a better and systematic manner by the theory of FDEs.

Most phenomena in physics, engineering, and other sciences can be modelled more accurately by using frac- tional calculus. In this direction of thought, FDEs have been developed as an interdisciplinary area of research in recent years. The non-local nature of fractional derivatives can be employed to simulate accurately varied natural phenomena having long memory [26–28].

0123456789().: V,-vol

In this work we mainly aim at presenting Caputo frac- tional derivatives and extend the application of SCM based on GPs for the following two samples:

1. 1D FCBE

D_t^α1u(x,t)=D^β_x¹u(x,t)+2u(x,t)Dxu(x,t)

−Dx(u(x,t)v(x,t)),

D^α_t²v(x,t)=D^β_x²v(x,t)+2v(x,t)D_xv(x,t)

−Dx(u(x,t)v(x,t)). (1)

2. 2D FCBE

D_t^α1u(x,y,t)+u(x,y,t)∂u(x,y,t)

∂x +v(x,y,t)∂u(x,y,t)

∂y

= 1 Re

∂²u(x,y,t)

∂x² + ∂²u(x,y,t)

∂y²

, D_t^α2v(x,y,t)+u(x,y,t)∂v(x,y,t)

∂x +v(x,y,t)∂v(x,y,t)

∂y

= 1 Re

∂²v(x,y,t)

∂x² +∂²v(x,y,t)

∂y²

, (2)

wheret >0. The factors αi, βi, 0≤α1 , α2 ≤1 and 1≤β1, β2 ≤2 denote the order of the fractional space and time derivatives, respectively. It is obvious that when at least one of the factors changes, different reaction sys- tems are obtained. Atα1 =α2 =1, β1 = β2 = 2 the fractional-order systems (1) and (2) define the conven- tional coupled Burgers’ equation in 1D and 2D spaces, respectively. In the proposed method, the shifted Gegen- bauer polynomials (SGPs) are applied as basis functions to approximate the solutions of the considered problems with the shifted Gegenbauer Gauss–Lobatto (SGGL) scheme for the spatial approximation and the shifted DE Gegenbauer Gauss (SGG) nodes scheme for tempo- ral discretisation. As a result of this approximation, sets of nonlinear algebraic equations (NLAEs) are obtained which are easy to solve by using any iterative technique.

To check out the accuracy of the proposed technique, numerical comparisons of several examples between our results and those attained by using existing meth- ods in the literature are presented. The current paper is planned as follows: Some mathematical notations which

are vital for our numerical results are introduced in §2.

The solution procedures of shifted Gegenbauer colloca- tion method (SGCM) for solving the 1D TSFCBEs and the 2D TFCBEs are discussed in §3. Numerical appli- cations are provided in §4. Conclusion is given in §5.

2. Preliminaries

DEFINITION 1

The Caputo time-fractional derivative operator of order β(left side) is defined as [26–28]

D_t^βu(x,t)=

⎧⎪

⎪⎨

⎪⎪

⎩ 1 (m−β)

t 0

(t−η)^m−β−1∂^mu(x, η)

∂η^m dη, m−1< β≤m, m ∈ N0.

∂^mu(x,t)

∂t^m , β =m, m ∈ N0.

(3)

The Caputo fractional derivative has the following prop- erties [26–28]:

(i) D^βC =0,Cis a constant.

(ii) D^βx^μ=

⎧⎪

⎪⎪

⎪⎪

⎨

⎪⎪

⎪⎪

⎪⎩

(μ+1)

(β+μ−1)x^μ−β, forμ∈ N0and μ≥βorμ /∈ N andμ >β .

0, forμ∈ N₀and

μ <β.

(iii) D_t^β(c₁h(x,t)+c₂ g(x,t)) = c₁D^β_t h(x,t) + c2D_t^β g(x,t).

βand β denote the ceiling and floor functions of β respectively. N and N0 stand for {1,2,3, . . .} and {0,1,2,3, . . .}respectively.

DEFINITION 2

The GPs,G^λ_n(x),n∈Z⁺[29,30] defined with the real parameterλ >−1/2, are a family of orthogonal polyno- mials and perform as the eigensolutions to the following singular Sturm–Liouville problem in[−1,1]:

d

dx 1−x²λ+¹₂ d

dx(G^λ_n(x)) +n(n+2λ)

1−x²_λ−(1/2)

G^λ_n(x)=0, with the first two terms

C₀^λ(x)=1, C₁^λ(x)=2λx

while the other polynomials are identified by (n+1)C_n+1^λ (x)

=2(λ+n)xC_n^λ(x)−(2λ+n−1)C_n^λ₋₁(x).

The GPs satisfy the orthogonality relation [30]

1

−1ω^λ(x)G^λ_n(x)G^λ_m(x)dz =γn^λδnm,

where the even functionω^λ(x)=(1−x²)^λ−(1/2)is the weight function for the GPs andγn^λis defined as γ_n^λ= π2^1−2λ(n+2λ)

(λ+n)((λ))²(n+1), andδnm is the Kronecker delta function.

It should be noted here that the GPs,G^λ_n(x)are nor- malised by Doha [31] such that

G^λ_n(1)=1, n =0,1,2, . . . in which the GPs are defined as G^λ_n(x)= (n+1)

λ+¹₂

λ+n+¹₂ P

λ−¹₂,λ−¹₂

n (x),

where P

λ−¹₂,λ−¹₂

n (x) is the Jacobi polynomial.

Throughout the rest of the paper, by the GPs we refer to those normalised by Doha [31]. This nor- malisation is characterised by an advantage that the polynomialsG⁰_n(x)are congruous with the Chebyshev polynomials of the first kind Tn(x), G¹_n^/2(x) are the Legendre polynomialsLn(x), andG¹_n(x)are equivalent to(1/(n+1))U_n(x), whereU_n(x)are the Chebyshev polynomials of the second kind. According to the new standardisation [31], the GPs satisfy the following orthogonality relation:

1

−1

ω^λ(x)G^λ_n(x)G^λ_m(x)dz =^λnδnm, where

^λn = 2^2λ−1(n+1)((¹₂ +λ))² (λ+n)(n+2λ) .

The approximation of function f(x)defined in the inter- val[−1,1]by GPs is given as [30]

f(x)= ∞ n=0

f_n^λG^λ_n(x), (4)

where the Gegenbauer coefficients are given by f_n^λ= 1

γ_n^λ ₁

−1ω^λ(x)f(x)G^λ_n(x)dx.

To use GPs for t ∈ [0,T], it is necessary to redefine the variable, x = _T²t−1. So the shifted Gegenbauer polynomials (SGPs), G^λ_n(_T²t−1),which is denoted by G^λ_T,n(t) are introduced. The SGPs satisfy the following orthogonality relation:

_T

0

ω^λ_T(x)G^λ_T,n(t)G^λ_T,m(t)dz =^λ_T,nδnm, where

ω^λT(x)=(T t −t²)^λ−12 and

^λT,n =(T/2)^2λ 2^2λ−1(n+1)((¹₂+λ))² (λ+n)(n+2λ) . An analytical form for the SGPs is known as G^λ_T,n(t)

= n m=0

(−1)^n−m n!(λ+¹₂)(n+m+2λ) T^mm!(n−m)!(m+λ+¹₂)(n+2λ)(t)^m,

0≤t ≤T. (5)

Note that the values of SGPs at the endpoints are deter- mined by

G^λ_T,n(0)=(−1)ⁿ(n+2λ) n!(2λ) and

G^λ_T,n(T)= (n+2λ) n!(2λ) .

The Caputo fractional derivative ofG^λ_T,n(t)is computed according to the following theorem.

Theorem 1. LetG^λ_T,n(t), 0<t<T denote the SGPs of degreenassociated with the parameterλ. Forβ >0, the Caputo fractional derivative of orderβ forG^λ_T,n(t) is

D^β(G^λ_T,n(t))= n

m=β

W_n^(λ,β)_,m ,

where W_n^(λ,β)_,m

= (−1)^n−mn!(λ+¹₂)(n+m+2λ) T^m(n−m)!(m+λ+¹₂)(n+2λ)(m−β+1). Proof. By operating the fractional derivative of orderβ to eq. (5), we have

D^β(G^λ_T,n(t))= n m=0

(−1)^n−mn!(λ+¹₂)(n+m+2λ) T^mm!(n−m)!

m+λ+¹₂

(n+2λ)D^β(t^m)

= n

m=β

(−1)^n−mn!(λ+¹₂)(n+m+2λ)(m+1)

T^m(n−m)!m!(m+λ+¹₂)(n+2λ)(m−β+1)t^m−β

= n

m=β

(−1)^n−mn!(λ+¹₂)(n+m+2λ)

T^m(n−m)!(m+λ+¹₂)(n+2λ)(m−β+1)t^m−β. We could easily complete the proof by noting that

D^β(t^m)=0 form =0,1,2, . . .β −1,β >0. The following theorem gives the estimate error for the function approximation.

Theorem 2. For a closed subspace Y in a Hilbert space H = L²[0,1]be such that Y = Span{G^λ_T,0(t), G^λ_T,1 (t), . . . , G^λ_T,N(t)}, let f(t)∈Cⁿ⁺¹[0,1], ifN

j=0 f_j^λ G^λ_T,j(t)is the best approximation of f(t)out ofY. Then

f(t)− N

j=0

f_j^λG^λ_T,j(t)

≤ h^(2n+3)/2 (n+1)!√

2n+3R t∈

tj,tj+1

⊆ [0,1], where

R= max

t∈[t_j,t_j+1]|f^(N+1)(τ)|

and

h =tj+1−t_j.

The proof of Theorem2is found in[32].

3. The proposed method

3.1 Procedure solution for 1D TSFCBE

In this subsection, we discuss the solution procedure for applying SGCM to numerically solve the 1D TSFCBEs of the following form:

D_t^αu(x,t)=D_x^βu(x,t)+2u(x,t)Dxu(x,t)

−Dx(u(x,t)v(x,t)), (x,t)∈[0,L]×[0,T],

D^α_t v(x,t)=D_x^βv(x,t)+2v(x,t)Dxv(x,t)

−Dx(u(x,t)v(x,t)),

(x,t)∈[0,L]× [0,T]. (6)

With initial and boundary conditions u(x,0)=h₁(x), v(x,0)=h₂(x), u(0,t)=h3(t),

u(L,t)=h4(t), v(0,t)=h5(t),

v(L,t)=h6(t), (7)

where h1(x),h2(x),h3(t),h4(t),h5(t) and h6(t) are given functions, the factorsα andβ;0 < α ≤ 1,1 <

β ≤ 2, point out the orders of the Caputo fractional space and time derivatives, respectively. In order to apply SGCM, the functions u(x,t) and v(x,t) are approximated as

u_N,M(x,t)= N i=0

M j=0

a_i,jG^λ_L,i(x)G^λ_T,j(t)

= N i=0

M j=0

ai,j f₀^i,j(x,t), (8)

v_N,M(x,t)= N i=0

M j=0

b_i,jG^λ_L,i(x)G^λ_T,j(t)

= N i=0

M j=0

bi,j f₀^i,j(x,t), (9) where f₀^i,j(x,t)= G^λ_L,i(x)G^λ_T,j(t).

Then the spatial partial derivatives(∂u(x,t)/∂x)and (∂v(x,t)/∂x)are calculated by

∂u(x,t)

∂x =

N i=0

M j=0

a_i,j ∂G^λ_L,i(x)

∂x G^λ_T,j(t)

= N i=0

M j=0

a_i,j f₁^i,j(x,t), (10)

∂v(x,t)

∂x =

N i=0

M j=0

b_i,j ∂G^λ_L,i(x)

∂x G^λ_T,j(t)

= N i=0

M j=0

b_i,j f₁^i,j(x,t), (11)

where

f₁^i,j(x,t)= ∂G^λ_L,i(x)

∂x G^λ_T,j(t).

Furthermore the Caputo fractional partial derivatives D_t^αu(x,t)andD^βxu(x,t)are assumed as

D_t^αu(x,t)= N i=0

M j=0

a_i,jG^λ_L,i(x)D_t^αG^λ_T,j(t)

= N i=0

M j=α

j k=α

a_i,jW_T^(λ,α)_,k G^λ_L,i(x)t^j−α

= N i=0

M j=α

a_i,j f₂^i,j(x,t), (12)

D^β_xu(x,t)= N

i=0

M j=0

a_i,jD_x^β G^λ_L,i(x)G^λ_L,j(t)

= N i=β

i k=β

M j=0

a_i,jW_L^(λ,β)_,k G^λ_T,j(t)x^i−β

= N i=β

M j=0

a_i,j f₃^i,j(x,t). (13) Similarly,

D^α_t v(x,t)= N i=0

M j=0

b_i,j f₂^i,j(x,t), (14)

D_t^βv(x,t)= N i=0

M j=0

b_i,j f₃^i,j(x,t), (15) where

f₂^i,j(x,t)= j k=α

W_T^(λ,α)_,k G^λ_L,i(x)t^j−α and

f₃^i,j(x,t)= i k=β

W_L^(λ,β)_,k G^λ_T,j(t)x^i−β.

Equations (8)–(15) enable one to write eq. (6) in the following form:

N i=0

M j=α

a_i,j f₂^i,j(x,t)= N i=β

M j=0

a_i,jf₃^i,j(x,t)

+2 N i=0

M j=0

a_i,j f₀^i,j(x,t) N i=0

M j=0

a_i,j f₁^i,j(x,t)

− N i=0

M j=0

a_i,j f₀^i,j(x,t) N i=0

M j=0

b_i,j fⁱ₁^,j(x,t)

− N i=0

M j=0

b_i,j f₀^i,j(x,t) N i=0

M j=0

ai,j f₁^i,j(x,t),

∀(x,t)∈[0,L]× [0,T] N

i=0

M j=α

b_i,j f₂^i,j(x,t)= N i=β

M j=0

b_i,j f₃^i,j(x,t)

+2 N i=0

M j=0

b_i,j f₀^i,j(x,t) N i=0

M j=0

b_i,j f₁^i,j(x,t)

− N i=0

M j=0

a_i,j f₀^i,j(x,t) N i=0

M j=0

b_i,j fⁱ₁^,j(x,t)

− N i=0

M j=0

b_i,j fⁱ₀^,j(x,t) N

i=0

M j=0

a_i,jf₁^i,j(x,t),

∀(x,t)∈[0,L]× [0,T]. (16) Also the initial and boundary conditions are approxi- mated as

u(x,0)= N i=0

M j=0

ai,j f₀^i,j(x,0)=h1(x),

v(x,0)= N i=0

M j=0

bi,j f₀^i,j(x,0)=h2(x) (17) and

u(0,t)= N i=0

M j=0

ai,j f₀^i,j(0,t)=h3(t),

u(L,t)= N i=0

M j=0

ai,j f₀^i,j(L,t)=h4(t),

v(0,t)= N i=0

M j=0

bi,j f₀^i,j(0,t)=h5(t),

v(L,t)= N i=0

M j=0

bi,j f₀^i,j(L,t)=h6(t), (18) respectively. In the SGCM, the residual of eq. (6) is set to zero at M(N −1) of SGGL and SGG points.

Accordingly, we have 2M(N −1) NLAEs of 2(M+ 1)(N +1)of the unknownsai,j andbi,j.

N i=0

M j=α

a_i,j f₂^i,j(x^λ_L,N,r,t_T^λ_,M,s)

− N i=β

M j=0

ai,j f₃^i,j(x^λ_L,N,r,t_T^λ_,M,s)

−2 N i=0

M j=0

ai,j f₀^i,j(x^λ_L,N,r,t_T^λ_,M,s)

× N i=0

M j=0

ai,j f₁^i,j(x^λ_L,N,r,t_T^λ_,M,s)

+ N i=0

M j=0

a_i,j f₀^i,j(x^λ_L,N,r,t_T^λ_,M,s)

× N i=0

M j=0

bi,j fⁱ₁^,j(x^λ_L,N,r,t_T^λ_,M,st)

+ N i=0

M j=0

bi,j f₀^i,j(x^λ_L,N,r,t_T^λ_,M,s)

× N i=0

M j=0

ai,j f₁^i,j(x^λ_L,N,r,t_T^λ_,M,s)=0, r =1, . . . ,N−1, s =1, . . . ,M N

i=0

M j=α

bi,j f₂^i,j(x^λ_L,N,r,t_T^λ_,M,s)

− N i=β

M j=0

bi,j f₃^i,j(x^λ_L,N,r,t_T^λ_,M,s)

−2 N i=0

M j=0

b_i,j f₀^i,j

x_L^λ_,N,r,t_T^λ_,M,s

× N i=0

M j=0

bi,j f₁^i,j

x^λ_L,N,r,t_T^λ_,M,s

+ N i=0

M j=0

ai,j f₀^i,j

x^λ_L,N,r,t_T^λ_,M,s

× N i=0

M j=0

bi,j f₁^i,j

x^λ_L,N,r,t_T,M,s^λ t

+ N i=0

M j=0

bi,j f₀^i,j

x^λ_L,N,r,t_T^λ_,M,s

× N i=0

M j=0

a_i,j f₁^i,j

x^λ_L,N,r,t_T^λ_,M,s

=0,

r =1, . . . ,N−1, s =1, . . . ,M. (19) For the initial and boundary conditions described by eqs (17) and (18), respectively, we have 2(N −1)and

4(M +1)NLAEs in 2(M +1)(N +1)of unknowns;

ai,j andbi,j

N i=0

M j=0

a_i,j f₀^i,j(x^λ_L,N,r,0)=h₁(x_L^λ_,N,r), r =1, . . . ,N−1

N i=0

M j=0

b_i,j f₀^i,j(x^λ_L,N,r,0)=h₂(x_L^λ_,N,r), r =1, . . . ,N−1

N i=0

M j=0

a_i,j f₀^i,j(0,t_T^λ_,M,s)=h₃(t_T^λ_,M,s), s =0, . . . ,M

N i=0

M j=0

ai,j f₀^i,j(L,t_T,M,s^λ )=h4(t_T,M,s^λ ), s =0, . . . ,M

N i=0

M j=0

bi,j f₀^i,j(0,t_T^λ_,M,s)=h5(t_T^λ_,M,s), s =0, . . . ,M

N i=0

M j=0

bi,j f₀^i,j(L,t_T^λ_,M,s)=h6(t_T^λ_,M,s),

s =0, . . . ,M. (20)

By combining eqs (19) and (20), we obtain a system of 2(M+1)(N +1)NLAEs which are easy to solve.

After the coefficientsai,j andbi,j are computed, it is straightforward to calculate the approximate solutions u_N,M(x,t) andvN,M(x,t) at any value of(x,t) in the defined domain of eqs (8) and (9).

3.2 2D TFCBEs

Consider the following 2D TFCBEs:

D_t^αu(x,y,t)+u(x,y,t)Dxu(x,y,t) +v(x,y,t)Dyu(x,y,t)

= 1

Re[Dx xu(x,y,t)+Dyyu(x,y,t)]

D_t^βv(x,y,t)+u(x,y,t)Dxv(x,y,t) +v(x,y,t)Dyv(x,y,t)

= 1

Re[Dx xv(x,y,t)+Dyyv(x,y,t)], (21)

∀(x,y,t)∈[0,L1]×[0,L2]× [0,T], where

D_x = ∂

∂x, D_y = ∂

∂y,