An Efficient Numerical Approach for Solving Two-Point Fractional Order Nonlinear Boundary Value Problems with Robin Boundary Conditions

This article proposes new strategies for solving two-point Fractional order Nonlinear Boundary Value Problems(FN-BVPs) with Robin Boundary Conditions(RBCs). In the new numerical schemes, a two-point FNBVP is transformed into a system of Fractional order Initial Value Problems(FIVPs) with unknown Initial Conditions(ICs). To approximate ICs in the system of FIVPs, we develop nonlinear shooting methods based on Newton's method and Halley's method using the RBC at the right end point. To deal with FIVPs in a system, we mainly employ High-order Predictor-Corrector Methods(HPCMs) with linear interpolation and quadratic interpolation into Volterra integral equations which are equivalent to FIVPs. The advantage of proposed schemes with HPCMs is that even though they are designed for solving two-point FNBVPs, they can handle both linear and nonlinear two-point Fractional order Boundary Value Problems(FBVPs) with RBCs and have uniform convergence rates of HPCMs, $O(h^2)$ and $O(h^3)$ for shooting techniques with Newton's method and Halley's method, respectively. A variety of numerical examples are demonstrated to confirm the effectiveness and performance of the proposed schemes. Also we compare the accuracy and performance of our schemes with another method.


Introduction
Fractional calculus has proven to describe many phenomena in science and engineering more accurately than integer order calculus because of the non-local property of the fractional derivative [2,3,4,5,6,7]. Many authors have introduced analytical and numerical methods of the Fractional order Initial Value Problems (FIVPs). Relatively, Fractional order Boundary Value Problems (FBVPs) have been paid less attention than FIVPs. Moreover, numerical methods of FNBVPs with Robin Boundary Conditions (RBCs) have hardly been studied. In this paper, we consider the two-point Fractional order Nonlinear Boundary Value Problem (FNBVP) with RBCs: where 0 < α 1 ≤ 1, 1 < α 2 < 2, α 1 , α 2 , γ t ∈ R. D α 1 a and D α 2 a are Caputo fractional differentiations defined as follows: For α = 0, J 0 a = I, the identity operator. Definition 1.2. let α ∈ R + . The operator D α a is defined by where ⌈ ⌉ is the ceiling function and ⌊ ⌋ is the floor function.
In this paper, we propose new schemes to deal with FNBVPs and the algorithms are summarized as follows: 1. In case of that 0 < α 1 < 1, we transform the FNBVP (1) with a = 0 into a system of FIVPs using Theorem 1.1.
In case of that α 1 = 1, i.e. the FNBVP (1) has a single term of fractional order α 2 , we substitute the integer order α 1 = 1 with the fractional order α 1 = 1 − ε, ε → 0+ so that the FNBVP satisfies the assumption, 0 < α 1 < 1 in Theorem 1.1. First Gronwall inequality for two-term equations in [8] guarantees that the difference between the solution of FNBVP with α 1 = 1 and with α 1 = 1 − ε approaches 0 as ε → 0+. The FNBVP with α 1 = 1 − ε is transformed into a system of FIVPs and then we reduce the number of equations in the system. We prove the reduced system is equivalent to the original system as ε → 0+ in Section 2. This article is organized as follows. In Section 2, we describe an idea about the transformation of FNBVP with RBCs (1) into a system of FIVPs. In Section 3, we describe nonlinear shooting methods based on Newton's method and Halley's method, to approximate unknown IC s := y(0) of FIVPs in the system. Also, we briefly mention how to apply the HPCMs into a system of FIVPs in Section 3. In Section 4, We demonstrate numerical examples verifying that proposed shooting techniques combined with HPCMs guarantee the global convergence rates of HPCMs. We also confirm the performance and effectiveness of the proposed methods by comparing with the modified integral discretization scheme in the paper [10]. A conclusion will be given in Section 5. Finally tables of numerical results and the linear explicit method which is an alternative method for solving FIVPs are described in the Appendix.
Applying Theorem 1.1 to (9), we obtain the following system of FIVPs: Now, we show that the system of FIVPs (10) is equivalent to the following system as ε → 0 using Lemmas 2.2 through 2.4 and Theorem 2.1.

Nonlinear Shooting Methods and High-Order Predictor-Corrector Methods
FBVPs have been transformed to systems of FIVPs in Section 2. Before we address how to deal with systems of FIVPs (6), (7), (11), and (21) using High-order Predictor-Corrector Methods (HPCMs), the unknown IC z(0) = s should be handled first. In this section, we describe two nonlinear shooting techniques based on Newton's method and Halley's method to approximate s. Both Newton's formula and Halley's formula are designed to determine the solution of a system of FIVPs satisfying the RBC at the right end point of an interval. Without loss of generality, we consider the system of FIVPs (7) that is equivalent to the FNBVP with RBCs (1).
In order that the RBC at the right end point We approximate the solution of the root-finding problem F(s) = 0 by using Newton's method and Halley's method, respectively. For convenience, we denote

Shooting with Newton's Method
The conventional Newton's formula for F(s) = 0 can be expressed as follows: where m is the maximum number of iterations and Observing y s (s k ) and y ts (s k ), it turns out that they are equal to ∂ ∂ s y(t)| s=s k ,t=b and ∂ ∂ s z(t)| s=s k ,t=b , respectively in the system of FIVPs (7). Thus we solve the following system obtained from system of FIVPs (7) by applying the operator ∂ ∂ s using HPCMs for each k: Since both t and s are independent variables, f s (t, y(t), w(t)) can be written as The detailed description of HPCMs dealing with a system of FIVPs is in Subsection 3.3. By solving the system (24), s k+1 in the Newton's formula (22) is computed. Using the updated approximate value of IC s, s k+1 , we update approximate solutions of systems of FIVPs (6), (7), (11), and (21). Repeating this process, we obtain an s k having an acceptable error of the root-finding problem F(s) = 0 at an appropriate number of iterations k.

Shooting with Halley's Method
The conventional Halley's formula for F(s) = 0 is as follows: where F s (s k ) is described in (23) and Similar to the way we found y s (s k ) and y ts (s k ) in the shooting with Newton's method, we find y ss (s k ) and z ss (s k ) by solving the following system of FIVPs obtained by applying the operator ∂ 2 ∂ s 2 using HPCMs for each k: Since t and s are independent variables, f ss (t, y(t), w(t)) can be written as

High-order Predictor-Corrector Methods for system of FIVPs
In order to find a s k with an acceptable accuracy, we iteratively solve systems of FIVPs (24) or (25). Once we find the s k , we solve systems of FIVPs (6), (7), (11), or (21). In this subsection, we describe how to deal with those systems of FIVPs using High-order Predictor-Corrector Methods (HPCMs) introduced in paper [1]. Without loss of generality, we consider the following FIVP: For convenience, let us denote y j as approximated value of y(t j ) except for y 0 = c 0 and let f j ≡ f (t j , y j ), y c j be a corrector of y j , y p j be a predictor of y j , and f p j ≡ f (t j , y p j ), j = 1, · · · , N. If j = 0 then, f 0 = f (0, c 0 ). We divide the domain Ω as follows: For simplicity, let step size be uniform which means t j+1 − t j = h, j = 0, 1 · · · , N − 1. Then (26) can be rewritten at time t n+1 as follows: We interpolate f (τ, y(τ)) using linear or quadratic Lagrange polynomials over each interval I j = [t j ,t j+1 ], j = 0, 1, · · · , N − 1. Then we obtain the following predictor-corrector schemes: 1. HPCM with linear Lagrange polynomial: 2. HPCM with quadratic Lagrange polynomial: and the predictor f p n+1 is found as follows: We implement HPCMs in a FIVP. For the purpose of the implementation to the system (6), as an example, we find predictors y p n+1 , w p n+1 , z p n+1 with s 0 individually. We compute f p n+1 ≡ f (t n+1 , y p n+1 , w p n+1 ), and we find correctors y c n+1 , w c n+1 , z c n+1 . Using proposed shooting techniques with HPCMs we find s 1 and then find predictors y p n+1 , w p n+1 , z p n+1 replacing s 0 with s 1 , compute f p n+1 , and update correctors y c n+1 , w c n+1 , z c n+1 . We repeat this process until the absolute value of the approximated error function is within a tolerance. Similarly, we apply the scheme with HPCMs into other systems of FIVPs (7), (11), and (21).
Remark 3.1. Alternatively, we only find predictors y p n+1 , w p n+1 . We then compute z n+1 using y p n+1 , w p n+1 , and then we compute w n+1 using z n+1 as a predictor, and y n+1 using w n+1 as a predictor. But it turns out that numerical results obtained by both ways of implementing HPCMs are nearly identical.
The following theorems [1] bound the Global error E n+1 of the HPCM with linear and quadratic interpolations, respectively.

Numerical Examples
In this section, we experimentally illustrate the performance of the proposed schemes. Numerically, we verify that our proposed schemes can deal with more complex FBVPs than the integral discretization schemes in [9,10]. For that purpose, the proposed schemes are implemented in FNBVPs with 0 < α 1 < 1 whose exact solutions are polynomial, exponential, and sine functions in Examples 4.1 through 4.3. We investigate absolute errors in maximum norm, convergence rates, and absolute values of the approximated error function |F(s k )| with various values of parameters. We discuss linear FBVPs with α 1 = 1 whose exact solutions have low regularity and high regularity in Examples 4.4 and 4.5, respectively. We compare numerical results obtained by our proposed schemes with the integral discretization schemes. But we emphasize that our proposed methods can deal with many different FBVPs unlike the another method in Examples 4.4 and 4.5. Regarding the numerical results shown in the Appendix, let us summarize the parameters used • h denotes the size of time sub-interval.
• s 0 denotes the initial approximation of the sequence {s k } in proposed shooting methods.
• k denotes the index of sequence {s k } generated by the proposed Newton's method or Halley's method. It can be considered as the number of iterations needed to meet a tolerance.
• m denotes the maximum number of iterations in Newton's and Halley's methods.
• Tol denotes the tolerance used to measure the error of the approximated error function |F(s k )| in Newton's method and Halley's method.
• N denotes the number of time sub-intervals. • E α,β (t) denotes the two-parameter function of Mittag-Leffler type [7].
In Examples 4.1 through 4.3, we transform the FNBVP into the system of FIVPs (7) and s 0 means an initial approximation to y ′ (0). In Examples 4.4 and 4.5, the linear FBVP is transformed into the system of FIVPs (21) with ε = 10 −10 and s 0 means an initial approximation to y(0 where the exact solution is y(t) = t 4 .

Example 4.3. Consider the following double-term FNBVP with RBCs
and the exact solution is y(t) = sin(λ t) − t + t 3 6 . In 4. Proposed methods are tested for a variety of values of α 1 , α 2 and numerical results for each pair of (α 1 , α 2 ) are shown in Tables A.3, A.6, A.9. For each pair of fractional orders (α 1 , α 2 ) pointwise absolute errors in the maximum norm, computed convergence rates, and number of iterations k such that |F(s k )|< Tol versus the number of subintervals N are listed in the tables. The initial approximation to s was set s 0 = 0.2 in all three tables. In order to minimize the number of iterations k, the tolerance was set Tol = 10 −5 for Newton's method and Tol = 10 −10 for Halley's method in Table A.3, Tol = 10 −10 for both shooting techniques in Table A.6, and Tol = 10 −15 for Newton's method and Tol = 10 −16 for Halley's method in Table A.9. Numerical results shown in the tables demonstrate that for all suggested pairs of fractional orders, rates of convergence approach 2 for Newton's method, 3 for Halley's method that are theoretical convergence rates of HPCMs. In Examples A.3 and A.9, we observe that the number of iterations k required to meet the tolerance at (0.9, 1.1) is relatively greater than other pairs of fractional orders for both Newton's and Halley's method.
Since D α 1 0 y(t), y ′ (t), D α 2 0 y(t) do not belong to C 3 [0, 1], global error estimates of HPCMs in Theorems 3.1 and 3.2 cannot be applied to Example 4.4. Alternatively, we adopt the linear explicit method described in Appendix B with proposed shooting techniques. In this example, we compare the accuracy and convergence rate of the approximate solution obtained by the proposed shooting technique based on Newton's method with the modified integral discretization scheme [10] for each α 2 = 1.1, 1.3, 1.5, 1.7, 1.9. In Table A.10, we can observe that our proposed method shows equal performance to the modified integral discretization scheme [10].
Example 4.5. Consider the following single-term Linear FBVP with RBCs whose 1 < α < 2 and the exact solution is In Example 4.5, we compare the performance of our proposed methods with the modified integral discretization scheme [10]. Table A.11 shows pointwise absolute errors and computed convergence profiles versus the number of subintervals for each α 2 = 1.1, 1.3, 1.5, 1.7, 1.9. In Table A.11, we can see that for all values of α 2 , the computed rates of convergence obtained by the proposed shooting technique based on Halley's method combined with third order HPCM are around 3.0 while the computed rates of convergence obtained by the modified integral discretization scheme [10] are around 2.0. The algorithm of the proposed shooting techniques with second-order HPCM requires less than the number of arithmetic operations needed by the modified integral discretization scheme to solve a FBVP with RBCs than the modified integral discretization scheme [10] because the predictor and corrector in HPCMs share the computation of the memory effect. In practice this results the proposed shooting technique based on Newton's method consume less CPU than the modified discretization scheme [10] and the CPU time executed by the proposed shooting technique based on Halley's method is approximately equal to the CPU time executed by the modified integral discretization scheme [10] as shown in Table A.11.

Conclusion
We introduced new numerical schemes for solving FNBVPs with any RBCs. The idea was to transform a FNBVP into a system of FIVPs. By doing that we could adopt a pre-existing numerical method for solving the system of FIVPs and we mainly employed HPCMs. The unknown IC s in the system was approximated by proposed shooting methods based on Newton's and Halley's method and this is the main algorithm of proposed schemes. Under the assumption that m is large enough so that |F(s m )| is small enough, theoretical convergence rates of proposed methods were O(h 2 ) for shooting with Newton's method and O(h 3 ) for shooting with Halley's method on account of global error estimates in HPCMs. In Examples 4.1 through 4.3, we verified that proposed schemes can handle double-term FNVBPs with RBCs whose exact solutions include polynomial, exponential, and sine function. Convergence profiles obtained by proposed schemes were computed as expected by the global error estimates. However, Tables A.3, A.6, and A.9 point out that the convergence rate of the sequence {s k } depends on fractional orders. We still need to address an error analysis of shooting techniques based on Newton's and Halley's methods for solving a system of FIVPs. This will be considered in a subsequent paper. Examples 4.4 and 4.5 demonstrated the performance of proposed methods for solving single-term linear FBVPs with exact solutions having low regularity and high regularity, respectively. Tables A.10 and A.11 showed that proposed methods can deal with not only nonlinear FBVPs but also linear FBVPs. In Example 4.4, we adopted the linear explicit method described in Appendix B and this shows that the proposed shooting techniques can be assembled with not only HPCMs but also other pre-existing numerical schemes for solving a system of FIVPs. In Example 4.5, we observed that computed convergence rates obtained by our proposed shooting technique based on Halley's method with third order HPCM are higher than the modified integral discretization scheme [10].