Abstract

The paper is devoted to the study of operational matrix method for approximating solution for nonlinear coupled system fractional differential equations. The main aim of this paper is to approximate solution for the problem under two different types of boundary conditions, -point nonlocal boundary conditions and mixed derivative boundary conditions. We develop some new operational matrices. These matrices are used along with some previously derived results to convert the problem under consideration into a system of easily solvable matrix equations. The convergence of the developed scheme is studied analytically and is conformed by solving some test problems.

1. Introduction

Fractional calculus is generalization of integer order integration and differentiation to its noninteger (fractional) order counterpart. Fractional calculus has proved to be a valuable tool in modeling many natural phenomena of physics, chemistry, engineering, aerodynamics, electrodynamics of complex medium, polymer rheology, and so forth [1]. It is well known that fractional order operator is a nonlocal operator. Unlike integer order models the fractional order models have the property that the current state of the system depends on all its historic states. This makes fractional models interesting and remains as an active and hot area of research.

Differential equations are used to model different physical and engineering phenomena. The aim is to know about the future state of the system under consideration. It is some time necessary to model differential equations with different kinds of initial and boundary conditions (depends on the nature of the problem). Various types of conditions which are used as boundary conditions are multipoint local boundary conditions, integral type boundary conditions, multipoint nonlocal boundary conditions, and mixed derivative boundary conditions. The multipoint boundary value problems appear in wave propagation and in elastic stability. For example, the vibrations of a guy wire of a uniform cross section composed of sections of different densities can be molded as a multipoint boundary value problem (see [2, 3] and the references therein).

In this paper, we are interested in finding approximate solutions for the following generalized class of nonlinear coupled systems of fractional order differential equations:where and are the unknown solutions to be determined. These solutions are defined on finite time domain; that is, . The superscripts represent the order of derivatives defined in Caputo sense and are defined as , , and . and are nonlinear functional of and and their fractional derivatives. The method developed in this paper is designed for approximating and under any of the following two types of boundary conditions. (i)M-point nonlocal boundary conditions are as follows:where and are some real constants and are intermediate points of the domain defined as (ii)Mixed derivative boundary conditions are as follows:where , , , and are some real constants.

We develop new operational matrices based on Jacobi polynomials and use these operational matrices to develop new formulation for finding approximate solutions for the problems. The technique basically lies in the domain of spectral methods. Spectral methods are widely used for solutions to differential equations, functions approximations, and variational problems. These methods demand approximation of solutions for a problem by truncated series of smooth global functions and provide very accurate approximations for a smooth solution with relatively few degrees of freedom. The main reason behind this accuracy is the behavior of spectral coefficients , which tend to zero faster than any algebraic power of their index .

Among others, some of the well known mathematicians who successfully applied spectral method are Dehghan [49], Bhrawy [10, 11], Doha et al. [12], and Saadatmandi [7, 13, 14]. In these papers the authors solved many scientific problems using spectral methods. Also Dehghan and Shakeri [8] used a seminumerical technique for solving the multipoint boundary value problems. In [4, 6] an efficient way is developed for the construction of operational matrices. In [5] the authors used operational matrices of Bernstein polynomials to solve nonlinear age-structured model. We are motivated by the work of Bhrawy and Al-Shomrani [11], who solved fractional order differential equations (both linear and nonlinear) with -point boundary conditions (local) via shifted Legendre polynomials and collocation technique.

To solve the nonlinear coupled system of boundary value problem we implement operational matrix method combined with quasilinearization method. Quasilinearization method was first introduced by Bellman and Kalaba [15] to solve nonlinear ordinary or partial differential equations as a generalization of the Newton-Raphson method. The origin of this method lies in the theory of dynamic programming. In this method, the nonlinear equations are expressed as a sequence of linear equations and these equations are solved recursively. The main advantage of this method is that it converges monotonically and quadratically to the exact solution to the original equations [16]. Also some other interesting works in which qasilinearization method is used for scientific problems are available in [1719].

In our previous work, we have constructed some efficient methods for the numerical simulations of couple systems of linear fractional differential equations and fractional order partial differential equations [2022]. Local and nonlocal boundary value problems can be found in [2325]. For the readers who are new to the field, we recommend studying our previous work in order to get a better understanding.

2. Preliminaries

In this section, we recall some basic definitions and known results from fractional calculus. More details can be found in [1].

Definition 1. Given an interval , the Riemann-Liouville fractional order integral of order of a function is defined byprovided the integral on right hand side exists.

Definition 2. For a given function , the Caputo fractional order derivative of order is defined as provided the right side is pointwise defined on , where in case is not an integer and in case is an integer.

Hence, it follows that , , , for a constant , and

2.1. The Shifted Jacobi Polynomials

The well known two parametric Jacobi polynomials defined on , with parameters and , are given by the following relation:whereThe Jacobi polynomials are orthogonal and the orthogonality condition iswhere is the weight function and is known as normalization constant.

Any squared integrable function on can be approximated by shifted Jacobi polynomials as follows: , where the coefficient can easily be calculated using (10). In practice we are interested in the truncated series, which can also be written in vector form aswhere , is the coefficient vector, and is terms function vector.

One important property of shifted Jacobi polynomials which will be frequently used in our analysis is given aswhere and

Lemma 3. The definite integral of the product of weight function with three Jacobi polynomials over the domain is constant, defined as , defined bywhere

Proof. In view of (8), we obtainUsing convolution theorem of Laplace transform, we have Using the value of in (15), we obtain

3. Operational Matrices

In this section we first recall some previous results and then derive some new operational matrices. The following results are known [22].

Lemma 4. Let be the function vector as defined in (11); then the -order integration of is given by , where is the operational matrix of integration of order . The entries of are defined aswhere

Proof. The detailed  proof of this lemma is available in [22].

Lemma 5. Let be the function vector as defined in (11); then the -order derivative of is given by , where is the operational matrix of fractional differentiation of order . The entries of are defined aswhere

Proof. The detailed proof of this lemma is available in [22].

In [22], these operational matrices are used for solutions to coupled systems of initial value problems for fractional order differential equations. These matrices do not have the ability to handle -point nonlocal boundary value problems. Hence the need to develop new operational matrices which have the ability to handle boundary conditions as well as initial conditions has been of great importance. The new operational matrices include the operational matrices studied in [22] as a factor and have ability to solve fractional order boundary value problems.

Lemma 6. Let and be any functions defined on . Thenwhere is the Jacobi coefficients vector of defined by (11) andThe matrix is the operational matrix of derivative as defined in Lemma 5 andThe entries are defined by the relationwhere are the Jacobi coefficients of and is as defined in Lemma 3.

Proof. Using (11) and Lemma 5, we obtainwhere . With rearranging, the above equation can be written aswhereApproximating with Jacobi polynomials, we haveHence (28) takes the formwhereNow approximating with Jacobi polynomials, we havewhereUsing (31) in (33) we obtainwhich in view of Lemma 3 takes the formRepeating the procedure for and , we obtainUsing (36) in (27), we get

Lemma 7. For , where and are real constants, , and, for , the following holds:where is operational matrix and is defined byAnd the entries are defined by the relation

Proof. Using the definition of fractional integral (5), on we obtainwhich in view of the definition of Jacobi polynomials yieldsEvaluating the integral and after simplification, we obtainFor simplicity, use the following notation:From (44) and (45), it follows thatWe may also writewhere are given byAfter simplification, we haveNow using relation (45) we can writeUsing (50) in (46), we obtainwhich can be written in matrix form as where the entries of the matrix are as defined in (50).

4. Application of the Operational Matrices

In this section we apply the new operational matrices to solve coupled system of nonlinear fractional differential equations. The following lemmas are of basic importance in our further investigation.

Lemma 8. If , where , andwhere and are as defined in (2), and if , then can be written asand the matrix is defined aswhere , , and .

Proof. ConsiderApplication of fractional integration of order implies thatwhere and are constants of integration. Using the initial condition implies that . Application of nonlocal boundary condition implies that Equating the above two equations we get Using the value of in (57) and making use of Lemma 6 we get the following estimates:where , , and After simplification we get which completes the proof.

Lemma 9. If , where andthenwhere the matrix is defined aswhere , , , and Here one assumes that and are defined as , , and

Proof. Considerwhich implies thatFrom the above equation we can also writeNow using and we getwhich can be easily solved for and , and we can write where , , and .
Now using the value of and in (67) we get the following estimates:where , , , and .
In view of Lemma 6 we can write the above relation aswhich can be written aswhich completes proof of the lemma.

4.1. Couple System of FDEs with Variable Coefficients

Consider the following coupled system of FDEs:subject to one of the above discussed boundary conditions ((2) or (4)). Assume that and Then in view of Lemma 8 or Lemma 9 (depends on the given set of boundary conditions) we can writeNote that the subscript is used to distinguish between different types of boundary condition. If we are given nonlocal boundary conditions we will use and if we are given mix derivative boundary condition we will use Superscripts and are used to distinguish between conditions for and , respectively. Now in view of Lemma 6, we can write (74) aswhich can be converted into easily solvable matrix equation (by following the same step as in [22] from (51) to (69)). The generalized form of the resulting matrix equation is given as where , , and are of compatible size. The vector is the unknown solution to be determined. The solution method of such type of equations is explained in [26] and a detailed algorithm is presented. command uses the same algorithm presented in [26] for solution to such type of matrix equations.

4.2. Convergence Analysis

Consider and as entries of , and , respectively. Also assume that are entries of . Now we can writeUsing the above estimates and using the same procedure we can writewherewhere and are the approximation coefficients of and , respectively. In view of maximum norm and (12) we can writeIn view of Lemma of [22] the expansion coefficients and decay faster than any of the algebraic order of index . We can write , , , and approach to zero as and as a result we conclude that We see that the convergence of the scheme depends on the decay of expansion coefficients, and the decay of expansion coefficients depends solely on the smoothness of solution for the problem. Therefore if the solution for the problem is smooth we get high accuracy at small scale level, and the accuracy will increase if we increase the scale level.

4.3. Nonlinear FDEs with -Point Boundary Conditions

To solve nonlinear FDEs we will implement operational matrix method combined with qasilinearization method. Consider the following nonlinear FDE:

In qasilinearization method the nonlinear differential equation is converted into linear differential equation with variable coefficients. The nonlinear part of the differential equation is linearized about some function. In most cases the initial function is selected according to the physical problem under consideration, although in some cases it is selected as the solution for the linear part [27]. We will first solve the linear part of (83):under given boundary conditions using the method developed in the previous section. and are the linear part of and , respectively. The solution for this problem will be labeled as and . The next step is to linearize the nonlinear part of and about and using Taylor series expansion. As a result the nonlinear differential equation will be converted into linear differential equation with variable coefficients, which can be solved by the proposed. The whole process can be seen as a recurrence relation aswhere and The above equation is linear fractional differential equations with variable coefficients and can be easily solved with the method developed in the previous section.

5. Illustrative Examples

We show the applicability of the method by solving some test problems. Where available we compare our results with results obtained with other methods.

Example 1. Consider the following linear fractional order three-point boundary value problem [28]: where the forcing term is defined as The exact solution for this problem is This problem is solved in [28] using Haar wavelets (HW). In [29] this problem is solved using improved reproducing kernel method (RKM). We compare the absolute error obtained with the proposed method with the error obtained with HW and RKM. The results are shown in Table 1.

Example 2. Consider the following fractional differential equations having variable coefficients and five-point boundary conditions [29, 30]:where The exact solution for the problem is We approximate solution for this problem with the proposed method and compare its absolute error with the error reported in [29, 30] (in these references reproducing kernel method and improved reproducing kernel method are used to solve this problem). We observe that the solution obtained with the proposed method is in good agreement with the exact solution. In Table 2 absolute error at different value of is compared with the error reported in [29, 30].

Example 3. Consider the following nonlinear fractional differential equations with 6-point nonlocal boundary value problem:where the forcing term is defined as The exact and unique solution for the problem is . We approximate the solution for this problem with the proposed iterative method. This example is solved using the parameters , , and We observe that the method provides good approximation to solution for the problem. The approximate solution approaches the exact solution as we make iteration. Absolute difference of exact and approximate solution at different stages of iteration is shown in Table 3. It can be easily seen that at fourth iteration the absolute difference is less than .

Example 4. Consider the following nonlinear differential equation with 6-point nonlocal boundary conditions:The forcing term is defined as The exact solution for the problem is . We carry out the same analysis as in the previous example. The same conclusion is made; the solution converges rapidly to the exact solution for the problem. In Table 4 absolute error at different stages of iteration is given. For this example the absolute error is less than at fourth iteration. For this example we set , , and

Example 5. Consider the Lane-Emden equations [31]: Systems of Lane-Emden equations arise in the modeling of several physical phenomena, such as pattern formation, population evolution, and chemical reactions. The parameters , , , and can be considered from the actual chemical reaction or dynamics under consideration. We solve this problem with the proposed method subject to the following type of boundary conditions:In [31], Adomian decomposition method is employed to get approximate solution to this problem. We compare our results with the results reported in [31]. We observe that the current method provides a very good approximation to the solutions for the problem. We calculate the error remainder term: where and are the approximate solution for the problem at different stages of iteration . We calculate the maximal error remainder term:The comparison of results of the proposed method with Adomain decomposition is presented in Table 5.

6. Conclusion and Future Work

Form the above analysis and observation we conclude that the method works very well and efficiently solve fractional order linear and nonlinear differential equation and coupled system under multipoint nonlocal boundary conditions and mixed derivative boundary conditions. When the absolute error is compared with reproducing kernel method and Adomain decomposition method we observe that the method yields very accurate solution. In this work shifted Jacobi polynomials are used; it may also be possible that the reader may get more accurate solution by using another class of orthogonal polynomials like Bernstein polynomials or Hermite polynomials and so forth. Our future work is related to investigating the best choice of orthogonal polynomials.

Competing Interests

The authors declare that they have no competing interests.