#### 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*.

#### 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:

The above model is subject to two-point nonlocal boundary conditions:

The functions and are given smooth functions. The parameters and are two positive constants. The parameter 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:

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 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.

#### 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 , the Riemann–Liouville fractional order integral of a function of order is defined by the following:provided that the integral on the right-hand side exists.

*Definition 2. *(Caputo derivative). For a given function , the Caputo fractional order derivative is defined as follows:where .

Hence, it follows that

##### 2.1. The Shifted Legendre Polynomials

The shifted Legendre polynomials [35] defined on are given by the following relation:

These polynomials are bounded by 1, and we have the following relation:

These polynomials are orthogonal on the domain , and the orthogonality condition is given as follows:

which implies that any can be approximated by Legendre polynomials as follows:

In vector notation, we write the following:where is the scale level of approximation. is the coefficient vector and is terms function vector. These notations can be easily extended to two-dimensional space [35] and two-dimensional Legendre polynomials of the order are defined as a product function of two Legendre polynomials

The orthogonality condition of is as follows:

Any can be approximated by the polynomials as follows:

For simplicity, use the notation where and rewrite (14) as follows:where is coefficient row vector and is column vector of functions defined by the following:where .

#### 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 or 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 be the function vector as defined in (16), then the fractional integral of order of w.r.t is given by the following:where is the operational matrix of the fractional integration of order and is defined as follows:where , , for and*

*Proof. *Taking the element defined by (12), then the fractional-order integration of w.r.t follows:Using the definition of fractional integration, we obtain the following:Approximating by terms of Legendre polynomials in two variables yieldswhere , which in view of the orthogonality conditions implies thatwhereand hence, it follows thatwhereUsing the notations, , and for , we get the desired result.

Lemma 2. *Let be as defined in (16), then the integration of order of w.r.t is given by the following:where is the operational matrix of derivative of order and is defined as follows:and , , for and*

*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 like and , where and are some positive constants and . 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. *Let , then for some constants and , the following relation holds:where, where , and*

*Proof. *LetFor some constant and , we can write the following expression:Now, considering the general term of , we can write the following:where we use the notation . By using the definition of Legendre polynomials, we can write the following:Approximating by Legendre polynomials in two variables yieldswhere , which in view of the orthogonality conditions implies thatwhereHence, it follows thatwhereUsing the notations, , and for , we get the desired result.

Lemma 4. *Let , then for some constants and , the following relation holds:where, where , and*

*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: 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 and , we get the following relation:

Using the values of and in (48), we get the following:

which in view of Lemma 3 can be written as follows:where .

Now approximating the source term , 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 we can easily ; however, the second constant is not known. We can use the two-point boundary conditions:

Using the equality , we get the following:

From which we can calculate the value of as follows:where . Using the value of and in (53), we can write the following:which in view of Lemma 4, can be written as follows:where . 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 , which can be used in (51) or (59) to get an approximate solution to the problem.

##### 4.1. Error Bound

Considering a sufficiently smooth function on , let is the space span by term Legendre polynomials. We assume that is its best approximation in . For this purpose, consider a polynomial is any polynomial of degree in variable and , respectively. Then, from the definition of best approximation,

The inequality in (62) also holds if is interpolating polynomial at point ; then by the similar arguments as in [36], the error of the approximation is given by the following:where

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 .

#### 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:The exact solution of the problem for fix is .

*Test Problem 2. *(see [22, 23]). Consider (45) with the following functions:The exact solution of the problem for the fix is .

*Test Problem 3. *(see [22, 23]). Consider (45) with the following functions.The exact solution of the problem for is .

#### 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:

If , 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:where is defined as follows:

The first problem is solved using Haar wavelets (HWCM), mesh-less method without splitting (MCTMQ), and with splitting (MCTSMQ). We fixed and obtained an approximate solution of Test Problem 1 for different values of scale level . We compared our results with HWCM, MCTMQ, and MCTSMQ. The comparison of of Test Problem 1 obtained with the proposed method and HWCM, MCTMQ, and MCTSMQ are shown in Table 1. We can easily note that obtained with HWCM at is ; note that at this level, HWCM converts the problem to a system of algebraic equations of 1152 unknowns. While the proposed method yields , while converting the problem to a system of algebraic equations of 81 unknowns. of this problem obtained with MCTMQ and MCTSMQ at is and , respectively. This shows the superiority of the proposed method over HWCM and meshless methods.

The parameters 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 for this problem obtained with the proposed method at is while HWCM yields at , and meshless method MCTMQ and MCTSMQ yield and .

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 of the proposed method with MCTMQ and MCTSMQ. It is observed that the proposed method is more robust than these methods. At , 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.

The accuracy of the present method is analyzed at different values of . We chose and calculate and for scale levels. We observed that the accuracy of the proposed method does not depend on the values of . The values of and obtained with the proposed method at different values of using scale level are compared with the and obtained with the HWCM at 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 and , the error norms are much less than the error norms obtained using MCTSMQ at . The at different values of and 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 , while the condition number of the proposed method is approximately equal to , 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.

**(a)**

**(b)**

**(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 when ). Note that if 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 and . We approximate solution for some fractional values of ; i.e., , for each value of we calculate the residual norms and 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 increases, the residual norms decreases, which guarantees the convergence of approximate solution for the fractional values of . The and for Test problem 1 are shown in Figure 3. From the simulation of fractional values of , we observe that the and condition number remain the same as for . These results, for some typical values of , are shown in Figure 4.

**(a)**

**(b)**

**(a)**

**(b)**

Test Problem 2 is also analyzed, and the same conclusion is made. We fix and and solve Test Problem 2 with the proposed method using different values of . We observe that the proposed method yields a more accurate solution as compared to HWCM, MCTMQ, and MCTSMQ. The error norm of the proposed method obtained using is , while for HWCM is , and for MCTMQ and MCTSMQ, the error norm is found to be and , respectively. The detailed results of the comparison of 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 is also compared with HWCM, MCTMQ, and MCTSMQ. The proposed method yields at , while the value of obtained with HWCM is and that for MCTMQ and MCTSMQ is and , respectively. It is a clear indication of the superiority of the proposed method over these methods. The detailed results of the comparison of are displayed in Table 6.

We fixed and , and simulated the algorithm using different choices of parameter ranges from to 8. For every value of , we record the value of and . It is found that maximum values of and are and , which are obtained at . We compared and , 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. and 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.

**(a)**

**(b)**

**(a)**

**(b)**

Test Problem 3 is analyzed using the proposed method. We simulate the algorithm using different scale levels and record the value of and 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 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 . Also, an increase of CPU time with an increase of scale level is observed.

**(a)**

**(b)**

**(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.

**(a)**

**(b)**

**(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.

#### 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.

#### Acknowledgments

This study was supported by Universiti Kebangsaan Malaysia (Grant no. GP-2020-K006388).