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 (FNBVPs) 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 (Nguyen and Jang in Fract. Calc. Appl. Anal. 20(2):447–476, 2017) into Volterra integral equations which are equivalent to FIVPs. The advantage of the 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(h2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal{O}(h^{2})$\end{document} and O(h3)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal{O}(h^{3})$\end{document} 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 nonlocal property of the fractional derivative [2][3][4][5][6][7]. Many authors have introduced numerical methods for solving fractional differential equations arising in science and engineering. The authors in Refs. [8,9] proposed a computational algorithm based on reproducing kernel Hilbert space for solving time-fractional partial differential equations in porous media and nonlinear homogeneous and nonhomogeneous time-fractional equations. In Ref. [10], a numerical method based on multiple fractional power series solution was introduced to deal with the Schrödinger equation. The authors in Refs. [11][12][13][14] proposed several numerical methods based on collocation method, finite difference method, and L1 approximation for solving time-fractional diffusion equations. Regarding Fractional order Boundary Value Problems (FBVPs), Ref. [15] (and the references therein) investigated a Caputo fractional hybrid twopoint boundary value problem describing the thermostat models. In Ref. [16], the authors studied a fractional-order nonlocal continuum model of a Euler-Bernoulli beam whose governing equation is described as a FBVP, using the fractional finite element model.
Recently, the authors in Ref. [17] developed a spectral collocation method to deal with two-point linear multi-term FBVPs with Caputo fractional operator. In Ref. [18], the authors reformulated two-point FBVPs with a Riemann-Liouville fractional operator to a Volterra integral equation of the second kind and then developed an integral discrete scheme based on finite difference method.
Definition 1.2 let α ∈ R + . The operator D α a is defined by is the ceiling function and is the floor function.
The multi-term Caputo sense FIVP can be transformed into the system of FIVPs by Theorem 1.1 [19].
In this paper, we propose new schemes to deal with FNBVPs and the algorithms are summarized as follows: 1 In the case that 0 < α 1 < 1, we transform the FNBVP (1) with a = 0 into a system of FIVPs using Theorem 1.1.
In the case of α 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, the Gronwall inequality for two-term equations in [19] 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 Sect. 2. 2 To deal with FIVPs, we adopt high-order predictor-corrector methods (HPCMs) with linear interpolation and quadratic interpolation [1] into Volterra integral equations which are equivalent to FIVPs. 3 ICs of the FIVPs in the system equivalent to (1) are obtained by RBC at t = 0. But ICs include s := y (0) and since s is unknown, we approximate s by means of nonlinear shooting techniques based on Newton's method and Halley's method. The error function |a 2 y(b, s) + b 2 y (b, s)γ 2 | is used to construct the root-finding problem in order to make the approximate solution to y(t) satisfy the RBC at t = b. 4 The algorithm of the proposed shooting technique is as follows: The system of FIVPs is solved with an initial approximation to s, s k at the kth iteration. Using the approximate solution to the system obtained by HPCMs with s 0 , we find s 1 by solving Newton's (Halley's) formula. We update the approximate solution to the system with s 1 and measure the norm of the error function. We repeat this process until the norm of the error function is within the tolerance. Similar to our proposed schemes, the authors in Refs. [20,21] introduced numerical methods for solving FBVPs with RBCs. In Refs. [20,21], the FBVP with RBCs is turned into the FIVP by using a shooting method with a guess for the unknown IC y(0) and then the FIVP is transformed into the Volterra integral equation. The integral-differential term in the Volterra integral equation is approximated by an integral discretization scheme with constant and first-order interpolating polynomials in paper [20] and [21], respectively. However, the integral discretization schemes can only handle linear FBVPs and the rate of convergence depends on the fractional order α. This is elaborately addressed in Sect. 4. The main advantages of our proposed schemes are as follows: 1 The proposed schemes can handle both linear and nonlinear FBVPs with general RBCs. 2 Our proposed schemes can deal with multi-term FBVPs where 0 < α 1 ≤ 1 and 1 < α 2 < 2.  [1]. 4 It is not required to solve a matrix system as Newton's method and Halley's method are applied into a system of FIVPs. This article is organized as follows. In Sect. 2, we describe an idea about the transformation of FNBVP with RBCs (1) into a system of FIVPs. In Sect. 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 Sect. 3. In Sect. 4, we demonstrate numerical examples verifying that the 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 Ref. [21]. A conclusion will be given in Sect. 5. Finally tables of numerical results and the linear explicit method which is an alternative method for solving FIVPs are described in the Appendix.
Case 1: 0 < α 1 < 1 First, we consider a FNBVP with Dirichlet boundary conditions as follows: Applying Theorem 1.1 with β 1 := α 1 , β 2 := 1α 1 , β 3 := α 2 -1, the FNBVP with Dirichlet boundary conditions (4) can be transformed as follows: From the system of fractional differential equations (5), we obtain the following system of FIVPs: where the IC s is unknown and so needs to be approximated. Similar to the case of Dirichlet boundary conditions, the FNBVP with RBCs (1) can be written as follows: Case 2: α 1 = 1 We consider the following FBVP with Dirichlet boundary conditions: where 1 < α 2 < 2, α 2 ∈ R. Since the fractional differential equation in (8) does not satisfy the assumption, 0 < α 1 < 1 in Theorem 1.1, we cannot apply the strategy used in (4) to (8). So we modify the equation in (8) to meet the assumption, with the same boundary conditions as in (8) as follows: where α 2 ∈ (1, 2), → 0+. By Lemma 2.1, solutions of the two FBVPs (8) and (9) are approximately equal and the absolute error depends on .
Lemma 2.1 (First Gronwall inequality for two-term equations in [19]) Let α 2 > 0 and α 1 ,α 1 ∈ (0, α 2 ) be chosen so that the equation subject to the initial conditions subject to the same initial conditions (where f satisfies a Lipschitz condition in its second and third arguments on a suitable domain) has unique continuous solutions y, z : [0, T] → R. We assume further that α 1 = α 1 . Then there exist constants K 1 and K 2 such that where E α n denotes the Mittag-Leffler function of order α n .

Nonlinear shooting methods and high-order predictor-corrector methods
FBVPs have been transformed to systems of FIVPs in Sect. 2. Before we address how to deal with systems of FIVPs (6), (7), (11), and (22) 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 a 2 y(b) + b 2 y (b) = γ 2 is involved in approximating s, we define y(s) := y(s, t)| t=b and let the error function be F(s) := a 2 y(s) + b 2 ∂ ∂t y(s)γ 2 . 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 throughout this section.

Shooting with Newton's method
The conventional Newton 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 the 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 Sect. 3.3. By solving the system (25), s k+1 in Newton's formula (23) 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 (22). 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 formula for F(s) = 0 is as follows: where F s (s k ) is described in (24) 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 (25) or (28). Once we find the s k , we solve systems of FIVPs (6), (7), (11), or (22). In this subsection, we describe how to deal with those systems of FIVPs using High-order Predictor-Corrector Methods (HPCMs) introduced in Ref. [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 . We divide the domain as follows: For simplicity, let the step size be uniform, which means t j+1t j = h, j = 0, 1, . . . , N -1. Then (30) can be rewritten at time t n+1 as follows: where 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: where and the predictor f p n+1 is found as follows: where 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 [20,21]. 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. max 1≤j≤N |y c jy(t j )|) • 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 (22) where the exact solution is y(t) = t 4 .
Example 4.3 Consider the following double-term FNBVP with RBCs: 6 (4-α 1 ) t 3-α 1 + D α 1 0 y(t),    Table 3, Tol = 10 -10 for both shooting techniques in Table 6, and Tol = 10 -15 for Newton's method and Tol = 10 -16 for Halley's method in Table 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 Tables 3 and 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.       Fig. 4) and Halley's method (in Fig. 5), respectively.
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 A 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 [21] for each α 2 = 1.1, 1.3, 1.5, 1.7, 1.9. In Table 10, we can observe that our proposed method consumes less CPU time than the modified integral discretization scheme [21] even though both methods shows the equal performance.
Example 4.5 Consider the following single-term linear FBVP with RBCs: 1-α 2 y (0) = γ 1 , y(1) + y (1) = γ 2 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 [21]. Table 11 shows pointwise absolute errors and computed convergence profiles versus the number of sub-intervals for each α 2 = 1.1, 1.3, 1.5, 1.7, 1.9. In Table 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 [21] are around 2.0.
Halley's method   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 [21] because the predictor and corrector in HPCMs share the computation of the memory effect. As a result the proposed shooting technique based on Newton's method consumes less CPU than the modified discretization scheme [21] 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 [21], as shown in Table 11.   Figure 6 illustrates convergence profiles obtained by the proposed methods and the modified integral discretization scheme [21] with the variety of fractional orders. From the graphs, we can see that computed rates of convergence are nearly O(h 2 ) for the proposed shooting method with Newton's technique and the modified integral discretization scheme [21], O(h 3 ) for the proposed shooting method with Halley's scheme, respectively. Plots exhibited in Figs. 7, 8, 9 display pointwise absolute errors in maximum norm versus the time step h and the approximate solution is computed by the modified integral discretization scheme [21] (in Fig. 7), the proposed shooting technique with Newton's  Fig. 8), the proposed shooting technique with Halley's method (in Fig. 9), respectively. The exact solution of Example 4.6 is unknown so we alter the measure to estimate the errors employing the uniform two-mesh differences and orders of convergence introduced in [22]: where y c 2j denotes the approximate solution at t 2j by the proposed methods with h = 1/2N . Using the errors e N , we estimate the convergence profiles as follows: In this example, we compare performances of our proposed shooting methods with the finite difference method proposed in [22] and the numerical results are shown in Table 12.
We observe the following from the results: 1 In Table 12, uniform two-mesh difference errors and convergence rates computed by (a) the finite difference method [22],  2 It is evident that the uniform two-mesh difference errors computed by our proposed methods are less than the finite difference method [22] for all N and α 2 .

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 the proposed shooting methods based on Newton's and Halley's method and this is the main algorithm of the 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 the proposed schemes can handle doubleterm FNVBPs with RBCs whose exact solutions include polynomial, exponential, and sine function. Convergence profiles obtained by the proposed schemes were computed as expected by the global error estimates. However, Tables 3, 6, and 9 suggest 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 FB-VPs with exact solutions having low regularity and high regularity, respectively. Tables 10 and 11 showed that the 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 A 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 [21].