A New Iterative Method for the Numerical Solution of High-Order Non-linear Fractional Boundary Value Problems

The boundary value problems (BVPs) have attracted the attention of many scientists from both practical and theoretical points of view, for these problems have remarkable applications in different branches of pure and applied sciences. Due to this important property, this research aims to develop an efficient numerical method for solving a class of non-linear fractional BVPs. The proposed method is free from perturbation, discretization, linearization, or restrictive assumptions, and provides the exact solution in the form of a uniformly convergent series. Moreover, the exact solution is determined by solving only a sequence of linear BVPs of fractional-order. Hence, from practical viewpoint, the suggested technique is efficient and easy to implement. To achieve an approximate solution with enough accuracy, we provide an iterative algorithm that is also computationally efficient. Finally, four illustrative examples are given verifying the superiority of the new technique compared to the other existing results.

1. Introduction

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.

2. Preliminaries

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

$\begin{array}{ll}{}_0{}^C{D}_t^{\alpha }(x(t))=\frac{1}{\Gamma (n-\alpha )}{\displaystyle {\int }_0^t{(t-\tau )}^{n-\alpha -1}}{x}^{(n)}(\tau )d\tau ,& (1)\end{array}$

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

$\begin{array}{ll}{}_0{}^C{I}_t^{\alpha }(x(t))=\frac{1}{\Gamma (\alpha )}{\displaystyle {\int }_0^t{(t-\tau )}^{\alpha -1}}x(\tau )d\tau .& (2)\end{array}$

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

$\begin{array}{ll}{}_0{}^C{D}_t^{\alpha }(a_1{x}_1(t)+a_2{x}_2(t))=a_1{}_0{}^C{D}_t^{\alpha }{x}_1(t)+a_2{}_0{}^C{D}_t^{\alpha }{x}_2(t).& (3)\end{array}$

Furthermore, the Caputo derivative of a constant function is zero, i.e., if x(t) ≡ k, then we have ${}_0{}^C{D}_t^{\alpha }k=0$. Additionally, the derivative and integral operators (1) and (2) satisfy the following anti-derivative property

$\begin{array}{ll}{}_0{}^C{I}_t^{\alpha }\left[{}_0{}^C{D}_t^{\alpha }(t)\right]=x(t)-x(0).& (4)\end{array}$

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

$\begin{array}{ll}\Vert {}_0{}^C{D}_t^{\alpha }{x}_1(t){-}_0^C{{D}}_t^{\alpha }{x}_2(t)\Vert \le L\Vert x_1(t)-x_2(t)\Vert ,& (5)\end{array}$

where L > 0 is the Lipschitz constant.

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

3. The Statement of the Problem

To formulate a fractional BVP, consider the following FDE

$\begin{array}{ll}{}_0{}^C{D}_t^{\alpha }(x(t))=f(x(t),t),n-1<\alpha \le n,t\in (0,T),& (6)\end{array}$

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

$\begin{array}{ll}{x}^{(2k)}(0)=a_{2k},x^{(2k)}(T)=b_{2k},k=0,1,\dots ,\frac{n}{2},& (7)\end{array}$

where a_2k, b_2k ($k=0,1,\dots ,\frac{n}{2}$) 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.

4. The Iterative Method

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,\dots ,\frac{n}{2}$.

Proof: Let x(·) be the solution of the BVP (6)-(7). Define ${\alpha }_i=x^{(i)}(0)$ and ${\beta }_j=x^{(j)}(T)$, i, j = 0, …, n − 1. Then x(·) is the solution of the following initial value problems (IVPs)

$\begin{array}{ll}\{\begin{array}{l}{}_0{}^C{D}_t^{\alpha }(x(t))=f(x(t),t),n-1<\alpha \le n,t\in (0,T),\hfill \\ x^{(i)}(0)={\alpha }_i,i=0,\dots ,n-1,\hfill \end{array}& (8)\end{array}$

$\begin{array}{ll}\{\begin{array}{l}{}_0{}^C{D}_t^{\alpha }(x(t))=f(x(t),t),n-1<\alpha \le n,t\in (0,T),\hfill \\ x^{(j)}(T)={\beta }_j,j=0,\dots ,n-1.\hfill \end{array}& (9)\end{array}$

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,\dots ,\frac{n}{2}$.

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)=\sum _{i=1}^{\infty }{\widehat{x}}_i(t)$, where ${\widehat{x}}_i(t)$ is attained by solving the sequence of linear fractional BVPs

$\begin{array}{ll}\{\begin{array}{l}{}_0{}^C{D}_t^{\alpha }({\widehat{x}}_1(t))={\lambda }_1(t){\widehat{x}}_1(t),\hfill \\ {\widehat{x}}_1^{(2k)}(0)=a_{2k},{\widehat{x}}_1^{(2k)}(T)=b_{2k},k=0,1,\dots ,\frac{n}{2},\hfill \end{array}& (10)\end{array}$

and for i = 2, 3, 4, …

$\begin{array}{ll}\{\begin{array}{l}{}_0{}^C{D}_t^{\alpha }({\widehat{x}}_i(t))={\lambda }_1(t){\widehat{x}}_i(t)+F_i(t,{\widehat{x}}_1(t),{\widehat{x}}_2(t),\dots ,{\widehat{x}}_{i-1}(t)),\hfill \\ {\widehat{x}}_i^{(2k)}(0)=0,{\widehat{x}}_i^{(2k)}(T)=0,k=0,1,\dots ,\frac{n}{2}.\hfill \end{array}& (11)\end{array}$

The non-homogeneous term F_i is determined by

$\begin{array}{l}{F}_i(t,{\widehat{x}}_1(t),{\widehat{x}}_2(t),\dots ,{\widehat{x}}_{i-1}(t))=\\ {\displaystyle \sum _{j=2}^i{\lambda }_j}(t){\displaystyle \sum _{k_1,\dots ,k_{i+1-j}}\frac{j!}{k_1!\cdots k_{i+1-j}!}}{\displaystyle \prod _{p=1}^{i+1-j}{\widehat{x}}_p^{k_p}}(t),& (12)\end{array}$

${\lambda }_j(t)=\frac{1}{j!}\frac{{\partial }^j}{\partial x^j}f(x,t){|}_{x=0}$, and the summation ${\displaystyle \sum _{k_1,\dots ,k_{i+1-j}}}$ is taken over all combinations of non-negative integer indices k₁ through k_i+1−j such that

$\begin{array}{ll}\{\begin{array}{l}{\displaystyle \sum _{p=1}^{i+1-j}{k}_p}=j,\hfill \\ {\displaystyle \sum _{p=1}^{i+1-j}p}{k}_p=i.\hfill \end{array}& (13)\end{array}$

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

$\begin{array}{ll}{}_0{}^C{D}_t^{\alpha }(x(t))={\lambda }_1(t)x(t)+{\lambda }_2(t)x^2(t)+{\lambda }_3(t)x^3(t)+\cdots ,& (14)\end{array}$

where ${\lambda }_j(t)=\frac{1}{j!}\frac{{\partial }^j}{\partial x^j}f(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

$\begin{array}{ll}x(t)=g(x_b,t),& (15)\end{array}$

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

$\begin{array}{l}x(t)=\underset{0}{\underbrace{{g(x_b,t)|}_{x_b=0}}}+\underset{{\widehat{x}}_1(t)}{\underbrace{{\frac{\partial }{\partial x_b}g(x_b,t)|}_{x_b=0}{x}_b}}\\ +\underset{{\widehat{x}}_2(t)}{\underbrace{x_b^T\left({\frac{1}{2!}\frac{{\partial }^2}{\partial x_b^2}g(x_b,t)|}_{x_b=0}\right)x_b}}+\cdots .& (16)\end{array}$

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

$\begin{array}{ll}x(t)=g(\epsilon x_b,t)=\epsilon {\widehat{x}}_1(t)+{\epsilon }^2{\widehat{x}}_2(t)+\cdots .& (17)\end{array}$

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

$\begin{array}{l}{}_0{}^C{D}_t^{\alpha }(\epsilon {\widehat{x}}_1(t)+{\epsilon }^2{\widehat{x}}_2(t)+\cdots )={\lambda }_1(t)\left(\epsilon {\widehat{x}}_1(t)+{\epsilon }^2{\widehat{x}}_2(t)+\cdots \right)\\ +{\lambda }_2(t){\left(\epsilon {\widehat{x}}_1(t)+{\epsilon }^2{\widehat{x}}_2(t)+\cdots \right)}^2+\cdots .& (18)\end{array}$

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

$\begin{array}{c}\epsilon {}_0{}^C{D}_t^{\alpha }({\widehat{x}}_1(t))+{\epsilon }^2{}_0{}^C{D}_t^{\alpha }({\widehat{x}}_2(t))+\cdots +{\epsilon }^i{}_0{}^C{D}_t^{\alpha }({\widehat{x}}_i(t))+\cdots =\\ \epsilon \left({\lambda }_1(t){\widehat{x}}_1(t)\right)+{\epsilon }^2\left({\lambda }_1(t){\widehat{x}}_2(t)+{\lambda }_2(t){\widehat{x}}_1^2(t)\right)+\cdots \\ +{\epsilon }^i\left({\lambda }_1(t){\widehat{x}}_i(t)+F_i(t,{\widehat{x}}_1(t),{\widehat{x}}_2(t),\dots ,{\widehat{x}}_{i-1}(t))\right)+\cdots ,& (19)\end{array}$

where

$\begin{array}{l}{F}_i(t,{\widehat{x}}_1(t),{\widehat{x}}_2(t),\dots ,{\widehat{x}}_{i-1}(t))=\\ {\displaystyle \sum _{j=2}^i{\lambda }_j}(t){\displaystyle \sum _{k_1,\dots ,k_{i+1-j}}\frac{j!}{k_1!\cdots k_{i+1-j}!}}{\displaystyle \prod _{p=1}^{i+1-j}{\widehat{x}}_p^{k_p}}(t),n& (20)\end{array}$

and the summation ${\displaystyle \sum _{k_1,\dots ,k_{i+1-j}}}$ is taken over all combinations of non-negative integer indices k₁ through k_i+1−j such that

$\begin{array}{ll}\{\begin{array}{l}{\displaystyle \sum _{p=1}^{i+1-j}{k}_p}=j,\hfill \\ {\displaystyle \sum _{p=1}^{i+1-j}p}{k}_p=i.\hfill \end{array}& (21)\end{array}$

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

$\begin{array}{ll}{\epsilon }^1:{}_0{}^C{D}_t^{\alpha }({\widehat{x}}_1(t))={\lambda }_1(t){\widehat{x}}_1(t),& (22)\end{array}$

$\begin{array}{ll}{\epsilon }^2:{}_0{}^C{D}_t^{\alpha }({\widehat{x}}_2(t))={\lambda }_1(t){\widehat{x}}_2(t)+{\lambda }_2(t){\widehat{x}}_1^2(t),& (23)\end{array}$

$\begin{array}{ll}\begin{array}{c}⋮\\ {\epsilon }^i:{}_0{}^C{D}_t^{\alpha }({\widehat{x}}_i(t))={\lambda }_1(t){\widehat{x}}_i(t)+F_i(t,{\widehat{x}}_1(t),{\widehat{x}}_2(t),\dots ,{\widehat{x}}_{i-1}(t)),\\ ⋮\end{array}& (24)\end{array}$

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

$\begin{array}{ll}{\epsilon }^1:{\widehat{x}}_1^{(2k)}(0)=a_{2k},{\widehat{x}}_1^{(2k)}(T)=b_{2k},k=0,1,\dots ,\frac{n}{2},& (25)\end{array}$

$\begin{array}{ll}{\epsilon }^i:{\widehat{x}}_i^{(2k)}(0)=0,{\widehat{x}}_i^{(2k)}(T)=0,k=0,1,\dots ,\frac{n}{2},i\ge 2,& (26)\end{array}$

and the proof is complete.

As can be seen, Equation (10) formulates a homogeneous linear BVP of fractional-order. By solving this problem, ${\widehat{x}}_1(t)$ is achieved in the first step. Following the proposed procedure in Theorem 4.1, we then obtain ${\widehat{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.

4.1. Approximate Solution

Although Theorem 4.1 suggests a closed-form expression for the solution of BVP (6)-(7), it is almost impossible to compute this solution in its present form since it is an infinite series. Hence, for the purpose of practical implementation, we need to truncate the series by considering its first M components where M is a positive integer number. Thus, the Mth-order approximate solution x_M(t) becomes

$\begin{array}{ll}{x}_M(t)={\displaystyle \sum _{i=1}^M}{\widehat{x}}_i(t).& (27)\end{array}$

To evaluate the value of M in Equation (27), the following criterion is considered according to the required accuracy. Indeed, the Mth-order approximate solution (27) has enough accuracy if for δ > 0, a given positive constant, the two consecutive solutions y_M−1(t) and y_M(t) satisfy

$\begin{array}{ll}{\Vert x_M(t)-x_{M-1}(t)\Vert }_{\infty }={\Vert {\widehat{x}}_M(t)\Vert }_{\infty }<\delta ,t\in (0,T).& (28)\end{array}$

Here, we present an iterative algorithm to design an approximate solution with enough accuracy.

Algorithm:

Step 1. Determine the first-order term ${\widehat{x}}_1(t)$ from Equation (10) and set i = 2.

Step 2. Determine the ith-order term ${\widehat{x}}_i(t)$ from Equation (11).

Step 3. Set M = i. By using the expression (27), compute x_M(t).

Step 4. If the condition (28) holds for a given small enough constant δ > 0, go to Step 5; else, replace i by i + 1 and go to Step 2.

Step 5. Consider x_M(t) as the appropriate approximate solution.

5. Numerical Simulations

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

$\begin{array}{ll}\{\begin{array}{l}{}_0{}^C{D}_t^{\alpha }x(t)=e^{-t}{x}^2(t),5<\alpha \le 6,t\in (0,1),\hfill \\ x^{(2k)}(0)=1,x^{(2k)}(1)=e,k=0,1,2,\hfill \end{array}& (29)\end{array}$

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., ${\Vert x_2(t)-x_1(t)\Vert }_{\infty }=2.2\times 10^{-3}<\delta$. 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.

TABLE 1

Table 1. The suggested technique at different iterations for Example 5.1.

FIGURE 1

Figure 1. Simulation curves of the exact solution and the second-order approximate solution for Example 5.1.

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.

TABLE 2

Table 2. Numerical comparison between the proposed iterative method and the other approximation techniques for Example 5.1.

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

$\begin{array}{ll}\{\begin{array}{l}{}_0{}^C{D}_t^{\alpha }x(t)=e^t{x}^2(t),5<\alpha \le 6,t\in (0,1),\hfill \\ x(0)=1,\dot{x}(0)=-1,\ddot{x}(0)=1,\hfill \\ x(1)=e^{-1},\dot{x}(1)=-e^{-1},\ddot{x}(1)=e^{-1},\hfill \end{array}& (30)\end{array}$

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., ${\Vert x_2(t)-x_1(t)\Vert }_{\infty }=1.3343\times 10^{-5}<\delta$. 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.

TABLE 3

Table 3. The suggested technique at different iterations for Example 5.2.

FIGURE 2

Figure 2. Simulation curves of the exact solution and the second-order approximate solution for Example 5.2.

FIGURE 3

Figure 3. Relation between the iteration time and the error for Examples 5.1 and 5.2.

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.

TABLE 4

Table 4. Numerical comparison between the proposed iterative method and the other approximation techniques for Example 5.2.

Example 5.3. Consider the following non-linear fractional BVP

$\begin{array}{ll}\{\begin{array}{l}{}_0{}^C{D}_t^{\alpha }x(t)=E_{\alpha }(-t^{\alpha })x^2(t),5<\alpha \le 6,t\in (0,1),\hfill \\ x^{(2k)}(0)={E_{\alpha }^{(2k)}(t^{\alpha })|}_{t=0},x^{(2k)}(1)={E_{\alpha }^{(2k)}(t^{\alpha })|}_{t=1},k=0,1,2,\hfill \end{array}& (31)\end{array}$

whose exact solution is $x(t)=E_{\alpha }(t^{\alpha })$ 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.

FIGURE 4

Figure 4. Simulation curves of the second-order approximate solution for α = 5.5, 5.6, 5.7, 5.8, 5.9, 6 [(A) Example 5.3 and (B) Example 5.4].

Example 5.4. Consider the non-linear fractional BVP

$\begin{array}{r}\{{}_0{}^C{D}_t^{\alpha }x(t)=E_{\alpha }(t^{\alpha })x^2(t),5<\alpha \le 6,t\in (0,1),\\ x^{(k)}(0)={E_{\alpha }^{(k)}(-t^{\alpha })|}_{t=0},x^{(k)}(1)={E_{\alpha }^{(k)}(-t^{\alpha })|}_{t=1},k=0,1,2,\text{}& (32)\end{array}$

whose exact solution is $x(t)=E_{\alpha }(-t^{\alpha })$.

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.

6. Conclusion

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.

Data Availability Statement

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

Author Contributions

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

Funding

This research was in part supported by a grant from University of Bojnord (No. 97/367/18077).

Conflict of Interest

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.

References

1. Toomore J, Zahn JP, Latour J, Spiegel EA. Stellar convection theory II: single-mode study of the second convection zone in A-type stars. Astrophys J. (1976) 207:545–63. doi: 10.1086/154522

CrossRef Full Text | Google Scholar

2. Glatzmaier GA. Numerical simulations of stellar convection dynamics III: at the base of the convection zone. Geophys Astrophys Fluid Dyn. (1985) 31:137–50. doi: 10.1080/03091928508219267

CrossRef Full Text | Google Scholar

3. Chandrasekhar S. Hydrodynamics and Hydromagnetic Stability. New York, NY: Dover (1981).

Google Scholar

4. Baldwin P. A localized instability in a Benard layer. Appl Anal. (1987) 24:1127–56. doi: 10.1080/00036818708839658

CrossRef Full Text | Google Scholar

5. Khalid A, Naeem MN, Ullah Z, Ghaffar A, Baleanu D, Nisar KS, et al. Numerical solution of the boundary value problems arising in magnetic fields and cylindrical shells. Mathematics. (2019) 7:508. doi: 10.3390/math7060508

CrossRef Full Text | Google Scholar

6. Akgul A, Akgul EK, Khan Y, Baleanu D. Solving the nonlinear system of third-order boundary value problems. In: Tas K, Baleanu D, Machado J, editors. Mathematical Methods in Engineering. Nonlinear Systems and Complexity. Cham: Springer (2019). p. 24. doi: 10.1007/978-3-319-90972-1_8

CrossRef Full Text | Google Scholar

7. Akgul A, Akgul EK, Baleanu D, Inc M. New numerical method for solving tenth order boundary value problems. Mathematics. (2018) 6:245. doi: 10.3390/math6110245

CrossRef Full Text | Google Scholar

8. Agarwal RP. Boundary Value Problems for Higher Order Differential Equations. Singapore: World Scientific (1986).

Google Scholar

9. Boutayeb A, Twizell EH. Numerical methods for the solution of special sixth-order boundary value problems. Int J Comput Math. (1992) 45:207–233. doi: 10.1080/00207169208804130

CrossRef Full Text | Google Scholar

10. Twizell EH, Boutayeb A. Numerical methods for the solution of special and general sixth-order boundary value problems, with applications to Benard layer eigen value problem. Proc R Soc A Math Phys Eng Sci. (1990) 431:433–50. doi: 10.1098/rspa.1990.0142

CrossRef Full Text | Google Scholar

11. Twizell EH. Numerical methods for sixth-order boundary value problems. In: Agarwal RP, Chow YM, Wilson SJ, editors. Numerical Mathematics Singapore 1988. International Series of Numerical Mathematics. Basel: Birkhauser (1988). p. 495–506. doi: 10.1007/978-3-0348-6303-2_40

CrossRef Full Text | Google Scholar

12. Al-Hayani W. Adomian decomposition method with Green's function for sixth-order boundary value problems. Comput Math Appl. (2011) 61:1567–75. doi: 10.1016/j.camwa.2011.01.025

CrossRef Full Text | Google Scholar

13. Noor MA, Mohyud-Din ST. Homotopy perturbation method for solving sixth order boundary value problems. Comput Math Appl. (2008) 55:2953–72. doi: 10.1016/j.camwa.2007.11.026

CrossRef Full Text | Google Scholar

14. Noor MA, Noor KI, Mohyud-Din ST. Variational iteration method for solving sixth-order boundary value problems. Commun Nonlinear Sci Num Simu. (2009) 14:2571–80. doi: 10.1016/j.cnsns.2008.10.013

CrossRef Full Text | Google Scholar

15. Kilbas AA, Srivastava HH, Trujillo JJ. Theory and Applications of Fractional Differential Equations. New York, NY: Elsevier (2006).

Google Scholar

16. Baleanu D, Asad JH, Jajarmi A. The fractional model of spring pendulum: new features within different kernels. Proc Roman Acad Ser A. (2018) 19:447–54.

Google Scholar

17. Baleanu D, Asad JH, Jajarmi A. New aspects of the motion of a particle in a circular cavity. Proc Roman Acad Ser A. (2018) 19:361–7.

Google Scholar

18. Baleanu D, Jajarmi A, Mohammadi H, Rezapour S. A new study on the mathematical modelling of human liver with Caputo-Fabrizio fractional derivative. Chaos Solitons Fractals. (2020) 134:109705. doi: 10.1016/j.chaos.2020.109705

CrossRef Full Text | Google Scholar

19. Jajarmi A, Yusuf A, Baleanu D, Inc M. A new fractional HRSV model and its optimal control: a non-singular operator approach. Phys A Stat Mech Appl. (2020) 547:123860. doi: 10.1016/j.physa.2019.123860

CrossRef Full Text | Google Scholar

20. Mohammadi F, Moradi L, Baleanu D, Jajarmi A. A hybrid functions numerical scheme for fractional optimal control problems: application to non-analytic dynamical systems. J Vibrat Control. (2018) 24:5030–43.

Google Scholar

21. Jajarmi A, Arshad S, Baleanu D. A new fractional modelling and control strategy for the outbreak of dengue fever. Phys A Stat Mech Appl. (2019) 535:122524. doi: 10.1016/j.physa.2019.122524

CrossRef Full Text | Google Scholar

22. Jajarmi A, Baleanu D, Sajjadi SS, Asad JH. A new feature of the fractional Euler-Lagrange equations for a coupled oscillator using a nonsingular operator approach. Front Phys. (2019) 7:196. doi: 10.3389/fphy.2019.00196

CrossRef Full Text | Google Scholar

23. Jajarmi A, Ghanbari B, Baleanu D. A new and efficient numerical method for the fractional modelling and optimal control of diabetes and tuberculosis co-existence. Chaos. (2019) 29:093111. doi: 10.1063/1.5112177

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Baleanu D, Jajarmi A, Sajjadi SS, Mozyrska D. A new fractional model and optimal control of a tumor-immune surveillance with non-singular derivative operator. Chaos. (2019) 29:083127. doi: 10.1063/1.5096159

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Baleanu D, Sajjadi SS, Jajarmi A, Asad JH. New features of the fractional Euler-Lagrange equations for a physical system within non-singular derivative operator. Eur Phys J Plus. (2019) 134:181. doi: 10.1140/epjp/i2019-12561-x

CrossRef Full Text | Google Scholar

26. Baleanu D, Jajarmi A, Sajjadi SS, Asad JH. The fractional features of a harmonic oscillator with position-dependent mass. Commun Theor Phys. (2020) 72:055002. doi: 10.1088/1572-9494/ab7700

CrossRef Full Text | Google Scholar

27. Yıldız TA, Jajarmi A, Yıldız B, Baleanu D. New aspects of time fractional optimal control problems within operators with nonsingular kernel. Discrete Continuous Dyn Syst S. (2020) 13:407–28. doi: 10.3934/dcdss.2020023

CrossRef Full Text | Google Scholar

28. Ali A, Shah K, Baleanu D. Ulam stability results to a class of nonlinear implicit boundary value problems of impulsive fractional differential equations. Adv Diff Equat. (2019) 2019:5. doi: 10.1186/s13662-018-1940-0

CrossRef Full Text | Google Scholar

29. Ugurlu E, Baleanu D, Tas K. On the solutions of a fractional boundary value problem. Turkish J Math. (2018) 42:1307–11. doi: 10.3906/mat-1609-64

CrossRef Full Text | Google Scholar

30. Baleanu D, Etemad S, Rezapour S. A hybrid Caputo fractional modeling for thermostat with hybrid boundary value conditions. Bound Value Probl. (2020) 2020:64. doi: 10.1186/s13661-020-01361-0

CrossRef Full Text | Google Scholar

31. Patnaik S, Sidhardh S, Semperlotti F. A Ritz-based finite element method for a fractional-order boundary value problem of nonlocal elasticity. arXiv preprint arXiv:200106885 (2020).

Google Scholar

32. Salem A, Alzahrani F, Alnegga M. Coupled system of nonlinear fractional Langevin equations with multipoint and nonlocal integral boundary conditions. Math Probl Eng. (2020) 2020:7345658. doi: 10.1155/2020/7345658

CrossRef Full Text | Google Scholar

33. Fazli H, Sun H, Aghchi S. Existence of extremal solutions of fractional Langevin equation involving nonlinear boundary conditions. Int J Comput Math. (2020) doi: 10.1080/00207160.2020.1720662

CrossRef Full Text | Google Scholar

34. Zhang S, Su X. The existence of a solution for a fractional differential equation with nonlinear boundary conditions considered using upper and lower solutions in reverse order. Comput Math Appl. (2011) 62:1269–1274. doi: 10.1016/j.camwa.2011.03.008

CrossRef Full Text | Google Scholar

35. Arqub OA, El-Ajou A, Zhour ZA, Momani S. Multiple solutions of nonlinear boundary value problems of fractional order: a new analytic iterative technique. Entropy. (2014) 16:471–93. doi: 10.3390/e16010471

CrossRef Full Text | Google Scholar

36. Khalil H, Al-Smadi M, Moaddy K, Khan RA, Hashim I. Toward the approximate solution for fractional order nonlinear mixed derivative and nonlocal boundary value problems. Discrete Dyn Nat Soc. (2016) 2016:5601821. doi: 10.1155/2016/5601821

CrossRef Full Text | Google Scholar

37. Cui Y, Sun Q, Su X. Monotone iterative technique for nonlinear boundary value problems of fractional order p ∈ (2, 3]. Adv Diff Equat. (2017) 2017:248. doi: 10.1186/s13662-017-1314-z

CrossRef Full Text | Google Scholar

38. Asaduzzaman M, Ali MZ. Existence of positive solution to the boundary value problems for coupled system of nonlinear fractional differential equations. AIMS Math. (2019) 4:880–95. doi: 10.3934/math.2019.3.880

CrossRef Full Text | Google Scholar

39. Diethelm K. (2010). Multi-term caputo fractional differential equations. In: The Analysis of Fractional Differential Equations. Lecture Notes in Mathematics, Vol. 2004. Berlin; Heidelberg: Springer (2010). doi: 10.1007/978-3-642-14574-2_8

CrossRef Full Text | Google Scholar