# Numerical Solutions of Coupled Systems of Fractional Order Partial Differential Equations

Authors: Yongjin Li; Kamal Shah
Journal: Advances in Mathematical Physics
 (2017)

---

## Abstract

We develop a numerical method by using operational matrices of fractional order integrations and differentiations to obtain approximate solutions to a class of coupled systems of fractional order partial differential equations (FPDEs). We use shifted Legendre polynomials in two variables. With the help of the aforesaid matrices, we convert the system under consideration to a system of easily solvable algebraic equation of Sylvester type. During this process, we need no discretization of the data. We also provide error analysis and some test problems to demonstrate the established technique.


---

## Body

## 1. Introduction

Recently fractional calculus has attracted the attention of many researchers as it has many applications in various disciplines of applied sciences and engineering (we refer to [1–4] and the references therein). In fact most of the engineering and physical processes can be modeled and well explained by using a system of fractional order differential and integral equations more realistically as compared to system of conventional differential and integral equations, (see [1, 5–15] and the references therein). The area of fractional calculus devoted to the existence and uniqueness of positive solution to (FDEs) and (FPDEs) is well studied and a lot of research work is available on it (we refer to [16–20]). In the last two decades, the area dealing with numerical solutions of FDEs and FPDEs has attracted the attention of many researchers. Researchers are taking keen interests in developing numerical procedures for FPDEs. Here, it is remarkable that most of the problems of applied sciences can be modeled by using PDEs and their systems. The study of coupled systems of PDEs can widely be found in engineering and biomechanics and other disciplines. For example, in biomechanics when modeling the phenomena of electrical activity in the heart, the resulting equations are coupled systems of PDEs (see, e.g., [21–23]). Coupled systems of PDEs also occur in modeling of some chemical and material engineering processes such as system containing a continuous stirred tank reactor (CSTR) and a plug flow reactor (PFR) in series (see [24, 25]). Various applications of coupled systems can be found in solid state physics and mechanics; for example, dynamics of multideformable bodies coupled by standard light fractional order discrete continuum layers is described by coupled systems of fractional order partial differential equations [7, 26–28]. The coupled systems of PDEs also appear in modeling of some important electromagnetic and gravitational problems, (see, e.g., [29, 30]).Most of the FPDEs do not have exact analytical solutions; therefore researchers need some appropriate numerical technique for the approximate solutions of such types (FPDEs). For the numerical solutions numerous techniques are available in the literature; for example, some of them are eigenvector expansion, Adomian decomposition method (ADM), fractional differential transform method (FDTM) [31, 32], generalized block pulse operational matrix method [33], and so on. In recent years many attentions have been devoted to develop operational matrices for fractional order differentiation and integrations. Such types of operational matrices are obtained by using Haar wavelets, Legendre wavelets, Sine Cosine, Chebyshev wavelets, and so on. Based on these operational matrices, some approximate methods for numerical solutions to FDEs have been developed (we refer to [34–36] and the references therein). In [37], the authors established a new numerical scheme based on the Haar wavelet for solving (FPDEs). More recently, the authors [38] developed some new results related to the Jacobi polynomials for operational matrices. All the operational matrices methods are used to solve FDEs and FPDEs. As much as we know they are not well studied and neither is properly applied to solve coupled systems of PFDEs. To the best of our knowledge very few articles are devoted to the numerical solutions of coupled systems of integral equations and ordinary and partial differential equations of fractional order; for details see [39–41].This manuscript is devoted to bring out the numerical solutions to a coupled system of FPDEs on an unbounded domain. The concerned coupled system contains mixed partial order derivative with respect tox, y as given by (1)DxpcUx,y+ω1Dxyp1cVx,y+ω2Dyq1cUx,y=θx,y,DxpcVx,y+ρ1Dcxyp2Ux,y+ρ2Dyq2cVx,y=ψx,y,corresponding to the conditions(2)Ui0,y=θiy,Vi0,y=ψiy,i=1⋯m,where Dzpc=∂p/∂zp represents partial derivative, p∈(k-1,k],k-2<pi,qi≤k-1, the real constants are denoted by ωi,ρi(i=1,2),θi,ψi, and (3)Ux,y,Vx,y,θx,y,ψx,y∈C0,1×0,1.It is remarkable that the proposed system contains the coupled system investigated in [38] as a special case. If we take θ(x,y)=ψ(x,y)=0 together with ω1=0,ω2=1 and ρ1=0,ρ2=1, the proposed system becomes a coupled system of Laplace equations of fractional order by taking p=q1=q2∈(1,2]. Further, if we involve the external (source) terms and consider ω1=ρ1=0, the proposed coupled system assumes the form of coupled system of Poisson’s equation of arbitrary order. Poisson’s equations have many applications in physics as well as in other various disciplines including electrostatics, mechanical engineering, and theoretical physics. For example, it has been applied to describe the potential energy field caused by a given charge or mass density distribution (see for details [42]).With the use of shifted Legendre polynomials in two variables, some operational matrices corresponding to fractional order differentiations and integrations are developed. Thanks to these operational matrices, the coupled system under consideration is transformed to a system of Sylvester type algebraic equations. Here, we remark that no discretization like for Tau-collocations method is required. By doing so, the proposed system (1) corresponding to the initial conditions is transformed to a system of Sylvester type algebraic equations. Our proposed method is simple and there is no computational complexity in the resulting algebraic system to solve. Further, we remark that our present method is a numerical method based on shifted Jacobi polynomials while the method discussed in [43] is an iterative method based on monotone iterative technique. Further we can easily establish a simple relationship for convergence of the propped method. The aforesaid procedure has now recently being applied for some nonlinear fractional order differential equations and some fruitful results were obtained; for details see [44]. In our future work we will use operational matrices method to compute approximate solutions of nonlinear FPDEs and their system.The manuscript is organized as follows: Section2 is concerning some basic definitions and preliminaries results related to fractional calculus and Legendre polynomials which are needed to form the proposed method. Section 3 is concerning the operational matrices corresponding to the fractional order derivatives and integration. Further, Section 4 is related to the application of the operational matrices to covert coupled system (1) to the corresponding system of algebraic equations. Section 5 is devoted to the numerical examples, and in the last section, we have provided a brief conclusion.


## 2. Preliminaries

Definition 1.
The Mittag-Leffler function subjected to one parameter is recalled by(4)Epx=∑n=0∞xnΓnp+1. The extension of (4) to two parameters is given as (5)Ep,qx=∑n=0∞xnΓnp+q.Definition 2 (see [45, 46]).
The fractional integral of Riemann-Liouville type of orderp∈(0,∞) of a function θ∈L1([0,a],R) is defined as (6)Ia+pθx=1Γp∫axx-ηp-1θηdη,provided that the right-hand side is pointwise defined on (0,∞).Definition 3.
For a given functionθ(x,y), the Caputo fractional derivative of order p is defined as (7)Dpcθx,y=1Γk-p∫axx-ηk-p-1θkη,ydη,p∈k-1,k,k∈N,provided that the right side is pointwise defined on (a,∞), where k=[p]+1.Hence, it follows that(8)Dpcxn=Γ1+nΓ1+n-pxn-p,Ipxn=Γ1+nΓ1+n+pxn+p,DpcC=0,foraconstantC.Theorem 4 (see [45]).
The solution of FDE(9)DxpcUx,y=hx,k-1<p≤k,is given by(10)Ux,y=Iphx+e0+e1x+⋯+ek-1xk-p,where each ei is a real number for i=0,1,2,…,k-1, where k=[p]+1, such that [p] represents integer part of p.Remark 5.
We will use Caputo derivative throughout this paper, because Caputo fractional derivative has clear geometrical representations like classical derivative. Further, it turns out that the Riemann-Liouville derivatives have certain disadvantages when trying to model real-world phenomena with fractional differential equations. For further details, see [45, 47].


### 2.1. The Shifted Legendre Polynomials

The well-known Legendre polynomials described over[-1,1] are representing as (11)Pi+1t=2i+1i+1tPit-ii+1Pi-1t,i=1,2,3,…,whereP0t=1,P1t=t.The transformation x=t+1/2 transforms the interval [-1,1] to J=[0,1] and the corresponding polynomials known as shifted Legendre polynomials are given as (12)Lix=∑n=0i-1i+nΓi+n+1Γi-n+1Γ2n+1xn,i=1,2,3,….The orthogonality condition is (13)∫01LixLjxdx=12i+1,ifi≠j,0,ifi=j.Therefore, any function θ(x)∈C(J) can be approximated in terms of shifted Legendre polynomials as (14)θx≈∑μ=0mCμLμx,whereCμ=θx,Lμx=2μ+1∫01θxLμxdx.In matrix form, we may write (14) as (15)θx=AKTΥK2x,y,whereK=m+1,where AK is the matrix of functions while ΥK2(x,y) is K term column matrix. The shifted Legendre polynomials of order K for two dimensions are recalled by (16)Lkx,y=LμxLνy,k=Kμ+ν+1,μ,ν=0,1,2,…,mand the corresponding orthognathic relation is provided as (17)∫01∫01Lμ1xLν1yLμ2xLν2ydxdy=12μ+12ν+1,ifμ1=μ2,ν1=ν2,0,ifμ1≠μ2,ν1≠ν2.Therefore any function U(x,y)∈C(J×J) can be approximated in terms of Lk(x,y) as (18)Ux,y≈∑μ=0m∑ν=0mCμνLμxLνy,where (19)Cμν=2μ+12ν+1∫01∫01Ux,yLμxLνydxdy.In matrix notation, we may write (18) as (20)Ux,y≈∑k=1K2CμνLnx,y=AK2TΥK2x,y,where the column vector of coefficients with order K2×1 is represented as AK2,k=Kμ+ν+1. Also column vector of order K2×1 is represented by ΥK2(x,y) and provided as (21)ΥK2x,y=φ11x,y⋯φ1Kx,yφ21x,y⋯φ2Kx,y⋯φKKx,yT,where φi+1,j+1(x,y)=(Li(x))(Lj(y)),i,j=0,1,2,…,m.


### 2.2. Maximum Absolute Error ([48])

Consider a sufficiently smooth functionU(x,y) on J×J; then the error in the numerical solutions is provided as (22)Ux,y-LK,Kx,y2≤C1+C2+C31KK+11KK+1,where(23)C1=14maxx,y∈J×J∂K+1∂xK+1Ux,y,C2=14maxx,y∈J×J∂K+1∂yK+1Ux,y,C3=116maxx,y∈J×J∂2K+2∂xK+1∂yK+1Ux,y.The above relation also holds at points (xi,yi) for interpolating polynomial L(K,K)(x,y).


## 2.1. The Shifted Legendre Polynomials

The well-known Legendre polynomials described over[-1,1] are representing as (11)Pi+1t=2i+1i+1tPit-ii+1Pi-1t,i=1,2,3,…,whereP0t=1,P1t=t.The transformation x=t+1/2 transforms the interval [-1,1] to J=[0,1] and the corresponding polynomials known as shifted Legendre polynomials are given as (12)Lix=∑n=0i-1i+nΓi+n+1Γi-n+1Γ2n+1xn,i=1,2,3,….The orthogonality condition is (13)∫01LixLjxdx=12i+1,ifi≠j,0,ifi=j.Therefore, any function θ(x)∈C(J) can be approximated in terms of shifted Legendre polynomials as (14)θx≈∑μ=0mCμLμx,whereCμ=θx,Lμx=2μ+1∫01θxLμxdx.In matrix form, we may write (14) as (15)θx=AKTΥK2x,y,whereK=m+1,where AK is the matrix of functions while ΥK2(x,y) is K term column matrix. The shifted Legendre polynomials of order K for two dimensions are recalled by (16)Lkx,y=LμxLνy,k=Kμ+ν+1,μ,ν=0,1,2,…,mand the corresponding orthognathic relation is provided as (17)∫01∫01Lμ1xLν1yLμ2xLν2ydxdy=12μ+12ν+1,ifμ1=μ2,ν1=ν2,0,ifμ1≠μ2,ν1≠ν2.Therefore any function U(x,y)∈C(J×J) can be approximated in terms of Lk(x,y) as (18)Ux,y≈∑μ=0m∑ν=0mCμνLμxLνy,where (19)Cμν=2μ+12ν+1∫01∫01Ux,yLμxLνydxdy.In matrix notation, we may write (18) as (20)Ux,y≈∑k=1K2CμνLnx,y=AK2TΥK2x,y,where the column vector of coefficients with order K2×1 is represented as AK2,k=Kμ+ν+1. Also column vector of order K2×1 is represented by ΥK2(x,y) and provided as (21)ΥK2x,y=φ11x,y⋯φ1Kx,yφ21x,y⋯φ2Kx,y⋯φKKx,yT,where φi+1,j+1(x,y)=(Li(x))(Lj(y)),i,j=0,1,2,…,m.


## 2.2. Maximum Absolute Error ([48])

Consider a sufficiently smooth functionU(x,y) on J×J; then the error in the numerical solutions is provided as (22)Ux,y-LK,Kx,y2≤C1+C2+C31KK+11KK+1,where(23)C1=14maxx,y∈J×J∂K+1∂xK+1Ux,y,C2=14maxx,y∈J×J∂K+1∂yK+1Ux,y,C3=116maxx,y∈J×J∂2K+2∂xK+1∂yK+1Ux,y.The above relation also holds at points (xi,yi) for interpolating polynomial L(K,K)(x,y).


## 3. Operational Matrices of Integrations and Differentiations

Extending the notion of the operational matrices of fractional order integration and differentiation from one dimension [36, 49] to two dimensions, the construction of the aforesaid operational matrices is provided as follows.Theorem 6.
The operational matrix corresponding to fractional integration of orderp of ΥK2(x,y) defined in (21) with respect to x is provided as (24)IxpΥK2x,y≃HK2×K2p,xΥK2x,y,where HK2×K2(p,x) is the operational matrix of fractional integration given as (25)HK2×K2p,x=¥1,1,n¥1,2,n⋯¥1,u,n⋯¥1,K2,n¥2,1,n¥2,2,n⋯Ω2,u,n⋯¥2,K2,n⋮⋮⋮⋮⋮⋮¥v,1,n¥v,2,n⋯¥v,u,n⋯¥v,K2,n⋮⋮⋮⋮⋮⋮¥K2,1,n¥K2,2,n⋯¥K2,u,n⋯¥K2,K2,n,u=Ki+j+1,v=Kμ+ν+1,¥v,u,n=∏i,j,μ,ν,n for i,j,ν,μ=0,1,2,…,m,(26)Ei,j,μ,ν,n=∑n=0μδi,ν2i+1∑l=0i-1i+l+μ+nΓi+l+1Γμ+n+1Γi-l+1Γ2l+1Γμ-n+1Γn+p+1n+l+p+1,δi,ν=0,ifi≠ν,1,ifi=ν.Proof.
ConsiderLk(x,y) as given in (16); then the arbitrary order integral with order p of Lk(x,y) corresponding to x is provided as (27)IxpLkx,y=IxpLμxLνy=∑n=0μ-1μ+nΓμ+n+1Γμ-n+1Γ2n+1IxpxnLνy,which gives that(28)IxpLμxLνy=∑n=0μ-1μ+nΓμ+n+1Γμ-n+1Γ2n+1Γn+p+1xn+pPνy.Further approximating xn+pLν(y) in terms of shifted Legendre polynomials, we have (29)xn+pLνy≈∑i=0m∑j=0mEi,jLixLjy,where Ei,j=(2i+1)(2j+1)∫01∫01xn+pLν(y)Li(x)Lj(y)dxdy. Thanks to orthogonality condition, we have (30)Ei,j=0,ifj≠ν,2i+1∑l=0j-1i+lΓi+l+1Γi-l+1Γ2l+1n+l+p+1,ifj=ν=δi,ν2i+1∑l=0i-1i+lΓi+l+1Γi-l+1Γ2l+1n+l+p+1,whereδi,ν=0,ifj≠ν,1,ifj=ν.Therefore (28) yields (31)IxpLμxLνy=∑n=0μ-1μ+nΓμ+n+1Γμ-n+1Γ2n+1Γn+p+1∑i=0m∑j=0mEi,j,νLixLjy=∑i=0m∑j=0m∑n=0μ-1μ+nΓμ+n+1Ei,j,νΓμ-n+1Γn+1Γn+p+1LixLjy=∑i=0m∑j=0mEi,j,μ,ν,nLixLjy,where(32)Ei,j,μ,ν,n=∑n=0μ-1μ+nΓμ+n+1Ei,j,νΓμ-n+1Γn+1Γn+p+1.We use these notations for easiness of u=Ki+j+1,v=Kμ+ν+1,¥v,u,n=Ei,j,ν,μ,nfori,j,μ,ν=0,1,2,3,…,m. Thus, we received the required resultTheorem 7.
The fractional order derivative ofΥK2(x,y) can be written as(33)DxpcΥK2x,y≃WK2×K2p,xΥK2x,y, such that WK2×K2(p,x) is the operational matrix corresponding to arbitrary order derivative of order p with respect to x provided as (34)WK2×K2p,x=¥1,1,n¥1,2,n⋯¥1,u,n⋯¥1,K2,n¥2,1,n¥2,2,n⋯¥2,u,n⋯¥2,K2,n⋮⋮⋮⋮⋮⋮¥v,1,n¥v,2,n⋯¥v,u,n⋯¥q,K2,n⋮⋮⋮⋮⋮⋮¥K2,1,n¥K2,2,n⋯¥K2,u,n⋯¥K2,K2,n,where v=Ki+j+1,u=Kμ+ν+1,¥v,u,n=Ωi,j,μ,ν,n  for  i,j,μ,ν=0,1,2,…,m, (35)℧i,j,μ,ν,n=∑l=pνδj,ν2i+1∑n=0i-1i+n+l+μΓi+n+1Γμ+l+1Γi-n+1Γ2n+1Γμ-n+1Γl+1Γl-p+1n+l-p+1with℧i,j,μ,ν,n=0,ifμ<p.Proof.
Takep the order fractional derivative of Ln(x,y) corresponding to x as(36)DxpcLkx,y=DxpcLμxLνy=∑n=0μ-1μ+nΓμ+n+1Γμ-p+1Γ2n+1DxpcxnLνy=∑l=pν-1μ+nΓμ+n+1Γμ-n+1Γ2n+1Γ1+n-pxn-pLνy.Approximate xn-pLν(y) in terms of shifted Legendre polynomials as(37)xn-pLνy≈∑i=0m∑j=0mDi,jLixLjy.Applying orthogonality condition, we may write (38)Di,j=0,ifj≠ν,2i+1∑n=0i-1i+nΓi+n+1Γi-n+1Γ2n+1n+l-p+1,ifj=ν.Therefore, we get(39)DxpcLkx,y=∑n=pμ-1μ+nΓμ+n+1Γμ-n+1Γn+12Γn-p+1∑i=0m∑j=0mDi,jLixLjy=∑i=0m∑j=0m∑n=pμ-1μ+nΓμ+n+1Γμ-n+1Γn+12Γn-p+1Di,jLixLjy=∑i=0m∑j=0m℧i,j,n,μ,νLixLjy,μ=p,…,K,where(40)℧i,j,n,μ,ν=∑n=pμ-1μ+nΓμ+n+1Γμ-n+1Γn+12Γn-p+1Di,j.Denoting u=Ki+j+1,v=Kμ+ν+1,¥v,u,n=℧i,j,ν,μ,n, such that i,j,μ,ν=0,1,2,3,…,m, for easiness, we get the required.Theorem 8.
The fractional derivative with orderp of ΥK2(x,y) with respect to y can be written as (41)DyαcΥK2x,y≃VK2×K2p,y¥K2x,y,such that VK2×K2(p,y) is the operational matrix given as (42)VK2×K2p,y=¥1,1,n¥1,2,n⋯¥1,u,n⋯¥1,K2,n¥2,1,n¥2,2,n⋯¥2,u,n⋯¥2,K2,n⋮⋮⋮⋮⋮⋮¥v,1,n¥v,2,n⋯¥v,u,n⋯¥v,K2,n⋮⋮⋮⋮⋮⋮¥K2,1,n¥K2,2,n⋯¥K2,u,n⋯¥K2,K2,n,u=Ki+j+1,v=Kμ+ν+1,¥v,u,n=℧i,j,μ,ν,n  for  i,j,μ,ν=0,1,2,…,m and(43)℧i,j,μ,ν,n=∑l=pνδi,μ∑l=0j-1j+l+ν+nΓj+l+1Γν+n+1Γj-l+1Γ2l+1Γν-n+1Γn+1Γn-p+1n+l-p+1,where℧i,j,μ,ν,n=0,ifν<p.Proof.
To prove this theorem it can be easily followed from the proof of Theorem7.Theorem 9.
The fractional derivative ofΥK2(x,y) with respect to x,y is provided as (44)DxypcΥK2x,y≃ZK2×K2p,x,yΥK2x,y,where ZK2×K2(p,x,y) represents the operational matrix of fractional derivative given as (45)ZK2×K2p,x,y=¥1,1,n¥1,2,n⋯¥1,u,n⋯¥1,K2,n¥2,1,n¥2,2,n⋯¥2,u,n⋯¥2,K2,n⋮⋮⋮⋮⋮⋮¥v,1,n¥v,2,n⋯¥v,u,n⋯¥v,K2,n⋮⋮⋮⋮⋮⋮¥K2,1,n¥K2,2,n⋯¥K2,u,n⋯¥K2,K2,n,and u=Ki+j+1,v=Kμ+b+1,¥u,v,n=℧i,j,μ,ν,n  for  i,j,μ,ν=0,1,2,…,m,(46)℧i,j,μ,ν,n=∑n=0μEμ,ν,n,l,pDi,j,μ,ν,Di,j,μ,ν=δi,j,μ,ν2i+12j+1∑n=0i∑l=0j-1i+j+n+lΓj+l+1Γi+n+1Γj-l+1Γi-n+1Γ2n+1Γ2l+1n+l-p/2+12,Eμ,ν,n,l,p=∑n=p/2μ∑l=p/2ν-1μ+nΓμ+n+1Γμ-n+1Γn+1-1ν+lΓν+l+1Γν-l+1Γl+1.Proof.
TakeLk(x,y) given in (17) and apply fractional order derivative of Lk(x,y) corresponding to y and x as (47)DxypcLkx,y=DxypcPμxLνy=∑n=0μ∑l=0ν-1μ+nΓμ+n+1Γμ-n+1Γn+1-1ν+lΓν+l+1Γν-1+1Γl+1Dxp/2xnDyp/2yl=∑n=p/2μ∑l=p/2ν-1μ+nΓμ+n+1Γμ-n+1Γn+1-1ν+lΓν+l+1Γν-l+1Γl+1xn-p/2yl-p/2Γn-p/2+1Γl-p/2+1.We approximate x(n-p/2)y(l-p/2) in terms of shifted Legendre polynomials as (48)xn-p/2yl-p/2≈∑i=0m∑j=0mDi,j,n,lLixLjy,such that (49)Di,j,n,l=2i+12j+1∫01∫01LμxLνyLixLjyxn-p/2yl-p/2dxdy.
Thanks to the orthogonality conditions and convolution theorem of Laplace transform, evaluating the integrals (49), we get (50)Di,j,n,l=δi,j,μ,ν2i+12j+1×∑n=0i∑l=0j-1i+j+n+lΓj+l+1Γi+n+1Γj-l+1Γi-n+1Γ2n+1Γ2l+1n+l-p/2+12,where (51)δi,j,μ,ν=1,ifi=μ,j=ν,0,ifi≠μ,j≠ν.Plugging (48) in (47), we may write(52)DxypcLkx,y=∑n=p/2μ∑l=p/2ν-1μ+nΓμ+n+1Γμ-n+1Γn+1-1ν+lΓν+l+1Γν-l+1Γl+1∑i=0m∑j=0mDi,j,n,lLixLjy=∑i=0m∑j=0m∑n=p/2a∑l=p/2ν-1μ+nΓμ+n+1Γμ-n+1Γn+1-1ν+lΓν+l+1Γν-l+1Γl+1Di,j,n,lLixLjy.Representing u=Kμ=ν+1,v=Ki+j+1, and ¥v,u,n=℧i,j,μ,ν,n,l for i,j,μ,ν=0,1,2,…,m, we obtain the required results.


## 4. Solutions of the Coupled Systems of Equations

Thanks to the operational matrices established in previous section, we now in position to obtain numerical solutions of the proposed coupled system (1) of FPDEs.With the help of operational matrices, we consider the approximations as(53)DxpcUx,y=AK2ΥK2x,y,DxpcVx,y=BK2ΥK2x,y.In view of Theorem 4, we can write (54)Ux,y-∑i=0neixi=AK2HK2×K2p,xΥK2x,y,Vx,y-∑i=0nbixi=BK2HK2×K2p,xΥK2x,y.Applying initial conditions, we get (55)Ux,y=AK2HK2×K2p,xΥK2x,y+XK21ΥK2x,y,Vx,y=BK2HK2×K2p,xΥK2x,y+XK22ΥK2x,y,where XK21ΥK2(x,y)=∑i=0nθi(y)xi and  XK22ΥK2(x,y)=∑i=0nψi(y)xi. In simple notation, we may write (56)AK2HK2×K2p,x+XK21=AK2^,BK2HK2×K2p,x+XK22=BK2^.Therefore (55) yields(57)Ux,y=AK2^ΥK2x,y,Vx,y=BK2^ΥK2x,y.Thanks to (57), the other terms of system (1) can be approximated as (58)Dyq1cUx,y=AK2^VK2×K2q1,yΥK2x,y,Dyq2cVx,y=BK2^VK2×K2q2,yΥK2x,y,Dxyp1cVx,y=BK2^ZK2×K2p1,x,yΥK2x,y,Dxyp2cUx,y=AK2^ZK2×K2p2,x,y,θx,y=XK23ΥK2x,y,ψx,y=XK24ΥK2x,y.Putting in (1), we get (59)AK2ΥK2x,y+ω1BK2^ZK2×K2p1,x,yΥK2x,y+ω2AK2^VK2×K2q1,yΥK2x,y=XK23ΥK2x,y,BK2ΥK2x,y+ρ1AK2^ZK2×K2p2,x,yΥK2x,y+ρ2BK2^VK2×K2q2,yΥK2x,y=XK24ΥK2x,y.Now we write matrix form of system (59) as(60)AK2ΥK2x,yBK2ΥK2x,y+ω1BK2^ZK2×K2p1,x,yΥK2x,yρ1AK2^ZK2×K2p2,x,yΥK2x,y+ω2AK2^VK2×K2γ2,yΥK2x,yρ2BK2^VK2×K2ρ2,yΥK2x,y-XK23ΥK2x,yXK24ΥK2x,y=O,which further gives (61)AK2BK2E^+AK2^BK2^OK2×K2ρ1ZK2×K2p2,x,yω1ZK2×K2p1,x,yOK2×K2E^+AK2^BK2^ω2VK2×K2q1,yOK2×K2OK2×K2ρ2VK2×K2q2,yE^-XK23XK24E^=OAK2BK2+AK2^BK2^ω2VK2×K2q1,yρ1ZK2×K2p2,x,yω1ZK2×K2p1,x,yρ2VK2×K2q2,y-XK23XK24=O,where(62)E^=ΥK2x,yOK2OK2ΥK2x,y.OK2 is a column null matrix and other null matrices are OK2×K2 and O.Simplifying (61), we get (63)AK2BK2-AK2^BK2^ω2VK2×K2q1,yρ1ZK2×K2p2,x,yω1ZK2×K2p1,x,yρ2VK2×K2q2,y-XK23XK24=O.Therefore by using AK2^ and BK2^, the aforesaid equation can be written as(64)AK2BK2-P^AK2BK2-XK21XK22P^-F1F2=O,where (65)P^=ω2HK2×K2p,xVK2×K2q1,yρ1HK2×K2p,xZK2×K2p2,x,yω2HK2×K2p,xZK2×K2p1,x,yρ2HK2×K2p,xVK2×K2q2,y,F1=ω2XK21VK2×K2q1,y+ω1XK22ZK2×K2p1,x,y-XK23,F2=ρ1XK21WK2×K2p2,x+ρ2XK22VK2×K2q2,y-XK24.Hence (64) is Sylvester type matrix equation, which on solving for unknown matrix AK2BK2 and using its value in (57), we get the required numerical solution of the proposed problem.Remark 10.
Here, we remark that we have used a machine type DESKTOP-V8125H8 with processor intel(R) Core (TM) i5-3210M CPU@ 2.50 GHz and installed memory (RAM) is 4 GB together with 64-bit operating system for computation and numerical simulations.


## 5. Numerical Test Problems

This section is concerning to the numerical test problems and their visualization.Problem 1.
Let the coupled system of FPDEs given as(66)Dx1.8cUx,y-Dxy2cVx,y-4Dy1.8cUx,y=θx,y,Dx1.8cVx,y-6Dxy2cUx,y+3Dy1.8cVx,y=ψx,y,U0,y=U′0,y=0,V0,y=V′0,y=0,such that the external functions θ(x,y) and ϕ(x,y) are given as(67)θx,y=27x2y29y-1x-127xy-2y-3-4x4y34+3y-0.016x2.5y4y-13125x2-175x+56,ϕx,y=-36x2y3y-1x-131-2x2y-1x-12-3x2y-12x-12+3y-12y-xy-1x-1-y-14x-4xy-1x-12-3xyy-1x-1-0.071x1.2y3y-121250x3-2625x2+1680x-308.The exact solution of the above system is (68)Ux,y=xy1-x21-y3,Vx,y=xy1-yxy-x2y-xy22.Evaluate the approximate solution of Problem 1 with our proposed method. After applying the method to the given problem, we see from Figure 1 that the considered scheme provides close agreement between numerical and exact solution. The comparison between exact and approximate solution is shown in Figure 1, while the absolute error corresponding to scale level K=10 is provided in Figure 2. Further, to check the efficiency of the method, we have also computed absolute error at various scale level and different points of the spaces as given in Table 1.
Further in Table2, we compute the CPU time for the computational of solutions for the test Problem 1.Table 1
Absolute error at various values of(x,y) for K=10,12 in U(x,y) and V(x,y) of Problem 1.
(x,y)
U
-
U
10
U
-
U
12
V
-
V
10
V
-
V
12
(
0.1,0.1
)
1.707
×
1
0
-
5
1.005
×
1
0
-
8
1.3250
×
1
0
-
6
8.006
×
1
0
-
9
(
0.1,0.5
)
1.956
×
1
0
-
5
6.006
×
1
0
-
8
1.596
×
1
0
-
6
7.508
×
1
0
-
9
(
0.1,0.9
)
2.236
×
1
0
-
5
5.803
×
1
0
-
8
4.773
×
1
0
-
6
6.009
×
1
0
-
9
(
0.5,0.1
)
1.873
×
1
0
-
5
4.600
×
1
0
-
8
5.385
×
1
0
-
6
7.036
×
1
0
-
9
(
0.5,0.5
)
1.887
×
1
0
-
5
3.008
×
1
0
-
8
7.858
×
1
0
-
6
5.005
×
1
0
-
9
(
0.5,0.9
)
2.092
×
1
0
-
5
2.519
×
1
0
-
8
6.619
×
1
0
-
6
4.001
×
1
0
-
9
(
0.9,0.1
)
1.973
×
1
0
-
5
2.363
×
1
0
-
8
5.363
×
1
0
-
6
2.506
×
1
0
-
9
(
0.9,0.5
)
2.907
×
1
0
-
5
5.009
×
1
0
-
8
6.009
×
1
0
-
6
3.005
×
1
0
-
9
(
0.9,0.9
)
2.306
×
1
0
-
5
6.000
×
1
0
-
8
5.000
×
1
0
-
6
2.500
×
1
0
-
9Table 2
Absolute error and CPU times at various values of(x,y) for K=12 in U(x,y) and V(x,y) of Problem 1.
(x,y)
CPU time
U
-
U
12
CPU time
V
-
V
12
(0.1,0.1)
55.78
s
1.005
×
1
0
-
8
57.88
s
8.006
×
1
0
-
9
(0.1,0.5)
56.71
s
6.006
×
1
0
-
8
57.98
s
7.508
×
1
0
-
9
(0.1,0.9)
56.91
s
5.803
×
1
0
-
8
57.99
s
6.009
×
1
0
-
9
(0.5,0.1)
56.98
s
4.600
×
1
0
-
8
58.70
s
7.036
×
1
0
-
9
(0.5,0.5)
57.31
s
3.008
×
1
0
-
8
58.72
s
5.005
×
1
0
-
9
(0.5,0.9)
57.62
s
2.519
×
1
0
-
8
58.88
s
4.001
×
1
0
-
9
(0.9,0.1)
58.11
s
2.363
×
1
0
-
8
59.23
s
2.506
×
1
0
-
9
(0.9,0.5)
59.41
s
5.009
×
1
0
-
8
59.38
s
3.005
×
1
0
-
9
(0.9,0.9)
59.88
s
6.000
×
1
0
-
8
60.08
s
2.500
×
1
0
-
9Figure 1
Evaluation of numerical and exact solutionsU(x,y), V(x,y) corresponding to scale level K=10 of Problem 1.Figure 2
Absolute error inU(x,y) and V(x,y) using K=10 of Problem 1.Problem 2.
To support the aforesaid established results, we consider the following problem:(69)Dx1.8cUx,y-Dxy2cVx,y+11Dy0.99cUx,y=θx,y,Dx1.8cVx,y+Dxy2cUx,y-11Dy0.99cVx,y=ψx,y,U0,y=y,U′0,y=y,V0,y=0,V′0,y=exp⁡y.The external source functions are provided as (70)θx,y=yx1.8E1,-0.8x-expy+11.046y0.1,ψx,y=4.59exp⁡yx0.8+exp⁡xx-11xy-0.99E1,0.1,where the exact solution of the Problem 2 is Ux,y=yexp⁡x,Vx,y=xexp⁡(y). We approximate the given test problem by using the proposed method. In Figure 3, we have shown the comparison between exact and approximate solutions at scale level K=10. Further, to verify the efficiency of the proposed method, also the absolute error at scale level K=10 in U(x,y) and V(x,y) has been shown in Figure 4. To demonstrate the maximum absolute error at different scale levels and various values of the space variables x,y has been provided in Table 3. From Table 3, it is clear that the absolute error is below 10-8, which is a very small quantity. The absolute error may further be decreased by enlarging scale level. CPU time for test Problem 2 is provided in Table 4.Table 3
Maximum absolute error atK=10,12 of U(x,y) and V(x,y) using various values of (x,y) for Problem 2.
(x,y)
U
-
U
10
U
-
U
12
V
-
V
10
V
-
V
12
(0.1,0.1)
1.000
×
1
0
-
5
2.005
×
1
0
-
7
1.006
×
1
0
-
5
7.908
×
1
0
-
8
(0.1,0.5)
7.000
×
1
0
-
5
1.596
×
1
0
-
7
2.000
×
1
0
-
5
7.009
×
1
0
-
8
(0.1,0.9)
006
×
1
0
-
5
1.003
×
1
0
-
7
2.126
×
1
0
-
5
5.019
×
1
0
-
8
(0.5,0.1)
6.003
×
1
0
-
5
5.000
×
1
0
-
7
2.451
×
1
0
-
5
5.008
×
1
0
-
8
(0.5,0.5)
1.448
×
1
0
-
5
4.000
×
1
0
-
7
1.081
×
1
0
-
5
4.052
×
1
0
-
8
(0.5,0.9)
2.500
×
1
0
-
5
8.619
×
1
0
-
7
5.003
×
1
0
-
6
4.051
×
1
0
-
8
(0.9,0.1)
1.527
×
1
0
-
5
7.063
×
1
0
-
7
6.050
×
1
0
-
6
3.094
×
1
0
-
8
(0.9,0.5)
2.5307
×
1
0
-
6
4.129
×
1
0
-
7
8.008
×
1
0
-
6
2.096
×
1
0
-
8
(0.9,0.9)
1.706
×
1
0
-
6
2.002
×
1
0
-
7
9.004
×
1
0
-
6
6.009
×
1
0
-
8Table 4
Absolute error and CPU times at various values of(x,y) for K=12 in U(x,y) and V(x,y) of Problem 2.
(x,y)
CPU time
U
-
U
12
CPU time
V
-
V
12
(0.1,0.1)
51.68
s
2.005
×
1
0
-
7
53.78
s
7.908
×
1
0
-
8
(0.1,0.5)
51.71
s
1.596
×
1
0
-
7
53.90
s
7.009
×
1
0
-
8
(0.1,0.9)
52.81
s
1.003
×
1
0
-
7
54.19
s
5.019
×
1
0
-
8
(0.5,0.1)
52.98
s
5.000
×
1
0
-
7
55.20
s
5.008
×
1
0
-
8
(0.5,0.5)
53.31
s
4.000
×
1
0
-
7
56.72
s
4.052
×
1
0
-
8
(0.5,0.9)
53.62
s
8.619
×
1
0
-
7
56.38
s
4.051
×
1
0
-
8
(0.9,0.1)
54.11
s
7.063
×
1
0
-
7
56.63
s
3.094
×
1
0
-
8
(0.9,0.5)
55.41
s
4.129
×
1
0
-
7
56.88
s
2.096
×
1
0
-
8
(0.9,0.9)
56.88
s
2.002
×
1
0
-
7
57.04
s
6.009
×
1
0
-
8Figure 3
Evaluation of exact and numerical solutionsU(x,y), V(x,y) corresponding to scale level K=10 for test Problem 2.Figure 4
The maximum absolute error inU(x,y) and V(x,y) at K=10 of Problem 2.Problem 3.
Consider another coupled system of FPDEs as(71)Dx1.8cUx,y+Dxy2cVx,y+4Dy0.8cUx,y=fx,y,Dx1.8cVx,y+9Dxy2cVx,y+3Dy0.8cVx,y=gx,y,U0,y=-y5,U′0,y=0,V0,y=y5,V′0,y=0.While the external source term is given by(72)θx,y=4.54x2+y2x2.8+9.18x0.2+2xy3x-2+58.08x2+y2x2.5+10.632x1.2x3-y3,ψx,y=5.508y2-x2x1.2-1.836x3+y3x0.2+54xyx-y+82.62y1.2x2-y2-5.50x3+y36y0.2.The exact solution of Problem 3 at p=2 is(73)Ux,y=x-yx2+y22+x2+y2xy,Vx,y=x+y2y-xx2+y2-xy.Apply the considered method to the given Problem 3 corresponding to a scale level K=8. Plot the exact and approximate solution as in Figure 5. We see from Figure 5 that the adapted techniques work excellently and there is a close agreement between approximate and exact solutions. Further, we have computed maximum absolute error for Problem 3 in Table 5 at scale level K=10,12. The maximum absolute error corresponding to the scale level K=8 has been provided in Figure 6, which is below 10-3. The absolute error can be further decreased by enlarging scale level with various values of x,y.
CPU time for test Problem3 is provided in Table 6.Table 5
Maximum absolute error inU(x,y) and V(x,y) of Problem 3 at K=10,12 corresponding to various values of (x,y).
(x,y)
U
-
U
10
U
-
U
12
V
-
V
10
V
-
V
12
(0.1,0.1)
1.707
×
1
0
-
6
4.325
×
1
0
-
8
6.826
×
1
0
-
6
1.000
×
1
0
-
8
(0.1,0.5)
2.956
×
1
0
-
6
3.196
×
1
0
-
8
5.638
×
1
0
-
6
2.009
×
1
0
-
8
(0.1,0.9)
4.236
×
1
0
-
6
2.703
×
1
0
-
8
3.126
×
1
0
-
6
2.900
×
1
0
-
8
(0.5,0.1)
8.873
×
1
0
-
6
7.038
×
1
0
-
8
5.451
×
1
0
-
6
2.945
×
1
0
-
8
(0.5,0.5)
7.000
×
1
0
-
5
8.058
×
1
0
-
7
4.781
×
1
0
-
5
1.520
×
1
0
-
8
(0.5,0.9)
8.002
×
1
0
-
5
1.019
×
1
0
-
7
4.933
×
1
0
-
5
3.051
×
1
0
-
8
(0.9,0.1)
6.027
×
1
0
-
5
3.063
×
1
0
-
7
4.565
×
1
0
-
5
3.194
×
1
0
-
7
(0.9,0.5)
9.007
×
1
0
-
5
2.029
×
1
0
-
7
5.808
×
1
0
-
5
1.096
×
1
0
-
7
(0.9,0.9)
9.006
×
1
0
-
5
2.008
×
1
0
-
7
2.746
×
1
0
-
5
3.209
×
1
0
-
7Table 6
Maximum absolute error and CPU times inU(x,y) and V(x,y) of Problem 3 at K=12 corresponding to various values of (x,y).
(x,y)
CPU time
U
-
U
12
CPU time
V
-
V
12
(0.1,0.1)
55.70
s
4.325
×
1
0
-
8
55.26
s
1.000
×
1
0
-
8
(0.1,0.5)
55.95
s
3.196
×
1
0
-
8
55.38
s
2.009
×
1
0
-
8
(0.1,0.9)
56.23
s
2.703
×
1
0
-
8
56.12
s
2.900
×
1
0
-
8
(0.5,0.1)
56.87
s
7.038
×
1
0
-
8
57.45
s
2.945
×
1
0
-
8
(0.5,0.5)
57.00
s
8.058
×
1
0
-
7
75.78
s
1.520
×
1
0
-
8
(0.5,0.9)
57.02
s
1.019
×
1
0
-
7
57.93
s
3.051
×
1
0
-
8
(0.9,0.1)
57.27
s
3.063
×
1
0
-
7
58.56
s
3.194
×
1
0
-
7
(0.9,0.5)
58.00
s
2.029
×
1
0
-
7
59.80
s
1.096
×
1
0
-
7
(0.9,0.9)
59.68
s
2.008
×
1
0
-
7
60.76
s
3.209
×
1
0
-
7Figure 5
Evaluation between exact and approximate solutionsU(x,y),V(x,y) for scale level K=8 of test Problem 3.Figure 6
Maximum absolute error inU(x,y) and V(x,y) at K=8 for different values of x,y for Problem 3.


## 6. Conclusion

In this article we have developed an efficient numerical technique by using shifted Legendre polynomials. A new operational matrix of mixed partial derivativeZK2×K2(p,x,y) is introduced to solve a multiterm coupled system of FPDEs. We have also tested our method by some test problems. The comparison between the exact and approximate solutions of the problems has been provided. We observed that the maximum absolute error is very small even at very small scale level. This minimum value of the absolute error indicates that the method is highly accurate and can be applied to a large numbers of FPDEs.



---

*Source: 1535826-2017-09-20.xml*