# A Novel Method for Solution of Fractional Order Two-Dimensional Nonlocal Heat Conduction Phenomena

Authors: Hammad Khalil; Ishak Hashim; Waqar Ahmad Khan; Abuzar Ghaffari
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this paper, we have extended the operational matrix method for approximating the solution of the fractional-order two-dimensional elliptic partial differential equations (FPDEs) under nonlocal boundary conditions. We use a general Legendre polynomials basis and construct some new operational matrices of fractional order operations. These matrices are used to convert a sample nonlocal heat conduction phenomenon of fractional order to a structure of easily solvable algebraic equations. The solution of the algebraic structure is then used to approximate a solution of the heat conduction phenomena. The proposed method is applied to some test problems. The obtained results are compared with the available data in the literature and are found in good agreement.
Dedicated to my father Mr. Sher Mumtaz, (1955-2021), who gave me the basic knowledege of mathematics.


---

## Body

## 1. Introduction

The first approach to study a physical experiment is to derive a mathematical expression, which formulates the dynamics of the experiment under certain assumption. Most of the physical phenomena are formulated in terms of ordinary differential equations. Some problems in which the quantity of interest also changes with respect to both space and time result in partial differential equations. A wide range of scientists devoted their very precise time to investigate various important aspects of partial differential equations. In [1], thermoelastic damping of in-plane vibration of a functionally graded material has been studied based on the Eringen nonlocal theory. In [2], a fractional sideways heat flow problem is investigated, in which the interior measurements at two interior points are given by continuous data with deterministic noises. The work in [3] deals with the exothermic reactions model having a constant heat source in the porous media with strong memory effects. This article explains the behavior of heat profile under the effect of different definitions of the derivative. In [4], the Newtonian liquid flow porous stretching/shrinking sheet utilizing a Brinkman mode is investigated.Nonlocal partial differential equations (PDEs) arise in the mathematical modeling of various problems in physics, engineering, ecology, and biological sciences [5–7]. The term nonlocal problems means that the solution of PDEs on the boundary is connected with the solution on some interior points of the domain. The case arises when the solution at the boundary is not known. Such formulation is placed in a separate class known as nonlocal boundary value problems. Some of the numerical investigations regarding PDEs with nonlocal boundary conditions reported in the literature can be found in [8–15]. Among others, some of the well-known methods that can be effectively applied to BVPs are finite difference methods, mesh-free methods, finite element methods, etc. For instance, the one-dimensional heat equation with nonlocal boundary conditions has been studied in [12, 16–18]. Two-dimensional diffusion problems with nonlocal boundary conditions have been discussed in [13, 19]. The numerical solution of the Laplace equation with integral boundary condition is explored in [8]. Similarly, the numerical solution of multidimensional linear elliptic equations with integral boundary conditions is explored in [20].To be specific, the fundamental problem of interest in this article is to find an approximate solution of the fractional-order two-dimensional Poisson equation, given as follows:(1)∂σu∂xσ+∂σu∂yσ=−fx,y.The above model is subject to two-point nonlocal boundary conditions:(2)ux,0=μ3x,ux,1=μ4x,0≤x≤1,u0,y=μ1y,u1,y=γuξ,y+μ2y,0≤y≤1,0≤ξ<1.The functionsfx,y,μ1y,μ2y,μ3x and μ4x are given smooth functions. The parameters ξ and γ are two positive constants. The parameter σ1<σ≤2 represents the order of derivative defined in Caputo sense.Recently, many authors devoted their studies to the approximate solution of integer order version of the above problem. In [21], Yang et al. approximated the solution of fractional order partial differential subjected to the simple initial condition of the form:(3)ux,0=gxu0,t=ht.The approach presented in [21], is interesting. However, it can not be utilized directly to the approximate solution of fractional order PDEs subject to nonlocal boundary condition (2). Islam et al. [22] presented a comprehensive text on the solution to the above problem. They implemented two different methods for the solution of integer order counterpart of (1). Their first approach is based on Haar wavelets. In the same paper, they also implemented a modified form of the mesh-less method to solve the integer order problem. Sajavicius [23] implemented the radial base function approach to integer order problems and studied some computational aspects of the proposed approach. The results presented in [22, 23] are the motivating factor of our interest to study the approximate solution of the fractional-order Poisson equation subject to nonlocal boundary conditions.Our approach is based on shifted Legendre polynomials and their operational matrices. We derived some new operational matrices to handle the problem. The new operational matrices can handle the nonlocal boundary conditions. The interesting readers may find useful results and some new strategies of this method in [24–27]. Application of orthogonal polynomials combined with Tau and Collection method can be found in [28–31].The rest of the paper is organized as follows: in Section2, we recall some primary results from the fractional calculus and approximation theory. In Section 3, we recall some previously derived operational matrices and develop some new operational matrices. In Section 4, a theoretical base is developed for the conversion of the nonlocal FPDEs to the matrix equation. Convergence analysis and error estimation are also developed in the same section. In Section 5, the proposed method is applied to some benchmark problems. In the same section, the obtained results are demonstrated and compared with other methods in the literature. Section 6 is devoted to the conclusions.


## 2. Preliminaries

In this section, we present some useful results and notations which are of primary importance in our further investigation.Definition 1.
(see [32–34]). Given an interval a,b⊂ℝ, the Riemann–Liouville fractional order integral of a function ϕ∈L1a,b,ℝ of order α∈ℝ+ is defined by the following:(4)ℐa+αϕt=1Γα∫att−sα−1ϕsds,provided that the integral on the right-hand side exists.Definition 2.
(Caputo derivative). For a given functionϕx∈Cna,b, the Caputo fractional order derivative is defined as follows:(5)Dαϕx=1Γn−α∫axϕntx−tα+1−ndt,n−1≤α<n,n∈N,where n=α+1.Hence, it follows that(6)Dαxk=Γ1+kΓ1+k−αxk−α,Iαxk=Γ1+kΓ1+k+αxk+α,DαC=0,for a constantC.


### 2.1. The Shifted Legendre Polynomials

The shifted Legendre polynomials [35] defined on 0,1 are given by the following relation:(7)Pix=∑k=0iΔi,kxk,whereΔi,k=−1i+ki+k!i−k!k!2.These polynomials are bounded by 1, and we have the following relation:(8)maxx∈0,1Pix=1.These polynomials are orthogonal on the domain0,1, and the orthogonality condition is given as follows:(9)∫01PixPjxdx=12i+1,ifi=j,0,ifi≠j.dydxwhich implies that anyfx∈C0,1 can be approximated by Legendre polynomials as follows: (10)fx≈∑a=0mCaPax,whereCa=2a+1∫01fxPaxdx.In vector notation, we write the following:(11)fx≈KMTP^Mx,where M=m+1 is the scale level of approximation. K is the coefficient vector and P^Mx is M terms function vector. These notations can be easily extended to two-dimensional space [35] and two-dimensional Legendre polynomials of the order M are defined as a product function of two Legendre polynomials (12)Pnx,y=PaxPby,n=Ma+b+1,a=0,1,2,…,m,b=0,1,2,…,m.The orthogonality condition ofPnx,y is as follows:(13)∫01∫01PaxPbyPcxPdydxdy=12a+12b+1,ifa=c,b=d;0,otherwise.Anyfx,y∈C0,1×0,1 can be approximated by the polynomials Pnx,y as follows:(14)fx,y≈∑a=0m∑b=0mCabPaxPby,whereCab=2a+12b+1∫01∫01fx,yPaxPbydxdy.For simplicity, use the notationCn=Cab where n=Ma+b+1 and rewrite (14) as follows:(15)fx,y≈∑n=1M2CnPnx,y=KM2Ψx,y,where KM2 is 1×M2 coefficient row vector and is M2×1 column vector of functions defined by the following:Ψx,y(16)Ψx,y=ψ11x,y,…,ψ1Mx,yψ21x,y,…,ψ2Mx,y,…,ψMMx,yT,where ψi+1,j+1x,y=PixPjy.


## 2.1. The Shifted Legendre Polynomials

The shifted Legendre polynomials [35] defined on 0,1 are given by the following relation:(7)Pix=∑k=0iΔi,kxk,whereΔi,k=−1i+ki+k!i−k!k!2.These polynomials are bounded by 1, and we have the following relation:(8)maxx∈0,1Pix=1.These polynomials are orthogonal on the domain0,1, and the orthogonality condition is given as follows:(9)∫01PixPjxdx=12i+1,ifi=j,0,ifi≠j.dydxwhich implies that anyfx∈C0,1 can be approximated by Legendre polynomials as follows: (10)fx≈∑a=0mCaPax,whereCa=2a+1∫01fxPaxdx.In vector notation, we write the following:(11)fx≈KMTP^Mx,where M=m+1 is the scale level of approximation. K is the coefficient vector and P^Mx is M terms function vector. These notations can be easily extended to two-dimensional space [35] and two-dimensional Legendre polynomials of the order M are defined as a product function of two Legendre polynomials (12)Pnx,y=PaxPby,n=Ma+b+1,a=0,1,2,…,m,b=0,1,2,…,m.The orthogonality condition ofPnx,y is as follows:(13)∫01∫01PaxPbyPcxPdydxdy=12a+12b+1,ifa=c,b=d;0,otherwise.Anyfx,y∈C0,1×0,1 can be approximated by the polynomials Pnx,y as follows:(14)fx,y≈∑a=0m∑b=0mCabPaxPby,whereCab=2a+12b+1∫01∫01fx,yPaxPbydxdy.For simplicity, use the notationCn=Cab where n=Ma+b+1 and rewrite (14) as follows:(15)fx,y≈∑n=1M2CnPnx,y=KM2Ψx,y,where KM2 is 1×M2 coefficient row vector and is M2×1 column vector of functions defined by the following:Ψx,y(16)Ψx,y=ψ11x,y,…,ψ1Mx,yψ21x,y,…,ψ2Mx,y,…,ψMMx,yT,where ψi+1,j+1x,y=PixPjy.


## 3. New Operational Matrices

The operational matrices of the fractional derivatives and integrals play a vital role in converting the FPDEs to the system of algebraic equations. The operational matrices of all derivatives are explicitly derived in our previous report [35]. We will need operational matrices in integration. The operational matrices of integration w.r.t x or y is not a difficult task and can be easily derived using the same procedure as in [35]. To make this study a self-contained material and for the ease and interest of our readers, we have provided detailed proof of deriving operational matrices of integration.Lemma 1.
LetΨx,y be the function vector as defined in (16), then the fractional integral of order σ of Ψx,y w.r.t y is given by the following:(17)IyσΨx,y≃GM2×M2σ,yΨx,y,where GM2×M2σ,y is the operational matrix of the fractional integration of order σ and is defined as follows:(18)GM2×M2σ,y=Δq,r′,where q=Mi+j+1, r=Ma+b+1, for i,j,a,b=0,1,2,…,m and(19)Δq,r′=Ci,j,b,a,k=∑k=0aδi,a2j+1∑l=0j−1j+l+b+kj+l!b+k!j−l!l!2k+l+σ+1b−k!k!Γk+σ+1.Proof.
Taking the elementPnx,y defined by (12), then the fractional-order integration of Pnx,y w.r.t y follows:(20)IyσPnx,y=Pax∑k=0b−1b+kb+k!b−k!k!2Iyσyk=∑k=0b−1b+kb+k!b−k!k!2PaxIyσyk.
Using the definition of fractional integration, we obtain the following:(21)IyσPaxPby=∑k=0b−1b+kb+k!b−k!k!Γk+σ+1Paxyk+σ,b=0,1,…,m,
ApproximatingPaxyk+σ by M terms of Legendre polynomials in two variables yields(22)Paxyk+σ≈∑i=0m∑j=0mCijPixPjy,where Cij=2i+12j+1∫01∫01Paxyk+σPixPjydxdy, which in view of the orthogonality conditions implies that(23)Cij,a=δi,a2j+1∑l=0j−1j+lj+l!j−l!l!2k+l+σ+1,where(24)δia=1,ifi=a,0,ifi≠a,and hence, it follows that(25)IyσPaxPby≈∑k=0b−1b+kb+k!b−k!k!Γk+σ+1∑i=0m∑j=0mCij,aPixPjy≈∑i=0m∑j=0m∑k=0b−1b+kb+k!b−k!k!Γk+σ+1Cij,aPixPjy≈∑i=0m∑j=0mCij,a,b,kPixPjy,b=0,1,…,M,where(26)Cij,a,b,k=∑k=0bδi,a2j+1∑l=0j−1j+l+b+kj+l!b+k!j−l!l!2k+l+σ+1b−k!k!Γk+σ+1.
Using the notations,q=Mi+j+1, r=Ma+b+1 and Δq,r′=Ci,j,b,a,k for i,j,a,b=0,1,2,3,…,m, we get the desired result.Lemma 2.
LetΨx,y be as defined in (16), then the integration of order σ of Ψx,y w.r.t x is given by the following:(27)IxσΨx,y≃GM2×M2σ,xΨx,y,where GM2×M2σ,x is the operational matrix of derivative of order σ and is defined as follows:(28)GM2×M2σ,x=Θq,r′,and r=Mi+j+1, q=Ma+b+1, Θq,r′=Si,j,b,a,k for i,j,a,b=0,1,2,…,m and(29)Si,j,b,a,k=∑k=0aδj,b2i+1∑l=0i−1i+l+a+ki+l!a+k!i−l!l!2k+l+σ+1a−k!k!Γk+σ+1.Proof.
The proof of this lemma is similar to the above lemma.The operational matrices derived in the above two lemmas are essential for our further analysis. However, these matrices are not enough to fulfill our requirements. In our analysis, we will face terms likecxn1/Γσ∫0ξξ−sσ−1us,yds and cyn1/Γσ∫0ξξ−sσ−1ux,sds, where c and n are some positive constants and 0<ξ≤1. Therefore to replace such terms with their equivalent matrix form, we need to derive two more operational matrices. The operational matrices used to replace such term by their equivalent matrix form are derived in the following lemmas.Lemma 3.
Letux,y=KM2Ψx,y, then for some constants c and n, the following relation holds:(30)cyn1Γσ∫0ξξ−sσ−1ux,sds=KM2PM2×M2c,n,σ,ξ,yΨx,y,where(31)PM2×M2c,n,σ,ξ,y=Δ1,1Δ1,2⋯Δ1,M2Δ1,2Δ2,2⋯Δ2,M2⋮⋮⋱⋮ΔM2,1ΔM2,2⋯ΔM2,M2,Δr,q=da,b,i,jc,n,σ,ξ, where r=Mi+j+1, q=Ma+b+1 and(32)da,b,i,jc,n,σ,ξ=δa,ic2b+1∑k=0j∑l=0bΔj,k−1b+lb+l!Γk+1ξσ+kb−l!l!2n+l+1Γσ+k.Proof.
Let(33)ux,y=KM2Ψx,y.
For some constantc and n, we can write the following expression:(34)cyn1Γσ∫0ξξ−sσ−1ux,sds=KM2cyn1Γσ∫0ξξ−sσ−1Ψx,sds.
Now, considering the general term ofΨx,s, we can write the following:(35)cyn1Γσ∫0ξξ−sσ−1Prx,sds=cynPix1Γσ∫0ξξ−sσ−1Pjsds,where we use the notation r=Mi+j+1. By using the definition of Legendre polynomials, we can write the following:(36)cyn1Γσ∫0ξξ−sσ−1Prx,sds=∑k=0jcΔj,kΓk+1ξσ+kΓσ+kynPix.
ApproximatingPixyn by Legendre polynomials in two variables yields(37)Pixyn≈∑a=0m∑b=0mdabPaxPby,where dab=2a+12b+1∫01∫01PixynPaxPbydxdy, which in view of the orthogonality conditions implies that(38)dab=δa,i2b+1∑l=0b−1b+lb+l!b−l!l!2n+l+1,where(39)δa,i=1,ifi=a,0,ifi≠a.
Hence, it follows that(40)cyn1Γσ∫0ξξ−sσ−1Prx,sds≈∑k=0jcΔj,kΓk+1ξσ+kΓσ+k∑a=0m∑b=0mdabPaxPby≈∑a=0m∑b=0m∑k=0jcΔj,kΓk+1ξσ+kΓσ+kdabPaxPby≈∑a=0m∑b=0mda,b,i,j,kc,n,σ,ξPaxPby,where(41)da,b,i,jc,n,σ,ξ=δa,ic2b+1∑k=0j∑l=0bΔj,k−1b+lb+l!Γk+1ξσ+kb−l!l!2n+l+1Γσ+k.
Using the notations,r=Mi+j+1, q=Ma+b+1 and Δr,q=da,b,i,jc,n,σ,ξ for i,j,a,b=0,1,2,3,…,m, we get the desired result.Lemma 4.
Letux,y=KM2Ψx,y, then for some constants c and n, the following relation holds:(42)cxn1Γσ∫0ξξ−sσ−1us,yds=KM2PM2×M2c,n,σ,ξ,xΨx,y,where(43)PM2×M2c,n,σ,ξ,x=Δ1,1Δ1,2⋯Δ1,M2Δ1,2Δ2,2⋯Δ2,M2⋮⋮⋱⋮ΔM2,1ΔM2,2⋯ΔM2,M2,Δr,q=da,b,i,jc,n,σ,ξ, where q=Mi+j+1, r=Ma+b+1 and(44)da,b,i,jc,n,σ,ξ=δb,jc2b+1∑k=0i∑l=0aΔi,k−1a+la+l!Γk+1ξσ+ka−l!l!2n+l+1Γσ+k.Proof.
This lemma can be easily proved by following similar steps as in the previous lemma. This is left as an exercise for interested readers.


## 4. Main Result: Application of Operational Matrices

The operational matrices developed in the previous section have a wide range of applications. As an application of the above matrices, we solve the following fractional order Poisson equation:(45)∂σu∂xσ+∂σu∂yσ=−fx,y,subject to nonlocal two-point boundary conditions(46)ux,0=μ3x,ux,1=μ4x,0≤x≤1,u0,y=μ1y,u1,y=γuξ,y+μ2y,0≤y≤1,0≤ξ<1.Readers may see how simple steps lead us to the approximate solution of such complicated problems with high accuracy. As usual, we seek the solution to the problem in terms of shifted Legendre polynomials given by the following:(47)∂σu∂yσ=KM2Ψx,y.On application of the fractional integral of orderσ and making use of Lemma 2, we get the following relation:(48)ux,y=KM2GM2×M2σ,yΨx,y+c0+c1y.Using the conditions aty=0 and y=1, we get the following relation:(49)c0=μ3x,c1=μ4x−KM2GM2×M2σ,yΨx,1+μ3x.Using the values ofc0 and c1 in (48), we get the following:(50)ux,y=KM2GM2×M2σ,yΨx,y+μ3x1−y+μ4xy−yKM2GM2×M2σ,yΨx,1.which in view of Lemma3 can be written as follows:(51)ux,y=KM2GM2×M2σ,y+P−1,1,σ,1,yΨx,y+F1M2Ψx,y,where F1M2Ψx,y=μ3x1−y+μ4xy.Now approximating the source termfx,y=FM2Ψx,y, and using (47) in (45), we can write the following:(52)∂σu∂xσ=−KM2Ψx,y−FM2Ψx,y.On application of fractional integral of orderσ w.r.t. x, we can write the following:(53)ux,y=−KM2GM2×M2σ,xΨx,y−FM2GM2×M2σ,xΨx,y+d0+d1x.Using the initial conditions atx=0 we can easily d0=μ1y; however, the second constant is not known. We can use the two-point boundary conditions:(54)u1,y=−KM2GM2×M2σ,xΨ1,y−FM2GM2×M2σ,xΨ1,y+μ1y+d1,uξ,y=−KM2GM2×M2σ,xΨξ,y−FM2GM2×M2σ,xΨξ,y+μ1y+d1ξ.Using the equalityu1,y=γuξ,y+μ2y, we get the following:(55)−KM2GM2×M2σ,xΨ1,y−FM2GM2×M2σ,xΨ1,y+μ1y+d1=−γKM2GM2×M2σ,xΨξ,y−γFM2GM2×M2σ,xΨξ,y+γμ1y+d1γξ+μ2y.From which we can calculate the value ofd1 as follows:(56)d1=r1KM2GM2×M2σ,xΨ1,y−r1γKM2GM2×M2σ,xΨξ,y+r1FM2GM2×M2σ,xΨ1,y−r1γFM2GM2×M2σ,xΨξ,y+γ−1r1μ1y+r1μ2y,where r1=1/1−γξ≠0. Using the value of d0 and d1 in (53), we can write the following:(57)ux,y=−KM2GM2×M2σ,xΨx,y−FM2GM2×M2σ,xΨx,y+μ1y+r1xKM2GM2×M2σ,xΨ1,y−r1xγKM2GM2×M2σ,xΨξ,y+r1xFM2GM2×M2σ,xΨ1,y−r1xγFM2GM2×M2σ,xΨξ,y+γ−1r1μ1y+r1μ2yx,which in view of Lemma 4, can be written as follows:(58)x,y=−KM2GM2×M2σ,xΨx,y−FM2GM2×M2σ,xΨx,y+KM2PM2×M2r1,1,σ,1,xΨ,y−KM2PM2×M2r1γ,1,σ,ξ,xΨx,y+FM2PM2×M2r1,1,σ,1,xΨx,y−FM2PM2×M2r1γ,1,σ,ξ,xΨx,y+F2M2Ψx,y,where F2M2Ψx,y=μ1y+γ−1r1μ1y+r1μ2yx. In simplified notation, we can write the following:(59)ux,y=KM2PM2×M2r1,1,σ,1,x−GM2×M2σ,x−PM2×M2r1γ,1,σ,ξ,xr1γ,1,σ,ξ,xΨx,y+FM2PM2×M2r1,1,σ,1,x−GM2×M2σ,x−PM2×M2r1γ,1,σ,ξ,xΨx,y+F2M2Ψx,y.On comparing equation (51) and (59), we can write the following:(60)KM2GM2×M2σ,y+PM2×M2−1,1,σ,1,yΨx,y=KM2PM2×M2r1,1,σ,1,x−GM2×M2σ,x−PM2×M2r1γ,1,σ,ξ,xΨx,y+FM2PM2×M2r1,1,σ,1,x−GM2×M2σ,x−PM2×M2r1γ,1,σ,ξ,xΨx,y+F2M2Ψx,y−F1M2Ψx,y.Canceling out the common term, we can write the following:(61)KM2GM2×M2σ,y+PM2×M2−1,1,σ,1,y=KM2PM2×M2r1,1,σ,1,x−GM2×M2σ,x−PM2×M2r1γ,1,σ,ξ,x+FM2PM2×M2r1,1,σ,1,x−GM2×M2σ,x−PM2×M2r1γ,1,σ,ξ,x+F2M2−F1M2.which is a linear system of equations and can be easily solved for the unknown vector KM2, which can be used in (51) or (59) to get an approximate solution to the problem.


### 4.1. Error Bound

Considering a sufficiently smooth functiongx,y on Δ=0,η×0,η, let ∏M,Mx,y is the space span by m term Legendre polynomials. We assume that gM,Mx,y is its best approximation in ∏M,Mx,y. For this purpose, consider a polynomial P^M,Mx,y is any polynomial of degree ≤M in variable x and y, respectively. Then, from the definition of best approximation,(62)gx,y−gM,Mx,y2≤gx,y−PM,Mx,y2.The inequality in (62) also holds if PM,Mx,y is interpolating polynomial at point xi,yj; then by the similar arguments as in [36], the error of the approximation is given by the following:(63)gx,y−PM,Mx,y2≤C1+C2+C31MM+11MM+1,where(64)C1=14maxx,y∈0,1×0,1∂M+1∂xM+1gx,y,C2=14maxx,y∈0,1×0,1∂M+1∂yM+1gx,y,C3=116maxx,y∈0,1×0,1∂2M+2∂xM+1∂yM+1gx,y.We refer the reader to [37] for the proof of the above result. From the above result, it is clear that the error of approximation of a function decreases with the increase of M.


## 4.1. Error Bound

Considering a sufficiently smooth functiongx,y on Δ=0,η×0,η, let ∏M,Mx,y is the space span by m term Legendre polynomials. We assume that gM,Mx,y is its best approximation in ∏M,Mx,y. For this purpose, consider a polynomial P^M,Mx,y is any polynomial of degree ≤M in variable x and y, respectively. Then, from the definition of best approximation,(62)gx,y−gM,Mx,y2≤gx,y−PM,Mx,y2.The inequality in (62) also holds if PM,Mx,y is interpolating polynomial at point xi,yj; then by the similar arguments as in [36], the error of the approximation is given by the following:(63)gx,y−PM,Mx,y2≤C1+C2+C31MM+11MM+1,where(64)C1=14maxx,y∈0,1×0,1∂M+1∂xM+1gx,y,C2=14maxx,y∈0,1×0,1∂M+1∂yM+1gx,y,C3=116maxx,y∈0,1×0,1∂2M+2∂xM+1∂yM+1gx,y.We refer the reader to [37] for the proof of the above result. From the above result, it is clear that the error of approximation of a function decreases with the increase of M.


## 5. Test Problems

We solve the fractional order generalization of some bench mark problems from [22, 23].Test Problem 1.
(see [22, 23]). Consider (45) with the following functions:(65)fx,y=−2ex+y,μ1y=ey,μ2y=e1+y−γeξ+y,μ3x=ex,μ4x=ex+1.
The exact solution of the problem for fixσ=2 is ex+y.Test Problem 2.
(see [22, 23]). Consider (45) with the following functions:(66)fx,y=−2π2sinπxsinπy,μ1y=0,μ2y=−sinπxsinπy,μ3x=0,μ4x=0.
The exact solution of the problem for the fixσ=2 is sinπxsinπy.Test Problem 3.
(see [22, 23]). Consider (45) with the following functions.(67)fx,y=−Γ5Γ5−σx4−σy−y4−σ−Γ4Γ4−σx3−σ3y−1−4y3−σ+3Γ3Γ3−σx2−σy−y2−σ,μ1y=y3−yy−13,μ2y=y16+78,μ3x=x3,μ4x=xx−13+x3+1.
The exact solution of the problem for1<σ≤2 is x3−yy−13+y3+xyx−13.


## 6. Results and Discussion

We solve the above problems with the proposed method. The first two problems are selected from [22, 23], while the third problem is a constructed problem. In [22, 23], these problems are studied and solved with two different methods, Haar wavelets and a family of the mesh-less method based on the radial base functions. We solved these problems with the operational matrices and compared our results with the results reported in these references. To measure accuracy, we calculate the following parameters:(68)L∞=maxj=1,2,…,Muxj,yj−u¯xj,yj,ΔyΔxLrms=∑j=0Muxj,yj−u¯xj,yj2M.Ifσ≠2, then the exact solution of the first two problems is not known. We use the residual error norms to measure the accuracy of the proposed method for the fractional values of σ. These residual norms are defined as follows:(69)R∞=maxj=1,2,…,Mrxj,yj,Rrms=∑j=0Mrxj,yjM,where rx,y is defined as follows:(70)rx,y=∂σu¯∂xσ+∂σu¯∂yσ+fx,y.The first problem is solved using Haar wavelets (HWCM), mesh-less method without splitting (MCTMQ), and with splitting (MCTSMQ). We fixedσ=2 and obtained an approximate solution of Test Problem 1 for different values of scale level M. We compared our results with HWCM, MCTMQ, and MCTSMQ. The comparison of Lrms of Test Problem 1 obtained with the proposed method and HWCM, MCTMQ, and MCTSMQ are shown in Table 1. We can easily note that Lrms obtained with HWCM at M=16 is 8.2022×10−3; note that at this level, HWCM converts the problem to a system of algebraic equations of 1152 unknowns. While the proposed method yields Lrms=1.3381×10−10, while converting the problem to a system of algebraic equations of 81 unknowns. Lrms of this problem obtained with MCTMQ and MCTSMQ at N=64 is 4.1256×10−4 and 3.2737×10−7, respectively. This shows the superiority of the proposed method over HWCM and meshless methods.Table 1
Lrms of Test Problem 1 obtained with the proposed method (PM) at different values of M and its comparison with HWCM [22], MCTMQ [22], and MCTSMQ [22].
HWCM [22]MCTMQ [22]MCTSMQ [22]Proposed methodM (n)LrmsNLrmsLrmsM (n)Lrms1124.0860×10−341.7629×10−41.2407×10−45251.0404×10−42321.7549×10−383.4559×10−65.5859×10−66365.0660×10−64986.4795×10−4162.5542×10−53.6689×10−67491.6107×10−783202.3142×10−4323.6732×10−46.6522×10−78644.3049×10−91611528.2022×10−3644.1256×10−43.2737×10−79811.3381×10−103243252.9017×10−3Not available101003.9221×10−12Here, we fixξ=0.5 and γ=1.The parametersL∞ for Test Problem 1 obtained with the proposed method are also compared with HWCM, MCTMQ, and MCTSMQ. The results are displayed in Table 2. It is observed that L∞ for this problem obtained with the proposed method at M=8 is 3.8661×10−10 while HWCM yields 2.5230×10−4 at M=16, and meshless method MCTMQ and MCTSMQ yield 1.4665×10−3 and 2.4120×10−6.Table 2
L∞ of Example 5.1 obtained with proposed method (PM) at different values of M and its comparison with HWCM [22], MCTMQ [22], and MCTSMQ [22].
HWCM [22]MCTMQ [22]MCTSMQ [22]Proposed methodM(n)L∞NL∞L∞M(n)L∞1123.8454×10−344.9876×10−42.8688×10−45252.9475×10−42321.4767×10−381.1276×10−51.5987×10−56361.4687×10−54983.9916×10−4167.7560×10−51.7883×10−57499.7986×10−783201.0055×10−4321.0744×10−33.4368×10−68649.9163×10−91611522.5230×10−4641.4665×10−32.4120 × 10−69813.8661×10−103243256.3147×10−4Not available101001.2467×10−11Here, we fixξ=0.5 and γ=1.The proposed method along with HWCM and the meshless method convert the problem to the system of linear algebraic equations. The computational cost and stability of the resulting algebraic equations are different for different methods. Often some method yields a very approximate solution, but the computational cost is much higher. We compared the condition numberκ and CPUtime of the proposed method with MCTMQ and MCTSMQ. It is observed that the proposed method is more robust than these methods. At M=64, MCTMQ solves the algebraic equations in 53.78 seconds, and MCTSMQ takes 51.72 seconds to solve the system, while the proposed method solves the system in 0.09516 seconds. The condition number of the proposed method is much less than MCTMQ and MCTSMQ. It means that the proposed method converts the problem to the system of algebraic equations, which is more stable as compared to MCTMQ and MCTSMQ. The comparison of CPU time and conditions number of the proposed method with MCTQM and MCTSMQ at different scale levels is shown in Table 3.Table 3
CPU time and condition numberκ of the proposed method and its comparison with MCTMQ [22] and MCTSMQ [22].
MCTMQ [22]MCTSMQ [22]Proposed methodNκCPU timeκCPU timeMκCPU time42.8541×10100.003.0574×10120.0061.3472×10170.0025589.3912×10180.004.7522×10160.0072.6632×10160.00257162.2347×10190.044.3691×10180.0284.4841×10160.00371325.4623×10210.865.7284×10190.7896.8584×10160.01689644.4325×102153.783.2699×102051.72101.3208×10170.09516The accuracy of the present method is analyzed at different values ofγ. We chose γ=−8,2,8 and calculate L∞ and Lrms for scale levels. We observed that the accuracy of the proposed method does not depend on the values of γ. The values of L∞ and Lrms obtained with the proposed method at different values of γ using scale level M=10 are compared with the L∞ and Lrms obtained with the HWCM at M=32 in Table 4. One can see that the accuracy obtained with the proposed method is very high as compared to HWCM. The error norms at different values of γ are also compared with MCTSMQ and the results are displayed in Figure 1. One can note that the accuracy remains the same for all values of γ, also at M=8 and M=9, the error norms are much less than the error norms obtained using MCTSMQ at N=64. The CPU time at different values of γ and M=10 are shown in Figure 2(a). While the condition number of the resulting matrix equation at different values of γ is shown in Figure 2(b), also in the same figure, we plot the condition number obtained with MCTSMQ. We see that the condition number for MCTSMQ is approximately equal to 1019, while the condition number of the proposed method is approximately equal to 1017, which guarantees the robustness and stability of the proposed method as compared to MCTSMQ. From the above observations, we see that the proposed method provides a very accurate estimate of the solution of the problem.Table 4
L∞ and Lrms of Example 5.1 for different values of γ at M=10 and its comparison with HWCM [22] M=32.
γHaar wavelets [22]Proposed methodL∞LrmsL∞Lrms−86.1545×10−61.5072×10−53.5506×10−118.5152×10−12−65.8644×10−61.4942×10−53.5660×10−113.4538×10−12−45.2613×10−61.4837×10−53.1595×10−116.7567×10−12−23.8832×10−51.5076×10−51.4811×10−113.6505×10−1204.3925×10−51.9345×10−57.4749×10−122.9916×10−1223.5678×10−59.5986×10−52.5216×10−118.7175×10−1243.3974×10−56.2947×10−53.3240×10−115.7501×10−1264.7081×10−56.8228×10−51.0725×10−101.6889×10−1183.4825×10−54.7649×10−51.4925×10−101.8421×10−11Figure 1
(a)L∞ of Test Problem 1 obtained with the proposed method at different values of γ and its comparison with MCTSMQ (black dots). (b) Lrms of Test Problem 1 obtained with the proposed method at different values of γ and its comparison with MCTSMQ (black dots). Here we fix ξ=0.5 and obtain results at different scale levels.
(a)(b)Figure 2
(a) CPU time of TestProblem 1 consumed by proposed method at different values of γ and different scale levels. (b) Condition number of proposed method at different scale level and different values of γ, and its comparison with condition number obtained with MCTSMQ(black dots).
(a)(b)HWCM, MCTMQ, and MCTSMQ can only handle integer order Poisson equations. Besides high accuracy, one of the significant advantages of the proposed method is that it can also solve the fractional order Poisson equation (the case when1<σ≤2). Note that if σ≠2 then the exact solution of the first two problems is not known. Therefore, to check the accuracy of the approximate solution, we use two parameters Rrms and R∞. We approximate solution for some fractional values of σ; i.e., σ=1.5,0.1,2, for each value of σ we calculate the residual norms Rrms and R∞ at scale level ranges from 5 to 25. We observe that the solution is not too accurate for low values of σ. As the value of σ increases, the solution becomes more and more accurate. Also, as the value of M increases, the residual norms decreases, which guarantees the convergence of approximate solution for the fractional values of σ. The Rrms and R∞ for Test problem 1 are shown in Figure 3. From the simulation of fractional values of σ, we observe that the CPU time and condition number κ remain the same as for σ=2. These results, for some typical values of σ, are shown in Figure 4.Figure 3
(a) Residual normRrms of Test Problem 1 at different scale levels for different values of σ . (b) Residual norm R∞ of Test Problem 1 at different scale levels for different values of σ.
(a)(b)Figure 4
(a) CPU time of TestProblem 1 at different scale levels for different values of σ . (b) Condition number κ of Problem 5.1 at different scale levels for different values of σ.
(a)(b)TestProblem 2 is also analyzed, and the same conclusion is made. We fix ξ=0.5 and σ=2 and solve Test Problem 2 with the proposed method using different values of M. We observe that the proposed method yields a more accurate solution as compared to HWCM, MCTMQ, and MCTSMQ. The error norm L∞ of the proposed method obtained using M=9 is 2.26363×10−7, while L∞ for HWCM is 5.2168×10−4, and for MCTMQ and MCTSMQ, the error norm is found to be 5.8808×10−3 and 1.1688×10−6, respectively. The detailed results of the comparison of L∞ of the proposed method and other methods are displayed in Table 5. One can easily see that the proposed method yields better results than other methods.Table 5
L∞ of Test Problem 2 obtained with the proposed method and its comparison with HWCM, MCTMQ, and MCTSMQ [22].
HWCMMCTMQMCTSMQProposed methodM (n)L∞NL∞L∞M (n)L∞1128.4747×10−245.3642×10−35.3050×10−35252.7635×10−32322.9680×10−283.1607×10−51.0564×10−46363.2254×10−34988.0890×10−3161.1891×10−46.2451×10−67497.6056×10−583202.0735×10−3321.4152×10−23.5764×10−68641.1291×10−41611525.2168×10−4645.8808×10−31.1688×10−69812.26363×10−73243251.3062×10−4Not available101001.9237×10−7These results are obtained usingξ=0.5 and γ=1.Similarly, the error normLrms is also compared with HWCM, MCTMQ, and MCTSMQ. The proposed method yields Lrms=9.0136×10−8 at M=9, while the value of Lrms obtained with HWCM is 1.7150×10−3 and that for MCTMQ and MCTSMQ is 1.7118×10−3 and 5.2915×10−7, respectively. It is a clear indication of the superiority of the proposed method over these methods. The detailed results of the comparison of Lrms are displayed in Table 6.Table 6
Lrms of Test Problem 2 obtained with the proposed method and its comparison with HWCM, MCTMQ, and MCTSMQ [22].
HWCMMCTMQMCTSMQProposed methodM (n)LrmsNLrmsLrmsM (n)Lrms1129.5549×10−242.8660×10−32.8212×10−35251.0591×10−32323.7649×10−281.6169×10−55.5406×10−56361.1128×10−34981.3624×10−2166.2722×10−53.3011×10−67491.0155×10−583204.8439×10−3323.6982×10−31.6629×10−68643.0552×10−51611521.7150×10−3641.7118×10−35.2915×10−79819.0136×10−83243256.0655×10−4Not available101008.5282×10−8These results are obtained usingξ=0.5 and γ=1.We fixedσ=2 and M=10, and simulated the algorithm using different choices of parameter γ ranges from −8 to 8. For every value of γ, we record the value of L∞ and Lrms. It is found that maximum values of L∞ and Lrms are 4.1572×10−6 and 7.4056×10−7, which are obtained at γ=8. We compared L∞ and Lrms, for every value of γ, with HWCM, and it is observed that for every value of γ, the proposed method yields a correct solution. A detailed comparison is presented in Table 7. L∞ and Lrms for different values of γ of the proposed method are also compared with MCTSMQ. The results are displayed in Figure 5. The CPU time and condition number κ of the proposed method for this problem are shown in Figure 6. In the same figure, the condition number is compared with the condition number obtained using MCTSMQ, and it is shown that the proposed method is more stable for the current problem.Table 7
L∞ and Lrms of Test Problem 3 for different values of γ at M=10 and its comparison with HWCM [22] M=32.
HWCM [22]Proposed methodγL∞LrmsL∞Lrms−81.5165×10−43.5243×10−42.6752×10−65.5132×10−7−61.4002×10−43.4396×10−41.3673×10−62.7300×10−7−41.2134×10−43.3437×10−48.6599×10−72.1811×10−8−28.6389×10−53.3085×10−42.1177×10−76.0997×10−801.0031×10−44.0150×10−41.3658×10−75.7502×10−826.8998×10−41.8333×10−32.2636×10−79.0136×10−745.7799×10−41.1488×10−32.2760×10−65.1555×10−763.5700×10−46.7976×10−43.5376×10−67.4677×10−782.9952×10−45.6825×10−44.1572×10−67.4056×10−7Figure 5
(a)L∞ of Test Problem 2 obtained with the proposed method at different values of γ and its comparison with MCTSMQ (black dots). (b) Lrms of Test Problem 2 obtained with the proposed method at different values of γ and its comparison with MCTSMQ (blue dots). Here we fix ξ=0.5 and obtain results at different scale levels.
(a)(b)Figure 6
(a) CPU time of TestProblem 2 consumed by the proposed method at different values of γ and different scale levels. (b) Condition number of the proposed method at different scale level and different values of γ, and its comparison with condition number obtained with MCTSMQ (red dots).
(a)(b)TestProblem 3 is analyzed using the proposed method. We simulate the algorithm using different scale levels and record the value of L∞ and Lrms at fractional values of σ ranges from 1.5 to 1.9. The results are displayed in Figure 7. It can be easily seen that the error norm decreases with the increase of scale level M and the rate of convergence is approximately the same for all values of σ. It is also observed that the error norm at low values of σ is relatively high as compared to the high value of σ. The CPU time and condition number κ for the fractional values of σ are shown in Figure 8. It can be easily noted that the condition number is approximately equal to 1017. Also, an increase of CPU time with an increase of scale level is observed.Figure 7
(a)L∞ of Test Problem 3 obtained with the proposed method at different values of σ and using different scale levels. (b) Lrms of Test Problem 3 obtained with the proposed method at different values of σ and using different scale levels. Here, we fix ξ=0.5 and obtain results at different scale levels.
(a)(b)Figure 8
(a) CPU time of TestProblem 3 obtained with the proposed method at different values of σ and using different scale levels. (b) Condition number κ of Test Problem 3 obtained with the proposed method at different values of σ and using different scale levels. Here, we fix ξ=0.5 and obtain results at different scale levels.
(a)(b)The error norms are also calculated at different values of the parameterγ and the fractional value of σ. We observed that the error norm for this problem is low at a low absolute value of γ; as the absolute value of γ increases, the error norms increases. These results are displayed in Figure 9. However, the CPU time and condition number do not show any considerable change with changing values of γ. These results are displayed in Figure 10.Figure 9
(a)L∞ of Test Problem 3 obtained with the proposed method at different values of γ and the fractional values of σ using scale level M=10. (b) Lrms of Problem 5.3 obtained with the proposed method at different values of γ and the fractional values of σ using scale level M=10. Here, we fix ξ=0.5 and obtain results at different scale levels.
(a)(b)Figure 10
(a) CPU time of TestProblem 3 obtained with the proposed method at different values of γ and the fractional values of σ using scale level M=10. (b) Condition number κ of Test Problem 3 obtained with the proposed method at different values of γ and the fractional values of σ using scale level M=10. Here, we fix ξ=0.5 and obtain results at different scale levels.
(a)(b)


## 7. Conclusion and Future Work

The main advantage of the proposed method is its applicability to the fractional order Poisson equations. The method can easily handle fractional order problems with two-point boundary conditions. The method converts the heat flow phenomena to an algebraic structure, whose condition number is independent of the order of derivative. The proposed method yields a very accurate approximation when applied to fractional order Poisson equations. The comparison of results of the proposed method with some recent methods, such as, HWCM, MCTMQ, and MCTSMQ, shows that the proposed method is more appropriate for integer order problems. One of the significant advantages of this method is the computational cost. The computational time is compared to the other mentioned methods, and it is observed that the proposed method solves the problem in a very short time. By measuring the condition number of the algebraic system, the proposed method shows that the condition number of the structure is very small compared to the other mentioned methods. The proposed method also solves the fractional order partial differential equations with two-point nonlocal boundary conditions. The convergence of the proposed method is shown with test problems. One of the main targets of our plan is to study the convergence and stability of the proposed method. The extension of this method to other applied problems also lies in the domain of our future work.



---

*Source: 1067582-2021-04-26.xml*