The application of boundary value problems (BVPs) can be found in different fields of pure and applied sciences; for instance, the narrow converting layers bounded by stable layers, which are believed to surround A-type stars, may be modeled by BVPs [1]. Also, these problems may model the dynamo action in some stars [2]. More discussions on the application of BVPs have also been provided in Chandrasekhar [3], Baldwin [4], and Khalid et al. [5]. More to the point, the approximation schemes to solve non-linear BVPs can be found in different sources of numerical analysis [6, 7]. In Agarwal [8], Agarwal discussed the existence of unique solution for these problems; however, no numerical method is contained therein. Boutayeb and Twizell [9] developed the finite difference methods to solve the above-mentioned problems effectively. They also improved a second-order method in Twizell and Boutayeb [10] to solve the general and special BVPs. Besides, Twizell [11] advanced a finite difference scheme of order two to investigate the solution of these problems. However, the existing methods suffer from enormous computational effort. To solve this difficulty, some alternative schemes have been presented including the Adomian decomposition method (ADM) with Green's function [12], homotopy perturbation method (HPM) [13], and variational iteration method (VIM) [14].

During the past decades, many scientists have frequently shown that the mathematical equations with fractional calculus architectures can describe the reality more precisely than the classic integer models with ordinary time-derivatives [15–19]. Recently, the advantages of this approach have been extensively investigated for various practical applications [20–27]. Concerning the fractional BVPs, some noticeable efforts have been done in Ali et al. [28] and Ugurlu et al. [29]. The aforesaid problems have also noteworthy real applications in different areas of science and technology. For instance, a hybrid Caputo fractional modeling was considered in Baleanu et al. [30] for thermostat with hybrid boundary conditions. In Patnaik et al. [31], the application of a fractional-order non-local continuum model was studied for a Euler-Bernoulli beam. The authors in Salem et al. [32] analyzed the coupled system of non-linear fractional Langevin equations with multi-point and non-local integral boundary conditions. The existence of extremal solutions of fractional Langevin equation involving non-linear boundary conditions was also investigated in Fazli et al. [33]. However, the properties of fractional BVPs should be studied deeply and approximation schemes should be continuously improved solving the above-mentioned problems appropriately. To this end, some valuable studies have been carried out, and a number of noteworthy results have been achieved. For instance, an existence theorem was discussed in Zhang and Su [34] for a linear fractional differential equation (FDE) with non-linear boundary conditions by using the method of upper and lower solutions in reverse order. In Arqub et al. [35], a new kind of analytical method was proposed to predict and represent the multiplicity of solutions to non-linear fractional BVPs. In Khalil et al. [36], the authors studied a coupled system of non-linear FDEs whose approximate solution was achieved under two different types of boundary conditions. In Cui et al. [37], a monotone iterative method was investigated for non-linear fractional BVPs while the fractional order was considered between 2 and 3. In Asaduzzaman and Ali [38], the existence of positive solution was investigated to the BVPs for coupled system of non-linear FDEs.

Motivated by the aforementioned statement, this manuscript aims to design a new iterative method to generate the approximate solution of non-linear fractional BVPs in the form of uniformly convergent series. The proposed method is free from perturbation, discretization, linearization, or restrictive assumptions. Moreover, contrary to the VIM [14] or the ADM [12], the suggested technique provides the exact solution without identifying the Lagrange multipliers or calculating the Adomian's polynomials. The new scheme just requires solving a sequence of linear fractional-order BVPs. Finally, four numerical examples are solved to verify the efficiency of the new technique.

The rest of paper is structured in the following way. Hereinafter, we review the fractional calculus approach and its main definitions. Section 3 describes the problem statement. A numerical technique is extended in section 4 solving non-linear fractional BVPs. Numerical and comparative results are reported in section 5, and finally, the paper is finished in section 6 by some concluding remarks.

This part is devoted to some preliminary results concerning the fractional operators. In the following, the Caputo derivative and the Riemann-Liouville integral are introduced, and their main properties are investigated as well [15].

Definition 2.1. For t ∈ (0, T) and n − 1 < α ≤ n, the αth-order Caputo derivative of a function x(t) is defined by

where Γ(·) is the gamma function. The corresponding Riemann-Liouville integral is also described as

With regard to the Caputo derivative (1), we can write

Furthermore, the Caputo derivative of a constant function is zero, i.e., if x(t) ≡ k, then we have Dtα0Ck=0. Additionally, the derivative and integral operators (1) and (2) satisfy the following anti-derivative property

More to the point, the Lipschitz condition is satisfied by the Caputo derivative (1)

where L > 0 is the Lipschitz constant.

For additional information, the interested readers can refer to Kilbas et al. [15].

To formulate a fractional BVP, consider the following FDE

where the function f(·) is analytic with regard to its arguments and f(0, t) = 0, ∀t ∈ (0, T). The expression Dtα0C denotes the αth-order Caputo derivative, and n is an even number. The boundary conditions for Equation (6) are given by

where a_2k, b_2k (k=0,1,…,n2) are real finite numbers. As is well-known, the exact solution of the fractional BVP (6)-(7) can hardly be achieved except in very special cases. Hence, an efficient iterative technique will be developed hereinafter in order to derive the corresponding approximate solution.

In this section, an efficient iterative method is improved to solve the fractional BVP (6), (7). To this end, first the following lemma is presented and proved.

Lemma 4.1. The solution of the fractional BVP (6)-(7) is analytic with respect to the boundary conditions a_2k, b_2k, k=0,1,…,n2.

Proof: Let x(·) be the solution of the BVP (6)-(7). Define αi=x(i)(0) and βj=x(j)(T), i, j = 0, …, n − 1. Then x(·) is the solution of the following initial value problems (IVPs)

Since f(x(t), t) is assumed to be analytic, x(·), as the solution of the IVPs (8) and (9), is analytic with respect to α_i and β_i, respectively [39]. Thus, x(·), as the solution of the BVP (6)-(7), is analytic with respect to a_2k, b_2k, k=0,1,…,n2.

Now, we state and prove the following theorem.

Theorem 4.1. The solution of the fractional BVP (6)-(7) is expressed by the uniformly convergent series x(t)=∑i=1∞x^i(t), where x^i(t) is attained by solving the sequence of linear fractional BVPs

and for i = 2, 3, 4, …

The non-homogeneous term F_i is determined by

λj(t)=1j!∂j∂xjf(x,t)|x=0, and the summation ∑k1,…,ki+1−j is taken over all combinations of non-negative integer indices k₁ through k_i+1−j such that

Proof: By using the Maclaurin series of f(x(t), t) with respect to x(t), we have

where λj(t)=1j!∂j∂xjf(x,t)|x=0. Besides, the solution of the fractional BVP (6)-(7) for an arbitrary vector x_b = (a₀, a₂, …, a_n, b₀, b₂, …, b_n) is expressed by

where the vector function g : ℝⁿ × (0, T) → ℝ is analytic based on Lemma 4.1. In addition, we have g(0, t) = 0, ∀t ∈ (0, T), since we have assumed that f(0, t) = 0 for all t ∈ (0, T). Therefore, by applying the Maclaurin series of g(x_b, t) with respect to x_b, from Equation (15) we derive

Since the function g(x_b, t) is analytic with respect to x_b, the Maclaurin series (16) exists and is uniformly convergent. Now, we perturb the boundary conditions by an arbitrary parameter ε > 0, i.e., x_b → εx_b. Then, Equation (16) is reformulated by

Substituting x(t) from Equation (17) into the expansion (14) yields

Rearranging Equation (18) with respect to the order of ε results

where

and the summation ∑k1,…,ki+1-j is taken over all combinations of non-negative integer indices k₁ through k_i+1−j such that

Since, Equation (19) must be satisfied for any ε > 0, we should equalize the coefficient of εⁱ on the left-hand side of Equation (19) with its corresponding coefficient on the right-hand side. This procedure yields

Now, we put t = 0 and t = T in Equation (17) and in its second- and fourth-order derivatives in order to achieve the boundary conditions for the sequence (22)-(24). Again, we should equalize the coefficients of εⁱ on the both sides of the resultant equations. Thus, we obtain

and the proof is complete.

As can be seen, Equation (10) formulates a homogeneous linear BVP of fractional-order. By solving this problem, x^1(t) is achieved in the first step. Following the proposed procedure in Theorem 4.1, we then obtain x^i(t) (i ≥ 2) by solving the non-homogeneous linear fractional BVP (11) in the ith step. Moreover, the non-homogeneous term in (11) is determined from Equation (12) by using the known functions provided in the previous steps. Thus, a recursive procedure should be employed here to solve the considered sequence.

In this part, four numerical examples are employed in order to verify the effectiveness of the new suggested technique. Here, we consider the examples form [13, 14] for the purpose of comparison with the other existing results.

Example 5.1. Consider a fractional BVP in the form below

whose exact solution is x(t) = e^t for α = 6.

Following the new technique as in section 4, we solve the presented sequence of fractional BVPs (10)-(11) in a recursive manner. Simulation results up to 10th iteration for α = 6 are reported in Table 1. As is shown, the error is reduced further by considering more components of x(t). To achieve an approximate solution with enough accuracy, the new algorithm is applied with δ = 0.01. From Table 1, we observe that the convergence is achieved just in the second step, i.e., ‖x2(t)-x1(t)‖∞=2.2×10-3<δ. Simulation curve of x₂(t) and the exact solution are plotted in Figure 1. This figure indicates that the second-order approximate solution is in good agreement with the exact solution.

The problem (31) for α = 6 has also been solved by using the HPM [13] and the VIM [14], respectively. Notice that the results of both methods are exactly the same as shown in Noor et al. [14]. Table 2 depicts the exact solution and the absolute errors achieved by applying two iterations of the HPM, VIM, and the proposed technique in this paper. Comparative results in this table verify the superiority of the suggested algorithm compared to the other approximation methods available in the literature.

Example 5.2. Consider the following non-linear BVP of fractional-order

whose exact solution is in the form x(t) = e^−t for α = 6.

Following the same procedure as in Example 5.1, we report the simulation results up to 10th iteration in Table 3. This table shows that considering more components of x(t) provides more precise results. From this table, it is also indicated that the proposed algorithm with δ = 10⁻⁴ converges after only two iterations, i.e., ‖x2(t)-x1(t)‖∞=1.3343×10-5<δ. In Figure 2, the simulation curve of x₂(t) is compared with the exact solution. Comparative results indicate that the second-order approximate solution is very close to the exact solution. Figure 3 shows the relation between the iteration time and the error given by the expression (28) using infinite norm for Examples 5.1 and 5.2. In this figure, the logarithmic scale is applied for the vertical axis. This figure verifies that the error decreases significantly as the iteration time increases.

The problem given by Equation (32) for α = 6 has also been solved by using the HPM and the VIM in Noor and Mohyud-Din [13] and Noor et al. [14], respectively. As can be seen in Noor et al. [14], the results of both methods are exactly the same. Table 4 exhibits the exact solution along with the absolute errors related to the HPM, VIM, and the proposed iterative algorithm. Comparing the results shows that the new approach is superior to the other existing methods.

Example 5.3. Consider the following non-linear fractional BVP

whose exact solution is x(t)=Eα(tα) where E_α(·) is known as the Mittag-Leffler function.

Simulation curve of x₂(t), i.e., the second-order approximate solution, for different values of α are plotted in Figure 4A. This figure indicates that the approximate solution tends to the classic integer solution for α = 6 when α → 6 as expected.

Example 5.4. Consider the non-linear fractional BVP

whose exact solution is x(t)=Eα(-tα).

In the same way as in Example 5.3, Figure 4B depicts the second-order approximate solution tending to the classic integer solution as α goes to 6.

This paper studied a new iterative scheme to provide the solution of non-linear fractional BVPs in terms of a uniformly convergent series. The proposed procedure was free from perturbation, discretization, linearization, or restrictive assumptions. Furthermore, contrary to the other approximation schemes such as ADM [12] and VIM [14], the suggested technique kept away from calculating the Adomian's polynomials or identifying the Lagrange multipliers, respectively. Hence, from practical viewpoint, the suggested technique is more efficient than the above-mentioned approximation methods. Simulation results, demonstrating the efficacy, high accuracy, and simplicity of the proposed method, were also included. In the following, we summarize the main aspects of our numerical findings. Tables 1, 3 provided the simulation results up to 10th iteration, and Figure 3 depicted the relation between the iteration time and the error given by the expression (28). From these results it is obvious that the error is reduced further by considering more components of x(t). Simulation curves in Figures 1, 2 also indicated that the second-order approximate solution is in good agreement with the exact solution. Tables 2, 4 exhibited the exact solution and the absolute error derived by employing two iterations of the HPM [13], VIM [14], and our new iterative algorithm. These tables clearly indicated the improvements made by employing the proposed method. The simulation curves for different values of α were given in Figures 4A,B verifying that the numerical approximate solution for α < 6 tends to the classic integer solution as α → 6. Future works can be focused on extending the suggested numerical technique to solve other types of BVPs.

All datasets generated for this study are included in the article/supplementary material.

All authors contributed equally to each part of this work. All authors read and approved the final manuscript.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.