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

In this paper, we have extended the operational matrix method for approximating the solution of the fractional-order twodimensional 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. (ese matrices are used to convert a sample nonlocal heat conduction phenomenon of fractional order to a structure of easily solvable algebraic equations. (e solution of the algebraic structure is then used to approximate a solution of the heat conduction phenomena. (e proposed method is applied to some test problems.(e obtained results are compared with the available data in the literature and are found in good agreement.


Introduction
e 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. e work in [3] deals with the exothermic reactions model having a constant heat source in the porous media with strong memory effects. is 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][6][7]. e term nonlocal problems means that the solution of PDEs on the boundary is connected with the solution on some interior points of the domain. e 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][9][10][11][12][13][14][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][17][18]. Two-dimensional diffusion problems with nonlocal boundary conditions have been discussed in [13,19]. e 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) e above model is subject to two-point nonlocal boundary conditions: u(1, y) � cu(ξ, y) + μ 2 (y), 0 ≤ y ≤ 1, 0 ≤ ξ < 1.
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: e 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. ey implemented two different methods for the solution of integer order counterpart of (1). eir 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. e results presented in [22,23] are the motivating factor of our interest to study the approximate solution of the fractionalorder 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. e new operational matrices can handle the nonlocal boundary conditions. e interesting readers may find useful results and some new strategies of this method in [24][25][26][27].
Application of orthogonal polynomials combined with Tau and Collection method can be found in [28][29][30][31]. e rest of the paper is organized as follows: in Section 2, 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.

Preliminaries
In this section, we present some useful results and notations which are of primary importance in our further investigation.
Definition 1 (see [32][33][34]). Given an interval [a, b] ⊂ R, the Riemann-Liouville fractional order integral of a function ϕ ∈ (L 1 [a, b], R) of order α ∈ R + is defined by the following: provided that the integral on the right-hand side exists.
Hence, it follows that 2 Mathematical Problems in Engineering ese polynomials are bounded by 1, and we have the following relation: ese polynomials are orthogonal on the domain [0, 1], and the orthogonality condition is given as follows: which implies that any f(x) ∈ C[0, 1] can be approximated by Legendre polynomials as follows: (10) In vector notation, we write the following: where M � m + 1 is the scale level of approximation. K is the coefficient vector and P M (x) is M terms function vector. ese 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 P n (x, y) � P a (x)P b (y), n � Ma + b + 1, a � 0, 1, 2, . . . , m, b � 0, 1, 2, . . . , m. (12) e orthogonality condition of P n (x, y) is as follows: can be approximated by the polynomials P n (x, y) as follows: where C ab � (2a + 1)(2b + 1) For simplicity, use the notation C n � C ab where n � Ma + b + 1 and rewrite (14) as follows: where K M 2 is 1 × M 2 coefficient row vector and is M 2 × 1 column vector of functions defined by the following:Ψ(x, y) where ψ i+1,j+1 (x, y) � P i (x)P j (y).

New Operational Matrices
e operational matrices of the fractional derivatives and integrals play a vital role in converting the FPDEs to the system of algebraic equations. e operational matrices of all derivatives are explicitly derived in our previous report [35]. We will need operational matrices in integration. e 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: where G σ,y M 2 ×M 2 is the operational matrix of the fractional integration of order σ and is defined as follows: where Proof. Taking the element P n (x, y) defined by (12), then the fractional-order integration of P n (x, y) w.r.t y follows: Mathematical Problems in Engineering Using the definition of fractional integration, we obtain the following: Approximating P a (x)y k+σ by M terms of Legendre polynomials in two variables yields where C ij � (2i + 1)(2j + 1) 1 0 1 0 P a (x)y k+σ P i (x)P j (y) dxdy, which in view of the orthogonality conditions implies that where and hence, it follows that where Using the notations, q � Mi + j + 1, r � Ma + b + 1 and Δ q,r ′ � C i,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: where G σ,x M 2 ×M 2 is the operational matrix of derivative of order σ and is defined as follows: and Proof.
e proof of this lemma is similar to the above lemma.
□ e 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 like cx n ( where c and n are some positive constants and 0 < ξ ≤ 1. erefore to replace such terms with their equivalent matrix form, we need to derive two more operational matrices. e operational matrices used to replace such term by their equivalent matrix form are derived in the following lemmas.
, then for some constants c and n, the following relation holds: where P (c,n,σ,ξ,y) Proof. Let For some constant c and n, we can write the following expression: Now, considering the general term of Ψ(x, s), we can write the following: where we use the notation r � Mi + j + 1. By using the definition of Legendre polynomials, we can write the following: Approximating P i (x)y n by Legendre polynomials in two variables yields where where Hence, it follows that where Using the notations, r � Mi + j + 1, q � Ma + b + 1 and Δ r,q � d c,n,σ,ξ a,b,i,j for i, j, a, b � 0, 1, 2, 3, . . . , m, we get the desired result. □ Lemma 4. Let u(x, y) � K M 2 Ψ(x, y), then for some constants c and n, the following relation holds:

Mathematical Problems in Engineering
where Proof. is lemma can be easily proved by following similar steps as in the previous lemma. is is left as an exercise for interested readers.

Main Result: Application of
Operational Matrices e 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: subject to nonlocal two-point boundary conditions 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: On application of the fractional integral of order σ and making use of Lemma 2, we get the following relation: Using the conditions at y � 0 and y � 1, we get the following relation: Using the values of c 0 and c 1 in (48), we get the following: which in view of Lemma 3 can be written as follows: where , and using (47) in (45), we can write the following: On application of fractional integral of order σ w.r.t. x, we can write the following: Using the initial conditions at x � 0 we can easily d 0 � μ 1 (y); however, the second constant is not known. We can use the two-point boundary conditions: Using the equality u(1, y) � cu(ξ, y) + μ 2 (y), we get the following: From which we can calculate the value of d 1 as follows: Mathematical Problems in Engineering where r 1 � (1/1 − cξ) ≠ 0. Using the value of d 0 and d 1 in (53), we can write the following: which in view of Lemma 4, can be written as follows: where F 2M 2 Ψ(x, y) � μ 1 (y) + (c − 1)r 1 μ 1 (y) + r 1 μ 2 (y))x.
In simplified notation, we can write the following: On comparing equation (51) and (59), we can write the following: Canceling out the common term, we can write the following: which is a linear system of equations and can be easily solved for the unknown vector K M 2 , which can be used in (51) e inequality in (62) also holds if P (M,M) (x, y) is interpolating polynomial at point (x i , y j ); then by the similar arguments as in [36], the error of the approximation is given by the following: where Mathematical Problems in Engineering zy M+1 g(x, y) , zy M+1 g(x, y) .

(64)
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.

Test Problems
We solve the fractional order generalization of some bench mark problems from [22,23].

Results and Discussion
We solve the above problems with the proposed method.
e 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: 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 σ. ese residual norms are defined as follows: where r(x, y) is defined as follows: 8 Mathematical Problems in Engineering e 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. e comparison of L rms of Test Problem 1 obtained with the proposed method and HWCM, MCTMQ, and MCTSMQ are shown in Table 1. We can easily note that L rms 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 L rms � 1.3381 × 10 − 10 , while converting the problem to a system of algebraic equations of 81 unknowns. L rms of this problem obtained with MCTMQ and MCTSMQ at N � 64 is 4.1256 × 10 − 4 and 3.2737 × 10 − 7 , respectively. is shows the superiority of the proposed method over HWCM and meshless methods. e parameters L ∞ for Test Problem 1 obtained with the proposed method are also compared with HWCM, MCTMQ, and MCTSMQ. e 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 . e proposed method along with HWCM and the meshless method convert the problem to the system of linear algebraic equations. e 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. e 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. e comparison of CPU time and conditions number of the proposed method with MCTQM and MCTSMQ at different scale levels is shown in Table 3.   [22], MCTMQ [22], and MCTSMQ [22].
Mathematical Problems in Engineering 9 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 10 19 , while the condition number of the proposed method is approximately equal to 10 17 , 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. HWCM, MCTMQ, and MCTSMQ can only handle integer order Poisson equations. Besides high accuracy, one Table 3: CPU time and condition number κ of the proposed method and its comparison with MCTMQ [22] and MCTSMQ [22].
MCTMQ [22] MCTSMQ [22]    of the significant advantages of the proposed method is that it can also solve the fractional order Poisson equation (the case when 1 < σ ≤ 2). Note that if σ ≠ 2 then the exact solution of the first two problems is not known. erefore, to check the accuracy of the approximate solution, we use two parameters R rms 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 R rms and R ∞ at scale level  Table 5: L ∞ of Test Problem 2 obtained with the proposed method and its comparison with HWCM, MCTMQ, and MCTSMQ [22].    Figure 3. From the simulation of fractional values of σ, we observe that the CPU time and condition number κ remain the same as for σ � 2. ese results, for some typical values of σ, are shown in Figure 4.  Test Problem 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. e 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. e 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. Similarly, the error norm L rms is also compared with HWCM, MCTMQ, and MCTSMQ. e proposed method yields L rms � 9.0136 × 10 − 8 at M � 9, while the value of L rms 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. e detailed results of the comparison of L rms are displayed in Table 6.
We fixed σ � 2 and M � 10, and simulated the algorithm using different choices of parameter c ranges from −8 to 8. For every value of c, we record the value of L ∞ and L rms . It is found that maximum values of L ∞ and L rms are 4.1572 × 10 − 6 and 7.4056 × 10 − 7 , which are obtained at c � 8. We compared L ∞ and L rms , for every value of c, with HWCM, and it is observed that for every value of c, the proposed method yields a correct solution. A detailed comparison is presented in Table 7. L ∞ and L rms for different values of c of the proposed method are also compared with MCTSMQ. e results are displayed in Figure 5. e 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.
Test Problem 3 is analyzed using the proposed method. We simulate the algorithm using different scale levels and record the value of L ∞ and L rms at fractional values of σ ranges from 1.5 to 1.9. e 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 σ. e 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 10 17 . Also, an increase of CPU time with an increase of scale level is observed. e error norms are also calculated at different values of the parameter c and the fractional value of σ. We observed that the error norm for this problem is low at a low absolute value of c; as the absolute value of c increases, the error norms increases.
ese results are displayed in Figure 9. However, the CPU time and condition number do not show any considerable change with changing values of c. ese results are displayed in Figure 10.

Conclusion and Future Work
e main advantage of the proposed method is its applicability to the fractional order Poisson equations. e method can easily handle fractional order problems with two-point boundary conditions. e method converts the heat flow phenomena to an algebraic structure, whose condition number is independent of the order of derivative. e proposed method yields a very accurate approximation when applied to fractional order Poisson equations. e 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. e 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. e proposed method also solves the fractional order partial differential equations with two-point nonlocal boundary conditions. e 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. e extension of this method to other applied problems also lies in the domain of our future work.

Data Availability
All the data are available in the manuscript.

Conflicts of Interest
All the authors declare that there are no conflicts of interest.

Authors' Contributions
All the authors have equally contributed to the preparation of the manuscript. All authors read and approved the final manuscript.