Multiple Solutions of Nonlinear Boundary Value Problems of Fractional Order: a New Analytic Iterative Technique

The purpose of this paper is to present a new kind of analytical method, the so-called residual power series, to predict and represent the multiplicity of solutions to nonlinear boundary value problems of fractional order. The present method is capable of calculating all branches of solutions simultaneously, even if these multiple solutions are very close and thus rather difficult to distinguish even by numerical techniques. To verify the computational efficiency of the designed proposed technique, two nonlinear models are performed, one of them arises in mixed convection flows and the other one arises in heat transfer, which both admits multiple solutions. Graphical results and tabulate data are presented and discussed quantitatively to illustrate the multiple solutions. The results reveal that the method is very effective, straightforward, and powerful for formulating these multiple solutions.


Introduction
Multiple or dual solutions of nonlinear boundary value problems (BVPs) of fractional order are an interesting subject in the area of mathematics, physics, and engineering.In fact, it is more consequential not to lose any solution of nonlinear BVPs of fractional order due to their wide application in scientific and engineering research.Based on this important fact, the present paper is going to present an analytical method, the so-called residual power series (RPS), that enables us to predict the multiplicity of solutions which nonlinear BVP of fractional order admits and furthermore to calculate the multiple solutions analytically at the same time.
BVPs of fractional order have received considerable attention in the recent years due to their wide applications in the areas of physics and engineering.Many important phenomena in electromagnetics, acoustics, viscoelasticity, electrochemistry, and material science are well described by fractional BVP [1][2][3][4].It is well known that the fractional order differential and integral operators are non-local operators.This is main reason why differential operators of fractional order provide an excellent instrument for description of memory and hereditary properties of various physical and engineering processes.For example, half-order derivatives and integrals proved to be more useful for the formulation of certain electrochemical problems than the classical models [5][6][7][8][9].Indeed, for example, applying fractional calculus theory to entropy theory has become a significant part and a hotspot research domain [10][11][12][13][14][15][16][17][18][19], since the fractional entropy could be used in the formulation of algorithms for image segmentation where traditional Shannon entropy has presented limitations [13] and in the analysis of anomalous diffusion processes and fractional diffusion equations [14][15][16][17][18][19].
In general, most BVPs of fractional order do not have exact solutions.Particularly, there is no known method for solving these types of equations in closed form solution.As a result, numerical and analytical techniques have been used to study such problems.The reader is referred to [20][21][22][23][24][25][26][27] in order to know more details about the fractional BVPs, including their history and kinds, their existence and uniqueness of solution, their applications and methods of solutions, etc.
Series expansions are a very important aid in numerical calculations, especially for quick estimates made in hand calculations, for example, in evaluating functions, integrals, or derivatives.Solutions to differential equations can often be expressed in terms of series expansions.Since, the advent of computers, it has, however, become more common to treat differential equations directly, using different approximation method instead of series expansions, but in connection with the development of automatic methods for formula manipulation, one can anticipate renewed interest in series methods.These methods have some advantages, especially in multidimensional and multiple solutions for BVPs of fractional order.
The RPS method was developed by the first author [28] as an efficient method for determining values of coefficients of the power series solution for first and the second-order fuzzy differential equations.It has been successfully applied in the numerical solution of the generalized Lane-Emden equation, which is a highly nonlinear singular differential equation [29] and in the numerical solution of higher-order regular differential equations [30].The RPS method is an effective and easy to construct power series solution for strongly linear and nonlinear equations without linearization, perturbation, or discretization [28][29][30].Different from the classical power series method, the RPS method does not need to compare the coefficients of the corresponding terms and a recursion relation is not required.This method computes the coefficients of the power series by a chain of equations of one or more variables.In fact, the RPS method is an alternative procedure for obtaining analytic solutions for BVPs of fractional order that admit multiple solutions.By using the residual error concept, we get a series solution, in practice a truncated series solution.Moreover, the multiple solutions and all their fractional derivatives are applicable for each arbitrary point in the given interval.On the other aspect as well, the RPS method does not require any conversion while switching from the low-order to the higher-order; as a result the method can be applied directly to a given problem by choosing an appropriate initial guess approximation.
In the present paper, the RPS method will investigate how to construct new algorithms for predicting and finding multiple solutions for those nonlinear BVPs of fractional order that admit multiple solutions.Furthermore, we will adapt a new generalization of Taylor's series formula that involves Caputo fractional derivatives in order to apply the RPS method.
The results dealing with multiple solutions of BVPs of fractional order are relatively scarce.Recently, many authors have discussed the multiple solutions to some problems using some of the well-known methods.However, the reader is referred to [31][32][33][34][35] in order to know more details about these methods, including their types and history, their motivation for use, their characteristics, and their applications.On the other hand, the numerical solvability of other version of differential equations and other related equations can be found in [36][37][38][39][40][41][42][43] and references therein.
The outline of the paper is as follows: in the next section, we utilize some necessary definitions and results from fractional calculus theory.In Section 3, the general form of generalized Taylor's formula is mentioned and proved.In Section 4, basic idea of the RPS method is presented in order to construct and predict multiple solutions for BVPs of fractional order.In Section 5, two nonlinear models are performed in order to illustrate the capability and simplicity of proposed method.Finally, conclusions are presented in Section 6.

Review of Fractional Calculus Theory
In this section, we present some necessary definitions and essentials results from fractional calculus theory.There are various definitions of fractional integration and differentiation, such as Grunwald-Letnikov's definition and Riemann-Liouville's definition [5,6,8].The Riemann-Liouville derivative has certain disadvantages when trying to model real-world phenomena with fractional differential equations (FDEs).Therefore, we shall introduce a modified fractional differential operator proposed by Caputo in his work on the theory of viscoelasticity [4].Throughout this paper, the set of integer numbers, the set of real numbers, and is the Gamma function.

General form of Generalized Taylor's Formula
In this section, we will introduce general form of generalized Taylor's formula that contains the Caputo definition for fractional derivatives.In fact, we need this generalization in the application of the RPS method in order to predict and find the multiplicity of solutions.
We will begin with the following definition which is needed throughout this work, especially, in the two succeeding sections.After that, we present a new and a fundamental theorem called general form of generalized Taylor's formula, which can formulate any function with certain properties in term of its fractional power series (FPS) representation.Definition 3.1: A power series of the form (3) is called FPS about , where is a variable and 's are constants called the coefficients of the series.As a special case, when , the expansion is called a fractional Maclaurin series.Notice that in writing out the term corresponding to and in Equation (3) we have adopted the convention that even when .Also, when each of terms of Equation (3) are vanished for and so.On the other hand, the FPS representation of Equation (3) always converges when .In the next lemma by we mean that ( -times).

Lemma 3.1: Suppose that and for
where .Then we have: with .
Proof: From the definition of the operator and by using the second mean value theorem for fractional integrals, one can find: (5) Theorem 3.1: Suppose that , , and can be differentiated -times on for , where .Then: with .
Proof: From the certain properties of the operator , we have: On the other direction, if we keep the repeating of this process, then after -times of computations, we will obtain: Entropy 2014, 16 (8) Thus, by using Lemma 3.1, the proof of the theorem will be complete.

Remark 3.1:
We mention here that, if we fixed , then the series representation of in Theorem 3.1 will leads to the following expansion of about : (9) with , which is the same as of Generalized Taylor's series that obtained in [44] for .As with any convergent series, this means that is the limit of the sequence of partial sums.In the case of general form of generalized Taylor's series, the partial sums are given as .In general, is the sum of its general form of generalized Taylor's series if .But on the other aspect as well, if we set , then is the remainder of general form of generalized Taylor's series.
That is, .

Corollary 3.1: If on
, where , then the reminder of general form of generalized Taylor's series satisfies the inequality: (10)

RPS Method for BVPs of Fractional Order
In this section, we predict and find multiple solutions for BVPs of fractional order that admit multiple solutions by substituting a FPS expansion with undetermined coefficients through the given equation.From the FDE a recursion formula for the computation of the coefficients was derived.On the other hand, the coefficients in the FPS expansion can be computed recursively by differentiating the FDEs.
For convenience, the reader is referred to [28][29][30] in order to know more details about the classical RPS methods, including their construction, their motivation for use, their characteristics compared to the conventional method, and their applications for solving different categories of linear and nonlinear differential equations of different types and orders.
In fact, the main goal of our work is to predict and find out multiple series solutions for nonlinear BVPs of fractional order.To illustrate the basic idea of the RPS method for solving fractional BVPs analytically, we first consider the following nonlinear fractional functional equation: (11) subject to the boundary conditions: Entropy 2014, 16 (12) where is a Caputo fractional derivative of order , is general nonlinear operator has Caputo fractional derivative term, is boundary operator, and is the boundary of the domain .The crucial step of the RPS method for solving Equations ( 11) and ( 12) analytically is based on the fact that the boundary conditions of Equation ( 12) should be transcribed into equivalent form, so that, the new version conditions involves an unknown parameter so-called prescribed parameter and are split to: (13) where is the forcing condition that comes from original conditions of Equation ( 12) and is the remainder boundary operator which contains the prescribed parameter .Now, we investigate and apply the RPS method on the following problem: (14) subject to the split conditions: (15) In order to apply the RPS method in the fractional sense, we assume firstly that we can apply the operator on the term of in Equation ( 14) and also we suppose that all solutions that satisfy Equations ( 14) and ( 15) can be expanded by FPS representation as follows: (16) where are coefficients to be determined and is the initial value of the independent variable of Equations ( 14) and (15).Assume that Equation ( 14) satisfies the initial conditions , with as a prescribed parameters, where the unknown parameters can be determined later by substituting and/or any other constraint conditions into the obtained solution form.It is worth noting that some of may be known from the given initial conditions in Equation (15).As the first step in the prediction of multiple solutions, we set be as an initial guess approximation of exact function solution of Equation ( 14).On the other hand, will be of the form which satisfies the known initial conditions in Equation (15) automatically.It was proved in [28][29][30] that the coefficients in Equation ( 16) will take the form . Therefore, we can consider the expansion formula: (17) Entropy 2014, 16 as an initial guess approximation for solution of Equations ( 14) and (15).But on the other aspect as well, depending on Equations ( 16) and ( 17), we can write: (18) Again, as the second step in the prediction of multiple solutions, we will let the double to denote the -truncated series approximation of .That is: (19) where the indices counter and whenever used mean that and .Prior to applying the RPS technique for finding the values of coefficients in the series expansion of Equation ( 19), we must define the residual function concept for the main nonlinear fractional functional Equation ( 14) as and the following truncated -resudial function: As in [28][29][30], it is clear that for each , where is the radius of convergence of Equation (18).In fact, this shows that for each and , since the fractional derivative of a constant function in the Caputo sense is zero.In the mean time, the fractional derivatives for each and of and are matching at ; it is obvious that: To obtain the value of coefficients in Equation ( 19) for and , we apply the following subroutine: substitute -truncated series approximation of into Equation (20), find the fractional derivative formula of at , and then finally solve the obtained algebraic equation to get the required coefficients.
To summarize the computation process of RPS method in numerical values, we apply the following: fixed and run the counter to find -truncated series expansion of suggested solution, next fixed and run the counter to obtain the -truncated series, and so on.In fact, to get -truncated series expansion for Equations ( 14) and ( 15), we use Equation ( 19) and write: (22) On the other hand, to determine the value of first unknown coefficient, , in Equation (22), we should substitute Equation ( 22) into both sides of the -residual function that obtained from Equation (20), to get the following result: (23) Now, depending on the result of Equation ( 21) for , Equation (23) gives .
Hence, using the -truncated series expansion of Equation ( 22), the -RPS approximation for Equations ( 14) and ( 15) can be expressed as: (24) Similarly, to find -truncated series expansion for Equations ( 14) and ( 15), we use Equation ( 19) and write: (25) Again, to find out the value of second unknown coefficient, c 11 in Equation ( 25), we must find and formulate -residual function based on Equation ( 20) and then substitute the form of in Equation ( 25) to find new discretized form of this residual function as follows: (26) while, on the other hand, by considering Equation (20) for and applying the operator to the both side of Equation ( 26), we get: Now, using the fact that , we can easily obtain: Hence, using the -truncated series expansion of Equation ( 25), the -RPS approximation for Equations ( 14) and ( 15) can be expressed as: (29) This procedure can be repeated till the arbitrary order coefficients of FPS solution for Equations ( 14) and ( 15) are obtained.Moreover, higher accuracy can be achieved by evaluating more components.

Remark 4.1:
It is worth indicating that there are still unknown prescribed parameters in the series expansion of Equation ( 19) (and simply in Equation ( 29)) that should be determined.It is essential that the existence of a unique or multiple solutions in terms of Equation (19) (and simply in Equation ( 29)) for the original BVP which is covered by Equations ( 11) and ( 12) depends on the fact that whether the forcing condition and/or any other constraint condition in Equation ( 13) admits a unique or multiple values for the formally introduced prescribed parameters .This stage is called rule of multiplicity of solutions that is a criterion in order to know how many solutions the BVP in Equations ( 11) and ( 12) admits.
Anyhow, as the final step in the construction, if we substitute and/or any other constraint conditions into the obtained solution form of Equation ( 19) (and simply in Equation ( 29)), then we obtain a system of nonlinear algebraic equations in the prescribed variables (here, we must recall that some of may be known from Equation ( 15)) which can be easy solved using symbolic computation software such as MAPLE or MATHEMATICA .In fact, if we substitute these values of prescribed parameters in the obtained solution form of Equation ( 19) (and simply in Equation ( 29)), then discretized form of the -truncated series approximation of of Equations ( 11) and ( 12) (and simply -truncated series approximation of as given by Equation ( 26)) will be obtained.

Applications and Numerical Discussions
The application problems are carried out using the proposed RPS method, which is one of the modern analytical techniques because of its iterative nature; it can handle any kind of boundary conditions and other constraints.The RPS method doesn't have mathematical requirements about the multiple solutions of fractional BVPs to be solved; the RPS method is also very effective in identifying global predicted solutions, and provides a great flexibility in choosing the initial guess approximations.However, in order to verify the computational efficiency of the designed RPS method, two nonlinear models are performed, one of them arises in mixed convection flows and the other one arises in heat transfer, which both admit multiple solutions.In the process of computation, all the symbolic and numerical computations were performed by using the MATHEMATICA software package.Throughout this section, we will try to give the results of the two applications; however, in some cases we will switch between the results obtained for the applications in order not to increase the length of the paper without the loss of generality for the remaining application and results.However, by easy calculations we can collect further results and discussion for the desire application.

Application 5.1:
The aim of this application is to apply the RPS method to analyze a kind of model in mixed convection flows namely, combined forced and free flow in the fully developed region of a vertical channel with isothermal walls kept at the same temperature.In this model, the fluid properties are assumed to be constant and the viscous dissipation effect is taken into account.The set of governing balance equations for the velocity field is reduced to the following [45,46]: (30) subject to boundary conditions: (31) where and are dimensionless velocity and transversal coordinate, respectively, and also , , , , , and in which , and are mean -component of the fluid velocity, fluid velocity, channel half-width, acceleration due to gravity, coefficient of thermal expansion, specific heat at constant pressure, dynamic viscosity, thermal conductivity, kinematic viscosity, Gebhart number, Prandtl number, and Reynolds number, respectively.
Next, we will show how one can find out the existence of multiple solutions for Equations ( 30) and ( 31) in aforesaid range for .To do so, we consider firstly Equations ( 30) and ( 31) and suppose that and .So, Equations ( 30) and ( 31) can be modified into the following form: (32) subject to the split conditions: (33) where is the additional forcing condition and is the additional constraint condition.Now, we apply the RPS method on Equations ( 32) and (33), where prescribed parameters and which are played an important and fundamental role to realize about multiplicity of solutions, will be obtained later by substituting the additional forcing condition and the additional constraint condition in resulting expansion formula that approximate Equations ( 32) and (33).
According to Equation ( 18), we assume that the series solution of Equations ( 32) and ( 33) can be written as: (34) where and which are hold from the conditions of Equation (33).Therefore, the initial guess approximation can be constructing as .Next, according to Equations ( 19) and ( 20) the -truncated series approximation of and the -residual function of Equation ( 32) can be defined and thus constructed, respectively, as: (35) (36) However, to determine the value of coefficient , we find out (1,0)-truncated series approximation of as and (1,0)-residual function as .On the other aspect as well, by using Equation (21) for and substituting , we obtain .Similarly, to find the value of coefficient , we evaluate -truncated series approximation of as and -residual function as . Thus, for , we conclude that , while the substitution of leads to .
To evaluate the value of coefficient , we need to write and .However, by considering the fact that , we can easily find .Similarly, the continuation in the same manner will leads also to .According to the initial guess approximation and the form of terms in Equation ( 32) taking into account the form of Equations ( 35) and ( 36), we can conclude that for .Therefore, according to Equation ( 34) the FPS solution of Equations ( 32) and ( 33) can be written in the form of the following expansion: (37) and hence the -truncated series approximation of can reformulated as: (38) Again, to determine the value of coefficient , we need solve the equation which gives .Similarly, we have , and so on.Consequently, based on the structure of Equation ( 38) the -truncated series approximation of generated from the RPS method can be written as: (39) It is clear that, Equation (39) contain two unknown certain parameters which are and .To determine their introductory values substituting the forcing condition and the constraint condition , and finally selecting some numerical values of .Now, to be specific, we consider two cases according to Equation ( 39) which consist of and .On the other hand, we generate and obtain the -truncated series approximation of using the same procedure.However, various values of and have been calculated and listed in Tables 1 and 2 when and when , respectively.For simplicity and not to conflict, we will let to denote the first approximate solution of and to denote the second approximate solution of .Consequently, we conclude that the RPS method furnishes multiple solutions for Equations (30) and (31).It is worth mentioning here that when , Table 1 indicates the existence of two solutions at , so that, , for the first branch solution and , for the second branch solution.In fact, these results answer the question how many solutions the nonlinear BVP in Equations ( 30) and ( 31) admits?The same procedure has been done at the case .As we see from Table 1, there exist multiple solutions namely for the first branch solution and for the second branch solution.Similar conclusion can be achieved when as shown in Table 2.The RPS technique has an advantage that it is possible to pick any point in the interval of integration and as well the approximate multiple solutions and all their fractional derivatives will be applicable.In other words, continuous approximate solutions can be obtained.Our next goal, is to show the mathematical behavior of the obtained multiple solutions geometrically.To do so, we plot the first and the second solutions obtained from the -truncated series approximation of at and when in Figure 1, while in Figure 2 we depict the first and the second approximate solutions at the same values when .The effective calculations of the two approximate branches solutions for Equations (30) and (31) with respect to some certain specific values of on and is explored next in which the obtained results are generated from the -truncated series approximation of .Table 3 gives the effect of the numeric value of when , while Table 4 gives the effect of the numeric value of when .The numeric value of lie within the range in step of .It is to be noted that, when the values of increasing gradually within mentioned range, the value of decreasing as well as the value of increasing for both branch solutions and for both order of derivatives.We mention here that, the case of correspond either to a very small viscous dissipation heating or to negligible buoyancy effects.However, Equations ( 30) and ( 31) are easily solved and admit the unique solution for both and .On the other direction, from the last tables, it can be seen that our results of the RPS method agree best with method of [34] when and method of [35] when .To measure the accuracy and the efficiency of the proposed RPS method for predicting and finding the multiple solutions for Equations ( 30) and ( 31), we report the residual error function at and when and when in Tables 5 and 6, respectively, in which the obtained results are generated from the -truncated series approximation of .The residual error function is defined using Equation (35) in which the grid points are building as , . For simplicity and not to conflict, we will let to denote the residual error function of the first approximate solution of and to denote the residual error function of the second approximate solution of .
In fact, the residual errors measure the extent of agreement between the th-order approximate RPS solutions and unknowns closed form solutions which are inapplicable in general for such nonlinear equations.However, from the tables, it can be seen that the RPS technique provides us with the accurate approximate solutions and explicates the rapid convergence in approximating the multiple solutions for Equations (30) and (31).In Tables 7 and 8 we tabulate the values of the approximate multiple solutions at the final grid node that generated from the -truncated series approximation of at and when and when .In fact, we do this to facilitate the calculations in order to show the validity and accuracy of the proposed RPS method in predicting and finding the multiple approximate solutions.In these tables, we can find that the values of , agree nicely and completely the forcing condition and the constraint condition , .In order to study the behavior of multiple approximate solutions in a better view, we plot the normalized of the two branch approximate solutions of Equations ( 30) and (31) with respect to , at some specific values of and in which the obtained results are generated from the -truncated series approximation of .However, Figure 3 shows the normalized function , at and when , while Figure 4 shows the normalized function , at and when .In these figures, we can see the almost similarity in the behavior of the two branches approximate solutions at the two mentioned specific value of when and when .
subject to the boundary conditions: The prediction and construction of multiple solutions for BVPs of fractional order is the fundamental target of this paper.Next, we will show in brief steps and calculations how we can predict and find out existence of multiple solutions for Equations (40) and (41).To do so, we consider firstly, Equations ( 40) and ( 41) and suppose that .So, a new discretized form of Equations ( 40) and ( 41) can be obtained as follows: (42) subject to the split conditions: (43) where is the additional forcing condition.Here, denotes temperature of the fin tip and will be determined later by the rule of multiplicity of solutions from the process of computations thought the RPS technique.
Similar to the previous procedure and discussions that used in Application 5.1, the FPS solution and the residual function for Equations ( 42) and ( 43) will take, respectively, the following form: (44) (45) while the -truncated series approximation of and the -resudial function that are derived from Equations ( 44) and ( 45) can be formulated, respectively, in form of: It is to be noted that the -truncated series solution of Equations ( 42) and ( 43) is and the -residual function is .Thus, using Equation ( 21) for , we get .Similarly, the -truncated series solution is and the -residual function is: More precisely, according to Equation ( 21) the solution of equation will gives . Thus, based on the initial guess approximation and the form of terms of Equation ( 42) taking into account the form of Equations ( 46) and (47), it easy to see that , .Therefore, according to Equation ( 44) a new discretized form of FPS solution for Equations ( 42) and ( 43) can be obtained and expressed as: (49) and hence the -truncated series approximation of can formulated as: (50) In the shape of shapes by continuing in this procedure and using Equations ( 46) and ( 47) taking into account Equation ( 21), we can easily obtained that , , and so on.Consequently, based on Equation (50) the -truncated series approximation of generated from the RPS method can be written as: (51) It is clear that, all the terms in Equation (51) contain an unknown certain parameter and to determine its introductory values we must substitute the boundary condition 1 back into Equation (51) to obtain a nonlinear algebraic equation in one variable, which can be easy solved using symbolic computation software.But on the other aspect as well, if we generate and obtain -truncated series approximation of by using the same procedure discussed, then two various values of have been calculated and listed in Table 9 when and when .
Table 9.The approximate numerical values of , when and when for .

0.447213595446
It is clear from the table that two -plateaus can be identified and consequently we conclude that the RPS method furnishes multiple solutions to Equations ( 40) and (41).It is worth mentioning here that Equation (46) (and simply Equation (51)) indicates existence of two solutions.On the other hand, the existence of a unique or multiple solutions in terms of Equation (46) (and simply Equation (51)) for the original BVP which is covered by Equations ( 40) and (41) depends on the fact that whether the forcing condition admits a unique or multiple values for the formally introduced prescribed parameters .
Finally, in Figure 5 we plot the first and the second approximate multiple solutions of Equations ( 40) and ( 41) that obtained from the -truncated series approximation of when and when .In fact, we do this for the same reasons that mentioned in the Application 5.1, where the same conclusion can be obtained too.On the other direction, as in the Entropy 2014, 16 previous application, we can see that the sketch of the two branches approximate RPS solutions that the problem admitted are agree best and nicely with method of [34] when and method of [35] .(a) (b)

Conclusions
It is very important not to lose any solution of nonlinear FDEs with boundary conditions in engineering and physical sciences.In this regard, the present paper has introduced a new methodology namely the RPS method to prevent this, so that the presented method is not only able to predict the existence of multiple solutions, but also to calculate all branches of solutions effectively at the same time by using an appropriate initial guess approximation.We also noted that the RPS solutions were computed via a simple algorithm without any need for perturbation techniques, special transformations, or discretization.The validity of this method has been checked by two nonlinear models, one of them arises in mixed convection flows and the other one arises in heat transfer, which both admit multiple or dual solutions.

Table 1 .
The

Table 2 .
The

Table 3 .
The effect values of on and , for the first and the second branch solutions when .

Table 4 .
The effect values of on and , for the first and the second branch solutions when .

Table 5 .
The values of absolute residual error function

Table 6 .
The values of absolute residual error function

Table 7 .
The

Table 8 .
The