Application of the Optimal Homotopy Asymptotic Approach for Solving Two-Point Fuzzy Ordinary Differential Equations of Fractional Order Arising in Physics

Abstract

This work focuses on solving and analyzing two-point fuzzy boundary value problems in the form of fractional ordinary differential equations (FFOBVPs) using a new version of the approximation analytical approach. FFOBVPs are useful in describing complex scientific phenomena that include heritable characteristics and uncertainty, and obtaining exact or close analytical solutions for these equations can be challenging, especially in the case of nonlinear problems. To address these difficulties, the optimal homotopy asymptotic method (OHAM) was studied and extended in a new form to solve FFOBVPs. The OHAM is known for its ability to solve both linear and nonlinear fractional models and provides a straightforward methodology that uses multiple convergence control parameters to optimally manage the convergence of approximate series solutions. The new form of the OHAM presented in this work incorporates the concepts of fuzzy sets theory and some fractional calculus principles to include fuzzy analysis in the method. The steps of fuzzification and defuzzification are used to transform the fuzzy problem into a crisp problem that can be solved using the OHAM. The method is demonstrated by solving and analyzing linear and nonlinear FFOBVPs at different values of fractional derivatives. The results obtained using the new form of the fuzzy OHAM are analyzed and compared to those found in the literature to demonstrate the method’s efficiency and high accuracy in the fuzzy domain. Overall, this work presents a feasible and efficient approach for solving FFOBVPs using a new form of the OHAM with fuzzy analysis.

1. Introduction

Fractional-order models are more accurate than integer-order models since there are more degrees of freedom in the fractional-order models. Fractional calculus apparently captures some of the hereditary properties of the system [1]. Fractional calculus is not modern; it is a generalization of traditional calculus theory, which deals with the integer order [2]. In fractional calculus, the derivative and integral found in classical calculus are generalized to the arbitrary real or complex order, that is, to non-integer order [3]. Fractional calculus is seen as an essential tool for managing such complicated problems that are reliant on long-term memory terms, even though classical calculus is a great tool for modeling many complex real-world phenomena [4]. Memory is the term used to describe the output or results that depend on the history of the variables from a previous period. Classical calculus cannot solve issues that depend on the memory of the variables [5].

During the past decade, fractional differential equations under the effect of uncertainty (FFDEs) have appeared more and more practically in different research areas, such as physics and engineering [6,7], in addition to many other fields [8]. FFDEs are characterized by a nonlocal derivative operator. This, in turn, contributes to modeling the complicated real-world problems that are based on the long memory term. Unfortunately, the stemming uncertainty caused by a lack of data or the difficulty of exactly determining the supplementary conditions will lead to errors in measurement, so using the nonlocal fractional derivative operators in the fuzzy environment will ensure a more accurate mathematical model that simulates human thinking.

Accurate modeling of complex real-world problems helps us provide a clear and explicit concept of complex dynamics by employing the definitions and theories of the fractional calculus theory and the fuzzy calculus theory. However, these models remain impractical until they are solved because the solutions provide a comprehensive view, in addition to the fact that the solutions aid in studying and understanding the physical and engineering properties of real-world problems [8]

In solving some of the FFODEs, the analytical approach aims to present a closed-form solution [9,10]. A closed-form solution is considered the exact solution to the problem [11]. The solution may be expressed as the sum of a finite number of elementary functions, such as polynomial, exponential, trigonometric, and hyperbolic functions. The advantage of a closed-form solution is that it provides an overall view of the solution to the problem. Moreover, in the analysis of results, using closed-form solutions generally does not require a huge amount of computation [12]. In many instances, analytical solutions cannot be found [13,14,15,16]. Nevertheless, the solutions to such equations are always in demand due to practical interests. Therefore, to deal with such instances in a more realistic manner, FFODEs are commonly solved using the approximation approach, which includes the numerical and approximate analytical methods.

Numerous methods were proposed for solving FFOBVPs; for instance, we refer the reader to explore [16,17,18,19]. These numerical methods demonstrated their ability to solve only linear cases of FFODEs. In the numerical approach, the aim is to obtain an approximate solution, where an open-form solution is sought instead of a closed-form. However, the numerical class of methods directly solves FFODEs of high orders; instead, they require a transformation into a system of the first order. Further, most studies employ numerical methods for linear first-order problems [20,21]. Unfortunately, most of the complicated real-world problems were modeled using nonlinear differential equations, which makes these methods inappropriate to deal with them—especially the problems governed by strong nonlinearity.

In addition to the optimal homotopy asymptotic method (OHAM) presented in this work, several other approximate analytical methods have been used to solve different types of FFOBVPs. These include the variational iteration method (VIM) [22,23], the reproducing kernel Hilbert space method (RKHSM) [24], the spectral collocation method (SCM) [7], the Adomian decomposition method (ADM) [25], the differential transform method [26], the residual power series method (RPSM) [27], and the fractional residual power series method (FRPSM) [28].

Hashim et al. [29] solved fuzzy IVPs with fractional derivative orders between 0 and 1 using the optimal homotopy asymptotic method (OHAM), and the paper presented the defuzzification of fuzzy fractional IVPs. The authors then introduced a framework for solving the considered problem using the OHAM. Upper and lower solutions were investigated in terms of the accuracy and convergence of the method by finding optimal values of the convergent parameters using a few terms of the series solution with higher accuracy than the fractional residual power method. The paper did not discuss the fuzzification of the boundary value problem or its solution. As one more section of this work, we will investigate the fuzzy theory and the OHAM to solve the fuzzy boundary value problem with a fractional derivative order between 1 and 2.

While these methods have shown promise in solving FFOBVPs, they often fail to provide a simple way to control and adjust the convergence area. This can be a significant obstacle in obtaining accurate solutions, especially for nonlinear problems. Therefore, the development of new approximate analytical methods is necessary to overcome these challenges and improve the accuracy of solutions.

The OHAM presented in this work addresses this issue by using multiple convergence control parameters to optimally manage the convergence of approximate series solutions. This allows for greater flexibility in controlling and adjusting the convergence area, resulting in more accurate solutions for FFOBVPs. This is due to the proposed methods’ inability to control the convergence region. Nowadays, homotopy methods are the most promising approaches for solving nonlinear real-world problems [30] due to their ability to simplify complicated problems, provide the freedom to choose the auxiliary functions, and provide a simple way to control the convergence, which helps us optimize the convergence series solutions for the strong nonlinearity problems. The OHAM has been used to solve various types of differential equations, including classical differential equations [31], fuzzy differential equations [32], and fractional differential equations [33]. The method’s effectiveness in controlling the convergence area has been demonstrated through numerical results for both linear and nonlinear problems. The OHAM provides multiple convergence control parameters that allow for greater flexibility in adjusting the convergence area and obtaining accurate solutions. This makes OHAM a promising method for solving complex differential equations encountered in science and engineering applications.

In order to solve FFOBVPs, this study intended to create novel approximative analytical techniques with convergence-control capabilities. The fundamental idea of the OHAM will be applied to the development of the new approach, which will be able to manage the significant challenge of managing the convergence of the approximative analytical solutions. This work also focuses on the development of the fuzzy OHAM’s fractional form, represented by the abbreviation FF-OHAM, on two different types of application problems that fall under the Caputo definitions of differentiability and involve linear and nonlinear application problems. This starts with the introduction of the basic tools of fuzzy fractional calculus in the second section, followed by providing the defuzzification procedure for the FFOBVPs in Section 3. Section 4 provides the new version of the FF-OHAM for solving FFOBVPs; then, the numerical simulation of the physical applications of the FF-OHAM will be provided in Section 5. Then, we will end with the conclusions regarding the effectiveness of the proposed method and the gained results.

2. Mathematical Background

In this section, we will present the basic concepts and definitions linked with fractional calculus theory in the fuzzy domain, which will help us comprehend the work in the next sections, such as the fuzzy fractional integral [34], which is a generalization of the classical fractional integral concept to the fuzzy-valued functions. It is a fuzzy operator that takes a fuzzy-valued function [35] as an input and produces another fuzzy-valued function as an output. The fuzzy fractional integral can be interpreted as a generalization of the fuzzy integral and the classical Riemann–Liouville integral of fractional order [36]. The Caputo derivative of fractional order is used to define the fuzzy fractional integral in the sense of Caputo [37]. It is worth noting that the concept of α-cut is also used in the fuzzy fractional integral theory to define the α-cut of a fuzzy fractional integral. The α-cut of a fuzzy fractional integral is a fuzzy number that corresponds to the fuzzy α-cut [38]. The fuzzy-valued function is obtained by taking the fuzzy fractional integral and the Riemann–Liouville integral of fractional order [24]. On the other hand, the following fundamental definition of the fuzzy fractional integral needs to be recalled:

Definition 1

([7]). For any continuous fuzzy valued function, $\tilde{g}\in C^{\mathcal{F}}\left[a,b\right]\cap L^{\mathcal{F}}\left[a,b\right]$, the fuzzy fractional Riemann–Liouville integration of $\tilde{g}\left(x\right)$ will be defined by the following form:

$${\tilde{\mathcal{J}}}^{\omega }\tilde{g}\left(x\right)=\frac{1}{\Gamma \left(\omega \right)}{\int }_0^x\tilde{g}\left(y\right){\left(x-y\right)}^{\omega -1}dy, \mathrm{for} \omega ,x\in \mathbb{R} \mathrm{and} \omega ,x>0$$

(1)

$\forall \alpha \in \left[\mathrm{0,1}\right], \alpha$-cuts for fuzzy-valued function, $\tilde{g}$, can be represented by 

$$\tilde{g}\left(x;\alpha \right)=\left[\underset{¯}{g}\left(x;\alpha \right),\overline{g}\left(x;\alpha \right)\right]$$

(2)

where ${\tilde{\mathcal{J}}}^{\omega }$ is the Riemann–Liouville integral operator of order $\omega$; $\Gamma \left(\omega \right)$ is the famous Gamma function; $C^{\mathcal{F}}\left[a,b\right]$ is the set of all fuzzy-valued measurable functions, $\tilde{g}$, on $[a,b]$; and $L^{\mathcal{F}}\left[a,b\right]$ is the space of fuzzy-valued functions, which are continuous on $[a,b]$.

Definition 2

([24]). Let $\omega \in \left(\mathrm{1,2}\right]$, and $\tilde{g}:\left[a,b\right]\to \tilde{U}$, such that $\tilde{g}$ and ${\tilde{g}}^{\prime }$ $\in C^{\mathcal{F}}\left[0,b\right]\cap L^{\mathcal{F}}\left[0,b\right]$. Then, $\mathcal{F}$ can define the fuzzy fractional derivative in the sense of the Caputo of the fuzzy function $\tilde{g}$ at $x\in \left(a,b\right)$, as follows:

$$\left({\tilde{D}}^{\omega }\tilde{g}\right)\left(x\right)=\frac{1}{\Gamma \left(2-\omega \right)}{\int }_0^x\frac{{\tilde{g}}^{″}\left(x\right)}{{\left(y-x\right)}^{\omega -1}}dx,x>0$$

(3)

where $D$ is the Housdorff metric of the fuzzy set $\tilde{U}$. Note that the fuzzy fractional Riemann–Liouville integration represents the left inverse operator of the fuzzy fractional Caputo derivative sense, such that $\forall$$\tilde{g}\left(x\right)\in C^{\mathcal{F}}\left[a,b\right]\cap L^{\mathcal{F}}\left[a,b\right]$. We have

$${\tilde{\mathcal{J}}}^{\omega }\left({\tilde{D}}^{\omega }\tilde{g}\right)\left(x\right)=\tilde{g}\left(x\right)-x{g}^{\prime }\left(0\right)-g\left(0\right),x\in \mathbb{R}, \mathrm{and} x>0.$$

(4)

3. Fuzzification and Defuzzification of FFODEs

The first step of the development of the proposed FF-OHAM for solving second-order FFOBVPs is the defuzzification step. This is a general step that applies to the general form of second-order FFOBVPs, as shown below.

Consider the second-order FFOBVP as follows:

$$\left\{\begin{array}{c}{\tilde{y}}^{\left(\omega \right)}\left(x\right)=\tilde{g}\left(x,\tilde{y}\left(x\right),{\tilde{y}}^{\left(1\right)}\left(x\right)\right),x\in \left[x_0,X\right],\\ 1<\omega \le 2,\end{array}\right.$$

(5)

subject to the following boundary conditions:

$$\left\{\begin{array}{c}\tilde{y}\left(x_0\right)={\tilde{a}}_0,{\tilde{y}}^{\left(1\right)}\left(x_0\right)={\tilde{a}}_1,\\ \tilde{y}\left(X\right)={\tilde{b}}_0,{\tilde{y}}^{\left(1\right)}\left(X\right)={\tilde{b}}_1,\end{array}\right.$$

(6)

where $\tilde{g}$ is the fuzzy function, while ${\tilde{y}}^{\left(\omega \right)}\left(x\right)$ is the fractional Caputo derivative of the fuzzy function $\tilde{y}\left(x\right)$; and the boundary conditions at the points $x_0$ and $X$ are fuzzy numbers.

For $x\in \left[x_0,X\right]$ and $\forall$$\alpha \in \left[\mathrm{0,1}\right]$, the fuzzy function will be defined by ${\left[\tilde{y}\right]}_{\alpha }={[\underset{¯}{y},\overline{y}]}_{\alpha }$ $\forall x\in \left[x_0,X\right]$ as follows:

$$\left\{\begin{array}{c}{\left[\tilde{y}\left(x_0\right)\right]}_{\alpha }=\left[\underset{¯}{y}\left(x_0;\alpha \right),\overline{y}\left(x_0;\alpha \right)\right],\\ {\left[{\tilde{y}}^{\left(1\right)}\left(x_0\right)\right]}_{\alpha }=\left[{\underset{¯}{y}}^{\left(1\right)}\left(x_0;\alpha \right),{\stackrel{-}{y}}^{\left(1\right)}\left(x_0;\alpha \right)\right],\\ {\left[\tilde{y}\left(X\right)\right]}_{\alpha }=\left[\underset{¯}{y}\left(X;\alpha \right),\overline{y}\left(X;\alpha \right)\right],\\ {\left[{\tilde{y}}^{\left(1\right)}\left(X\right)\right]}_{\alpha }=\left[{\underset{¯}{y}}^{\left(1\right)}\left(X;\alpha \right),{\stackrel{-}{y}}^{\left(1\right)}\left(X;\alpha \right)\right].\end{array}\right.$$

(7)

Now, by assuming $\widetilde{\mathrm{Ŷ}}\left(x\right)=\left\{\tilde{y}\left(x\right),{\tilde{y}}^{\left(1\right)}\left(x\right)\right\}$, for defuzzification we have:

$$\widetilde{\mathrm{Ŷ}}\left(x,\alpha \right)=\left[\underset{¯}{\mathrm{Ŷ}}\left(x,\alpha \right),\stackrel{-}{\mathrm{Ŷ}}\left(x,\alpha \right)\right]=[\underset{¯}{y}\left(x,\alpha \right),{\underset{¯}{y}}^{\left(1\right)}\left(x,\alpha \right),\overline{y}\left(x;\alpha \right),{\stackrel{-}{y}}^{\left(1\right)}\left(x;\alpha \right)].$$

(8)

In addition, by utilizing the concepts of the extension principle theory, we can write the $\alpha$-cut of the fuzzy function, as shown below:

$${\left[\tilde{g}\left(x,\widetilde{\mathrm{Ŷ}}\right)\right]}_{\alpha }=\tilde{g}\left(x,\widetilde{\mathrm{Ŷ}}\left(x;\alpha \right)\right)=\left[\underset{¯}{g}\left(x,\widetilde{\mathrm{Ŷ}};\alpha \right),\overline{g}\left(x,\widetilde{\mathrm{Ŷ}};\alpha \right)\right],$$

(9)

where

$$\left\{\begin{array}{c}\underset{¯}{g}\left(x,\widetilde{\mathrm{Ŷ}}\left(x;\alpha \right)\right)=g_l\left(x,\underset{¯}{\mathrm{Ŷ}}\left(t;\alpha \right),\overline{\mathrm{\Upsilon }}\left(x;\alpha \right)\right)=g_l\left(x,\overline{\mathrm{Ŷ}}\left(x;\alpha \right)\right),\\ \overline{g}\left(x,\widetilde{\mathrm{Ŷ}}\left(x;\alpha \right)\right)=g_u\left(x,\underset{¯}{\mathrm{Ŷ}}\left(x;\alpha \right),\overline{\mathrm{\Upsilon }}\left(x;\alpha \right)\right)=g_u\left(x,\overline{\mathrm{Ŷ}}\left(x;\alpha \right)\right).\end{array}\right.$$

(10)

which means that $\forall \alpha \in \left[\mathrm{0,1}\right]$. We have

$$\left\{\begin{array}{c}{\underset{¯}{y}}^{\left(\omega \right)}\left(x;\alpha \right)=g_l\left(x,\widetilde{\mathrm{Ŷ}}\left(x;\alpha \right)\right),\\ {\overline{y}}^{\left(\omega \right)}\left(x;\alpha \right)=g_u\left(x,\widetilde{\mathrm{Ŷ}}\left(x;\alpha \right)\right).\end{array}\right.$$

(11)

where

$$\left\{\begin{array}{c}{g}_l\left(x,\widetilde{\mathrm{Ŷ}}\left(x;\alpha \right)\right)=\min \left\{{\tilde{y}}^{\left(\omega \right)}\left(x,\tilde{\mu }\left(\alpha \right)\right)|\tilde{\mu }\left(\alpha \right)\in {\left[\widetilde{\mathrm{Ŷ}}\left(x;\alpha \right)\right]}_{\alpha }\right\},\\ g_u\left(x,\widetilde{\mathrm{Ŷ}}\left(x;\alpha \right)\right)=\max \left\{{\tilde{y}}^{\left(\omega \right)}\left(x,\tilde{\mu }\left(\alpha \right)\right)|\tilde{\mu }\left(\alpha \right)\in {\left[\widetilde{\mathrm{Ŷ}}\left(x;\alpha \right)\right]}_{\alpha }\right\}.\end{array}\right.$$

(12)

4. FF-OHAM for FFTBVPs

In this section, the F-OHAM presented by [21] for solving the integer order of ODEs is fuzzified and then defuzzified using some concepts of the fuzzy set theory in Section 2 to create a new form of the method denoted by the FF-OHAM for solving linear and nonlinear second-order FFOBVPs approximately.

$$\left\{\begin{array}{c}{\tilde{y}}^{\left(\omega \right)}\left(x\right)=\tilde{g}\left(x,\tilde{y}\left(x\right),{\tilde{y}}^{\left(1\right)}\left(x\right)\right)+\tilde{\mathrm{G}}\left(x\right)\hspace{1em}x\in \left[x_0,X\right],\\ \tilde{y}\left(x_0\right)={\tilde{a}}_0,\hspace{1em}{\tilde{y}}^{\prime }\left(x_0\right)={\tilde{a}}_1,\\ \tilde{y}\left(X\right)={\tilde{b}}_0,\hspace{1em}{\tilde{y}}^{\prime }\left(X\right)={\tilde{b}}_1,\\ \omega \in \left(\mathrm{1,2}\right],\end{array}\right.$$

(13)

Followed by the defuzzification of Equation (5), such that for all $\alpha \in \left[\mathrm{0,1}\right]$, we have the following lower bound:

$$\left\{\begin{array}{c}{\underset{¯}{\mathcal{L}}}_{{\beta }_1}\left(\tilde{y}\left(x;\alpha \right)\right)-g_l\left(x,\widetilde{\mathrm{Ŷ}}\left(x\right)\right)-\underset{¯}{\mathrm{G}}\left(x\right)=0,\hspace{1em}x\in \left[x_0,X\right],\\ \mathrm{ß}\left(\underset{¯}{y}\left(x;\alpha \right),\frac{\partial {\left[\underset{¯}{y}\right]}_{\alpha }}{\partial x}\right)=0.\end{array}\right.$$

(14)

and the following upper bound:

$$\left\{\begin{array}{c}{\overline{\mathcal{L}}}_{{\beta }_1}\left(\tilde{y}\left(x;\alpha \right)\right)-g_u\left(x,\widetilde{\mathrm{Ŷ}}\left(x\right)\right)-\overline{\mathrm{G}}\left(x\right)=0,x\in \left[x_0,X\right],\\ \mathrm{ß}\left(\overline{y}\left(x;\alpha \right),\frac{\partial {\left[\underset{¯}{\overline{y}}\right]}_{\alpha }}{\partial x}\right)=0.\end{array}\right.$$

(15)

According to [3], Equations (14) and (15) can be written as the following lower and upper zeroth order deformation homotopy equation:

$$\left\{\begin{array}{c}\left(1-q\right)\left[{\underset{¯}{\mathcal{L}}}_{\omega }\left({\left[\underset{¯}{y}\left(x;q\right)\right]}_{\alpha }\right)-\underset{¯}{\mathrm{G}}\left(x;\alpha \right)\right]=\underset{¯}{\mathcal{H}}\left(q;\alpha \right)\left[{\underset{¯}{\mathcal{L}}}_{\omega }\left({\left[\underset{¯}{y}\left(x;q\right)\right]}_{\alpha }\right)\right]\\ -\underset{¯}{\mathcal{H}}\left(q;\alpha \right)\left[\underset{¯}{\mathrm{G}}\left(x;\alpha \right)\right]-\underset{¯}{\mathcal{H}}\left(q;\alpha \right)\left[g_l\left({\left[\tilde{y}\left(x;q\right)\right]}_{\alpha }\right)\right],\\ \left(1-q\right)\left[{\overline{\mathcal{L}}}_{\omega }\left({\left[\overline{y}\left(x;q\right)\right]}_{\alpha }\right)-\overline{\mathrm{G}}\left(x;\alpha \right)\right]=\overline{\mathcal{H}}\left(q;\alpha \right)\left[{\overline{\mathcal{L}}}_{\omega }\left({\left[\overline{y}\left(x;q\right)\right]}_{\alpha }\right)\right]\\ -\overline{\mathcal{H}}\left(q;\alpha \right)\left[\overline{\mathrm{G}}\left(x;\alpha \right)\right]-\overline{\mathcal{H}}\left(q;\alpha \right)\left[g_l\left({\left[\tilde{y}\left(x;q\right)\right]}_{\alpha }\right)\right],\end{array}\right.$$

(16)

subject to the following fuzzy boundary conditions

$$\mathrm{ß}\left({\left[\tilde{y}\left(x;q\right)\right]}_{\alpha },\frac{\partial {\left[\tilde{y}\left(x;q\right)\right]}_{\alpha }}{\partial x}\right)=0,$$

(17)

where ${\tilde{\mathcal{L}}}_{\omega }=\left[{\underset{¯}{\mathcal{L}}}_{\omega },{\overline{\mathcal{L}}}_{\omega }\right]=\left[\frac{{\partial }^{\left(\omega \right)}{\left[\underset{¯}{y}\left(x;q\right)\right]}_{\alpha }}{\partial x^{\left(\omega \right)}},\frac{{\partial }^{\left(\omega \right)}{\left[\overline{y}\left(x;q\right)\right]}_{\alpha }}{\partial x^{\left(\omega \right)}}\right]$ are the linear operators and $q\in [0,1]$ is an embedding parameter. Here, $\tilde{\mathcal{H}}\left(q;\alpha \right)={\left[\underset{¯}{\mathcal{H}}\left(q\right),\overline{\mathcal{H}}\left(q\right)\right]}_{\alpha }$ is a nonzero auxiliary fuzzy function for $q\ne 0$, and ${\left[\tilde{y}\left(x;q\right)\right]}_{\alpha }$ is an unknown fuzzy function.

Obviously, for $q=0$ and $q=1$, we obtain the initial approximation, and the exact solution, respectively, as follows:

$$\left\{\begin{array}{c}{\left[\underset{¯}{y}\left(x;0\right)\right]}_{\alpha }={\underset{¯}{y}}_0\left(x;\alpha \right),{\left[\underset{¯}{y}\left(x;1\right)\right]}_{\alpha }=\underset{¯}{Y}\left(x;\alpha \right),\\ {\left[\overline{y}\left(x;0\right)\right]}_{\alpha }={\overline{y}}_0\left(t;\alpha \right).{\left[\overline{y}\left(x;1\right)\right]}_{\alpha }=\overline{Y}\left(x;\alpha \right).\end{array}\right.$$

(18)

Thus, as $q$ increases from 0 to 1, the series solution, ${\left[\tilde{y}\left(x;q\right)\right]}_{\alpha }$, changes from ${\tilde{y}}_0\left(x;\alpha \right)$ to the solution of Equations (14) and (15), $\tilde{Y}\left(x;\alpha \right)$, where ${\tilde{y}}_0\left(x;\alpha \right)$ is obtained from Equation (16) for $q=0$ as follows:

$$\left\{\begin{array}{c}{\underset{¯}{y}}_0\left(x;\alpha \right)={\tilde{\mathcal{J}}}^{\left({\beta }_1\right)}\underset{¯}{\mathrm{G}}\left(x;\alpha \right),\\ {\overline{y}}_0\left(x;\alpha \right)={\tilde{\mathcal{J}}}^{\left({\beta }_1\right)}\overline{\mathrm{G}}\left(x;\alpha \right),\end{array}\right.$$

(19)

subject to the following fuzzy boundary condition

$$\mathrm{ß}\left({\tilde{y}}_0\left(x;\alpha \right),\frac{\partial {\left[{\tilde{y}}_0\right]}_{\alpha }}{\partial x}\right)=0$$

(20)

We choose the auxiliary function $\tilde{\mathcal{H}}(q;\alpha )$ for Equation (16) in the following form:

$$\left\{\begin{array}{c}\underset{¯}{\mathcal{H}}\left(q;\alpha \right)={\underset{¯}{S}}_1(\alpha )q+{\underset{¯}{S}}_2(\alpha )q^2+\dots ={\displaystyle \sum _{j=1}^k}{\underset{¯}{S}}_j{\left(\alpha \right)q}^j,\\ \underset{¯}{\mathcal{H}}\left(q;\alpha \right)={\overline{S}}_1(\alpha )q+{\overline{S}}_2{\left(\alpha \right)q}^2+\dots ={\displaystyle \sum _{j=1}^k}{\overline{S}}_j\left(\alpha \right)q^j,\end{array}\right.$$

(21)

where ${\tilde{S}}_1\left(\alpha \right)=\left[{\underset{¯}{S}}_1\left(\alpha \right),{\overline{S}}_1\left(\alpha \right)\right],{\tilde{S}}_2\left(\alpha \right)=\left[{\underset{¯}{S}}_2\left(\alpha \right),{\overline{S}}_2\left(\alpha \right)\right],\dots$ are the constants to be found for all $\alpha \in \left[\mathrm{0,1}\right]$. Now, by expanding ${\left[\tilde{y}\left(x;q,S_j\left(\alpha \right)\right)\right]}_{\alpha }$ into Taylor’s series about $q$, we obtain the following approximate series solution:

$$\left\{\begin{array}{c}{\left[\underset{¯}{y}\left(x;q,{\underset{¯}{S}}_j\left(\alpha \right)\right)\right]}_{\alpha }={\underset{¯}{y}}_0\left(x;\alpha \right)+{\displaystyle \sum _{j=1}^k}{\left[{\underset{¯}{y}}_j\left(x{,\underset{¯}{S}}_j\left(\alpha \right)\right)\right]}_{\alpha }{q}^j,\\ {\left[\overline{y}\left(x;q,{\overline{S}}_j\left(\alpha \right)\right)\right]}_{\alpha }={\overline{y}}_0\left(x;\alpha \right)+{\displaystyle \sum _{j=1}^k}{\left[{\overline{y}}_j\left(x{,\overline{S}}_j\left(\alpha \right)\right)\right]}_{\alpha }{q}^j.\end{array}\right.$$

(22)

According to [3], by substituting Equation (22) into Equation (16) and then collecting the coefficient of like powers of $q$, we will obtain the following system of linear equations—where the zeroth-order problem is given by Equation (19), while the first to $k^{th}$ -order problems are given as in the general $k^{th}$-order formula with respect to ${\tilde{y}}_k\left(x;\alpha \right),$ for $k\ge 1$:

$$\left\{\begin{array}{c}{\underset{¯}{y}}_k\left(x;\alpha \right)={\underset{¯}{y}}_{k-1}\left(x;\alpha \right)+{\displaystyle \sum _{j=1}^{k-1}}{\underset{¯}{S}}_j(\alpha )\left({\underset{¯}{y}}_{k-j}\left(x;\alpha \right)\right)\\ {\underset{¯}{\mathcal{J}}}^{\left(\omega \right)}\left({\underset{¯}{S}}_k\left(\alpha \right){g_l}_0\left({\tilde{y}}_0\left(x;\alpha \right)\right)+{\sum }_{j=1}^{k-1}{\underset{¯}{S}}_j(\alpha ){g_l}_{k-j}\left({\displaystyle \sum _{i=0}^{k-1}}{\tilde{y}}_i\left(x;\alpha \right)\right)\right)\\ {\overline{y}}_k\left(x;\alpha \right)={\overline{y}}_{k-1}\left(x;\alpha \right)+{\displaystyle \sum _{j=1}^{k-1}}{\overline{S}}_j(\alpha )\left({\overline{y}}_{k-j}\left(x;\alpha \right)\right)\\ {\overline{\mathcal{J}}}^{\left(\omega \right)}\left({\overline{S}}_k\left(\alpha \right){g_l}_0\left({\tilde{y}}_0\left(x;\alpha \right)\right)+{\sum }_{j=1}^{k-1}{\overline{S}}_j(\alpha ){g_u}_{k-j}\left({\displaystyle \sum _{i=0}^{k-1}}{\tilde{y}}_i\left(x;\alpha \right)\right)\right)\end{array}\right.$$

(23)

$$\mathrm{ß}\left({\tilde{y}}_k\left(x;\alpha \right),\frac{\partial {\left[{\tilde{y}}_k\right]}_{\alpha }}{\partial x}\right)=0$$

(24)

where ${g_l}_{k-j}\left({\displaystyle {\sum }_{i=0}^{k-1}}{\tilde{y}}_i\left(x;\alpha \right)\right)$ and ${g_u}_{k-j}\left({\displaystyle {\sum }_{i=0}^{k-1}}{\tilde{y}}_i\left(x;\alpha \right)\right)$ are the coefficients of $q^j$ in the expansion of $g_l{\left[\tilde{y}\left(x;q\right)\right]}_{\alpha }$ and $g_u{\left[\tilde{y}\left(x;q\right)\right]}_{\alpha }$ about the embedding parameter $q$. We have the lower and upper bounds as follows:

$$\left\{\begin{array}{c}{g}_l\left({\left[\tilde{y}\left(x;q,{\displaystyle \sum _{j=1}^k}{\tilde{S}}_j\left(\alpha \right)\right)\right]}_{\alpha }\right)={g_l}_0({\tilde{y}}_0\left(x;\alpha \right))+{\displaystyle \sum _{j=1}^k}{g_l}_j\left({\displaystyle \sum _{j=0}^k}{\left[{\tilde{y}}_j\right]}_{\alpha }\right)q^j,\\ g_u\left({\left[\tilde{y}\left(x;q,{\displaystyle \sum _{j=1}^k}{\tilde{S}}_j\left(\alpha \right)\right)\right]}_{\alpha }\right)={g_u}_0({\tilde{y}}_0\left(x;\alpha \right))+{\displaystyle \sum _{j=1}^k}{g_u}_j\left({\displaystyle \sum _{j=0}^k}{\left[{\tilde{y}}_j\right]}_{\alpha }\right)q^j.\end{array}\right.$$

(25)

It has been observed that the convergence of the series in Equation (22) depends upon the auxiliary constants ${\tilde{S}}_1(\alpha ),{\tilde{S}}_2(\alpha ),\dots {\tilde{S}}_k(\alpha )$, then, at $q=1$, we obtain the exact solution shown below:

$$\left\{\begin{array}{c}{\left[\underset{¯}{Y}\left(x,{\displaystyle \sum _{j=1}^{\infty }}{\underset{¯}{S}}_j\left(\alpha \right)\right)\right]}_{\alpha }={\underset{¯}{y}}_0\left(x;\alpha \right)+{\displaystyle \sum _{j=1}^{\infty }}{\left[{\underset{¯}{y}}_j\left(x;{\displaystyle \sum _{j=1}^{\infty }}{\underset{¯}{S}}_j\left(\alpha \right)\right)\right]}_{\alpha },\\ {\left[\overline{Y}\left(x,{\displaystyle \sum _{j=1}^{\infty }}{\overline{S}}_j\left(\alpha \right)\right)\right]}_{\alpha }={\overline{y}}_0\left(x;\alpha \right)+{\displaystyle \sum _{j=1}^{\infty }}{\left[{\overline{y}}_j\left(x;{\displaystyle \sum _{j=1}^{\infty }}{\overline{S}}_j\left(\alpha \right)\right)\right]}_{\alpha }.\end{array}\right.$$

(26)

5. Convergence Dynamic of the FF-OHAM

Substituting Equation (22) into Equations (14) and (15) yields the following residual:

$$\left\{\begin{array}{c}\underset{¯}{RE}\left(x,{\displaystyle \sum _{j=1}^k}{\underset{¯}{S}}_j\left(\alpha \right);\alpha \right)={\underset{¯}{\mathcal{L}}}_{\omega }\left(\underset{¯}{y}\left(x,{\displaystyle \sum _{j=1}^k}{\underset{¯}{S}}_j\left(\alpha \right);\alpha \right)\right)-\underset{¯}{\mathrm{G}}\left(x;\alpha \right)\\ -g_l\left(\tilde{y}\left(x,{\displaystyle \sum _{j=1}^k}{\tilde{S}}_j\left(\alpha \right);\alpha \right)\right),\\ \overline{RE}\left(x,{\displaystyle \sum _{j=1}^k}{\overline{S}}_j\left(\alpha \right);\alpha \right)={\stackrel{-}{\mathcal{L}}}_{\omega }\left(\overline{y}\left(x,{\displaystyle \sum _{j=1}^k}{\overline{S}}_j\left(\alpha \right);\alpha \right)\right)-\overline{\mathrm{G}}\left(x;\alpha \right)\\ -g_u\left(\tilde{y}\left(x,{\displaystyle \sum _{j=1}^k}{\tilde{S}}_j\left(\alpha \right);\alpha \right)\right).\end{array}\right.$$

(27)

As mentioned in [22], if $\widetilde{RE}=0$, then $\tilde{y}$ yields the exact solution $\tilde{Y}$, although, generally, it does not happen, especially in nonlinear FFOBVPs. To identify the auxiliary constants of ${\tilde{S}}_j\left(\alpha \right)$, $j=\mathrm{1,2},\dots k$, we choose $x_0$ and $X,$ such that the optimum values of ${\tilde{S}}_j\left(\alpha \right)$ for the convergent solution of the desired problem is obtained. To find the optimal values of ${\tilde{S}}_j\left(\alpha \right)$ here, we apply the method of least squares as follows:

$$\widetilde{SRE}\left(x,{\displaystyle \sum _{j=1}^k}{\tilde{S}}_j\left(\alpha \right);\alpha \right)={\displaystyle \underset{x_0}{\overset{X}{\int }}}{\widetilde{RE}}^2\left(x,{\displaystyle \sum _{j=1}^k}{\tilde{S}}_j(\alpha );\alpha \right)dx,$$

(28)

where $\widetilde{RE}$ is the residual,

$$\left\{\begin{array}{c}{\left[\underset{¯}{RE}\right]}_{\alpha }={\underset{¯}{\mathcal{L}}}_{\omega }\left({\left[\underset{¯}{y}\right]}_{\alpha }\right)-\underset{¯}{\mathrm{G}}\left(x;\alpha \right)-g_l\left({\left[\tilde{y}\right]}_{\alpha }\right)\\ {\left[\overline{RE}\right]}_{\alpha }={\stackrel{-}{\mathcal{L}}}_{\omega }\left({\left[\overline{y}\right]}_{\alpha }\right)-\overline{\mathrm{G}}\left(x;\alpha \right)-g_u\left({\left[\tilde{y}\right]}_{\alpha }\right)\end{array}\right.$$

(29)

and

$$\frac{\partial \widetilde{SRE}}{\partial {\tilde{S}}_1(\alpha )}=\frac{\partial \widetilde{SRE}}{\partial {\tilde{S}}_2(\alpha )}=\dots \frac{\partial \widetilde{SRE}}{\partial {\tilde{S}}_k\left(\alpha \right)}=0.$$

(30)

It should be noted that our process included the fuzzy level set $\alpha$, so the best values of ${\tilde{S}}_k\left(\alpha \right)$ are determined from Equation (30) for each $\alpha \in \left[\mathrm{0,1}\right]$, which provides us with an easy way to set and optimally control the convergent area and the rate of the solution series.

6. Numerical Simulation of the Physical Applications via FF-OHAM

This section reflects the use of the FF-OHAM from Section 3 and Section 4 for some fuzzy models in physics. The method’s performance is tested in two linear and nonlinear FFOBVPs applications.

• Mechanical Application: Fuzzy Fractional Bagley–Torvik Equation

Consider the fuzzy fractional Bagley–Torvik equation [7]:

$$D^{\left(1.5\right)}\tilde{y}\left(x\right)+\tilde{y}\left(x\right)=\tilde{F}\left(x;\alpha \right),\hspace{1em}x\in \left[\mathrm{0,1}\right],$$

(31)

such that

$$\tilde{F}\left(x;\alpha \right)=\left\{\begin{array}{c}\underset{¯}{F}\left(x;\alpha \right)\\ \overline{F}\left(x;\alpha \right)\end{array}\right.=\left\{\begin{array}{c}\alpha \left(x^2-x\right)+4\alpha \frac{\sqrt{x}}{\sqrt{\mathrm{╥}}},\\ \left(2-\alpha \right)\left(x^2-x\right)+4\left(2-\alpha \right)\frac{\sqrt{x}}{\sqrt{\mathrm{╥}}}.\end{array}\right.$$

(32)

subject to the following fuzzy boundary condition

$$\left\{\begin{array}{c}\underset{¯}{y}\left(0;\alpha \right)=\underset{¯}{y}\left(1;\alpha \right)=\left(\alpha -1\right),\\ \overline{y}\left(0;\alpha \right)=\overline{y}\left(1;\alpha \right)=\left(1-\alpha \right).\end{array}\right.$$

(33)

with the following fuzzy exact solution

$$\left\{\begin{array}{c}\underset{¯}{Y}\left(x;\alpha \right)=\alpha \left(x^2-x\right),\\ \overline{Y}\left(x;\alpha \right)=\left(2-\alpha \right)\left(x^2-x\right).\end{array}\right.$$

(34)

we can contract the FF-OHAM series solution for all $\alpha \in \left[\mathrm{0,1}\right]$ of Equation (31) as follows:

$$\left\{\begin{array}{c}\left(1-q\right)\left[{\tilde{D}}^{\left(1.5\right)}\left(\tilde{y}\left(x;\alpha \right)\right)-\tilde{F}\left(x;\alpha \right)\right]={\displaystyle \sum _{j=1}^5}{\tilde{S}}_j{\left(\alpha \right)q}^j\\ {\tilde{S}}_j{\left(\alpha \right)q}^j\left[{\tilde{D}}^{\left(1.5\right)}\left(\tilde{y}\left(x;\alpha \right)\right)+\tilde{y}\left(x;\alpha \right)-\left(\tilde{F}\left(x;\alpha \right)\right)\right]\end{array}\right.$$

(35)

such that

$$\tilde{y}\left(x;\alpha \right)={\tilde{y}}_0(x;\alpha )+{\displaystyle \sum _{j=1}^k}{\tilde{y}}_j(x,S_1,\dots ,S_j;\alpha )q^j$$

(36)

Zeroth-order problem:

$$\left\{\begin{array}{c}{\tilde{y}}_0\left(x,\alpha \right)={\tilde{\mathcal{J}}}^{\left(1.5\right)}\left[\alpha \left(\left(x^2-x\right)+4\frac{\sqrt{x}}{\sqrt{\mathrm{╥}}}\right),\left(2-\alpha \right)\left(\left(x^2-x\right)+4\frac{\sqrt{x}}{\sqrt{\mathrm{╥}}}\right)\right]\\ \mathrm{ß}\left({\tilde{y}}_0\left(x;\alpha \right),\frac{\partial {\left[{\tilde{y}}_0\right]}_{\alpha }}{\partial x}\right)=0\end{array}\right.$$

(37)

First-order problem:

$$\left\{\begin{array}{c}{\tilde{y}}_1\left(x,{\tilde{S}}_1\left(\alpha \right);\alpha \right)={\left(1+{\tilde{S}}_1\left(\alpha \right)\right)\tilde{y}}_0\left(x;\alpha \right)+{\tilde{S}}_1\left(\alpha \right){\tilde{\mathcal{J}}}^{\left(1.5\right)}{\tilde{y}}_0\left(x;\alpha \right)\\ -\left(1+{\tilde{S}}_1\left(\alpha \right)\right){\tilde{\mathcal{J}}}^{\left(1.5\right)}\left(\tilde{F}\left(x;\alpha \right)\right),\\ \mathrm{ß}\left({\tilde{y}}_1\left(x;\alpha \right),\frac{\partial {\left[{\tilde{y}}_1\right]}_{\alpha }}{\partial x}\right)=0.\end{array}\right.$$

(38)

Second-order problem:

$$\left\{\begin{array}{c}{\tilde{y}}_2\left(x,{\tilde{S}}_1\left(\alpha \right),{\tilde{S}}_2\left(\alpha \right);\alpha \right)={\left(1+{\tilde{S}}_1\left(\alpha \right)\right)\tilde{y}}_1\left(x,{\tilde{S}}_1\left(\alpha \right);\alpha \right)+{\tilde{S}}_2\left(\alpha \right){\tilde{y}}_0\left(x;\alpha \right)+\\ {\tilde{S}}_1\left(\alpha \right){\tilde{\mathcal{J}}}^{\left(1.5\right)}{\tilde{y}}_1\left(x,{\tilde{S}}_1\left(\alpha \right);\alpha \right)+{\tilde{S}}_2\left(\alpha \right){\tilde{\mathcal{J}}}^{\left(1.5\right)}{\tilde{y}}_0\left(x;\alpha \right)-{\tilde{S}}_2\left(\alpha \right){\tilde{\mathcal{J}}}^{\left(1.5\right)}\left(\tilde{F}\left(x;\alpha \right)\right)\\ \mathrm{ß}\left({\tilde{y}}_2\left(x;\alpha \right),\frac{\partial {\left[{\tilde{y}}_2\right]}_{\alpha }}{\partial x}\right)=0.\end{array}\right.$$

(39)

Third-order problem:

$$\left\{\begin{array}{c}{\tilde{y}}_3\left(x,{\tilde{S}}_1\left(\alpha \right),{\tilde{S}}_2\left(\alpha \right),{\tilde{S}}_3\left(\alpha \right);\alpha \right)={\left(1+{\tilde{S}}_1\left(\alpha \right)\right)\tilde{y}}_2\left(x,{\tilde{S}}_1\left(\alpha \right),{\tilde{S}}_2\left(\alpha \right);\alpha \right)+\\ {\tilde{S}}_2\left(\alpha \right){\tilde{y}}_1\left(x,{\tilde{S}}_1\left(\alpha \right);\alpha \right)+{\tilde{\mathcal{J}}}^{\left(1.5\right)}{\displaystyle \sum _{j=1}^3}{\tilde{S}}_j\left(\alpha \right){\tilde{y}}_{3-j}\left(x,{\tilde{S}}_1\left(\alpha \right),\dots ,{\tilde{S}}_{3-j}\left(\alpha \right);\alpha \right)\\ -{\tilde{S}}_3\left(\alpha \right){\tilde{\mathcal{J}}}^{\left(1.5\right)}\left(\tilde{F}\left(x;\alpha \right)\right)+{\tilde{S}}_3\left(\alpha \right){\tilde{y}}_0\left(x;\alpha \right)\\ \mathrm{ß}\left({\tilde{y}}_3\left(x;\alpha \right),\frac{\partial {\left[{\tilde{y}}_3\right]}_{\alpha }}{\partial x}\right)=0.\end{array}\right.$$

(40)

Fourth-order problem:

$$\left\{\begin{array}{c}{\tilde{y}}_4\left(x,{\tilde{S}}_1\left(\alpha \right),\dots ,{\tilde{S}}_4\left(\alpha \right);\alpha \right)={\left(1+{\tilde{S}}_1\left(\alpha \right)\right)\tilde{y}}_3\left(x,{\tilde{S}}_1\left(\alpha \right),{\tilde{S}}_2\left(\alpha \right),{\tilde{S}}_3\left(\alpha \right);\alpha \right)+\\ {\tilde{S}}_2\left(\alpha \right){\tilde{y}}_2\left(x,{\tilde{S}}_1\left(\alpha \right),{\tilde{S}}_2\left(\alpha \right);\alpha \right)+{\tilde{S}}_3\left(\alpha \right){\tilde{y}}_1\left(x,{\tilde{S}}_1\left(\alpha \right);\alpha \right)+{\tilde{S}}_4\left(\alpha \right){\tilde{y}}_0\left(x;\alpha \right)+\\ {\tilde{\mathcal{J}}}^{\left(1.5\right)}{\displaystyle \sum _{j=1}^4}{\tilde{S}}_j\left(\alpha \right){\tilde{y}}_{4-j}\left(x,{\tilde{S}}_1\left(\alpha \right),\dots ,{\tilde{S}}_{4-i}\left(\alpha \right);\alpha \right)-{\tilde{S}}_4\left(\alpha \right){\tilde{\mathcal{J}}}^{\left(1.5\right)}\left(\tilde{F}\left(x;\alpha \right)\right)\\ \mathrm{ß}\left({\tilde{y}}_4\left(x;\alpha \right),\frac{\partial {\left[{\tilde{y}}_4\right]}_{\alpha }}{\partial x}\right)=0\end{array}\right.$$

(41)

Fifth-order problem:

$$\left\{\begin{array}{c}{\tilde{y}}_5\left(x,{\tilde{S}}_1\left(\alpha \right),\dots ,{\tilde{S}}_5\left(\alpha \right);\alpha \right)={\left(1+{\tilde{S}}_1\left(\alpha \right)\right)\tilde{y}}_4\left(x,{\tilde{S}}_1\left(\alpha \right),\dots ,{\tilde{S}}_4\left(\alpha \right);\alpha \right)+\\ {\tilde{S}}_2\left(\alpha \right){\tilde{y}}_3\left(x,{\tilde{S}}_1\left(\alpha \right),{\tilde{S}}_2\left(\alpha \right),{\tilde{S}}_3\left(\alpha \right);\alpha \right)+{\tilde{S}}_3\left(\alpha \right){\tilde{y}}_2\left(x,{\tilde{S}}_1\left(\alpha \right),{\tilde{S}}_2\left(\alpha \right);\alpha \right)+\\ {\tilde{S}}_4\left(\alpha \right){\tilde{y}}_1\left(x,{\tilde{S}}_1\left(\alpha \right);\alpha \right)+{\tilde{S}}_5\left(\alpha \right){\tilde{y}}_0\left(x;\alpha \right)-{\tilde{S}}_5\left(\alpha \right){\tilde{\mathcal{J}}}^{\left(1.5\right)}\left(\tilde{F}\left(x;\alpha \right)\right)+\\ {\tilde{\mathcal{J}}}^{\left(1.5\right)}{\displaystyle \sum _{j=1}^5}{\tilde{S}}_j\left(\alpha \right){\tilde{y}}_{5-j}\left(x,{\tilde{S}}_1\left(\alpha \right),\dots ,{\tilde{S}}_{5-i}\left(\alpha \right);\alpha \right)\\ \mathrm{ß}\left({\tilde{y}}_5\left(x;\alpha \right),\frac{\partial {\left[{\tilde{y}}_5\right]}_{\alpha }}{\partial x}\right)=0.\end{array}\right.$$

(42)

Next, we will solve Equation (31) with a third-order series FF-OHAM using the Mathematica 13 Dsolve package

$$\tilde{y}\left(x;\alpha \right)={\tilde{y}}_0\left(x;\alpha \right)+{\sum }_{j=1}^3{\tilde{y}}_j\left(x{,\tilde{S}}_1\left(0.5\right),\dots {,\tilde{S}}_j\left(0.5\right);0.5.\right)$$

(43)

For this linear application, we found that the fuzzy convergence parameters at $\alpha =0.5$ provide an appropriate and accurate series solution at each $\alpha \in \left[\mathrm{0,1}\right]$, such that ${\tilde{S}}_1\left(0.5\right)=-1.065291064957493$, ${\tilde{S}}_2\left(0.5\right)=-0.00004338985509905724$, and ${\tilde{S}}_3\left(0.5\right)=-0.0020739269082325523$.

Next, we will employ the fuzzy convergence parameters ${\tilde{S}}_1\left(0.5\right),{\tilde{S}}_2\left(0.5\right),$ and ${\tilde{S}}_3\left(0.5\right)$ in Equation (43) to find the third-order FF-OHAM approximate series solution for Equation (31), as shown in

Table 1. The approximate solutions and errors for Equation (31) using the third-order FF-OHAM at $\mathrm{x}=0.5$ for all $\mathsf{\alpha }\in \left[\mathrm{0,1}\right]$.

$\mathit{\alpha }$	${\left[\underset{¯}{\mathit{E}\mathit{R}}\right]}_{\mathit{\alpha }},{\tilde{\mathit{S}}}_{\mathit{j}}$	${\left[\overline{\mathit{E}\mathit{R}}\right]}_{\mathit{\alpha }},{\tilde{\mathit{S}}}_{\mathit{j}}$	${\left[\underset{¯}{\mathit{y}}\right]}_{\mathit{\alpha }},{\tilde{\mathit{S}}}_{\mathit{j}}$	${\left[\overline{\mathit{y}}\right]}_{\mathit{\alpha }},{\tilde{\mathit{S}}}_{\mathit{j}}$
0	$-1.36061\times {10}^{-5}$	$0$	$-0.49998$	$0$
0.2	$-1.22455\times {10}^{-5}$	$-1.36061\times {10}^{-6}$	$-0.44998$	$-0.04999$
0.4	$-1.08849\times {10}^{-5}$	$-2.72123\times {10}^{-6}$	$-0.39999$	$-0.09999$
0.6	$-9.52429\times {10}^{-6}$	$-4.08184\times {10}^{-6}$	$-0.34999$	$-0.14999$
0.8	$-8.16368\times {10}^{-6}$	$-5.44245\times {10}^{-6}$	$-0.29999$	$-0.19999$
1	$-6.80306\times {10}^{-6}$	$-6.80306\times {10}^{-6}$	$-0.24999$	$-0.24999$

, as follows.

Using a three-dimensional graph, we summarize the solutions using the third-order FF-OHAM over all $x\in \left[\mathrm{0,0.5}\right]$ and $\alpha \in \left[\mathrm{0,1}\right]$ corresponding with the best optimal convergence control values—${\tilde{S}}_1\left(0.5\right),{\tilde{S}}_2\left(0.5\right),$ and ${\tilde{S}}_3\left(0.5\right)$—of Equation (31) in Figure 1.

To analyze the behavior of FF-OHAM for solving second-order FFOBVPs, we shall proceed to solve Equation (31) using the same data, $x\in \left[\mathrm{0,0.5}\right]$ and $\omega =1.5$, and the fifth-order FF-OHAM instead of the third-order FF-OHAM to illustrate the convergence dynamic of FF-OHAM for different terms of the approximate series solution; therefore, the series solution will take the following form:

$$\tilde{y}\left(x;\alpha \right)={\tilde{y}}_0\left(x;\alpha \right)+{\sum }_{j=1}^5{\tilde{y}}_j\left(x{,\tilde{S}}_1\left(0.5\right),\dots {,\tilde{S}}_j\left(0.5\right);0.5.\right)$$

(44)

such that the optimal convergence control parameters calculated using the Mathematica 13 Dsolve package are

${\tilde{S}}_1\left(0.5\right)=-1.0270653590282228$	${\tilde{S}}_2\left(0.5\right)=6.734609572815013\times {10}^{-7}$	${\tilde{S}}_3\left(0.5\right)=-0.000025083072228706767$
${\tilde{S}}_4\left(0.5\right)=0.000131089732250741$	${\tilde{S}}_5\left(0.5\right)=0.0001858737255737089$

The above convergence parameters will be employed in Equation (44) to find the fifth-order FF-OHAM approximate series solution for Equation (31), as shown in

Table 2. The approximate solutions and errors for Equation (31) using the fifth-order FF-OHAM at $\mathrm{x}=0.5$ for all $\mathsf{\alpha }\in \left[\mathrm{0,1}\right]$.

$\mathit{\alpha }$	${\left[\underset{¯}{\mathit{E}\mathit{R}}\right]}_{\mathit{\alpha }},{\tilde{\mathit{S}}}_{\mathit{j}}\left(0.5\right)$	${\left[\overline{\mathit{E}\mathit{R}}\right]}_{\mathit{\alpha }},{\tilde{\mathit{S}}}_{\mathit{j}}\left(0.5\right)$	${\left[\underset{¯}{\mathit{y}}\right]}_{\mathit{\alpha }},{\tilde{\mathit{S}}}_{\mathit{j}}\left(0.5\right)$	${\left[\overline{\mathit{y}}\right]}_{\mathit{\alpha }},{\tilde{\mathit{S}}}_{\mathit{j}}\left(0.5\right)$
0	$-2.79607\times {10}^{-8}$	$0$	$-0.50000$	$0$
0.2	$-2.51647\times {10}^{-8}$	$-2.79607\times {10}^{-9}$	$-0.45000$	$-0.05000$
0.4	$-2.23686\times {10}^{-8}$	$-5.59214\times {10}^{-9}$	$-0.40000$	$-0.10000$
0.6	$-1.95725\times {10}^{-8}$	$-8.38822\times {10}^{-9}$	$-0.35000$	$-0.15000$
0.8	$-1.67764\times {10}^{-8}$	$-1.11843\times {10}^{-8}$	$-0.30000$	$-0.20000$
1	$-1.39804\times {10}^{-8}$	$-1.39804\times {10}^{-8}$	$-0.25000$	$-0.25000$

, as follows.

Figure 2 illustrates the summary of the solutions using the fifth-order FF-OHAM over all $x\in \left[\mathrm{0,0.5}\right]$ and $\alpha \in \left[\mathrm{0,1}\right]$ corresponding with the best optimal convergence control values—${\tilde{S}}_1\left(0.5\right),{\tilde{S}}_2\left(0.5\right),{\tilde{S}}_3\left(0.5\right),{\tilde{S}}_4\left(0.5\right),$ and ${\tilde{S}}_5\left(0.5\right)$—of Equation (43) in a three-dimensional figure.

Table 1 and Table 2 and Figure 1 and Figure 2 illustrate that the third- and fifth-order FF-OHAMs satisfy the triangular solution of the fuzzy differential equations for Equation (31) [2]. On the other hand, we can conclude that the series solution of the linear physical application involving FFOBVP using the FF-OHAM will approach the exact solutions whenever the order of the FF-OHAM increases. The developed FF-OHAM is compared with the spectral collection method (SCM) for solving the mechanical application described in Equation (31). Figure 3 and Figure 4 illustrate the lower and upper accuracy of the fifth-order FF-OHAM in comparison to the fifth-order SCM for solving the mechanical pplication $\forall x\in \left[\mathrm{0,1}\right]$ at $\alpha =0.5$ based on the absolute error defined below in Equation (45).

$$\left\{\begin{array}{c}\underset{¯}{ERR}=\left|\underset{¯}{Y}\left(x;\alpha \right)-\underset{¯}{y}\left(x;\alpha \right)\right|\\ \overline{ERR}=\left|\overline{Y}\left(x;\alpha \right)-\overline{y}\left(x;\alpha \right)\right|\end{array},\forall x\in \left[\mathrm{0,1}\right]\right.,\mathrm{and} \forall \alpha \in \left[\mathrm{0,1}\right]$$

(45)

We can conclude from Figure 3 and Figure 4 that the accuracy of the approximate solution solved for using the fifth-order FF-OHAM series provides better accuracy compared to SCM for all $x\in \left[\mathrm{0,1}\right]$.

• Thermal Conductivity of a Material: Nonlinear Fractional Temperature Distribution Equation

Consider the mathematical model, a nonlinear fractional temperature distribution equation of order $\omega \in \left(\mathrm{1,2}\right]$, which explains the distribution of the temperature in the lumped convection system in a layer comprised of materials with varying thermal conductivity [39]:

$$\left\{\begin{array}{c}{D}^{\left(\omega \right)}y\left(x\right)-\eta {\left(y\left(x\right)\right)}^4=0,\hspace{1em}x\in \left[\mathrm{0,1}\right]\\ y^{\prime }\left(0\right)=0,\hspace{1em}y\left(1\right)=1.\end{array}\right.$$

(46)

where $x$ is the time-independent variable, and $y\left(x\right)$ is the dimensionless temperature.

The following is the fuzzy version of Equation (46):

$$\left\{\begin{array}{c}{\tilde{D}}^{\left(\omega \right)}\tilde{y}\left(x;\alpha \right)-\eta {\left(y\left(x;\alpha \right)\right)}^4=0,\hspace{1em}x\in \left[\mathrm{0,1}\right]\\ {\underset{¯}{y}}^{\prime }\left(0;\alpha \right)=\left(0.1\alpha -0.1\right),\hspace{1em}\underset{¯}{y}\left(1;\alpha \right)=\left(0.1\alpha +0.9\right),\\ {\overline{y}}^{\prime }\left(0;\alpha \right)=\left(0.1-0.1\alpha \right),\hspace{1em}\overline{y}\left(1;\alpha \right)=\left(1.1-0.1\alpha \right).\end{array}\right.$$

(47)

To solve the fuzzy fractional model of the thermal conductivity using the FF-OHAM, sccording to Section 3, we can build the approximate series solution for Equation (47) of order $\omega \in \left(\mathrm{1,2}\right]$ for all $\alpha \in \left[\mathrm{0,1}\right]$ as follows:

For $k\ge 0$, we can construct the following FF-OHAM form

$$\left(1-q\right)\left[{\tilde{D}}^{\left(\omega \right)}\left(\tilde{y}\left(x;\alpha \right)\right)\right]={\displaystyle \sum _{j=1}^k}{\underset{¯}{S}}_j{\left(\alpha \right)q}^j\left[{\tilde{D}}^{\left(\omega \right)}\left(\tilde{y}\left(x;\alpha \right)\right)-\eta {\left(\tilde{y}\left(x;\alpha \right)\right)}^4\right],$$

(48)

Then, the approximate series solution is introduced in Equation (49) below:

$$\tilde{y}\left(x;\alpha \right)={\tilde{y}}_0(x;\alpha )+{\displaystyle \sum _{j=1}^k}{\tilde{y}}_j(x,S_1,\dots ,S_j;\alpha )q^j$$

(49)

For a zeroth-order problem:

$${\tilde{y}}_0\left(x;\alpha \right)=\left[\tilde{0}\right].$$

(50)

For first- to tenth-order problems:

$$\left\{\begin{array}{c}\left(1-q\right)\left[{\tilde{D}}^{\left(\omega \right)}\left(\tilde{y}\left(x;\alpha \right)\right)\right]={\displaystyle \sum _{j=0}^9}{\tilde{S}}_j{\left(\alpha \right)q}^j[{\tilde{D}}^{\left(\omega \right)}\left(\tilde{y}\left(x;\alpha \right)\right)-\\ \eta \left({\displaystyle \sum _{i=0}^9}{\tilde{y}}_{9-i}\left(x;\alpha \right){\displaystyle \sum _{j=0}^i}{\tilde{y}}_{i-j}\left(x;\alpha \right){\displaystyle \sum _{s=0}^j}{\tilde{y}}_s\left(x;\alpha \right){\tilde{y}}_{j-s}\left(x;\alpha \right)\right)],\\ {{\tilde{y}}_k}^{\prime }\left(0;\alpha \right)={\tilde{y}}_k\left(1;\alpha \right)=\tilde{0}.\end{array}\right.$$

(51)

Next, using the Mathematica 13 Dsolve package to find the series solutions for the lower and the upper bounds of Equation (47), for $j=\mathrm{1,2},\dots ,10$, we obtain

$$\tilde{y}\left(x;\alpha \right)={\tilde{y}}_0(x;\alpha )+{\sum }_{j=1}^{10}{\tilde{y}}_j(x,S_1,\dots ,S_j;\alpha )q^j$$

(52)

such that the optimal lower and upper convergence control parameters calculated and coded using Mathematica 13 to find the most accurate solution for Equation (47) via the tenth-order FF-OHAM are listed in

Table 3. Lower auxiliary convergence parameters of the tenth-order FF-OHAM for solving Equation (47) at $\mathsf{\omega }=1.9$, $\mathrm{x}=0.1,$ for all $\mathsf{\alpha }\in \left[\mathrm{0,1}\right]$.

$\mathit{\alpha }$	${{\displaystyle \underset{¯}{\mathit{S}}}}_{\mathit{j}}$
0	${\underset{¯}{S}}_1=-0.42754243116098856$	${\underset{¯}{S}}_2=0.10029988098228566$	${\underset{¯}{S}}_3=-0.13497958972180282$
	${\underset{¯}{S}}_4=0.057027658640051423$	${\underset{¯}{S}}_5=-0.026736473824494192$	${\underset{¯}{S}}_6=-0.09334881995323582$
	${\underset{¯}{S}}_7=0.14476539458261192$	${\underset{¯}{S}}_8=2.6782513149178886$	${\underset{¯}{S}}_9=-3.1547021023298236$
		${\underset{¯}{S}}_{10}=0$
0.5	${\underset{¯}{S}}_1=-0.4275424758687012$	${\underset{¯}{S}}_2=0.10581949295995502$	${\underset{¯}{S}}_3=-0.1375566084081214$
	${\underset{¯}{S}}_4=0.07545735284696302$	${\underset{¯}{S}}_5=-0.049825680474098175$	${\underset{¯}{S}}_6=-0.21083114235974942$
	${\underset{¯}{S}}_7=0.33623100353843177$	${\underset{¯}{S}}_8=-1.780319356305217$	${\underset{¯}{S}}_9=1.8279961729200187$
		${\underset{¯}{S}}_{10}=0$
1	${\underset{¯}{S}}_1=-0.4275424759052609$	${\underset{¯}{S}}_2=0.10838719398757574$	${\underset{¯}{S}}_3=-0.1393691078519383$
	${\underset{¯}{S}}_4=0.0819534530845722$	${\underset{¯}{S}}_5=-0.0540$4354434171497<	${\underset{¯}{S}}_6=-0.06431130239382181$
	${\underset{¯}{S}}_7=0.12107016205095011$	${\underset{¯}{S}}_8=-0.9412597479377718$	${\underset{¯}{S}}_9=1.0018880129871244$
		${\underset{¯}{S}}_{10}=0$

and

Table 4. Upper auxiliary convergence parameters of the tenth-order FF-OHAM for solving Equation (47) at $\mathsf{\omega }=1.9$, $\mathrm{x}=0.1,$ for all $\mathsf{\alpha }\in \left[\mathrm{0,1}\right]$.

$\mathit{\alpha }$	${\overline{\mathit{S}}}_{\mathit{j}}$
0	${\overline{S}}_1=-0.4499999715464292$	${\overline{S}}_2=0.07033774093204971$	${\overline{S}}_3=-0.08588077852784994$
	${\overline{S}}_4=0.02200367504397761$	${\overline{S}}_5=-0.013738000249538687$	${\overline{S}}_6=-0.13177749373901043$
	${\overline{S}}_7=0.20247748961630552$	${\overline{S}}_8=-1.167757698543094$	${\overline{S}}_9=1.1715997608744646$
		${\overline{S}}_{10}=0$
0.5	${\overline{S}}_1=-0.4500000000118868$	${\overline{S}}_2=0.11085869384502318$	${\overline{S}}_3=-0.1322901831399701$
	${\overline{S}}_4=0.053720614964490584$	${\overline{S}}_5=-0.019748829952898384$	${\overline{S}}_6=-0.12141795332965101$
	${\overline{S}}_7=0.1700603758070842$	${\overline{S}}_8=-0.7863871032395353$	${\overline{S}}_9=0.7686499426055967$
		${\overline{S}}_{10}=0$
1	${\overline{S}}_1=-0.4499999999893346$	${\overline{S}}_2=0.1095056263140715$	${\overline{S}}_3=-0.1302513333365527$
	${\overline{S}}_4=0.0473415579627094$	${\overline{S}}_5=-0.0185338348210904$	${\overline{S}}_6=-0.14589007991813402$
	${\overline{S}}_7=0.2076338115632854$	${\overline{S}}_8=-1.0720018208884141$	${\overline{S}}_9=1.0605759385659808$
		${\overline{S}}_{10}=0$

below.

The above lower and upper convergence parameters in Table 3 and Table 4 bare employed in Equation (52) to find the tenth-order FF-OHAM approximate series solution for Equation (47), as shown in

Table 5. The approximate solutions and errors for Equation (47) using the tenth-order FF-OHAM when $\mathsf{\omega }=1.9$ at $\mathrm{x}=0.1$ for all $\mathsf{\alpha }\in \left[\mathrm{0,1}\right]$.

$\mathit{\alpha }$	${\left[\underset{¯}{\mathit{E}\mathit{R}}\right]}_{\mathit{\alpha }},{\underset{¯}{\mathit{S}}}_{\mathit{j}}$	${\left[\overline{\mathit{E}\mathit{R}}\right]}_{\mathit{\alpha }},{\overline{\mathit{S}}}_{\mathit{j}}$	${\left[\underset{¯}{\mathit{y}}\right]}_{\mathit{\alpha }},{\underset{¯}{\mathit{S}}}_{\mathit{j}}$	${\left[\overline{\mathit{y}}\right]}_{\mathit{\alpha }},{\overline{\mathit{S}}}_{\mathit{j}}$
0	$-7.73663\times {10}^{-5}$	$-1.78658\times {10}^{-5}$	$0.81730$	$0.83216$
0.5	$-6.65123\times {10}^{-6}$	$1.68215\times {10}^{-6}$	$0.82167$	$0.82909$
1	$6.42966\times {10}^{-5}$	$-2.62733\times {10}^{-5}$	$0.82560$	$0.82560$

and summarized Figure 5 below.

Figure 5 illustrates the summary of the solutions using the tenth-order FF-OHAM over all $x\in \left[\mathrm{0,0.1}\right]$ and $\alpha \in \left[\mathrm{0,1}\right]$ corresponding with the best optimal convergence control values, ${\tilde{S}}_j$, of Equation (47) in the three-dimensional graph.

Morever, the fuzzy solutions, shown in Table 5 and Figure 5, clarify that the new construction of the FF-OHAM satisfies the fuzzy solution of the new fuzzy version of the distribution of the model of the temperature in the lumped convection system in a layer comprised of materials with varying thermal conductivity. Furthermore, the FF-OHAM provides an appropriate approximate series solution for the strong nonlinearity fractional differential equation with the presence of the uncertainty, which makes this approach applicable and suitable for solving the most complicated, nonlinear real-world problems.

7. Conclusions

The present study focused on developing an approximate analytical method, called FF-OHAM, for solving linear and nonlinear fractional order boundary value problems (FFOBVPs). The FF-OHAM method has the ability to control the convergence of the series solution by selecting the optimal convergence parameter for each method. The Bagley–Torvik equation, which is an inhomogeneous linear FFOBVP, was used as a case study to demonstrate the accuracy of the FF-OHAM method in solving linear cases. The method was found to provide an accurate series solution as the series order increases and the obtained solution converges to the exact solution. The solutions obtained using the FF-OHAM were found to be more accurate than those obtained using the SCM. Furthermore, the study also introduced a new fuzzy version of the fractional temperature distribution equation and utilized the FF-OHAM to find the series solution for this nonlinear problem. The FF-OHAM method was found to provide an accurate series solution for solving nonlinear cases without needing an exact solution. The convergence dynamic of the FF-OHAM was also used to obtain optimal convergence parameters for this problem. Finally, it is noted that all the fuzzy fractional solutions obtained using the FF-HAM and the FF-OHAM satisfy the triangular fuzzy solution, which is a desirable property for fuzzy systems. It is a good idea to explore techniques to improve the computational efficiency of the developed FF-OHAM method, as this can lead to faster and more efficient solutions for FFOBVPs, such as parallelization. The FF-OHAM method can be parallelized to run on multiple processors or cores simultaneously. This can help to reduce the computational time required to obtain a solution, especially for large and complex problems.