Approximate Solutions of Fractional Differential Equations Using Optimal q-Homotopy Analysis Method: A Case Study of Abel Differential Equations

Abstract

In this study, the optimal q-Homotopy Analysis Method (optimal q-HAM) has been used to investigate fractional Abel differential equations. This article is designed as a case study, where several forms of Abel equations, containing Bernoulli and Riccati equations, are given with ordinary derivatives and fractional derivatives in the Caputo sense to present the application of the method. The optimal q-HAM is an improved version of the Homotopy Analysis Method (HAM) and its modification q-HAM and focuses on finding the optimal value of the convergence parameters for a better approximation. Numerical applications are given where optimal values of the convergence control parameters are found. Additionally, the correspondence of the approximate solutions obtained for these optimal values and the exact or numerical solutions are shown with figures and tables. The results show that the optimal q-HAM improves the convergence of the approximate solutions obtained with the q-HAM. Approximate solutions obtained with the fractional Differential Transform Method, q-HAM and predictor–corrector method are also used to highlight the superiority of the optimal q-HAM. Analysis of the results from various methods points out that optimal q-HAM is a strong tool for the analysis of the approximate analytical solution in Abel-type differential equations. This approach can be used to analyze other fractional differential equations arising in mathematical investigations.

1. Introduction

Mathematical modeling has been instrumental for people in a variety of ways for centuries. The roots of today’s modeling studies in economy, medicine and engineering can be traced back to the first applications of algebraic equations for analyzing inheritance and trade problems. Recent studies in modeling focus on the use of differential equations using the tools of calculus to analyze changes in certain components of various systems. In most of these studies, fractional calculus has become the center of attention.

Certain events or systems in nature can be modeled more efficiently using non-integer order derivation. Dating back to the letter of Leibniz in the late 17th century, studies on fractional calculus improved constantly with a certain increase in the numbers in the last few decades. Riemann–Liouville, Grünwald–Letnikov and Caputo fractional derivative definitions have been extensively studied in the literature. More recent definitions, such as the Atangana–Baleanu fractional derivative and Caputo–Fabrizio fractional derivative, are also popular tools in modeling studies. Katugampola, Hilfer and many other derivatives have been given in the literature (see reference [1]) along with conformable fractional derivative which has become the basis of conformable calculus [2,3].

A certain number of studies on fractional calculus focus on numerical or approximate analytical approaches rather than analytical techniques. Some of the well-known methods such as the Differential Transformation Method (DTM) and the Variational Iteration Method (VIM) have fractional counterparts that are frequently used in applications [4,5]. He’s Homotopy Perturbation Method [6] and the Homotopy Analysis Method proposed by Liao [7] share a common relation with the homotopy concept in topology [8] and are both among the most studied methods in the literature. Liao [8] showed that the Homotopy Perturbation Method (HPM) is a special case of the Homotopy Analysis Method (HAM), and the former has many applications in fractional calculus. HAM and its modifications such as the optimal q-HAM [9], which is based on the improvement of the convergence by using auxiliary convergence parameters, have been used in the literature for several applications. Some of these applications include the investigation of heat transfer [10], option pricing [11], the convection–diffusion equation [12], the fractional order logistic equation [13] and non-Newtonian fluid flow [14] in their ordinary and fractional forms. A more general form of HAM, the q-Homotopy Analysis Method (q-HAM), has been proposed by El-Tawil and Huseen and the method uses two parameters to achieve faster convergence in comparison to HAM [15,16]. The q-HAM and the q-HATM, which is a combination of q-HAM and the Laplace transform method, have also been applied to many modeling problems such as Burger’s equation [17], non-Darcy Flow problem [18], the fractional vibration equation [19] and the fractional heat equation [20]. Recent studies on HAM, q-HAM and their modifications focus mostly on fractional applications such as the fractional KdV equation [21], the fractional Sawada-Kotera equation [22] and the fractional Fisher’s equation [23]. Some other examples of the most recent studies for the method can be exemplified as follows. Biswas and Ghosh have used the q-HAM to analyze the time-fractional Harry Dym equation [24]. Hussein et al. have analyzed the Cahn-Hilliard equation using the q-HAM [25]. Cheng et al. have used the q-HAM to solve the time-fractional Keller–Segel-type equations. In the mentioned study, the authors have also performed a symmetry analysis alongside the implementation of the q-HAM [26]. Sunita et al. have used the q-HAM with Elzaki transform to investigate the two-dimensional solute transform problem [27]. One modification of the q-HAM, called the optimal q-HAM (or Oq-HAM), focuses on finding the optimal values for the auxiliary parameter. This method uses the optimal value for the parameter, controlling the rate and region of convergence, and hence offers a better approximation to the solution. This method has been recently used to investigate the Kaup–Kupershmidt equation [28] and various fractional partial differential equations [29]. Another problem that can be investigated using the fractional optimal q-HAM is the Abel differential equation.

The Abel equation is named after the Norwegian mathematician Niels Henrik Abel who lived in the 19th century. The Abel differential equation of the first kind is a generalization of the Riccati equation, and the Abel differential equation of the second kind is a further generalization [30]. Abel equations have been used to describe the relativistic evolution of a causal dissipative cosmological fluid in conformally flat space–time [31], magnetostatic problems [32] or inflationary dynamics [33]. Abel equations are ordinary differential equations and several studies on their general or periodic solutions have been given in the literature (see references [34,35]). These equations have also been studied recently in the fractional sense using various fractional derivatives. The use of the short memory principle for the solution of the fractional equation [36], a numerical analysis of the equation with a Caputo–Fabrizio derivative [37], a numerical approach using Genocchi polynomials for the fractional model with a Caputo derivative [38] and a numerical investigation of the fractional Abel equation using generalized Bessel functions [39] are some of the recent fractional studies.

In this study, optimal q-HAM will be used to investigate the fractional Abel equation in the Caputo sense. The motivation of this study is the existing literature on the HAM method and the improvements in the convergence achieved by the Oq-HAM modification of the scheme. In accordance with the focus of the recent literature on mathematical modeling using tools of fractional calculus, it is aimed at achieving such improvement for Abel differential equations within the fractional framework. Although the Caputo fractional derivative has been used for this study, other definitions of fractional derivation could also be employed, whenever it is checked that they satisfy the necessary conditions to apply these methods. This approach enables an investigation taking advantage of non-integer order derivation, such as a more generalized investigation and the ability to provide a more accurate modeling of real phenomena. This article adds to the current literature on fractional Abel equations, such as the study by Jafari et al. which uses HAM to demonstrate numerical results [40]. Several numerical examples are given with a variety of coefficients for Abel equations of the first type or equations that are reduced to Bernoulli or Riccati equations in the ordinary and fractional forms. Exact solutions of the equations are, whenever available, compared with the results from optimal q-HAM to demonstrate that the method is a valuable tool to analyze the fractional Abel differential equation. It should be noted that, as mentioned above, optimal q-HAM is a well-established method, and this manuscript is structured as a case study for the analysis of fractional Abel equations using the method. This case study concentrates on the improvement of the approximation obtained with q-HAM through determining the optimal values of the auxiliary parameter. A comprehensive comparison of the improved approximation is presented using results from various methods to verify the benefits of the application of optimal q-HAM in the case of fractional Abel equations. Key contributions of this study can be given as the following bullet points for a clearer presentation:

-

A fractional optimal q-HAM has been presented and the approximate solutions for several ordinary and fractional Abel differential equations have been given.

-

The approximate solutions have been graphically presented for examples and comparisons of the approximations have been made with exact or numerical solutions.

-

The optimal values of h and the convergence regions have been given for various selections of the auxiliary parameters.

-

Relative errors, h-curves and numerical comparisons of solutions have been given to underline the suitability of the method and the improvements achieved by using Oq-HAM.

-

Approximate analytical solutions from fractional DTM and q-HAM have been compared with solutions obtained with Oq-HAM to underline the superiority of the method.

-

Results from q-HAM have been used to point out that a better approximation is achieved as the value of the auxiliary parameter h approaches its optimal value, verifying the improvement through Oq-HAM.

2. Abel Differential Equation of the First Kind and Its Analysis with Optimal q-HAM

The Abel differential equation of the first kind is given by the following:

$$y^{\prime }\left(x\right)=f_3\left(x\right)y^3\left(x\right)+f_2\left(x\right)y^2\left(x\right)+f_1\left(x\right)y\left(x\right)+f_0\left(x\right)$$

(1)

where $f_3\left(x\right)\ne 0$, $f_2\left(x\right)$, $f_1\left(x\right)$ and $f_0\left(x\right)$ are meromorphic functions. A Riccati equation will be obtained in the case where $f_3\left(x\right)=0,$ and a Bernoulli equation can also be obtained in the case where $f_0\left(x\right)=0$ and $f_2\left(x\right)=0$ or $f_3\left(x\right)=0$. An Abel equation of the second kind can be given by the following:

$$\left[g_0\left(x\right)+g_1\left(x\right)y\left(x\right)\right]y^{\prime }\left(x\right)=f_3\left(x\right)y^3\left(x\right)+f_2\left(x\right)y^2\left(x\right)+f_1\left(x\right)y\left(x\right)+f_0\left(x\right)$$

which is a generalization of the equation of the first kind. The equation of the second kind reduces to the equation of the first kind for $g_0\left(x\right)=1 \mathrm{a}\mathrm{n}\mathrm{d} g_1\left(x\right)=0$ [30].

In order to apply q-HAM to analyze the Abel equation of the first kind, Equation (1) is rewritten as follows:

$$y^{\prime }\left(x\right)-f_3\left(x\right)y^3\left(x\right)-f_2\left(x\right)y^2\left(x\right)-f_1\left(x\right)y\left(x\right)-f_0\left(x\right)=0$$

(2)

or

$$N\left[y\left(x\right)\right]=0$$

(3)

where $N\left[y\left(x\right)\right]={\displaystyle \frac{d}{dx}}y\left(x\right)-f_3\left(x\right)y^3\left(x\right)-f_2\left(x\right)y^2\left(x\right)-f_1\left(x\right)y\left(x\right)-f_0\left(x\right)$ is the nonlinear operator. The zero-order deformation of q-HAM is given as follows:

$$\left(1-nq\right)L\left[\varphi \left(x;q\right)-y\left(0\right)\right]=qhH\left(x\right)N\left[\varphi \left(x;q\right)\right]$$

(4)

where $0\le q\le {\displaystyle \frac{1}{n}},n\ge 1$ is the embedded parameter, $L\left[\varphi \left(x;q\right)\right]={\displaystyle \frac{d}{dx}}\varphi (x)$ is the linear operator, $h\ne 0$ is an auxiliary parameter and $H\left(x\right)\ne 0$ is an auxiliary function [15]. Equation (4) becomes the following:

$$\varphi \left(x;0\right)=y\left(0\right)$$

for $q=0$ and similarly the following:

$$\varphi \left(x;{\displaystyle \frac{1}{n}}\right)=y\left(x\right)$$

for $q={\displaystyle \frac{1}{n}}$. Hence, the solution $\varphi \left(x;q\right)$ approaches the solution $y\left(x\right)$ from the initial condition $y\left(0\right)$ as $q$ changes from $0$ to ${\displaystyle \frac{1}{n}}$ [12].

Using the Taylor series expansion for $\varphi \left(x;q\right)$, we obtain the following:

$$\varphi \left(x;q\right)=y\left(0\right)+\sum _{m=1}^{\infty }{u}_m\left(x\right)q^m$$

(5)

where

$$u_m\left(x\right)={\left.{\displaystyle \frac{1}{m!}}{\displaystyle \frac{{\mathit{\partial }}^m\varphi \left(x;q\right)}{\partial q^m}}\right|}_{q=0}.$$

(6)

Assuming the auxiliary linear parameter $h$, the auxiliary function $H\left(x\right)$ and the initial condition $y\left(0\right)$ are properly selected such that the series (5) converges as $q\to {\displaystyle \frac{1}{n}}$, the approximate solution is given as follows:

$$y\left(x\right)=\varphi \left(x;{\displaystyle \frac{1}{n}}\right)=y\left(0\right)+\sum _{m=1}^{\infty }{u}_m\left(x\right){\left({\displaystyle \frac{1}{n}}\right)}^m.$$

(7)

Define the vector as follows:

$${\overrightarrow{u}}_r\left(x\right)=\left\{u_0\left(x\right),u_1\left(x\right),\dots ,u_r\left(x\right)\right\}.$$

(8)

If Equation (4) is differentiated $m$ times with respect to $q$ and divided by $m!$ and $q$ is set to zero, the $m$ -th order deformation equation is obtained as follows:

$$L\left[u_m\left(x\right)-{\chi }_m{u}_{m-1}\left(x\right)\right]=hH\left(x\right)R_m\left({\overrightarrow{u}}_{m-1}\left(x\right)\right)$$

(9)

with the initial condition as follows:

$$u_m^{\left(k\right)}\left(x\right)=0,k=\mathrm{0,1},2,\dots ,m-1$$

(10)

where

$$R_m\left({\overrightarrow{u}}_{m-1}\left(x\right)\right)={\displaystyle \frac{1}{\left(m-1\right)!}}{\left.{\displaystyle \frac{{\partial }^{m-1}N\left[\varphi \left(x;q\right)\right]}{\partial q^{m-1}}}\right|}_{q=0}$$

(11)

and

$${\chi }_m=\left\{\begin{array}{c}0, if m\le 1\\ n, if m>1\end{array}\right..$$

(12)

Note that, for $n=1$ in Equation (7), the method reduces to the standard HAM. However, if $n\ne 1$, the existence of the factor ${\left({\displaystyle \frac{1}{n}}\right)}^m$ in Equation (7) enables a much faster convergence for q-HAM compared to the standard HAM.

Considering the Abel equation of the first kind (1), differentiation of Equation (4) with respect to $q=0$ gives the following first order deformation equation:

$${\displaystyle \frac{d{u}_1}{dx}}\left(x\right)=hH\left(x\right)\left[{\displaystyle \frac{d{u}_0}{dx}}\left(x\right)-f_3\left(x\right)u_0^3\left(x\right)-f_2\left(x\right)u_0^2\left(x\right)-f_1\left(x\right)u_0\left(x\right)-f_0\left(x\right)\right].$$

(13)

The general form of the $m$ -th order deformation equation can be given as follows:

$${\displaystyle \frac{d{u}_m}{dx}}\left(x\right)=n{\displaystyle \frac{d}{dx}}\left[u_{m-1}\left(x\right)\right]+hH\left(x\right)\left[{\displaystyle \frac{d{u}_{m-1}}{dx}}\left(x\right)-f_3\left(x\right)\sum _{i=0}^{m-1}{u}_{m-1-i}\sum _{j=0}^i{u}_j{u}_{i-j}-f_2\left(x\right)\sum _{k=0}^{m-1}{u}_k{u}_{m-1-k}-f_1\left(x\right)u_{m-1}\left(x\right)\right].$$

(14)

The solution of Equation (14) considering the initial condition (10) gives $u_m\left(x\right)$ and the numerical solution is obtained through Equation (7).

Improving the Convergence through Minimizing the Residual Error

Several studies by Liao and others [9,41] have presented a methodology for obtaining the optimal convergence control parameters by minimizing the square residual errors over the whole region [42]. The method, up to the stage where the solution $u_m\left(x\right)$ is obtained, is called q-HAM. The naming “optimal q-HAM” is based on the idea of selecting the optimal $h$ value for better convergence. A limited number of studies, such as [40], focus on the use of the analysis of Abel equations with HAM or other methods. The main goal of this study is to improve the convergence of solutions obtained by HAM and q-HAM through optimizing the auxiliary parameter. The methodology used to obtain this improvement works in a similar manner to the least-squares approach. If the square residual error is denoted as follows:

$${\mathsf{\Delta }}_m\left(h\right)={\int }_{\mathsf{\Omega }}{\left(N\left(u_m\left(x\right)\right)\right)}^{2}d\mathsf{\Omega },$$

(15)

then the optimal value of the auxiliary parameter $h$ will be the value minimizing the square residual error, which can be obtained by solving the equation below:

$${\displaystyle \frac{d}{dh}}{\mathsf{\Delta }}_m\left(h\right)=0.$$

(16)

This study presents results for a case study of the application of optimal q-HAM for the analysis of Abel equations and focuses on the determination of the optimal values of the auxiliary parameter for analyzing fractional Abel differential equations defined in the Caputo sense. Hence, the definitions of Riemann–Liouville and Caputo fractional derivative [43] operators needed for the fractional analysis are given in this section.

Definition 1. 

[44] Let $n\in {\mathbb{R}}_{+}$. The operator $J_a^n$, defined on $L_1[a,b]$ by

$$J_a^{n}f\left(x\right)={\displaystyle \frac{1}{\Gamma \left(n\right)}}{\int }_a^x{\left(x-t\right)}^{n-1}f\left(t\right)dt$$

for $a\le x\le b$, is called the Riemann–Liouville fractional integral operator of the order $n$. It is also known that when $n=0$, $J_a^0=I$ is obtained (where $I$ is the identity operator).

Definition 2. 

[44] Let $n\in {\mathbb{R}}_{+}$ and $m=⌈n⌉$. The operator $D_a^n$, defined by following:

$$D_a^{n}f\left(x\right)=D^m{J}_a^{m-n}f\left(x\right)$$

is called the Riemann–Liouville fractional derivative operator of the order $n$ (where $⌈n⌉$ is the ceiling function and denotes the smallest integer that is larger than $n$). Once again, for $n=0$, $D_a^0=I$ is obtained.

Definition 3. 

[44] Let $n=⌈\alpha ⌉$. Then, the Caputo fractional derivative operator ${}_{ }{}^{C}D_t^{\alpha }$ is given as follows:

$${}_{ }{}^{C}D_t^{\alpha }f(t)={\displaystyle \frac{1}{\Gamma (n-\alpha )}}{\int }_0^t{\displaystyle \frac{f^{\left(n\right)}\left(\tau \right)}{{\left(t-\tau \right)}^{\alpha +1-n}}}d\tau .$$

Properties of the Riemann–Liouville fractional integral and derivative operators and the Caputo fractional derivative operator can be found in the literature. We will be denoting the Caputo fractional derivative as $D_{\mathrm{x}}^{\alpha }$, where $\alpha$ is the order or derivation, in correspondence with the general notation in the literature.

3. Numerical Examples

In this section, various forms of the (ordinary) Abel differential equation of the first kind, including Bernoulli and Riccati equations with varying coefficients, are given as numerical examples in addition to fractional Abel differential equations with a fractional Caputo derivative.

Example 1 (Riccati Differential Equation).

In Equation (1), if we take $f_3\left(x\right)=0$, $f_2\left(x\right)=1-x$, $f_1\left(x\right)=2x-1$ and $f_0\left(x\right)=-x$, the following Riccati differential equation is obtained, where $0\le x\le 1$.

$$y^{\prime }\left(x\right)=\left(1-x\right)y^2\left(x\right)+\left(2x-1\right)y\left(x\right)-x$$

Putting the initial condition $\mathrm{y}\left(0\right)=2$, the analytical solution can be formulated as follows:

$$y\left(x\right)={\displaystyle \frac{3-e^x+x{e}^x}{3-2e^x+x{e}^x}}.$$

Using $\mathrm{H}\left(\mathrm{x}\right)=1$ in q-HAM, the first three terms are obtained as follows:

$$\begin{array}{c}{u}_1\left(x\right)=h\left(-2x+{\displaystyle \frac{x^2}{2}}\right),\hfill \\ u_2\left(x\right)=hn\left(-2x+{\displaystyle \frac{x^2}{2}}\right)+h\left(-2hx+{\displaystyle \frac{7h{x}^2}{2}}-{\displaystyle \frac{11h{x}^3}{6}}+{\displaystyle \frac{h{x}^4}{4}}\right),\hfill \\ u_3\left(x\right)=h\left(-2h^{2}x-2hnx+{\displaystyle \frac{13h^2{x}^2}{2}}+{\displaystyle \frac{7}{2}}hn{x}^2-8h^2{x}^3-{\displaystyle \frac{11}{6}}hn{x}^3+{\displaystyle \frac{39h^2{x}^4}{8}}+{\displaystyle \frac{1}{4}}hn{x}^4-{\displaystyle \frac{4h^2{x}^5}{3}}+{\displaystyle \frac{h^2{x}^6}{8}}\right)\hfill \\ +n\left(hn\left(-2x+{\displaystyle \frac{x^2}{2}}\right)+h\left(-2hx+{\displaystyle \frac{7h{x}^2}{2}}-{\displaystyle \frac{11h{x}^3}{6}}+{\displaystyle \frac{h{x}^4}{4}}\right)\right).\end{array}$$

Solution (7) from the q-HAM can be given as a finite series as follows:

$$y_k\left(x\right)=\sum _{m=1}^{\infty }{u}_m\left(x\right){\left({\displaystyle \frac{1}{n}}\right)}^m\cong \sum _{m=1}^k{u}_m\left(x\right){\left({\displaystyle \frac{1}{n}}\right)}^m.$$

(17)

This solution depends on the value of the auxiliary parameter $h$. Some of the relative errors have been given below in

Table 1. Relative errors obtained for the Riccati differential equation with q-HAM.

$\mathit{x}$	$\mathit{n}=1, \mathit{h}=-0.5$	$\mathit{n}=1, \mathit{h}=-1$	$\mathit{n}=1, \mathit{h}=-1.25$	$\mathit{n}=20, \mathit{h}=-1$
$0.1$	$0.00037533$	$2.74631072\times {10}^{-10}$	$2.54314884\times {10}^{-9}$	$0.06599027$
$0.2$	$0.00220035$	$3.68633053\times {10}^{-7}$	$1.16212931\times {10}^{-9}$	$0.14343628$
$0.3$	$0.00791511$	$0.00002018$	$1.08409111\times {10}^{-7}$	$0.23032080$
$0.4$	$0.02150092$	$0.00029051$	$0.00001266$	$0.32344195$
$0.5$	$0.04762474$	$0.00196197$	$0.00024490$	$0.41849102$
$0.6$	$0.08942998$	$0.00800457$	$0.00184491$	$0.51032883$
$0.7$	$0.14554139$	$0.02243673$	$0.00758219$	$0.59345618$
$0.8$	$0.20780964$	$0.04627964$	$0.01987730$	$0.66261677$
$0.9$	$0.26161939$	$0.07256799$	$0.03572155$	$0.71341999$
$1.0$	$0.28971143$	$0.08673558$	$0.04471525$	$0.74285215$

with the corresponding $n$ and $h$ values used for obtaining the approximate solutions (Table 1).

The values in Table 1 represent approximate solutions with growing relative errors. It is also seen that the amount of error in these solutions, obtained with q-HAM, change with respect to the changes in the value of the auxiliary parameter $h$. The use of optimal q-HAM offers a methodology for determining the approximate solution from q-HAM that has the least error, through the analysis of the optimal $h$.

The square residual error for the Riccati differential equation is written below to obtain the optimum value of the auxiliary parameter $h$.

$${\Delta }_m={\int }_{\mathsf{\Omega }}{\left[y_k^{\prime }\left(x\right)-\left(1-x\right)y_k^2\left(x\right)-\left(2x-1\right)y\left(x\right)+x\right]}^{2}dx.$$

(18)

The optimal value of $h$ is obtained by minimizing the square residual error given in Equation (18). The $h$-curves for various values of the parameter $n$ at $x=1$ for $k=10$ have been shown in Figure 1. The convergence regions have been determined from these graphs using the intervals where the lines are parallel to the $x$-axis. The optimal values of $h$ are selected in these regions through the roots of the nonlinear equation as follows:

$${\displaystyle \frac{d{\Delta }_m}{dh}}=0.$$

(19)

The convergence regions, the optimal $h$ values and the minimum values of Equation (18) for $k=10$ are given in

Table 2. Optimal $h$ and the convergence region of $h$ for the Riccati differential equation.

$\mathit{n}$	Convergence Region	$\mathbf{Optimal} \mathit{h}$	${\Delta }_{\mathit{m}}$
$1$	$-2.1793\le h\le -1.4460$	$-1.87312216$	$0.00946045$
$5$	$-11.0632\le h\le -5.8422$	$-9.36561082$	$0.01025391$
$10$	$-21.7929\le h\le -14.4595$	$-18.73122164$	$0.01025391$
$15$	$-32.6281\le h\le -22.3157$	$-28.09683247$	$0.00982666$
$20$	$-43.2345\le h\le -34.4072$	$-37.46244329$	$0.01025391$

for the Riccati differential equation.

The approximation solution and the analytical solution with the optimal $h$ values are shown in Figure 2. The figure contains the exact solution, the solution obtained with q-HAM and the solution obtained by using Optimal q-HAM with $n=1$. Note that the solution curve for $n=20$ is inseparable from the curve for $n=1$ in this example and hence, not plotted additionally.

Figure 2 shows that the analytical solution and the approximate solution (17) are similar for optimal $h$ values. The figure also contains the solution curve obtained with q-HAM for $h=-1$. Note that all of the solution curves obtained with optimal q-HAM are plotted similarly in all figures for a consistent presentation. In addition, these two solutions have been numerically compared in the following table for the various values of the independent variable (

Table 3. The exact and optimal q-HAM approximate solutions compared for $k=10$.

$\mathit{x}$	$\mathbf{Relative} \mathbf{Error} (\mathit{n}=1$)	$\mathbf{Relative} \mathbf{Error} (\mathit{n}=20$)
$0.1$	$0.00011073$	$0.00011073$
$0.2$	$0.00142025$	$0.00142025$
$0.3$	$0.00143402$	$0.00143402$
$0.4$	$0.00096602$	$0.00096602$
$0.5$	$0.00021452$	$0.00021452$
$0.6$	$0.00057228$	$0.00057228$
$0.7$	$0.00138279$	$0.00138279$
$0.8$	$0.00299686$	$0.00299686$
$0.9$	$0.00593229$	$0.00593229$
$1.0$	$0.00767123$	$0.00767123$

).

Table 3 shows that, considering the relative errors shown in the last column, which all satisfy ${\epsilon }_{Relative}<0.8\%$ in the analyzed interval, it can be concluded that the approximate solution is a suitable approximation of the analytical solution. It can also be seen that, for $n=1$ and $n=20$, the relative errors in the solutions obtained with q-HAM decrease as the value of the auxiliary parameter $h$ approaches the optimal value used for Oq-HAM.

Example 2 (Abel Differential Equation).

In Equation (1), if we take $f_3\left(x\right)=-x^2$, $f_2\left(x\right)=5x$, $f_1\left(x\right)=-2x$ and $f_0\left(x\right)=x^3$, the following Abel differential equation of the first kind is obtained, where $0\le x\le 0.7$.

$$y^{\prime }\left(x\right)=-x^2{y}^3\left(x\right)+5x{y}^2\left(x\right)-2xy\left(x\right)+x^3.$$

The initial condition will be used as $\mathrm{y}\left(0\right)=1$. The approximate solution from q-HAM will be compared to the numerical solution obtained with Mathematica using NDSolve.

Using $\mathrm{H}\left(\mathrm{x}\right)=1$ in the q-HAM, the first three terms are obtained as follows:

$$\begin{gathered} \begin{array}{c}{u}_1\left(x\right)=h\left(-{\displaystyle \frac{3x^2}{2}}+{\displaystyle \frac{x^3}{3}}-{\displaystyle \frac{x^4}{4}}\right),\hfill \\ u_2\left(x\right)=hn\left(-{\displaystyle \frac{3x^2}{2}}+{\displaystyle \frac{x^3}{3}}-{\displaystyle \frac{x^4}{4}}\right)+h\left(-{\displaystyle \frac{3h{x}^2}{2}}+{\displaystyle \frac{h{x}^3}{3}}+{\displaystyle \frac{11h{x}^4}{4}}-{\displaystyle \frac{43h{x}^5}{30}}+{\displaystyle \frac{h{x}^6}{2}}-{\displaystyle \frac{3h{x}^7}{28}}\right),\hfill \\ u_3\left(x\right)=h\left(-{\displaystyle \frac{3}{2}}{h}^2{x}^2-{\displaystyle \frac{3}{2}}hn{x}^2+{\displaystyle \frac{h^2{x}^3}{3}}+{\displaystyle \frac{1}{3}}hn{x}^3+{\displaystyle \frac{23h^2{x}^4}{4}}+{\displaystyle \frac{11}{4}}hn{x}^4-{\displaystyle \frac{43h^2{x}^5}{15}}-{\displaystyle \frac{43}{30}}hn{x}^5-{\displaystyle \frac{39h^2{x}^6}{8}}+{\displaystyle \frac{1}{2}}hn{x}^6 \\ +{\displaystyle \frac{1843h^2{x}^7}{420}}-{\displaystyle \frac{3}{28}}hn{x}^7-{\displaystyle \frac{2809h^2{x}^8}{1440}}+{\displaystyle \frac{485h^2{x}^9}{756}}-{\displaystyle \frac{127h^2{x}^{10}}{1120}}+{\displaystyle \frac{3h^2{x}^{11}}{176}}\right) \\ +n\left(hn\left(-{\displaystyle \frac{3x^2}{2}}+{\displaystyle \frac{x^3}{3}}-{\displaystyle \frac{x^4}{4}}\right)+h\left(-{\displaystyle \frac{3h{x}^2}{2}}+{\displaystyle \frac{h{x}^3}{3}}+{\displaystyle \frac{11h{x}^4}{4}}-{\displaystyle \frac{43h{x}^5}{30}}+{\displaystyle \frac{h{x}^6}{2}}-{\displaystyle \frac{3h{x}^7}{28}}\right)\right).\hfill \end{array} \end{gathered}$$

Similarly, the q-HAM solution for this problem is given in the form of Equation (17). To calculate the optimal $h$ values, the derivative of ${\Delta }_m$ in the equation is calculated with respect to $h$ and set to zero. The $h$-curves for different $n$ values and $k=6$ are obtained at $x=0.6$ as shown in Figure 3.

The convergence regions, the optimal $h$ values and the minimum values of Equation (1) for $k=6$ for the Abel differential equation are given in

Table 4. Optimal $h$ and the convergence region of $h$ for the Abel differential equation.

n	Convergence Region	$\mathbf{Optimal} \mathit{h}$	${\Delta }_{\mathit{m}}$
$1$	$-2.4872\le h\le -1.0394$	$-1.17183989$	$0.00427890$
$5$	$-12.4361\le h\le -5.1972$	$-5.83815009$	$0.00468601$
$10$	$-22.1298\le h\le -10.9513$	$-11.68006881$	$0.00465427$
$15$	$-24.8946\le h\le -15.5917$	$-17.60030257$	$0.00399791$
$20$	$-30.9369\le h\le -22.6137$	$-23.36013762$	$0.00465427$

. Also, the optimal q-HAM solution and the NDSolve solution with the optimal $h$ values are shown in Figure 4 for $n=20$. Once again, the solution curve for $n=1$ has been omitted to prevent repeated use of an inseparable plot. Figure 4 shows that the approximate solution is in accordance with the solution obtained with NDSolve for $k=6$.

Example 3 (Fractional Bernoulli Differential Equation).

A fractional Bernoulli differential equation can be obtained from the ordinary Abel differential equation using fractional derivatives and suitable coefficient functions. The following fractional Bernoulli initial value problem will be considered for $x\in \left(\mathrm{0,1}\right)$.

$$D_x^{\alpha }y=x^{4}y-x^2{y}^2, y\left(0\right)=1$$

(20)

The $\mathrm{m}$-th order deformation can be obtained for fractional differential equations in similarity to the q-HAM approximation approach for ordinary differential equations. In this case, the ordinary derivation in Equation (14) needs to be replaced with fractional derivatives. The Caputo derivative (Definition 3) is used in this regard to obtain a fractional differential equation.

Using $\mathrm{H}\left(\mathrm{x}\right)=1,\mathsf{\alpha }=0.5$ in the q-HAM, the first three terms are obtained as

$$\begin{gathered} \begin{array}{c}{u}_1\left(x\right)=0.5641895835477563h\left(1.0666666666666667x^{2.5}-0.8126984126984127x^{4.5}\right),\hfill \\ u_2\left(x\right)=0.5641895835477563h^2\left(0.9305382717253959x^5-0.9995109922628649x^7 \\ +0.2671648878690047x^9\right) \\ +\left(h+n\right)\left(0.5641895835477563h\left(1.0666666666666667x^{2.5}-0.8126984126984127x^{4.5}\right)\right),\hfill \\ u_3\left(x\right)=0.5641895835477563h^2\left(n\right(0.930538271725396x^5-0.9995109922628648x^7 \\ +0.2671648878690047x^9)+h(0.930538271725396x^5-0.9995109922628648x^7 \\ +0.8988704875870189x^{7.5}+0.2671648878690047x^9-1.2512690912236168x^{9.5} \\ +0.5561133343471595x^{11.5}-0.07204306295221984x^{13.5}\left)\right)+(h \\ +n)\left(0.5641895835477563h^2\right(0.9305382717253959x^5-0.9995109922628649x^7 \\ +0.2671648878690047x^9)+(h+n\left)\right(0.5641895835477563h(1.0666666666666667x^{2.5} \\ -0.8126984126984127x^{4.5})\left)\right).\hfill \end{array} \end{gathered}$$

The approximate results from the q-HAM will be compared to the numerical results obtained with the Predictor–Corrector method for fractional differential equations given in [45]. The numerical method is an improved version of the Adams–Bashforth–Moulton algorithm and will be referred to as the Predictor–Corrector (PC) method throughout the study.

To calculate the optimal value of the auxiliary parameter $h$, the $h$-curve graphs are drawn for $\alpha =0.5$ and $\alpha =1$ with $k=4$, which is the number of terms taken in the q-HAM solution at $x=1$ in Figure 5a,b.

For the fractional Bernoulli differential equation, the exact square residual error is obtained as follows:

$${\Delta }_m={\int }_{\mathsf{\Omega }}{\left[D_x^{\alpha }{y}_k\left(x\right)-x^4{y}_k\left(x\right)+x^2{y}_k^2\left(x\right)\right]}^{2}dx.$$

If the derivative of ${\Delta }_m$ with respect to $h$ is calculated and solved for zero, the tables of optimal $h$ values obtained for $\alpha =0.5$ and $\alpha =1$ are given in

Table 5. The optimal $h$ and the convergence region of $h$ for $\alpha =0.5$.

$\mathit{n}$	Convergence Region	$\mathbf{Optimal} \mathit{h}$	${\Delta }_{\mathit{m}}$
$1$	$-0.7635\le h\le -0.5590$	$-0.71310376$	$3.37203607\times {10}^{-8}$
$5$	$-3.8175\le h\le -2.7950$	$-3.56551974$	$3.37178963\times {10}^{-8}$
$10$	$-8.0043\le h\le -5.2082$	$-7.13103948$	$3.37178963\times {10}^{-8}$
$15$	$-12.0063\le h\le -7.8123$	$-10.69655936$	$3.37170178\times {10}^{-8}$
$20$	$-16.5847\le h\le -9.8154$	$-14.26207897$	$3.37178963\times {10}^{-8}$

and

Table 6. The optimal $h$ and the convergence region of $h$ for $\alpha =1$.

n	Convergence Region	$\mathbf{Optimal} \mathit{h}$	${\Delta }_{\mathit{m}}$
$1$	$-0.9087\le h\le -0.7178$	$-0.86488950$	$3.40577788\times {10}^{-10}$
$5$	$-4.5790\le h\le -3.5435$	$-4.30383827$	$3.40577233\times {10}^{-10}$
$10$	$-8.6973\le h\le -7.6468$	$-8.60767655$	$3.40577233\times {10}^{-10}$
$15$	$-12.9723\le h\le -11.5535$	$-12.91151482$	$3.40577955\times {10}^{-10}$
$20$	$-18.0723\le h\le -15.8497$	$-17.21535309$	$3.40577233\times {10}^{-10}$

, respectively.

From Table 5 and Table 6, it is seen that, for the optimal values of $h$, the square residual errors are minimal. In addition, the solutions for the optimal $h$ values obtained with optimal q-HAM have been found to be similar to the solutions from the PC method (for $\alpha =0.5$), as shown in Figure 6, and the exact solutions (for $\alpha =1$), as shown in Figure 7. Hence, it can be concluded that the optimal $h$-values given in Table 5 and Table 6 are suitable for the fractional Bernoulli differential equation. Note that $n=1$ and $n=20$ produce inseparable curves for Figure 6 and Figure 7, hence, only the figures for $n=1$ have been given.

In order to give a comparison of the solutions obtained with the optimal q-HAM and the solutions from another approximation method for fractional differential equations, we use the fractional Differential Transform Method [46]. The fractional counterpart of the Differential Transform Method (DTM) has been another widely used tool for analyzing equations with fractional order derivatives and has been applied to study projectile motion [47] and COVID-19 transmission [48]. The reader is advised to refer to the study of Arikoglu and Ozkol for the theorems and properties of fractional DTM [46]. The solution of the fractional Bernoulli differential Equation (20) using fractional DTM has been shown below (Figure 8).

The values for the numerical solution (PC method) and the approximate solutions have been compared for $x\in \left[\mathrm{0,1}\right]$ in the table below (

Table 7. Numerical and approximate solutions of (20) compared for $x\in \left[\mathrm{0,1}\right]$ and $n=1$.

$\mathit{x}$	PC Method	DTM	Oq-HAM	q-HAM $\left(\mathit{h}=-0.5\right)$	q-HAM $\left(\mathit{h}=-1\right)$	q-HAM $\left(\mathit{h}=-1.25\right)$
$0.1$	$0.99809477$	$0.99811661$	$0.99812903$	$0.99823303$	$0.99811661$	$0.99812438$
$0.2$	$0.98968409$	$0.98972089$	$0.98978053$	$0.99032479$	$0.98972089$	$0.98977507$
$0.3$	$0.97341252$	$0.97347280$	$0.97358992$	$0.97487314$	$0.97347292$	$0.97369769$
$0.4$	$0.95054829$	$0.95063802$	$0.95077595$	$0.95282813$	$0.95064119$	$0.95145505$
$0.5$	$0.92440765$	$0.92453475$	$0.92462529$	$0.92713087$	$0.92456096$	$0.92725681$
$0.6$	$0.89960938$	$0.90008460$	$0.89976706$	$0.90220525$	$0.90002374$	$0.90798397$
$0.7$	$0.88137849$	$0.88596782$	$0.88143340$	$0.88332823$	$0.88287250$	$0.90323390$
$0.8$	$0.87523525$	$0.91332540$	$0.87512733$	$0.87617499$	$0.87987603$	$0.92411447$
$0.9$	$0.88720130$	$1.13112410$	$0.88689226$	$0.88683630$	$0.89823968$	$0.97849497$
$1.0$	$0.92444239$	$2.18855686$	$0.92399200$	$0.92244773$	$0.94414273$	$1.06294033$

).

Note that the approximate solutions for fractional DTM have been obtained with $30$ terms in the approximation. Although a detailed approximation has been obtained for (20) with fractional DTM, the results show that the optimal q-HAM provides a better approximation with optimal $h$ values. Hence, it can be concluded that the convergence is improved in comparison to the classical q-HAM and the DTM method as well by optimizing the value of the auxiliary parameter. The relative errors in

Table 8. Relative errors in approximate solutions of (20) compared.

$\mathit{x}$	DTM	q-HAM $\left(\mathit{h}=-0.5\right)$	q-HAM $\left(\mathit{h}=-1\right)$	q-HAM $\left(\mathit{h}=-1.25\right)$	Oq-HAM
$0.1$	$0.00002188$	$0.00013853$	$0.00002188$	$0.00002967$	$0.03431887\times {10}^{-3}$
$0.2$	$0.00003718$	$0.00064738$	$0.00003718$	$0.00009193$	$0.09744808\times {10}^{-3}$
$0.3$	$0.00006192$	$0.00150051$	$0.00006204$	$0.00029295$	$0.18224420\times {10}^{-3}$
$0.4$	$0.00009440$	$0.00239845$	$0.00009774$	$0.00095394$	$0.23951227\times {10}^{-3}$
$0.5$	$0.00013749$	$0.00294591$	$0.00016584$	$0.00308215$	$0.23543438\times {10}^{-3}$
$0.6$	$0.00052825$	$0.00288554$	$0.00046059$	$0.00930913$	$0.17527668\times {10}^{-3}$
$0.7$	$0.00520699$	$0.00221215$	$0.00169508$	$0.02479686$	$0.06229900\times {10}^{-3}$
$0.8$	$0.04351990$	$0.00107370$	$0.00530232$	$0.05584695$	$0.12330404\times {10}^{-3}$
$0.9$	$0.27493513$	$0.00041140$	$0.01244180$	$0.10290074$	$0.34832931\times {10}^{-3}$
$1.0$	$1.36743456$	$0.00215768$	$0.02131052$	$0.14981782$	$0.48719902\times {10}^{-3}$

show how the approximation obtained with optimal q-HAM has incomparably smaller relative errors. For instance, at $x=0.9$, the relative error for DTM ($0.274935128211391$) is more than 700 times larger than the relative error for optimal q-HAM ($0.348329311959152\times {10}^{-3}$).

The solutions obtained with q-HAM have varying relative errors, depending on the choice of the value of the auxiliary parameter $h$. Figure 6 presents the solution curve obtained by using the optimal $h$ value in comparison to the solution curve obtained with q-HAM using $h=-1.$ As the value of the parameter $h$ approaches its optimal value, the relative errors are expected to decrease, which also means the curves from Oq-HAM and q-HAM become closer. The results in Table 7 and Table 8 and the curves in Figure 6 are a clear indicator of the importance of using the optimal $h$ value. This value, obtained by analyzing the square errors to obtain the minimal error, clearly results in a better approximation to the exact solution.

Example 4 (Fractional Abel Differential Equation).

Consider the following fractional Abel differential equation of the first kind.

$$D_x^{\alpha }y=x^{4}y-x^2{y}^2-x{y}^3+x, y\left(0\right)=1.$$

Using $\mathrm{H}\left(\mathrm{x}\right)=1$ and $\mathsf{\alpha }=0.5$ in the q-HAM, the first three terms are obtained as

$$\begin{gathered} \begin{array}{c}{u}_1\left(x\right)=0.5641895835477563h\left(1.0666666666666667x^{2.5}-0.8126984126984127x^{4.5}\right),\hfill \\ u_2\left(x\right)=0.5641895835477563h^2(1.550897119542327x^4+0.930538271725396x^5 \\ -0.974849617998034x^6-0.999510992262865x^7+0.2671648878690047x^9) \\ +\left(h \\ +n\right)\left(0.5641895835477563h\left(1.0666666666666667x^{2.5}-0.8126984126984127x^{4.5}\right)\right),\hfill \\ u_3\left(x\right)={\displaystyle \frac{1}{x^{0.5}}}0.5641895835477563h^2\left(n\left(1.550897119542327x^{4.5}+0.930538271725396x^{5.5} \\ -0.974849617998034x^{6.5}-0.999510992262865x^{7.5}+0.2671648878690048x^{9.5}\right) \\ +h\left(1.550897119542327x^{4.5}+0.930538271725396x^{5.5}+1.9393939393939406x^6 \\ -0.974849617998034x^{6.5}+3.0085738865913085x^7-0.999510992262865x^{7.5} \\ -0.15138592266939213x^8-3.1885000918855044x^9+0.2671648878690048x^{9.5} \\ -1.2512690912236175x^{10}+0.882618761181763x^{11}+0.5561133343471595x^{12} \\ -0.07204306295221984x^{14}\right)\right)+(h+n)\left(0.5641895835477563h^2\right(1.550897119542327x^4 \\ +0.930538271725396x^5-0.974849617998034x^6-0.999510992262865x^7 \\ +0.2671648878690047x^9)+(h+n\left)\right(0.5641895835477563h(1.0666666666666667x^{2.5} \\ -0.8126984126984127x^{4.5})\left)\right).\hfill \end{array} \end{gathered}$$

The $h$-curve graphs for $\alpha =0.5$ and $\alpha =1$ at $x=1$ are given in Figure 9a,b. Here, four terms $(k=4)$ are used for the q-HAM solution. The optimal $h$ values can be determined in the convergence region where $u_m^{″}\left(x\right)=0$.

The convergence regions, the optimal $h$-values and square residual errors for $\alpha =0.5$ and $\alpha =1$ are given, in

Table 9. The optimal $h$ and the convergence region of $h$ for $\alpha =0.5$.

n	Convergence Region	$\mathbf{Optimal} \mathit{h}$	${\Delta }_{\mathit{m}}$
$1$	$-0.4733<h<-0.3412$	$-0.44457431$	$0.00001744$
$5$	$-1.7703<h<-3.0819$	$-2.22287156$	$0.00001744$
$10$	$-4.7330<h<-3.5605$	$-4.44574312$	$0.00001744$
$15$	$-7.8677<h<-5.2408$	$-6.66861468$	$0.00001744$
$20$	$-7.0409<h<-11.4363$	$-8.89148624$	$0.00001744$

and

Table 10. The optimal $h$ and the convergence region of $h$ for $\alpha =1$.

$\mathit{n}$	Convergence Region	$\mathbf{Optimal} \mathit{h}$	${\Delta }_{\mathit{m}}$
$1$	$-0.8361<h<-0.5755$	$-0.66406054$	$0.00000287$
$5$	$-3.9465<h<-2.9461$	$-3.32030272$	$0.00000287$
$10$	$-8.3606<h<-6.0067$	$-6.64060545$	$0.00000287$
$15$	$-11.8395<h<-9.0364$	$-9.96090817$	$0.00000287$
$20$	$-14.7504<h<-12.0772$	$-13.28121090$	$0.00000287$

, respectively.

The numerical solution graphs of the optimal q-HAM and the PC method for the $h$ values obtained are shown in Figure 10 and Figure 11. Only one plot has been given for $n=1$ and $n=20$ in both figures, since these values of the parameter produce inseparable solution curves.

It is seen from Figure 10 that the Oq-HAM solution and the PC method exhibit similar behaviors. To analyze the graphs numerically, the relative errors at several points can be analyzed. The relative errors are given in the following table for $\alpha =0.5$ and $\alpha =1$ (

Table 11. Relative errors from Oq-HAM compared for various values of $\alpha$ and $n$.

$\mathit{x}$	$\mathit{n}=1,\mathit{\alpha }=0.5$	$\mathit{n}=20,\mathit{\alpha }=0.5$	$\mathit{n}=1,\mathit{\alpha }=1$
$0.1$	$0.00018036$	$0.00018036$	$0.00000398$
$0.2$	$0.00061307$	$0.00061307$	$0.00002582$
$0.3$	$0.00084274$	$0.00084274$	$0.00005754$
$0.4$	$0.00058553$	$0.00058553$	$0.00006133$
$0.5$	$0.00003805$	$0.00003805$	$0.00000251$
$0.6$	$0.00040642$	$0.00040642$	$0.00010365$
$0.7$	$0.00055971$	$0.00055971$	$0.00018170$
$0.8$	$0.00046573$	$0.00046573$	$0.00014956$
$0.9$	$0.00008034$	$0.00008034$	$0.00003481$
$1.0$	$0.00077584$	$0.00077584$	$0.00044462$

).

Once again, the tables and figures show that for the optimal $h$ values, minimal square residual errors are obtained and the solutions from the PC and Oq-HAM methods are almost identical within the interval of interest.

The approximate solution of the fractional Abel differential equation obtained with fractional DTM using $30$ terms in the approximation is shown below (Figure 12).

Several values of $x\in \left[\mathrm{0,1}\right]$ have been compared below for the numerical solutions and the approximate solutions of the fractional Abel differential equation (

Table 12. Solutions of the fractional Abel equation compared.

$\mathit{x}$	PC Method	DTM	Optimal q-HAM $(\mathit{n}=1)$
$0.1$	$0.99816551$	$0.99819974$	$0.99834554$
$0.2$	$0.99081774$	$0.99090882$	$0.99142519$
$0.3$	$0.97841889$	$0.97858007$	$0.97924344$
$0.4$	$0.96353307$	$0.96378185$	$0.96409724$
$0.5$	$0.94931777$	$0.94939298$	$0.94935389$
$0.6$	$0.93876790$	$0.93435200$	$0.93838637$
$0.7$	$0.93448129$	$0.89709713$	$0.93395826$
$0.8$	$0.93869199$	$0.74365559$	$0.93825481$
$0.9$	$0.95338380$	$0.17671416$	$0.95330720$
$1.0$	$0.98037835$	$-1.58723431$	$0.98113897$

).

The relative errors (relative to the numerical solution) for the approximate solutions obtained with the fractional DTM and optimal q-HAM methods have been compared below (

Table 13. Relative errors for the fractional Abel equation compared.

$\mathit{x}$	DTM	Oq-HAM
$0.1$	$0.00003430$	$0.00018036$
$0.2$	$0.00009192$	$0.00061307$
$0.3$	$0.00016474$	$0.00084274$
$0.4$	$0.00025820$	$0.00058553$
$0.5$	$0.00007923$	$0.00003805$
$0.6$	$0.00470394$	$0.00040642$
$0.7$	$0.04000525$	$0.00055971$
$0.8$	$0.20777465$	$0.00046573$
$0.9$	$0.81464531$	$0.00008034$
$1.0$	$2.61900180$	$0.00077584$

).

Note that the improved convergence obtained through the use of the optimal value for the auxiliary parameter in the optimal q-HAM can be observed in the figures. The relative error for the optimal q-HAM is more than $1000$ times smaller compared to the relative error for the fractional DTM method at $x=0.9$. The results in Table 12 and Table 13 and Figure 10, Figure 11 and Figure 12 also underline that optimal q-HAM provides a better approximation for this equation.

4. Discussion

There are only a limited number of studies on the approximate solutions of fractional Abel differential equations in the literature. One study that might be considered related, by Jafari et al. [40], presents the application of HAM for the analysis of Abel equations of fractional order. This study improves the mentioned study and the following related studies in a way that adds to the current literature with a detailed presentation of Oq-HAM for Abel-type differential equations. The convergence analysis and the resulting optimality investigation for the auxiliary parameter $h$, comparison of the results with q-HAM, DTM, PC methods and exact solutions, etc. provide an in-depth analysis that extends the existing methodology. Some of the numerical results found in the literature are given below to compare the findings of this article with the existing literature numerically as well. The study by Al-Smadi et al. presents numerical values for the residual errors in the numerical solution of Abel differential equations analyzed fractionally with a Caputo–Fabrizio derivative using the reproducing kernel method [37]. The results of this study show that, for $t\in \left[\mathrm{0,1}\right]$, the residual errors have a minimal value of $2.675\times {10}^{-3}$ for $\alpha =0.9$ and $2.679\times {10}^{-3}$ for $\alpha =0.85$. Although the error decreases for non-integer orders of derivation that are closer to $1$, considering that the relative errors in Table 8 are given for $\alpha =0.5$, it can be said that the results in our study are better in some regions. For instance, the Oq-HAM relative error at $x=1.0$ (the highest amount in the inspected interval) is $0.487\times {10}^{-3}$. This relative error is also better than the error obtained with fractional DTM as well. The similarity between the approximate solutions of the optimal q-HAM and the exact/numerical solutions for the examples is also in accordance with the similarity of the results from the reproducing kernel method and the exact solutions given in Al-Smadi et al. [37]. Another study by Rigi and Tajadodi presents numerical results obtained for fractional Abel differential equations in the Caputo sense using Genocchi polynomials [38]. Results are given for $x\in \left[\mathrm{0,1}\right]$ for two examples, and it can be said that the similarity of the approximate solutions to the exact/numerical solutions in this study matches the similarity between the approximate and exact solutions presented in the study by Rigi and Tajadodi [38]. For instance, the relative error in Example 4, for $n=1,\alpha =0.5$ is around $0.776\times {10}^{-3}$ at $x=1$ in our study, whereas the absolute error for $\alpha =1$ with $N=4$ in the Genocchi polynomials method is obtained as $8.56552\times {10}^{-3}$ at $x=1$ for the first example in the referred study. Results from existing numerical approximation methods are also found to be similar to the findings presented in our study. However, considering that Oq-HAM is an approximate analytical method, it can be said that the method provides better results in comparison to the other approximate analytical solution methods available for fractional Abel equations in the literature.

5. Conclusions

In this study, the optimal q-Homotopy Analysis Method (Oq-HAM) has been used to analyze Abel differential equations with various coefficient functions. Bernoulli and Riccati equations, which are special forms of Abel differential equations, are used, along with the general case in ordinary and fractional frameworks, to demonstrate the application of q-HAM. The Caputo fractional derivative is used to obtain fractional Abel differential equations from ordinary equations, and applications of optimal q-HAM are given for ordinary Abel equations and fractional Abel equations. This flow has been specifically designed to start from the application of q-HAM for ordinary Abel equations and finally reach the application of the optimal method for fractional equations. Four numerical examples have been given to present the method for obtaining the approximate solutions of the equations using q-HAM. The solutions are then analyzed to determine the optimal value of the convergence control parameter by minimizing the square residual errors. A theoretical analysis of the improvement in the convergence obtained using the auxiliary parameter in optimal q-HAM has been given by Liao [9] for the case with ordinary derivatives. Hence, this study is structured to present the improvement gained for the analysis of fractional Abel equations by focusing on comparisons of the application of optimal q-HAM with various other schemes. The optimal values of the auxiliary parameter $h$ have been given to show the convergence regions, and the approximate solutions corresponding to these optimal values have been compared with the exact or numerical solutions of the problems. The results show that the use of Oq-HAM enables finding the optimal values of $h$ and, thus, finding better approximations than those obtained by using q-HAM. The exact solutions or numerical solutions obtained with the NDSolve function are shown to be almost identical to the approximate solutions from Oq-HAM when the optimal values of the auxiliary parameter are used. The specific aim of this study is to analyze the approximations with q-HAM for Abel equations and to improve the convergence through optimizing auxiliary parameters. This is verified through the presentation of the results showing that Oq-HAM enhances the convergence of the solutions that would be obtained by using HAM, q-HAM or similar other methods.

The first example, the ordinary Riccati equation, has been analyzed with q-HAM and Oq-HAM and the solutions have been compared with the exact solution. Relative error analysis shows that, for $n=1$, the optimal value of the auxiliary parameter $h$ is obtained as $h=-1.87312216$. This is verified by the analysis of approximate solutions obtained with q-HAM, where the solutions were analyzed for $h=-0.5,-1$ and $-1.25$ when $n=1$. Results show that the relative errors (relative to the exact solution) are growing along with the increase in $x$ and for $x=1$, the relative errors are obtained as $0.28971143$ $\left(h=-0.5\right)$, $0.08673558$ $\left(h=-1\right)$, $0.04471525$ $\left(h=-1.25\right)$ and $0.00767123$ (optimal $h$). This shows that the relative error for optimal $h$ is at least 170 times smaller for Oq-HAM. The improvement in the approximation using the optimal $h$ value can also be seen for the case where $n=20$. Figure 2 shows the similarity of the solution curves for the exact solution and the Oq-HAM solution. The improvement from q-HAM to Oq-HAM is further analyzed for the fractional Bernoulli equation in the third example. The relative errors in Table 8, shown for various selections of the parameters, are better than the relative errors given for the solutions from q-HAM in all cases analyzed for this example.

Approximate solutions of fractional Abel equations have also been obtained with fractional DTM to highlight the improved convergence obtained with optimal q-HAM. A similar growth in relative errors, along with the increase in $x$, is seen for both methods. However, as Figure 12 suggests, the relative error of DTM at $x=1$ $\left(2.61900180\right)$ is almost 3375 times more than the relative error obtained with Oq-HAM $\left(0.00077584\right)$. Further discussion of Oq-HAM is given in comparison to the exact solutions and numerical solutions obtained with the predictor–corrector method. Figures and tables provide results that verify that optimal q-HAM is a reliable tool for the investigation of fractional Abel differential equations.

This manuscript presents a case study for the application of the well-established optimal q-HAM to achieve improved approximations for fractional Abel differential equations. A detailed comparison of the approximate results obtained by using the optimal values of the auxiliary parameter determined with Oq-HAM justifies the advantages of the method. This methodology can be generalized to other fractional equations containing different fractional derivatives for a wider range of applications. Many differential equations and mathematical models can be evaluated through the use of optimal q-HAM to obtain better approximations of the solutions. Possible future studies can be made for the comparison of the solutions obtained with the q-HAM and optimal q-HAM methods with other approximation methods such as the VIM method for fractional Abel differential equations. Comparison of the approximate solutions with numerical schemes, such as the fractional Euler method, to investigate the convergence of solutions could also form the basis for new studies. The use of random components and other definitions of fractional derivation are among the possible ideas for future investigations.