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

#### 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 to , as given by corresponding to the conditionswhere represents partial derivative, , the real constants are denoted by , and It is remarkable that the proposed system contains the coupled system investigated in [38] as a special case. If we take together with and , the proposed system becomes a coupled system of Laplace equations of fractional order by taking Further, if we involve the external (source) terms and consider , 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: Section 2 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 The extension of (4) to two parameters is given as

*Definition 2 (see [45, 46]). *The fractional integral of Riemann-Liouville type of order of a function is defined as provided that the right-hand side is pointwise defined on .

*Definition 3. *For a given function , the Caputo fractional derivative of order is defined as provided that the right side is pointwise defined on , where .

Hence, it follows that

Theorem 4 (see [45]). *The solution of FDEis given bywhere each is a real number for , where , such that represents integer part of .*

*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 are representing as The transformation transforms the interval to and the corresponding polynomials known as shifted Legendre polynomials are given as The orthogonality condition is Therefore, any function can be approximated in terms of shifted Legendre polynomials as In matrix form, we may write (14) as where is the matrix of functions while is term column matrix. The shifted Legendre polynomials of order for two dimensions are recalled by and the corresponding orthognathic relation is provided as Therefore any function can be approximated in terms of as where In matrix notation, we may write (18) as where the column vector of coefficients with order is represented as Also column vector of order is represented by and provided as where .

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

Consider a sufficiently smooth function on ; then the error in the numerical solutions is provided as whereThe above relation also holds at points for interpolating polynomial

#### 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 order of defined in (21) with respect to is provided as where is the operational matrix of fractional integration given as for ,*

*Proof. *Consider as given in (16); then the arbitrary order integral with order of corresponding to is provided as which gives thatFurther approximating in terms of shifted Legendre polynomials, we have where Thanks to orthogonality condition, we have Therefore (28) yields whereWe use these notations for easiness of Thus, we received the required result

Theorem 7. *The fractional order derivative of can be written as such that is the operational matrix corresponding to arbitrary order derivative of order with respect to provided as where for , *

*Proof. *Take the order fractional derivative of corresponding to asApproximate in terms of shifted Legendre polynomials asApplying orthogonality condition, we may write Therefore, we getwhereDenoting , such that , for easiness, we get the required.

Theorem 8. *The fractional derivative with order of with respect to can be written as such that is the operational matrix given as for and*

*Proof. *To prove this theorem it can be easily followed from the proof of Theorem 7.

Theorem 9. *The fractional derivative of with respect to is provided as where represents the operational matrix of fractional derivative given as and for ,*

*Proof. *Take given in (17) and apply fractional order derivative of corresponding to and as We approximate in terms of shifted Legendre polynomials as such that Thanks to the orthogonality conditions and convolution theorem of Laplace transform, evaluating the integrals (49), we get where Plugging (48) in (47), we may writeRepresenting , and for , 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 In view of Theorem 4, we can write Applying initial conditions, we get where and In simple notation, we may write Therefore (55) yieldsThanks to (57), the other terms of system (1) can be approximated as Putting in (1), we get Now we write matrix form of system (59) aswhich further gives where is a column null matrix and other null matrices are and .

Simplifying (61), we get Therefore by using and , the aforesaid equation can be written aswhere Hence (64) is Sylvester type matrix equation, which on solving for unknown matrix 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 assuch that the external functions and are given asThe exact solution of the above system is 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 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 Table 2, we compute the CPU time for the computational of solutions for the test Problem 1.

*Problem 2. *To support the aforesaid established results, we consider the following problem:The external source functions are provided as where the exact solution of the Problem 2 is