Toward the Approximate Solution for Fractional Order Nonlinear Mixed Derivative and Nonlocal Boundary Value Problems

1Department of Mathematics, The University of Poonch Rawalakot, Azad Jammu and Kashmir, Pakistan 2Department of Applied Science, Ajloun College, Al-Balqa Applied University, Ajloun 26816, Jordan 3Department of Mathematics, Faculty of Science and Arts, Shaqra University, Shaqra 11691, Saudi Arabia 4Department of Mathematics, University of Malakand, Chakdara, Lower Dir, Khyber Pakhtunkhwa, Pakistan 5School of Mathematical Sciences, Universiti Kebangsaan Malaysia, 43600 Bangi, Selangor, Malaysia 6Research Institute, Center for Modeling and Computer Simulation, King Fahd University of Petroleum and Minerals, Dharan 31261, Saudi Arabia


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 V = V() are the unknown solutions to be determined.These solutions are defined on finite time domain; that is,  ∈ [0, ].The superscripts represent the order of derivatives defined in Caputo sense and are defined as 1 <  1 ,  2 ≤ 2, and 0 <  1 ≤ 1.  and  are nonlinear functional of  and V and their fractional derivatives.The method developed in this paper is designed for approximating  and V under any of the following two types of boundary conditions.
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 [4][5][6][7][8][9], 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 [17][18][19].
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 [20][21][22].Local and nonlocal boundary value problems can be found in [23][24][25].For the readers who are new to the field, we recommend studying our previous work in order to get a better understanding.

Preliminaries
In this section, we recall some basic definitions and known results from fractional calculus.More details can be found in [1].
provided 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  = [] + 1 in case  is not an integer and  =  in case  is an integer.
Hence, it follows that     = (Γ(1 + )/Γ(1 +  + )) + for  > 0,  ≥ 0,    = 0, for a constant , and 2.1.The Shifted Jacobi Polynomials.The well known two parametric Jacobi polynomials defined on [0, ], with parameters  and , are given by the following relation: where where  (,)  () = (−)    is the weight function and Any squared integrable function V() on [0, ] 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 as where  =  + 1, H   is the coefficient vector, and Ψ  () is  terms function vector.
Using convolution theorem of Laplace transform, we have Using the value of Υ (,,) in (15), we obtain

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 × is the operational matrix of integration of order .The entries Θ ,, of H , × are defined as where ( 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 × is the operational matrix of fractional differentiation of order .The entries Φ ,, of G , × are defined as where 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 [0, ].Then where W   is the Jacobi coefficients vector of () defined by ( 11) and The matrix G , × is the operational matrix of derivative as defined in Lemma 5 and The entries are defined by the relation where   are the Jacobi coefficients of   () and ϝ (,,) (,,) is as defined in Lemma 3.
Proof.Using (11) and Lemma 5, we obtain where () =    Ψ  ().With rearranging, the above equation can be written as where Approximating   () with Jacobi polynomials, we have Hence (28) takes the form where Now approximating ℵ  () with Jacobi polynomials, we have where Using ( 31) in (33) we obtain which in view of Lemma 3 takes the form Repeating the procedure for  = 0, 1, . . .,  and  = 0, 1, . . ., , we obtain Using ( 36) in ( 27), we get Lemma 7.For    =   , where  and  are real constants, () = W   Ψ  (), and, for 0 <  ≤ , the following holds: where × is  ×  operational matrix and is defined by = 0, 1, . . ., ,  = 0, 1, . . ., . (39) And the entries  (,) are defined by the relation Proof.Using the definition of fractional integral (5), on we obtain which in view of the definition of Jacobi polynomials yields Evaluating the integral and after simplification, we obtain For simplicity, use the following notation: From ( 44) and (45), it follows that We may also write ⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞ ℧ (,,,) (,) where  (,) are given by After simplification, we have (49) Now using relation (45) we can write Using (50) in (46), we obtain which can be written in matrix form as where the entries of the matrix  (,   ,) × are as defined in (50).

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   () = W   Ψ  (), where 1 <  ≤ 2, and where   and   are as defined in (2), and if ∑ −2 =0     −  ̸ = 0, then () can be written as and the matrix V 1 is defined as where Application of fractional integration of order  implies that where  0 and  1 are constants of integration.Using the initial condition (0) =  0 implies that  0 =  0 .Application of nonlocal boundary condition implies that Equating the above two equations we get Using the value of  1 in (57) and making use of Lemma 6 we get the following estimates: where  = 1/(∑ −2 =0     − ),   =   /(∑ −2 =0     − ), and After simplification we get which completes the proof.
, where 1 <  ≤ 2 and then where the matrix V 2 is defined as where Here one assumes that  2 ,  3 ̸ = 0 and are defined as which implies that From the above equation we can also write Now using   (0) =  1 (0)+ 2 () and   () =  1 (0)+ 2 () we get which can be easily solved for  0 and  1 , and we can write where ), and  3 = ( 1 +  2 ).Now using the value of  0 and  1 in (67) we get the following estimates: where In view of Lemma 6 we can write the above relation as which can be written as which completes proof of the lemma.

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  ( 1 ) () = X   Ψ  () and V ( 1 ) () = Y   Ψ  ().Then in view of Lemma 8 or Lemma 9 (depends on the given set of boundary conditions) we can write Note that the subscript  is used to distinguish between different types of boundary condition.If we are given nonlocal boundary conditions we will use  = 1 and if we are given mix derivative boundary condition we will use  = 2. Superscripts  and V are used to distinguish between conditions for () and V(), respectively.Now in view of Lemma 6, we can write (74) as which 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.

Convergence Analysis.
Consider Ω  , and Ω V , as entries of V   , and V V  , respectively.Also assume that Λ 24 .Now we can write Using the above estimates and using the same procedure we can write where where   and   are the approximation coefficients of () and (), respectively.In view of maximum norm and ( 12 In view of Lemma 2.1 of [22] the expansion coefficients   and V  decay faster than any of the algebraic order of index .We can write   , V  ,   , and   approach to zero as  → ∞ and as a result we conclude that       ()     → 0 as  → ∞,      V ()     → 0 as  → ∞.
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.

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 Table 1: Comparison of absolute error of Example 1 obtained with the proposed method and its comparison with Haar wavelets [28] and RKM [29]. = 8 (HW) [28]  = 8 (RKM) [ 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. 1 and  2 are the linear part of  and , respectively.The solution for this problem will be labeled as  0 () and V 0 ().The next step is to linearize the nonlinear part of  and  about  0 () and V 0 () 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 as where   = (,   1  , V  1  ,   , V  ) and   = (,   1  , V  1  ,   , V  ).The above equation is linear fractional differential equations with variable coefficients and can be easily solved with the method developed in the previous section.

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 () =  5 − 37 3 /20 + 33 2 /40.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 () = 2 2 + 2 cos() + Γ(3)/Γ(1.7) 0.7 .The exact solution for the problem is () =  2 .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 where the forcing term () is defined as The exact and unique solution for the problem is () =  5 − 2 + 1.We approximate the solution for this problem with the proposed iterative method.This example is solved using the parameters  = 0,  = 2, and  = 8.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 =  (1) .
The forcing term () is defined as (92) The exact solution for the problem is () = cos(/3).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 10 −16 at fourth iteration.For this example we set  = 1,  = 1, and  = 10.(94) 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: ER   = The comparison of results of the proposed method with Adomain decomposition is presented in Table 5.

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.

Table 2 :
[29,30]son of absolute error of Example 2 obtained with the proposed method and its comparison with RKM[30]and improved reproducing kernel method[29].theproposedmethod is in good agreement with the exact solution.In Table2absolute error at different value of  is compared with the error reported in[29,30].

Table 3 :
Absolute difference of exact and approximate solution of Example 3 at different stages of the iteration.

Table 4 :
Absolute difference of exact and approximate solution of Example 4 at different stages of the iteration.

Table 5 :
[31]mal error of Example 5 obtained with the proposed method and its comparison with error reported in[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  11 ,  12 ,  21 , and  22 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: ()  () +  1 ()    () −  11  2  () −  12   () V  () , ER V  = V ()  () +  2 () V   () −  21  2  () −  22   () V  () ,(95)where   and V  are the approximate solution for the problem at different stages of iteration .We calculate the maximal error remainder term: ∈[0,1]