Mona Alsulami^*, Mariam Al-Mazmumy, Bashair Banweer
Department of Mathematics and Statistics, College of Science, University of Jeddah, Jeddah, Saudi Arabia.
DOI: 10.4236/am.2026.176020   PDF    HTML   XML   8 Downloads   49 Views   Citations

Abstract

The current paper numerically tackles the class of fractional ordinary boundary value problems using the modification method of the classical Adomian decomposition method and through a fractional inverse operator. Among the advantages of the adopted fractional inverse operator is that the operator directly incorporates all boundary conditions without introducing unknown constants, thereby making the approximate solution valid and preferred, with the satisfaction of the involved boundary conditions. Moreover, the results of the study show that the proposed method provides exact solutions for problems that have available exact solutions, and on the other hand, gives highly accurate approximate solutions in comparison with other computational techniques, confirming its efficiency and reliability. Indeed, several comparison tables and figures are provided, assessing the effectiveness of the proposed method over existing approaches in the literature.

Keywords

Adomian Decomposition Method, Fractional Boundary Value Problems, Fractional Inverse Operator, Dirichlet Boundary Conditions

Share and Cite:

1. Introduction

Non-integer order calculus is a generalization of classical calculus by extending the notion of integration and differentiation to fractional orders. This field, in recent decades has been of increasing concern owing to its efficiency in perfectly describing a wide range of applications, where classical calculus fails or gives half-backed information. Thus, applications of the new concept have emerged in a wide range of fields in engineering and science, including electromagnetism [1], control engineering [2] [3], bio-sciences [4], fluid dynamics and water waves [5], electrochemistry [6], thermodynamical processes [7], continuum and statistical mechanics [8], as well as in viscoelasticity [9]. Building on this foundation, fractional differential equations (FDEs) have emerged as a powerful tool in mathematical modeling. Their significance lies in their ability to accurately represent a variety of physical, scientific, and engineering processes, which makes them highly valuable in many practical and industrial applications, including physics, dynamics, mechanics, and other applied sciences [10]-[13]. In this study, we consider fractional ordinary differential equations with Dirichlet boundary conditions, commonly known as fractional boundary value problems (FBVPs) have been significantly explored in the recent times due to their relevance in modelling different phenomena in many fields such as the modeling of viscoelastic materials, wave propagation in porous media, image processing, various problems in fluid mechanics, control systems, and electrical circuits. Also applied in chemical and biological systems, FBVPs have been used in chemical graph models to represent molecular structures such as cyclopentasilane [14] and enable modeling of nonlocal and multi-dimensional interactions on molecular graphs [15] [16], demonstrating their practical significance and the need for effective analytical and numerical solution methods. To solve such problems, scientists and mathematicians developed several methods for obtaining both the approximate closed-form and numerical solutions, including the homotopy analysis method (HAM) [17], homotopy perturbation method (HPM) [18], variational iteration method (VIM) [19], the generalized differential transform method [20], shifted Legendre approximation method [21], shooting method [22], the Chebyshev spectral method [23], the collocation-shooting method [24], the monotone iterative sequences method [25], and the piecewise polynomial collocation approach [26]. Moreover, among these approaches, the Adomian decomposition method (ADM) and its variants [27]-[30], have attracted significant interest. Recently, considerable attention has been devoted to the ADM for solving a wide range of deterministic and stochastic differential equations [31]-[36]. The strength of the ADM lies in its ability to provide analytical approximations for broad classes of differential equations, without resorting to either of discretization, perturbation, linearization, or closure approximations, which often lead to extensive numerical computations.

However, the purpose of the present study is to focus on the ADM and its variants for treating the class of FBVPs for ordinary differential equations, by adapting a fractional inverse operator that directly incorporates all boundary conditions, eliminating the need to introduce unknown constants, in contrast to earlier modified ADM approaches. On the other hand, unlike in [37], where the inverse operator is employed to construct two separate formulations for the solution using the initial and boundary conditions independently, which are then combined through an averaging process to derive the recursive relations, the present work focuses only on the boundary conditions in constructing the recursive scheme. This yields a simpler decomposition procedure with fewer computational steps. The convergence of the Adomian decomposition series and the error analysis for fractional differential equations have been investigated in [38]-[41]. The significance of this work lies in its ability to provide exact solutions for certain problems while producing high-quality approximate solutions for such problems that have no known exact solutions. Additionally, the range of studies applying the classical Adomian method directly to such FBVPs remains limited, let alone the incorporation of the unexplored fractional inversion operator, highlighting the need for further research in this area. The current study is restricted to Caputo fractional ordinary boundary value problems with Dirichlet boundary conditions, and the illustrative examples mainly involve orders $0<\alpha \le 1$ and $1<\alpha \le 2$ . Thus, the present work is arranged in the following way: Section 2 presents some preliminaries on non-integer order calculus. Section 3 describes the ADM and its modified form suitable for solving FBVPs. Section 4 gives several illustrative examples, highlighting the efficiency and accuracy of the proposed method. Finally, Section 5 summarizes the main results and conclusions.

2. Fractional Derivatives and Fractional Integrals

This section presents some preliminaries on the non-integer order calculus, which includes the definitions of the Riemann-Liouville (RL) and Caputo fractional integrals and derivatives, and other definitions and properties; see the book by Kilbas et al. [42] for more on the theory of non-integer calculus.

Definition 2.1. Assume $\left[a,b\right]$ to be a finite interval in $ℝ$ . Then, the RL fractional integral for the function $f\left(x\right)$ of order $\alpha (>0)\in ℝ$ for $x>a,$ and $x<b$ are respectively defined as

$\left(I_{a+}^{\alpha }f\right)\left(x\right)=\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle {\int }_a^x{\left(x-t\right)}^{\alpha -1}f\left(t\right)\text{d}t},$ (1)

and

$\left(I_{b-}^{\alpha }f\right)\left(x\right)=\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle {\int }_x^b{\left(t-x\right)}^{\alpha -1}f\left(t\right)\text{d}t}.$ (2)

Definition 2.2. The RL fractional derivative for the function $f\left(x\right)$ of order $\alpha \left(\ge 0\right)\in ℝ$ for $x>a,$ and $x<b$ when $n=\left[\alpha \right]+1$ are respectively defined as

$\left(D_{a+}^{\alpha }f\right)\left(x\right)={\left(\frac{\text{d}}{\text{d}x}\right)}^n\left(I_{a+}^{n-\alpha }f\right)\left(x\right)=\frac{1}{\Gamma \left(n-\alpha \right)}{\left(\frac{\text{d}}{\text{d}x}\right)}^n{\displaystyle {\int }_a^x{\left(x-t\right)}^{n-\alpha -1}f\left(t\right)\text{d}t},$ (3)

and

$\left(D_{b-}^{\alpha }f\right)\left(x\right)={\left(-\frac{\text{d}}{\text{d}x}\right)}^n\left(I_{b-}^{n-\alpha }f\right)\left(x\right)=\frac{1}{\Gamma \left(n-\alpha \right)}{\left(-\frac{\text{d}}{\text{d}x}\right)}^n{\displaystyle {\int }_x^b{\left(t-x\right)}^{n-\alpha -1}f\left(t\right)\text{d}t},$ (4)

where $\left[\alpha \right]$ denotes the greatest integer value for the fractional-order $\alpha$ .

Definition 2.3. Assume $\left[a,b\right]$ is a finite interval in $ℝ$ . The fractional derivatives $\left({}^c{D}_{a+}^{\alpha }f\right)\left(x\right)$ and $\left({}^c{D}_{b-}^{\alpha }f\right)\left(x\right)$ of order $\alpha \in {ℝ}^{+}$ in $\left[a,b\right]$ via (3) and (4) are the left-sided, and right-sided Caputo derivatives of order $\alpha$ , are re-expressed respectively defined as

$\left({}^c{D}_{a+}^{\alpha }f\right)\left(x\right)=\left(D_{a+}^{\alpha }\left[f\left(t\right)-{\displaystyle \sum }_{j=0}^{n-1}\frac{f^{\left(j\right)}\left(a\right)}{j!}{\left(t-a\right)}^j\right]\right)\left(x\right),$ (5)

and

$\left({}^c{D}_{b-}^{\alpha }f\right)\left(x\right)=\left(D_{b-}^{\alpha }\left[f\left(t\right)-{\displaystyle \sum }_{j=0}^{n-1}\frac{f^{\left(j\right)}\left(b\right)}{j!}{\left(b-t\right)}^j\right]\right)\left(x\right),$ (6)

where $n=\left[\alpha \right]+1$ for $\alpha \notin {\mathbb{N}}_0;n=\alpha$ for $\alpha \in {\mathbb{N}}_0$ .

Lemma 2.1. Assume $\alpha \in {ℝ}^{+},$ and $n=\left[\alpha \right]+1$ for $\alpha \notin {\mathbb{N}}_0;n=\alpha$ for $\alpha \in {\mathbb{N}}_0$ . Then, if $f\left(x\right)\in C^n\left[a,b\right]$ or $f\left(x\right)\in A{C}^n\left[a,b\right]$ , one gets

$\left(I_{a+}^{\alpha }{}^c{D}_{a+}^{\alpha }f\right)\left(x\right)=f\left(x\right)-{\displaystyle \sum }_{j=0}^{n-1}\frac{f^{\left(j\right)}\left(a\right)}{j!}{\left(x-a\right)}^j,$ (7)

$\left(I_{b-}^{\alpha }{}^c{D}_{b-}^{\alpha }f\right)\left(x\right)=f\left(x\right)-{\displaystyle \sum }_{j=0}^{n-1}\frac{{\left(-1\right)}^j{f}^{\left(j\right)}\left(b\right)}{j!}{\left(b-x\right)}^j.$ (8)

Theorem 2.1. [43] For $\alpha \in {ℝ}^{+}$ $\left(n-1<\alpha \le n\right)$ , and $f\left(x\right)\in C_{-1}^n,n\in N$ , then

$\left(I^{\alpha }{D}^{\alpha }f\right)\left(x\right)=f\left(x\right)-{\displaystyle \sum }_{k=0}^{n-1}{f}^{\left(j\right)}\left(0^{+}\right)\frac{x^j}{j!}.$ (9)

Property 2.1. For the power function $f\left(x\right)=x^m$ , $m>-1$ , $x>0$ and $\alpha >0$ , one gets

$I^{\alpha }{x}^m=\frac{\Gamma \left(1+m\right)}{\Gamma \left(1+m+\alpha \right)}{x}^{m+\alpha }.$ (10)

Property 2.2. The Caputo fractional derivative of the function $f\left(x\right)=x^m$ for $m>-1$ and $n-1<\alpha \le n$ satisfies

${}^c{D}_{a+}^{\alpha }{x}^m=\{\begin{array}{ll}0,\hfill & \text{for}m<\lceil \alpha \rceil ,\hfill \\ \frac{\Gamma \left(1+m\right)}{\Gamma \left(1+m-\alpha \right)}{x}^{m-\alpha },\hfill & \text{for}m\ge \lceil \alpha \rceil .\hfill \end{array}$ (11)

3. Description of the Proposed Method

In proposing the new method for solving the class of FBVPs, the following nonlinear FDE is considered

$D^{\alpha }u\left(x\right)+Ru\left(x\right)+Nu\left(x\right)=g\left(x\right),a<x<b,$ (12)

subject to the following Dirichlet boundary conditions

$u\left(a\right)={\gamma }_1,u\left(b\right)={\gamma }_2,$ (13)

where $D^{\alpha }$ is the fractional derivative operator in the sense of Caputo, $R$ is a linear differential operator, which may equally feature other fractional operators of order less than $\alpha$ . At the same time, $N$ is an operator that is nonlinear, and $g\left(x\right)$ is the nonhomogeneous term, also referred to as the source function.

What is more, the governing fractional integral operator $I^{\alpha }$ of much concern in this study has the form [37]

$\begin{array}{c}{I}^{\alpha }(.)={\displaystyle {\int }_a^x\frac{{\left(x-s\right)}^{\alpha -1}}{\Gamma \left(\alpha \right)}(.)\text{d}s}-\left(\frac{x-a}{b-a}\right){\displaystyle {\int }_a^b\frac{{\left(b-s\right)}^{\alpha -1}}{\Gamma \left(\alpha \right)}(.)\text{d}s},\\ =I^{\alpha }(.)-\left(\frac{x-a}{b-a}\right)I^{\alpha }{(.)}_{x=b}.\end{array}$ (14)

Accordingly, upon applying the inverse operator (14) to both sides of (12), one gets

$I^{\alpha }\left[D^{\alpha }u\left(x\right)\right]+I^{\alpha }\left[Ru\left(x\right)\right]+I^{\alpha }\left[Nu\left(x\right)\right]=I^{\alpha }\left[g\left(x\right)\right],$

such that upon using Lemma 2.1 and the boundary conditions (13) the latter equation becomes

$\begin{array}{l}u\left(x\right)-{\displaystyle \sum }_{k=0}^{n-1}{u}^{\left(k\right)}\left(a\right)\frac{{\left(x-a\right)}^k}{k!}-\left(\frac{x-a}{b-a}\right)\left[u\left(b\right)-{\displaystyle \sum }_{k=0}^{n-1}{u}^{\left(k\right)}\left(a\right)\frac{{\left(b-a\right)}^k}{k!}\right]\\ +I^{\alpha }\left[Ru\left(x\right)\right]+I^{\alpha }\left[Nu\left(x\right)\right]=I^{\alpha }\left[g\left(x\right)\right],\end{array}$

or equally,

$u\left(x\right)=u\left(a\right)+\left(\frac{x-a}{b-a}\right)\left[u\left(b\right)-u\left(a\right)\right]-I^{\alpha }\left[Ru\left(x\right)\right]-I^{\alpha }\left[Nu\left(x\right)\right]+I^{\alpha }\left[g\left(x\right)\right].$ (15)

Next, expressing the solution $u\left(x\right)$ using an infinite series of the form:

$u\left(x\right)={\displaystyle \sum }_{n=0}^{\infty }{u}_n\left(x\right),$ (16)

and expanding the nonlinear term $Nu\left(x\right)$ via the use of Adomian polynomials [44]:

$N\left(u\right)={\displaystyle \sum }_{n=0}^{\infty }{A}_n\left(u_0,u_1,\cdots ,u_n\right),$ (17)

where the polynomials are plainly determined using the relation

$A_n=\frac{1}{n!}\frac{{\text{d}}^n}{\text{d}{\lambda }^n}{\left[N\left({\displaystyle \sum }_{i=0}^n{\lambda }^i{u}_i\left(x\right)\right)\right]}_{\lambda =0},n=0,1,2,\cdots$ (18)

one thus substitutes the latter series expressions in (16) and (17) into (15) to get

$\begin{array}{c}{\displaystyle \sum }_{n=0}^{\infty }{u}_n\left(x\right)=u\left(a\right)+\left(\frac{x-a}{b-a}\right)\left[u\left(b\right)-u\left(a\right)\right]+I^{\alpha }\left[g\left(x\right)\right]\\ -I^{\alpha }R\left[{\displaystyle \sum }_{n=0}^{\infty }{u}_n\left(x\right)\right]-I^{\alpha }\left[{\displaystyle \sum }_{n=0}^{\infty }{A}_n\right],\end{array}$

which thus leads to the resulting recursive scheme as follows

$\{\begin{array}{l}{u}_0\left(x\right)=h\left(x\right),\\ u_{n+1}\left(x\right)=-I^{\alpha }\left[R{u}_n\left(x\right)\right]-I^{\alpha }\left[A_n\right],n\ge 0.\end{array}$ (19)

where

$h\left(x\right)=u\left(a\right)+\left(\frac{x-a}{b-a}\right)\left[u\left(b\right)-u\left(a\right)\right]+I^{\alpha }\left[g\left(x\right)\right].$

Notably, on spitting the latter expression for $h$ , that is, letting $h=h_1+h_2$ , the reliable modification of the ADM [29] thus recasts the recursive scheme in (19) to the following

$\{\begin{array}{l}{u}_0\left(x\right)=h_1,\\ u_1\left(x\right)=h_2-I^{\alpha }\left[R{u}_0\left(x\right)\right]-I^{\alpha }\left[A_0\right],\\ u_{n+1}\left(x\right)=-I^{\alpha }\left[R{u}_n\left(x\right)\right]-I^{\alpha }\left[A_n\right],n\ge 1.\end{array}$ (20)

Hence, the overall solution is obtained via the infinite series in (16) by substituting the determined components $u_0,u_1,u_2,\cdots$ . In fact, truncation yields an approximate solution, while summation - when possible - leads to a closed form.

4. Applications

The present section shows the efficiency of the proposed new decomposition method through the application of an inverse fractional operator. Several illustrative examples will be extensively examined, portraying Dirichlet boundary conditions and endowed with the Caputo fractional derivative. For each example, the obtained numerical solution will be benchmarked against the exact solution, if available, by computing the consequent absolute error. Moreover, a comparative analysis will be established with the existing numerical methods available in the literature.

Example 4.1. Consider the FBVP through the Bagley-Torvik equation

$D^{3/2}u\left(x\right)+u^{″}\left(x\right)+u\left(x\right)=1+x,x\in \left[0,1\right],$ (21)

$u\left(0\right)=1,u\left(1\right)=2,$ (22)

that satisfies the exact solution $u\left(x\right)=1+x$ .

First, with the help of an operator notation, we write (21) as follows

$D^{3/2}u\left(x\right)=1+x-u^{″}\left(x\right)-u\left(x\right).$ (23)

Then, with the application of the inversion operator (14) on the both side of (23) using Theorem 2.1 and boundary conditions (22), one gets

$u\left(x\right)=1+x+I^{3/2}\left[1+x\right]-I^{3/2}\left[u^{″}\left(x\right)+u\left(x\right)\right],$

where, from the latter equation, the proposed method yields the recursive scheme as follows

$\{\begin{array}{l}{u}_0\left(x\right)=1+x,\hfill \\ u_1\left(x\right)=I^{3/2}\left[1+x\right]-I^{3/2}\left[{u^{″}}_0\left(x\right)+u_0\left(x\right)\right],\hfill \\ u_{n+1}\left(x\right)=I^{3/2}\left[{u^{″}}_n\left(x\right)+u_n\left(x\right)\right],n\ge 1.\hfill \end{array}$

In particular, from the above scheme, computing the value of $u_1\left(x\right)$ by using $I^{3/2}$ from (14), one gets

$u_1\left(x\right)=I^{3/2}\left[1+x\right]-I^{3/2}\left[{u^{″}}_0\left(x\right)+u_0\left(x\right)\right]=I^{3/2}\left[1+x\right]-I^{3/2}\left[1+x\right]=0.$

Hence, $u_n\left(x\right)=0$ $\forall n\ge 1$ . Accordingly, the closed-form solution is obtained as $u\left(x\right)=u_0\left(x\right)=1+x$ , which coincides with the reported exact solution. Moreover, Table 1 highlights the superior accuracy of the present method in comparison with the Bessel collocation method (BCM) presented in [45]. As shown, the present method reaches the exact solution with two components, demonstrating its high precision and rapid convergence.

Table 1. Comparison of absolute error for Example 4.1.

$x$	BCM [45] ($N=6$ )	Present method ($N=2$ )
0.0	0	0
0.1	$6.1919\times {10}^{-16}$	0
0.2	$1.0292\times {10}^{-15}$	0
0.3	$1.4779\times {10}^{-15}$	0
0.4	$1.9697\times {10}^{-15}$	0
0.5	$2.4941\times {10}^{-15}$	0
0.6	$3.0365\times {10}^{-15}$	0
0.7	$3.5805\times {10}^{-15}$	0
0.8	$4.1090\times {10}^{-15}$	0
0.9	$4.6047\times {10}^{-8}$	0
1.0	$5.0500\times {10}^{-15}$	0

Example 4.2. Consider the FBVP

$D^{3/2}u\left(x\right)+u\left(x\right)=\frac{2x^{1/2}}{\Gamma \left(3/2\right)}+x^2-x,x\in \left[0,1\right],$ (24)

$u\left(0\right)=0,u\left(1\right)=0,$ (25)

which satisfies the exact solution $u\left(x\right)=x^2-x$ .

Now, one begins by re-writing (24) as follows

$D^{3/2}u\left(x\right)=\frac{2x^{1/2}}{\Gamma \left(3/2\right)}+x^2-x-u\left(x\right).$ (26)

Next, with the application of the inversion operator $I^{3/2}$ from (14) on (26), one gets after using the property in (7) and the conditions (25) the following

$u\left(x\right)=I^{3/2}\left[\frac{2x^{1/2}}{\Gamma \left(3/2\right)}+x^2-x\right]-I^{3/2}\left[u\left(x\right)\right],$

upon which the application of (14) and the property (10) subsequently yields

$u\left(x\right)=x^2+\frac{2}{\Gamma \left(9/2\right)}{x}^{7/2}-\frac{1}{\Gamma \left(7/2\right)}{x}^{5/2}-x-\left[\frac{2}{\Gamma \left(9/2\right)}-\frac{1}{\Gamma \left(7/2\right)}\right]x-I^{3/2}\left[u\left(x\right)\right].$

Consequently, employing the proposed decomposition method reveals the overall recursive scheme as follows

$\{\begin{array}{l}{u}_0\left(x\right)=x^2-x,\\ u_1\left(x\right)=\frac{2}{\Gamma \left(9/2\right)}{x}^{7/2}-\frac{1}{\Gamma \left(7/2\right)}{x}^{5/2}-\left[\frac{2}{\Gamma \left(9/2\right)}-\frac{1}{\Gamma \left(7/2\right)}\right]x-I^{3/2}\left[u_0\left(x\right)\right],\\ u_{n+1}\left(x\right)=-I^{3/2}\left[u_n\left(x\right)\right],n\ge 1.\end{array}$

Computing the value of $u_1\left(x\right)$ by using $I^{3/2}$ from (14) and the property in (10), one gets

$\begin{array}{c}{u}_1\left(x\right)=\frac{2}{\Gamma \left(9/2\right)}{x}^{7/2}-\frac{1}{\Gamma \left(7/2\right)}{x}^{5/2}-\left[\frac{2}{\Gamma \left(9/2\right)}-\frac{1}{\Gamma \left(7/2\right)}\right]x-I^{3/2}\left[x^2-x\right],\\ =\frac{2}{\Gamma \left(9/2\right)}{x}^{7/2}-\frac{1}{\Gamma \left(7/2\right)}{x}^{5/2}-\left[\frac{2}{\Gamma \left(9/2\right)}-\frac{1}{\Gamma \left(7/2\right)}\right]x-\frac{2}{\Gamma \left(9/2\right)}{x}^{7/2}\\ +\frac{1}{\Gamma \left(7/2\right)}{x}^{5/2}+\left[\frac{2}{\Gamma \left(9/2\right)}-\frac{1}{\Gamma \left(7/2\right)}\right]x,\\ =0.\end{array}$

Hence, $u_n\left(x\right)=0$ $\forall n\ge 1$ . Moreover, one gets the solution by summing the components as $u\left(x\right)=u_0\left(x\right)=x^2-x$ , coinciding with the already reported exact solution. What is more, Table 2 shows that the present technique with two components produces zero absolute error for all values of $x$ , demonstrating that the method attains the exact solution. Compared with Chelyshkov operational matrix method (COMM) [46], Pseudo spectral method (PSM) [47], and Generalized Adams-Bashforth-Moulton method (GABMM) [48], which show small numerical errors. Indeed, the present method exhibits superior accuracy and fast convergence. For GABMM, $h$ denotes the step size used in the simulation.

Table 2. Comparison of absolute error for Example 4.2.

$x$	COMM [46] ($N=10$ )	PSM [47] ($N=10$ )	GABMM [48] ($N=6$ ) $h=1/20$	GABMM [48] ($N=6$ ) $h=1/320$	Present method ($N=2$ )
0.0	$6.204\times {10}^{-5}$	-	-	-	0
0.1	$9.450\times {10}^{-6}$	-	-	-	0
0.2	$8.600\times {10}^{-7}$	-	-	-	0
0.3	$1.261\times {10}^{-5}$	-	-	-	0
0.4	$1.273\times {10}^{-5}$	-	-	-	0
0.5	$5.170\times {10}^{-6}$	-	-	-	0
0.6	$1.224\times {10}^{-6}$	-	-	-	0
0.7	$6.27\times {10}^{-6}$	-	-	-	0
0.8	$8.13\times {10}^{-6}$	-	-	-	0
0.9	$9.30\times {10}^{-6}$	-	-	-	0
1.0	$2.379\times {10}^{-5}$	$6.83\times {10}^{-4}$	$3.42\times {10}^{-3}$	$5.71\times {10}^{-5}$	0

Example 4.3. Consider the FBVP

$D^{2\alpha }u\left(x\right)=u\left(x\right)+x,0<\alpha \le 1,x\in \left[0,1\right],$ (27)

$u\left(0\right)=0,u\left(1\right)=0,$ (28)

that satisfies the exact solution $u\left(x\right)=\frac{\text{sinh}\left(x\right)}{\text{sinh}\left(1\right)}-x$ while assuming the integer order unity.

In the same way, applying the inversion operator $I^{2\alpha }$ from (14) to (27), and upon using properties in (7) and (10) through the conditions (28), one gets

$u\left(x\right)=\frac{1}{\Gamma \left(2+2\alpha \right)}{x}^{1+2\alpha }-\frac{1}{\Gamma \left(2+2\alpha \right)}x+I^{2\alpha }\left[u\left(x\right)\right],$

that yields the consequent ADM recursive scheme as follows

$\{\begin{array}{l}{u}_0\left(x\right)=\frac{1}{\Gamma \left(2+2\alpha \right)}{x}^{1+2\alpha }-\frac{1}{\Gamma \left(2+2\alpha \right)}x,\hfill \\ u_{n+1}\left(x\right)=I^{2\alpha }\left[u_n\left(x\right)\right],n\ge 0.\hfill \end{array}$

Subsequently, one computes some of the successive iterates as

$\begin{array}{c}{u}_1\left(x\right)=I^{2\alpha }\left[u_0\left(x\right)\right]\\ =\frac{1}{\Gamma \left(2+4\alpha \right)}{x}^{1+4\alpha }-\frac{1}{{\left[\Gamma \left(2+2\alpha \right)\right]}^2}{x}^{1+2\alpha }-\left[\frac{1}{\Gamma \left(2+4\alpha \right)}-\frac{1}{{\left[\Gamma \left(2+2\alpha \right)\right]}^2}\right]x,\end{array}$

$\begin{array}{c}{u}_2\left(x\right)=I^{2\alpha }\left[u_1\left(x\right)\right]\\ =\frac{1}{\Gamma \left(2+6\alpha \right)}{x}^{1+6\alpha }-\frac{1}{\Gamma \left(2+2\alpha \right)\Gamma \left(2+4\alpha \right)}{x}^{1+4\alpha }\\ -\left[\frac{1}{\Gamma \left(2+2\alpha \right)\Gamma \left(2+4\alpha \right)}-\frac{1}{{\left[\Gamma \left(2+2\alpha \right)\right]}^3}\right]x^{1+2\alpha }\\ -\left[\frac{1}{\Gamma \left(2+6\alpha \right)}-\frac{2}{\Gamma \left(2+2\alpha \right)\Gamma \left(2+4\alpha \right)}+\frac{1}{{\left[\Gamma \left(2+2\alpha \right)\right]}^3}\right]x,\end{array}$

$\begin{array}{c}{u}_3\left(x\right)=I^{2\alpha }\left[u_2\left(x\right)\right]\\ =\frac{1}{\Gamma \left(2+8\alpha \right)}{x}^{1+8\alpha }-\frac{1}{\Gamma \left(2+2\alpha \right)\Gamma \left(2+8\alpha \right)}{x}^{1+6\alpha }\\ -\left[\frac{1}{{\left[\Gamma \left(2+4\alpha \right)\right]}^2}-\frac{1}{{\left[\Gamma \left(2+2\alpha \right)\right]}^2\Gamma \left(2+4\alpha \right)}\right]x^{1+4\alpha }\end{array}$

$\begin{array}{c}-\left[\frac{1}{\Gamma \left(2+6\alpha \right)\Gamma \left(2+2\alpha \right)}-\frac{2}{{\left[\Gamma \left(2+2\alpha \right)\right]}^2\Gamma \left(2+4\alpha \right)}+\frac{1}{{\left[\Gamma \left(2+2\alpha \right)\right]}^4}\right]x^{1+2\alpha }\\ -[\frac{1}{\Gamma \left(2+8\alpha \right)}-\frac{1}{\Gamma \left(2+2\alpha \right)\Gamma \left(2+8\alpha \right)}-\frac{1}{{\left[\Gamma \left(2+4\alpha \right)\right]}^2}\\ +\frac{3}{{\left[\Gamma \left(2+2\alpha \right)\right]}^2\Gamma \left(2+4\alpha \right)}-\frac{1}{\Gamma \left(2+6\alpha \right)\Gamma \left(2+2\alpha \right)}+\frac{1}{{\left[\Gamma \left(2+2\alpha \right)\right]}^4}]x,\end{array}$

$⋮$

Accordingly, the approximate solution follows from the summation in (16). What is more, Figure 1 illustrates the nature of the solution for different fractional orders $\alpha$ . It is observed that when $\alpha =1$ , the reported exact solution is attained. In addition, Table 3 presents the absolute error values, where the present method consistently yields smaller errors than the method in [49], proving its superior accuracy and efficiency in solving FBVPs.

Figure 1. Exact and ADM solutions of Example 4.3 for various values of $\alpha =0.6,0.8,1$ .

Table 3. Comparison of absolute error for Example 4.3.

$x$	FCDM [49]			Present method
	$\alpha =0.6$	$\alpha =0.8$	$\alpha =1$	$\alpha =0.6$	$\alpha =0.8$	$\alpha =1$
0.0	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
0.1	$2.810\times {10}^{-2}$	$9.600\times {10}^{-3}$	$2.000\times {10}^{-4}$	$1.813\times {10}^{-2}$	$7.899\times {10}^{-3}$	$1.932\times {10}^{-7}$
0.2	$4.880\times {10}^{-2}$	$1.730\times {10}^{-2}$	$1.200\times {10}^{-3}$	$3.248\times {10}^{-2}$	$1.452\times {10}^{-2}$	$3.674\times {10}^{-7}$
0.3	$6.170\times {10}^{-2}$	$2.260\times {10}^{-2}$	$3.000\times {10}^{-4}$	$4.281\times {10}^{-2}$	$1.949\times {10}^{-2}$	$5.058\times {10}^{-7}$
0.4	$6.750\times {10}^{-2}$	$2.530\times {10}^{-2}$	$3.000\times {10}^{-4}$	$4.907\times {10}^{-2}$	$2.267\times {10}^{-2}$	$5.946\times {10}^{-7}$
0.5	$6.670\times {10}^{-2}$	$2.530\times {10}^{-2}$	$2.000\times {10}^{-4}$	$5.128\times {10}^{-2}$	$2.394\times {10}^{-2}$	$6.253\times {10}^{-7}$
0.6	$6.020\times {10}^{-2}$	$2.300\times {10}^{-2}$	$6.000\times {10}^{-4}$	$4.940\times {10}^{-2}$	$2.325\times {10}^{-2}$	$5.947\times {10}^{-7}$
0.7	$4.910\times {10}^{-2}$	$1.870\times {10}^{-2}$	$5.115\times {10}^{-1}$	$4.341\times {10}^{-2}$	$2.056\times {10}^{-2}$	$5.059\times {10}^{-7}$
0.8	$3.450\times {10}^{-2}$	$1.290\times {10}^{-2}$	$3.300\times {10}^{-4}$	$3.325\times {10}^{-2}$	$1.583\times {10}^{-2}$	$3.674\times {10}^{-7}$
0.9	$1.770\times {10}^{-2}$	$6.400\times {10}^{-3}$	$4.000\times {10}^{-4}$	$1.883\times {10}^{-2}$	$8.992\times {10}^{-3}$	$1.933\times {10}^{-7}$
1.0	0.0000	0.0000	0.0000	0.0000	$9.738\times {10}^{-15}$	0.0000

Example 4.4. Consider the FBVP for the nonlinear Bratu equation

$D^{\alpha }u\left(x\right)+\mu {\text{e}}^{u\left(x\right)}=0,1<\alpha \le 2,x\in \left[0,1\right],$ (29)

$u\left(0\right)=0,u\left(1\right)=0.$ (30)

The Bratu problem, introduced by Bratu in 1914 [50] has vast applications, including the fuel ignition [50], radiative heat transfer [51], and electrospinning of polymer nanofibers [52].

In this regard, the exact solution for the governing Bratu model when $\alpha =2$ takes the expression

$u\left(x\right)=-2\text{ln}\left[\frac{\text{cosh}\left(\left(x-\frac{1}{2}\right)\frac{\theta }{2}\right)}{\text{cosh}\left(\frac{\theta }{4}\right)}\right],$ (31)

with $\theta$ satisfying $\theta =\sqrt{2\mu }\cosh \left(\frac{\theta }{4}\right)$ . Moreover, with the discovery of the critical value ${\mu }_c$ as ${\mu }_c=3.513830719$ , and satisfying $1=\frac{1}{4}\sqrt{2{\mu }_c}\sinh \left(\frac{{\theta }_c}{4}\right)$ , $\mu$ is thus classified as:

1) When $\mu <{\mu }_c$ , the Bratu equation has two solutions.

2) When $\mu ={\mu }_c$ , the Bratu equation has one solution.

3) When $\mu >{\mu }_c$ , the Bratu equation has a zero solution.

Accordingly, to solve the Bratu problem in (29)-(30), one re-writes the equation using the operator notation when $\mu =1$ as $D^{\alpha }u\left(x\right)=-{\text{e}}^{u\left(x\right)}$ . Next, applying the inversion operator $I^{\alpha }$ in (??) to the equation and via the property (7) and the conditions (30) to get

$u\left(x\right)=-I^{\alpha }\left[{\text{e}}^{u\left(x\right)}\right].$

We decompose the nonlinear term $Nu={\text{e}}^u$ by the Adomian polynomials from (17)-(18) as follows:

$A_0={\text{e}}^{u_0},$

$A_1={\text{e}}^{u_0}{u}_1,$

$A_2=\frac{1}{2}{\text{e}}^{u_0}\left(u_1^2+2u_2\right),$

$⋮$

Hence, through the classical ADM procedure, the recursive scheme is obtained as follows

$\{\begin{array}{l}{u}_0\left(x\right)=0,\hfill \\ u_{n+1}\left(x\right)=-I^{\alpha }\left[A_n\right],n\ge 0.\hfill \end{array}$

Consequently, one gets some of the solution components from the latter scheme as

$u_1\left(x\right)=-I^{\alpha }\left[{\text{e}}^{u_0}\right]=\frac{1}{\Gamma \left(1+\alpha \right)}\left(x-x^{\alpha }\right),$

$u_2\left(x\right)=-I^{\alpha }\left[{\text{e}}^{u_0}{u}_1\right]=\frac{1}{\Gamma \left(1+2\alpha \right)}\left(x^{2\alpha }-x\right)+\frac{1}{\Gamma \left(1+\alpha \right)\Gamma \left(2+\alpha \right)}\left(x-x^{1+\alpha }\right),$

$\begin{array}{c}{u}_3\left(x\right)=-I^{\alpha }\left[\frac{1}{2}{\text{e}}^{u_0}\left(u_1^2+2u_2\right)\right]\\ =-\left[\frac{\Gamma \left(1+2\alpha \right)}{2{\left[\Gamma \left(1+\alpha \right)\right]}^2}+1\right]\frac{x^{3\alpha }}{\Gamma \left(1+3\alpha \right)}+\left[\frac{\Gamma \left(2+\alpha \right)}{\Gamma \left(1+\alpha \right)}+1\right]\frac{x^{1+2\alpha }}{\Gamma \left(1+\alpha \right)\Gamma \left(2+2\alpha \right)}\\ +\left[\frac{1}{\Gamma \left(1+2\alpha \right)}-\frac{1}{\Gamma \left(1+\alpha \right)\Gamma \left(2+\alpha \right)}\right]\frac{x^{1+\alpha }}{\Gamma \left(2+\alpha \right)}\\ -\left[\frac{1}{{\left[\Gamma \left(1+\alpha \right)\right]}^2\Gamma \left(3+\alpha \right)}\right]x^{2+\alpha }\\ -[\frac{2{\left[\Gamma \left(1+\alpha \right)\right]}^2+\Gamma \left(1+3\alpha \right)}{2{\left[\Gamma \left(1+\alpha \right)\right]}^2\Gamma \left(1+3\alpha \right)}+\frac{\Gamma \left(2+\alpha \right)+\Gamma \left(1+\alpha \right)}{{\left[\Gamma \left(1+\alpha \right)\right]}^2\Gamma \left(2+2\alpha \right)}\\ +\frac{\Gamma \left(1+2\alpha \right)-\Gamma \left(2+\alpha \right)\Gamma \left(1+\alpha \right)}{{\left[\Gamma \left(2+\alpha \right)\right]}^2\Gamma \left(1+2\alpha \right)\Gamma \left(1+\alpha \right)}+\frac{1}{{\left[\Gamma \left(1+\alpha \right)\right]}^2\Gamma \left(3+\alpha \right)}]x,\end{array}$

$⋮$

Moreover, one gets the approximate solution from (16); see Figure 2 that gives the graphical illustration of the solution while responding to variation of the fractional order $\alpha$ when fixing $\mu =1$ . The solution curves for the different values of $\alpha$ are very close to each other, with almost perfect overlap when $\alpha =2$ . Additionally, Table4 reports absolute errors comparison, showing that the present method attains higher accuracy for all values of $\alpha$ using only 4 components, while the reproducing kernel method (RKM) [53] employs 20 components, yet yields less accurate results. This clearly reaffirms the superiority of the proposed method over the contending method for solving nonlinear FBVPs.

Figure 2. Exact and ADM solutions of Example 4.4 for various values of $\alpha =1.8,1.9,2$ .

Table 4. Comparison of absolute error for Example 4.4 when $\mu =1$ .

$x$	RKM [53] $\left(N=20\right)$			Present method $\left(N=4\right)$
	Error			Error
	$\alpha =1.8$	$\alpha =1.9$	$\alpha =2$	$\alpha =1.8$	$\alpha =1.9$	$\alpha =2$
0.0	-	-	-	0.0000	0.0000	0.0000
0.1	$6.3195\times {10}^{-3}$	$3.1488\times {10}^{-3}$	$2.0159\times {10}^{-4}$	$6.3662\times {10}^{-3}$	$3.4466\times {10}^{-3}$	$2.9236\times {10}^{-5}$
0.2	$1.5522\times {10}^{-2}$	$7.5193\times {10}^{-3}$	$2.8370\times {10}^{-4}$	$8.2881\times {10}^{-3}$	$4.7708\times {10}^{-3}$	$5.6468\times {10}^{-5}$
0.3	$2.1884\times {10}^{-2}$	$1.0915\times {10}^{-2}$	$2.9549\times {10}^{-4}$	$7.7752\times {10}^{-3}$	$4.8318\times {10}^{-3}$	$7.9000\times {10}^{-5}$
0.4	$2.6817\times {10}^{-2}$	$1.3327\times {10}^{-2}$	$2.5664\times {10}^{-4}$	$5.8428\times {10}^{-3}$	$4.0924\times {10}^{-3}$	$9.3988\times {10}^{-5}$
0.5	$2.9159\times {10}^{-2}$	$1.4594\times {10}^{-2}$	$1.9630\times {10}^{-4}$	$3.2486\times {10}^{-3}$	$2.9089\times {10}^{-3}$	$9.9257\times {10}^{-5}$
0.6	$2.9046\times {10}^{-2}$	$1.4586\times {10}^{-2}$	$1.3779\times {10}^{-4}$	$6.2844\times {10}^{-4}$	$1.5876\times {10}^{-3}$	$9.3988\times {10}^{-5}$
0.7	$2.6148\times {10}^{-2}$	$1.3186\times {10}^{-2}$	$9.1660\times {10}^{-5}$	$1.4621\times {10}^{-3}$	$4.0252\times {10}^{-4}$	$7.9000\times {10}^{-5}$
0.8	$2.0365\times {10}^{-2}$	$1.0310\times {10}^{-2}$	$5.7160\times {10}^{-5}$	$2.5348\times {10}^{-3}$	$3.9939\times {10}^{-4}$	$5.6468\times {10}^{-5}$
0.9	$1.1642\times {10}^{-2}$	$5.9153\times {10}^{-3}$	$2.8320\times {10}^{-5}$	$2.1665\times {10}^{-3}$	$3.9939\times {10}^{-4}$	$2.9236\times {10}^{-5}$
1.0	-	-	-	$2.0950\times {10}^{-13}$	$9.5133\times {10}^{-13}$	$2.6886\times {10}^{-13}$

Example 4.5. Consider the system of nonhomogeneous FBVPs for $x\in \left[0,1\right]$ as follows

$\{\begin{array}{ll}{u}^{″}\left(x\right)-v^{\prime }\left(x\right)+D^{0.4}v\left(x\right)+v\left(x\right)=g_1\left(x\right),\hfill & (32)\hfill \\ v^{″}\left(x\right)-u^{\prime }\left(x\right)+D^{0.6}u\left(x\right)+u\left(x\right)=g_2\left(x\right),\hfill & (33)\hfill \\ u\left(0\right)=0=u\left(1\right),v\left(0\right)=0=v\left(1\right),\hfill & (34)\hfill \end{array}$

where

$g_1\left(x\right)=11x^2+3x-1-\frac{2x^{1.6}}{\Gamma \left(2.6\right)}+\frac{2x^{0.6}}{\Gamma \left(1.6\right)},$ (35)

$g_2\left(x\right)=x^4-4x^3-x-1+\frac{24x^{3.4}}{\Gamma \left(4.4\right)}-\frac{x^{0.4}}{\Gamma \left(1.4\right)},$ (36)

that satisfies the exact solution $v\left(x\right)=x-x^2$ and $u\left(x\right)=x^4-x$ .

Accordingly, rewriting the fractional system in the operator form, and thereafter applying the inverses $I^{0.4}$ and $I^{0.6}$ on (32) and (33), respectively, via the application of (14) and (7), together with the conditions (34), one gets

$v\left(x\right)=I^{0.4}\left[g_1\left(x\right)\right]+I^{0.4}\left[v^{\prime }\left(x\right)-v\left(x\right)\right]-I^{0.4}\left[u^{″}\left(x\right)\right],$ (37)

$u\left(x\right)=I^{0.6}\left[g_2\left(x\right)\right]+I^{0.6}\left[u^{\prime }\left(x\right)-u\left(x\right)\right]-I^{0.6}\left[v^{″}\left(x\right)\right].$ (38)

Next, by applying (14) and the property (10) on $I^{0.4}\left[g_1\left(x\right)\right]$ and $I^{0.6}\left[g_2\left(x\right)\right]$ , one obtains

$\begin{array}{c}v\left(x\right)=\frac{22}{\Gamma \left(3.4\right)}{x}^{2.4}+\frac{3}{\Gamma \left(2.4\right)}{x}^{1.4}-x^2+x-\left[\frac{22}{\Gamma \left(3.4\right)}+\frac{3}{\Gamma \left(2.4\right)}\right]x\\ +I^{0.4}\left[v^{\prime }\left(x\right)-v\left(x\right)\right]-I^{0.4}\left[u^{″}\left(x\right)\right],\end{array}$

$\begin{array}{c}u\left(x\right)=\frac{\Gamma \left(5\right)}{\Gamma \left(5.6\right)}{x}^{4.6}-\frac{4\Gamma \left(4.4\right)}{\Gamma \left(4.6\right)}{x}^{3.6}-\frac{1}{\Gamma \left(2.6\right)}{x}^{1.6}+x^4-x\\ -\left[\frac{\Gamma \left(5\right)}{\Gamma \left(5.6\right)}-\frac{4\Gamma \left(4.4\right)}{\Gamma \left(4.6\right)}-\frac{1}{\Gamma \left(2.6\right)}\right]x+I^{0.6}\left[u^{\prime }\left(x\right)-u\left(x\right)\right]-I^{0.6}\left[v^{″}\left(x\right)\right].\end{array}$

Lastly, through the use of the reliable modification of the ADM, the proposed new method reveals the recursive relations for $v\left(x\right)$ and $u\left(x\right)$ , respectively, as

$\{\begin{array}{l}{v}_0\left(x\right)=x-x^2,\\ v_1\left(x\right)=\frac{22}{\Gamma \left(3.4\right)}{x}^{2.4}+\frac{3}{\Gamma \left(2.4\right)}{x}^{1.4}-\left[\frac{22}{\Gamma \left(3.4\right)}+\frac{3}{\Gamma \left(2.4\right)}\right]x\\ +I^{0.4}\left[{v^{\prime }}_0\left(x\right)-v_0\left(x\right)\right]-I^{0.4}\left[{u^{″}}_0\left(x\right)\right],\\ v_{n+1}\left(x\right)=I^{0.4}\left[{v^{\prime }}_n\left(x\right)-v_n\left(x\right)\right]-I^{0.4}\left[{u^{″}}_n\left(x\right)\right],n\ge 1,\end{array}$ (39)

and

$\{\begin{array}{l}{u}_0\left(x\right)=x^4-x,\\ u_1\left(x\right)=\frac{\Gamma \left(5\right)}{\Gamma \left(5.6\right)}{x}^{4.6}-\frac{4\Gamma \left(4.4\right)}{\Gamma \left(4.6\right)}{x}^{3.6}-\frac{1}{\Gamma \left(2.6\right)}{x}^{1.6}\\ -\left[\frac{\Gamma \left(5\right)}{\Gamma \left(5.6\right)}-\frac{4\Gamma \left(4.4\right)}{\Gamma \left(4.6\right)}-\frac{1}{\Gamma \left(2.6\right)}\right]x+I^{0.6}\left[{u^{\prime }}_0\left(x\right)\right]\\ -I^{0.6}\left[u_0\left(x\right)\right]-I^{0.6}\left[{v^{″}}_0\left(x\right)\right],\\ u_{n+1}\left(x\right)=I^{0.6}\left[{u^{\prime }}_n\left(x\right)+u_n\left(x\right)\right]-I^{0.6}\left[{v^{″}}_n\left(x\right)\right],n\ge 1.\end{array}$ (40)

Moreover, computing the values of $v_1\left(x\right)$ and $u_1\left(x\right)$ by using $I^{0.4}$ and $I^{0.6}$ , respectively using (14) and the property (10) yields $v_1\left(x\right)=u_1\left(x\right)=0$ . Thus, one obtains $v_{n+1}\left(x\right)=0$ $\forall n\ge 1$ and $u_{n+1}\left(x\right)=0$ $\forall n\ge 1$ , which subsequently leads to the solution of the fractional system as follows

$\{\begin{array}{l}v\left(x\right)=v_0\left(x\right)=x-x^2,\hfill \\ u\left(x\right)=u_0\left(x\right)=x^4-x,\hfill \end{array}$ (41)

which exactly matches the referenced analytical solution.

In the same fashion, Table 5 demonstrates the superior accuracy of the proposed iterative technique over the Sinc-collocation method (SCM) deployed in [54]. Notably, the proposed method attains the exact solution at all considered points using only two components, yielding zero absolute error for both fractional solutions $v$ and $u$ . In contrast, the SCM [54], despite using 40 components, produces nonzero errors of order 10⁻⁵ and 10⁻⁶. These results affirm the high correctness, competence, and swift convergence of the present techniques for treating systems of FBVPs.

Table 5. Comparison of absolute error for Example 4.5.

$x$	SCM [54] $\left(N=40\right)$		Present method $\left(N=2\right)$
	Error in $v$	Error in $u$	Error in $v$	Error in $u$
0.0	0	0	0	0
0.1	$1.27\times {10}^{-5}$	$7.02\times {10}^{-7}$	0	0
0.2	$2.07\times {10}^{-5}$	$2.10\times {10}^{-6}$	0	0
0.3	$1.58\times {10}^{-5}$	$3.15\times {10}^{-6}$	0	0
0.4	$1.28\times {10}^{-5}$	$3.71\times {10}^{-6}$	0	0
0.5	$1.14\times {10}^{-5}$	$3.62\times {10}^{-6}$	0	0
0.6	$1.11\times {10}^{-5}$	$3.30\times {10}^{-6}$	0	0
0.7	$1.14\times {10}^{-5}$	$3.06\times {10}^{-6}$	0	0
0.8	$1.01\times {10}^{-5}$	$2.22\times {10}^{-6}$	0	0
0.9	$8.14\times {10}^{-5}$	$1.24\times {10}^{-6}$	0	0
1.0	0	0	0	0

5. Conclusion

In this study, the variant of modified ADM is combined with the fractional inverse operator and successfully applied to treat a wide class of FBVPs. The proposed modification allowed the direct incorporation of boundary conditions into the operator, thereby eliminating the need for unknown constants and simplifying the computational process. Through illustrative examples, the efficacy and accuracy of the devised iterative scheme were demonstrated. The numerical results showed that the new scheme not only produced exact solutions for certain problems but also provided highly accurate approximations for others, outperforming several existing numerical techniques in terms of precision and convergence speed. These results confirm the efficiency and reliability of the modified ADM as a powerful semi-analytical tool for tackling FDEs with boundary conditions. Future work may focus on extending the proposed approach to other types of boundary conditions, such as Neumann or mixed cases, to evaluate its flexibility and robustness.