A New Hybrid Optimal Auxiliary Function Method for Approximate Solutions of Non-Linear Fractional Partial Differential Equations

Abstract

This study uses the optimal auxiliary function method to approximate solutions for fractional-order non-linear partial differential equations, utilizing Riemann–Liouville’s fractional integral and the Caputo derivative. This approach eliminates the need for assumptions about parameter magnitudes, offering a significant advantage. We validate our approach using the time-fractional Cahn–Hilliard, fractional Burgers–Poisson, and Benjamin–Bona–Mahony–Burger equations. Comparative testing shows that our method outperforms new iterative, homotopy perturbation, homotopy analysis, and residual power series methods. These examples highlight our method’s effectiveness in obtaining precise solutions for non-linear fractional differential equations, showcasing its superiority in accuracy and consistency. We underscore its potential for revealing elusive exact solutions by demonstrating success across various examples. Our methodology advances fractional differential equation research and equips practitioners with a tool for solving non-linear equations. A key feature is its ability to avoid parameter assumptions, enhancing its applicability to a broader range of problems and expanding the scope of problems addressable using fractional calculus techniques.

1. Introduction

Many challenges that emerge in natural phenomena, including fluid mechanics, biology, and thermodynamics, can be effectively represented through mathematical models. These models can be translated into a mathematical framework using differential equations. Classifying these differential equations as either linear or non-linear is contingent upon the specific characteristics of the issues that manifest across various scientific domains. In contemporary times, a significant research thrust has been directed towards fractional-order differential equations. Fractional calculus stands as a refined adaptation of classical calculus. Addressing the challenge of solving non-linear fractional-order partial differential equations (FPDE) entails the utilization of a diverse array of both numerical and analytical methods.

Some non-linear fractional partial differential equations exist among these PDEs, such as time-fractional Cahn–Hilliard (TFCH) equations, fractional Burgers–Poisson equations, and BBM Burger equations. The Cahn–Hillard equation was named after Cahn and Hilliard in 1958 [1]. This equation is critical in order to comprehend a variety of interesting physical phenomena, such as spinodal decomposition, phase ordering dynamics, and phase separation processes. It also explains a critical qualitative distinguishing feature of two-phase systems (see [1,2,3,4] for a detailed discussion); furthermore, long memory processes and fractional integration in econometrics [5] and kinetics of phase decomposition processes [6] have also been discussed. Researchers have investigated mathematical and numerical solutions of TFCH equations [7,8,9,10,11] since the appearances of their real-life applications in the abovementioned fields. The fractional Burgers–Poisson (FBP) equation was proposed for the first time in 2004 to describe the propagation of long waves in dispersive media in a single direction [12]. The FBP equation models a unidirectional water wave with weaker dispersive effects than the KdV equation in shallow water. The BBM–Burger equation defines the mathematical model of the propagation of small-amplitude long waves in non-linear dispersive media.

The BBM equation is a refinement of the KdV equation, as is well known. The wave-breaking models are affected by the BBM–Burger equation and the KdV equation [13]. Water waves inspired the KdV equation, which was later used as a basis for long waves in a variety of other physical systems. However, the KdV equation was not valid in certain long-wave physical systems. As a result, the BBM–Burger model was proposed, representing unidirectional long-wave propagation in a non-linear dispersive system [13,14,15]. Solving fractional differential equations and fractional partial differential equations has been a significant focus for researchers. Since most fractional differential equations do not have exact analytic solutions, approximate and numerical methods are commonly used [16,17,18], Baleanu et al. discussed the planar system-masses in an equilateral triangle [19], Ghanim et al. discussed certain implementations in fractional calculus operators, some new extensions on fractional differential and integral properties, and some analytical merits of the Kummer-type function [20,21,22], Almalahi et al. discussed qualitative analysis of the Langevin integro-fractional differential equation [23], Jajarmi et al. and Sajjadi et al. discussed a new iterative method for the numerical solution of high-order non-linear fractional boundary value problems and fractional optimal control problems with a general derivative operator [24,25,26], and Mohammadi et al. discussed a hybrid functions numerical scheme for fractional optimal control problems [27]. On the well-posedness of the sub-diffusion equation with the conformable derivative model, mathematical modeling for the adsorption process of the dye removal Laplace–Carson integral transformation for exact solutions [28,29,30,31] has been discussed in the literature in order to find the numerical and approximate solutions. Many analytical methods have been attempted in order to solve non-linear problems, including the new iterative method (NIM), homotopy perturbation method (HPM), homotopy analysis method (q-HAM), residual power series method (RPSM), and FHATM.

Similarly, using Riemann–Liouville’s (R-L) fractional integral and the Caputo derivative, we present another approach known as the optimal auxiliary function method (OAFM) for higher dimensional equations, which was first introduced by Vasile Marinca and colleagues for thin film flow problem [32]. Laiq Zada et al. later expanded the approach to partial differential and generalized seventh-order KDV equations [33]. Several other approaches and models have been developed recently to deal with different types of fractional-order equations and PDEs in general; see [34,35,36,37]. Auxiliary convergence control parameters and auxiliary functions are included in this approach to control and accelerate the method’s convergence. The OAFM operates without the need for presuming any parameter to be small or large. This proposed technique possesses the benefit of adeptly handling both linear and non-linear challenges while preserving a broad scope of applicability and effectiveness. The rest of this work is structured as follows. We begin with an initial section that provides essential definitions aimed at aiding the readers’ comprehension. Following this, the subsequent section outlines the fundamental concepts underlying the OAFM. Moving forward, in the third section, we apply the OAFM to tackle fractional-order problems, leading to the derivation of valuable solutions that are accurately presented through tables and graphs. Exploring deeper, the fourth and fifth sections provide a comprehensive analysis and discussion of the presented tables and graphs, respectively. Finally, we conclude our study in the sixth section, summarizing the key findings and suggestions drawn from our work.

2. Preliminaries

We need some basic definitions to investigate our problems with the help of OAFM.

Definition 1 

([38]). Riemann–Liouville fractional integral left-sided:

$${}^{RL}\mathsf{I}{}_{a^{+}}^{\alpha }\left[f\left(r\right)\right]=\frac{1}{\mathsf{\Gamma }\left(\alpha \right)}\hspace{0.17em}{\displaystyle \underset{a}{\overset{\eta }{\int }}\left(\eta -r\right)}{\hspace{0.17em}}^{\alpha -1}f\left(r\right)dr\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\eta \ge a,$$

(1)

Definition 2 

([38]). Riemann–Liouville right-sided integral:

$${}^{RL}\mathsf{I}{}_{b^{-}}^{\alpha }\left[f\left(\eta \right)\right]=\frac{1}{\mathsf{\Gamma }\left(\alpha \right)}\hspace{0.17em}{\displaystyle \underset{\eta }{\overset{b}{\int }}\left(r-\eta \right)}{\hspace{0.17em}}^{\alpha -1}f\left(r\right)dr\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\eta \le b,$$

(2)

Definition 3 

([38]). Modified Riemann–Liouville fractional derivative:

$$D\left[f\left(r\right)\right]=\frac{1}{\mathsf{\Gamma }\left(1-\alpha \right)}\hspace{0.17em}\frac{d}{d\eta }{\displaystyle \underset{0}{\overset{\eta }{\int }}\left(\eta -r\right)}{\hspace{0.17em}}^{-\alpha }\hspace{0.17em}[f\left(r\right)-\hspace{0.17em}f\left(0\right)]\hspace{0.17em}dr.$$

(3)

Definition 4 

([38]). Liouville derivative:

$$D^{\alpha }\left[f(\eta )\right]=\frac{1}{\mathsf{\Gamma }\left(1-\alpha \right)}\frac{d}{d\eta }{\displaystyle \underset{-\infty }{\overset{\eta }{\int }}{\left(\eta -r\right)}^{-\alpha }f\left(r\right)dr,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-\infty \prec }\hspace{0.17em}\eta \hspace{0.17em}\hspace{0.17em}\prec \hspace{0.17em}\infty$$

(4)

Definition 5 

([38]). Liouville left-sided derivative:

$$D_{0^{+}}^{\alpha }\left[f(\eta )\right]=\frac{1}{\mathsf{\Gamma }\left(\eta -\alpha \right)}\frac{d^n}{d{\eta }^n}{\displaystyle \underset{0}{\overset{\eta }{\int }}{\left(\eta -r\right)}^{-\alpha +n+1}f\left(r\right)dr,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\eta \succ 0.$$

(5)

Definition 6 

([38]). Liouville right-sided derivative:

$$D_{-}^{\alpha }\left[f(\eta )\right]=\frac{{\left(-1\right)}^n}{\mathsf{\Gamma }\left(\eta -\alpha \right)}\frac{d^n}{d{\eta }^n}{\displaystyle \underset{\eta }{\overset{\infty }{\int }}{\left(\eta -r\right)}^{-\alpha +n-1}f\left(r\right)dr,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\eta \prec \infty .\hspace{0.17em}\hspace{0.17em}$$

(6)

Definition 7 

([38]). Right Riemann–Liouville integral of variable fractional order:

$${}_{\eta }\mathsf{I}{}_b^{\alpha }\left[f\left(r\right)\right]=\hspace{0.17em}{\displaystyle \underset{\eta }{\overset{b}{\int }}\left(\eta -r\right)}{\hspace{0.17em}}^{\alpha \left(r,\eta \right)-1}f\left(r\right)\frac{dr}{\mathsf{\Gamma }\hspace{0.17em}\left[\alpha \left(r,\eta \right)\right]}\hspace{0.17em}.$$

(7)

3. The Basic Idea of the Optimal Auxiliary Function Method

In this paper, we successfully applied the OAFM to solve the analytical approximate non-linear fractional solution of partial differential equations. We consider the most general form of a non-linear differential equation.

The following are the general PDEs:

$$\frac{{\partial }^{\alpha }\theta (\eta ,r)}{\partial \hspace{0.17em}{r}^{\alpha }}=\wp (\eta ,r)+N\left(\theta (\eta ,r)\right)\hspace{0.17em}=0,$$

(8)

They are subject to boundary conditions:

$$\begin{array}{l}\frac{{\partial }^{\alpha -k}}{r}\theta (\eta ,r)=h_k(\eta ).(k=0,1,\dots ,n-1)\frac{{\partial }^{\alpha -n}}{r}\theta (\eta ,0)=0,\hspace{0.17em}n=[\alpha ].\\ \frac{{\partial }^k}{r}\theta (\eta ,0)=g_k(\eta ).(k=0,1,\dots ,n-1)\frac{{\partial }^n}{r}\theta (\eta ,0)=0,\hspace{0.17em}n=[\alpha ].\end{array}$$

(9)

In Equation (8), $\frac{{\partial }^{\alpha }}{\partial \hspace{0.17em}{r}^{\alpha }}$ shows the Caputo or R-L operator $\theta (\eta ,r)$ unknown function while $\wp (\eta ,r)$ is a known analytic function.

Step 1: To find the approximate solution of Equation (9), we have to consider the approximate solution in the form of two components shown in Equation (10).

$$\widehat{\theta }\hspace{0.17em}(\eta ,r)={\theta }_0(\eta ,r)+{\theta }_1(\eta ,r,c_{\iota }),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\iota =1,2,3,4,5\dots \delta .$$

(10)

Step 2: We arrange an equation to find the zero- and first-order solution (10) into Equation (8). Its results are

$$\frac{{\partial }^{\alpha }{\theta }_0(\eta ,r)}{\partial r^{\alpha }}+\frac{{\partial }^{\alpha }{\theta }_1(\eta ,r)}{\partial r^{\alpha }}+\wp (\eta ,r)+{\rm N}\left[\frac{{\partial }^{\alpha }{\theta }_0(\eta ,r)}{\partial r^{\alpha }}+\frac{{\partial }^{\alpha }{\theta }_1(\eta ,r),C_i}{\partial r^{\alpha }}\right]=0.$$

(11)

Step 3: The initial approximation ${\theta }_0(\eta ,r)$ can be obtained from the linear equation as

$$\frac{{\partial }^{\alpha }{\theta }_0(\eta ,r)}{\partial r^{\alpha }}+\wp (\eta ,r)=0,$$

(12)

Applying the inverse operator, we obtain ${\theta }_0(\eta ,r)$ as follows:

$${\theta }_0(\eta ,r)={\wp }^{\prime }\left(\eta ,r\right)$$

(13)

Step 4: Expand the non-linear term from Equation (11) in the form of

$${\rm N}\left[\frac{{\partial }^{\alpha }{\theta }_0(\eta ,r)}{\partial r^{\alpha }}+\frac{{\partial }^{\alpha }{\theta }_1(\eta ,r,C_i)}{\partial r^{\alpha }}\right]={\rm N}\left[{\theta }_0(\eta ,r)\right]+{\displaystyle \sum _{k=1}^{\infty }\frac{{\theta }_1^k}{k!}}{{\rm N}}^{(k)}\left[{\theta }_0(\eta ,r)\right].\hspace{0.17em}$$

(14)

Step 5: To solve Equation (14) easily and accelerate the convergence of the first-order approximation, we introduce another expression that can be written as follows:

$$\frac{{\partial }^{\alpha }{\theta }_1(\eta ,r,C_i)}{\partial r^{\alpha }}=-E_1\left[{\theta }_0(\eta ,r)\right]{\rm N}\left[{\theta }_0(\eta ,r)\right]-E_2\left[{\theta }_0(\eta ,r),C_j\right],$$

(15)

Remark 1. 

Where $E_1$ and $E_2$ are two auxiliary functions depending upon ${\theta }_0(\eta ,r)$ and the convergence control parameter $C_{\iota }$ and $C_J\hspace{0.17em},$$\iota =1,2,3,4\dots ,\hspace{0.17em}J=s+1,s+2,\dots \delta$.

Remark 2. 

$E_1$ and $E_2$ are of the form ${\theta }_0(\eta ,r),$$N\left[{\theta }_0(\eta ,r)\right]$, or the combination of both ${\theta }_0(\eta ,r),$ and $N\left[{\theta }_0(\eta ,r)\right]$, but they are not unique.

Remark 3. 

If ${\theta }_0(\eta ,r),$ or $N\left[{\theta }_0(\eta ,r)\right]$ are the functions of a polynomial, then $E_1\left[{\theta }_0(\eta ,r),C_{\iota }\right]$ and $E_2\left[{\theta }_0(\eta ,r),C_j\right]$ are taken as the summation. If ${\theta }_0(\eta ,r),$ or $N\left[{\theta }_0(\eta ,r)\right]$ are the exponential functions, then $E_1\left[{\theta }_0(\eta ,r),C_{\iota }\right]$ and $E_2\left[{\theta }_0(\eta ,r),C_j\right]$ are taken as the functions of addition of an exponential function. If ${\theta }_0(\eta ,r),$ or $N\left[{\theta }_0(\eta ,r)\right]$ are in the form of a trigonometric function, then $E_1\left[{\theta }_0(\eta ,r),C_{\iota }\right]$ and $E_2\left[{\theta }_0(\eta ,r),C_j\right]$ are taken as the addition of a trigonometric function. A special case for If $N\left[{\theta }_0(\eta ,r)\right]=0$ then ${\theta }_0(\eta ,r),$ is the exact solution of Equation (10).

Step 6: To find the square of the residual error to obtain the values of $C_{\iota }$ and $C_j$ we use either the collocation method, the Galerkin method, the Ritz method, or the least square method.

$$H\left(C_{\iota },C_j\right)={\displaystyle \underset{0}{\overset{r}{\int }}{\displaystyle \underset{\mathsf{\Omega }}{\int }{R}^2\left(\eta ,r;C_{\iota },C_j\right)d\eta dr}},$$

(16)

where R is the residual,

$$\begin{array}{l}R\left(\eta ,r,C_i,C_j\right)=\frac{{\partial }^{\alpha }\widehat{\theta }\hspace{0.17em}(\eta ,r),C_i,C_J}{\partial r}+\wp (\eta ,r)+{\rm N}\left[\widehat{\theta }\hspace{0.17em}\left(\eta ,r,C_i,C_j\right)\right],\hspace{0.17em}\hspace{0.17em}i=1,2,3\dots s,\hspace{0.17em}\hspace{0.17em}j=s+1,s+2,\hspace{0.17em}s+3\dots \rho ,\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{\partial J}{\partial \hspace{0.17em}{C}_1}=\frac{\partial J}{\partial \hspace{0.17em}{C}_2}=\dots \frac{\partial J}{\partial \hspace{0.17em}{C}_q}=0.\end{array}$$

(17)

In solving the above equations simultaneously, we obtain the values of the constants $C_{\iota }$.

Remark 4. 

This powerful tool does not depend on small or large parameters. Our procedure consists of auxiliary functions $E_1$ and $E_2$ which control the convergence of the approximate solution after only one iteration.

4. Numerical Experiments and Results

In this section, we present three problems and demonstrate how their numerical results compare with other methods from the literature.

Problem 1. 

Consider the time-fractional equation of Cahn–Hilliard [39]:

$$\begin{array}{l}\frac{{\partial }^{\alpha }\theta (\eta ,r)}{\partial r^{\alpha }}=\mu \frac{\partial \theta (\eta ,r)}{\partial \eta }+6\theta (\eta ,r)\frac{\partial {\theta }^2(\eta ,r)}{\partial \eta }+3{\theta }^2(\eta ,r)\frac{{\partial }^2\theta (\eta ,r)}{\partial {\eta }^2}-\frac{{\partial }^2\theta (\eta ,r)}{\partial {\eta }^2}-\frac{{\partial }^4\theta (\eta ,r)}{\partial {\eta }^4},\hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}0<\alpha \le 1\hspace{0.17em},\end{array}$$

(18)

Then, it is subject to the initial condition:

$${\theta }_0\left(\eta ,0\right)=\tanh \left(\frac{\eta }{\sqrt{2}}\right)$$

(19)

The exact solution of Equation (18) when α, µ = 1 is

$$\theta (\eta ,r)=\tanh \left(\frac{\eta +r}{\left(\sqrt{2}\right)}\right)$$

(20)

So, Equation (18) can be written as

$$\frac{{\partial }^{\alpha }\theta (\eta ,r)}{\partial r^{\alpha }}-\mu \frac{\partial \theta (\eta ,r)}{\partial \eta }-6\theta (\eta ,r)\frac{\partial {\theta }^2(\eta ,r)}{\partial \eta }-3{\theta }^2(\eta ,r)\frac{{\partial }^2\theta (\eta ,r)}{\partial {\eta }^2}+\frac{{\partial }^2\theta (\eta ,r)}{\partial {\eta }^2}+\frac{{\partial }^4\theta (\eta ,r)}{\partial {\eta }^4}=0.$$

(21)

In Equation (21), we take the linear and non-linear parts as

$$\left\{\begin{array}{l}L\left(\theta (\eta ,r)\right)=\frac{{\partial }^{\alpha }\theta (\eta ,r)}{\partial r^{\alpha }}\\ N\left(\theta (\eta ,r)\right)=-\mu \frac{\partial \theta (\eta ,r)}{\partial \eta }-6\theta (\eta ,r)\frac{\partial {\theta }^2(\eta ,r)}{\partial \eta }-3{\theta }^2(\eta ,r)\frac{{\partial }^2\theta (\eta ,r)}{\partial {\eta }^2}+\frac{{\partial }^2\theta (\eta ,r)}{\partial {\eta }^2}+\frac{{\partial }^4\theta (\eta ,r)}{\partial {\eta }^4}.\end{array}\right.$$

(22)

The initial approximation ${\theta }_0(\eta ,r)$ is obtained from Equation (12):

$$\frac{{\partial }^{\alpha }\theta (\eta ,r)}{\partial r^{\alpha }}=0,$$

(23)

In applying the inverse operator as mentioned in Equation (13), we obtain the following solution,

$${\theta }_0\left(\eta ,r\right)=\tanh \left(\frac{\eta }{\sqrt{2}}\right)$$

(24)

By using Equation (24) in Equation (22), the non-linear operator becomes

$$N\left({\theta }_0(\eta ,r)\right)=-\mu \frac{\partial {\theta }_0(\eta ,r)}{\partial \eta }-6{\theta }_0(\eta ,r)\frac{\partial {\theta }_0^2(\eta ,r)}{\partial \eta }-3{\theta }_0^2(\eta ,r)\frac{{\partial }^2{\theta }_0(\eta ,r)}{\partial {\eta }^2}+\frac{{\partial }^2{\theta }_0(\eta ,r)}{\partial {\eta }^2}+\frac{{\partial }^4{\theta }_0(\eta ,r)}{\partial {\eta }^4}.$$

(25)

The first approximation is given by (15):

$$\frac{{\partial }^{\alpha }{\theta }_1(\eta ,r,C_i)}{\partial r^{\alpha }}=-E_1\left[{\theta }_0(\eta ,r)\right]{\rm N}\left[{\theta }_0(\eta ,r)\right]-E_2\left[{\theta }_0(\eta ,r),C_j\right],$$

(26)

Here, we select $E_1$ and $E_2$ according to the non-linear operator,

$$\begin{array}{l}{E}_1=-\left\{\begin{array}{l}{C}_1\left(\frac{1}{\sqrt{2}}-\frac{{\tanh }^2\left(\frac{\eta }{\sqrt{2}}\right)}{\sqrt{2}}\right)+C_1\left(\frac{1}{2}\left(-\tanh \left(\frac{\eta }{\sqrt{2}}\right)+{\tanh }^3\left(\frac{\eta }{\sqrt{2}}\right)\right)\right)r\\ +C_1\left(\frac{1}{12}\left(-\sqrt{2}+4\sqrt{2}{\tanh }^2\left(\frac{\eta }{\sqrt{2}}\right)-3\sqrt{2}{\tanh }^4\left(\frac{\eta }{\sqrt{2}}\right)\right)\right)r^2\\ +C_1\left(\frac{1}{12}\left(2\tanh \left(\frac{\eta }{\sqrt{2}}\right)-5{\tanh }^3\left(\frac{\eta }{\sqrt{2}}\right)+3{\tanh }^5\left(\frac{\eta }{\sqrt{2}}\right)\right)\right)r^3\\ +C_1\left(\frac{1}{120}\left(2\sqrt{2}-17\sqrt{2}{\tanh }^2\left(\frac{\eta }{\sqrt{2}}\right)+30\sqrt{2}{\tanh }^4\left(\frac{\eta }{\sqrt{2}}\right)-15\sqrt{2}{\tanh }^6\left(\frac{\eta }{\sqrt{2}}\right)\right)r^4\right)\end{array}\right\}\\ E_2=0\end{array}$$

(27)

Using Equations (24) and (25) in Equation (26), and apply the inverse operator, we obtain the first approximation as

$${\theta }_1\left(\eta ,r\right)=\frac{1}{120\alpha \hspace{0.17em}\mathsf{\Gamma }\left(\alpha \right)}{C}_1\hspace{0.17em}{r}^{\alpha }Sec{h}^2\left(\frac{\eta }{\sqrt{2}}\right)\hspace{0.17em}\left\{\begin{array}{l}2\sqrt{2}\left(30-5r^2+r^4\right)+20r\left(-3+r^2\right)\tanh \left(\frac{\eta }{\sqrt{2}}\right)\\ -15\sqrt{2}{r}^2\left(-2+r^2\right){\tanh }^2\left(\frac{\eta }{\sqrt{2}}\right)-30r^3{\tanh }^3\left(\frac{\eta }{\sqrt{2}}\right)\\ +15\sqrt{2}{r}^4{\tanh }^4\left(\frac{\eta }{\sqrt{2}}\right)\end{array}\right\}$$

(28)

in adding Equations (24) and (28), we obtain first-order approximate solutions as

$$\begin{array}{c}\widehat{\theta }\hspace{0.17em}(\eta ,r)=\hspace{0.17em}\hspace{0.17em}{\theta }_0(\eta ,r)+{\theta }_0(\eta ,r,C_1).\\ \widehat{\theta }\hspace{0.17em}(\eta ,r)=\tanh \left(\frac{\eta }{\sqrt{2}}\right)+\frac{1}{120\hspace{0.17em}\alpha \hspace{0.17em}\mathsf{\Gamma }\left(\alpha \right)}{C}_1{r}^{\alpha }Sec{h}^2\left(\frac{\eta }{\sqrt{2}}\right)\left\{\begin{array}{l}2\sqrt{2}\left(30-5r^2+r^4\right)+20r\left(-3+r^2\right)\tanh \left(\frac{\eta }{\sqrt{2}}\right)\\ -15\sqrt{2}{r}^2\left(-2+r^2\right){\tanh }^2\left(\frac{\eta }{\sqrt{2}}\right)-30r^3{\tanh }^3\left(\frac{\eta }{\sqrt{2}}\right)\\ +15\sqrt{2}{r}^4{\tanh }^4\left(\frac{\eta }{\sqrt{2}}\right)\end{array}\right\}\end{array}$$

(29)

Below, the results of Problem 1 are presented in

Table 1. The numerical values of the control parameters, with varying values of $\alpha$ for Equation (18).

α = 1		α = 0.8	α = 0.7	α = 0.6
$C_1$	0.9999999999999998	1.0689593321155948	1.1649666232352796	1.3213063996776493

and

Table 2. Error analysis of the OAFM, with NIM and q-HAM using the exact solution (α = µ = 1) for Equation (18).

$\mathit{r}$	$\mathit{\eta }$	OAFM Solution				Exact Solution	NIM Error [39]	q-HAM Error [39]	OAFM Error
		α = 0.6	α = 0.7	α = 0.8	α = 1
0.01	0.0	0.065975	0.0361	0.020385	0.007070	0.007070	1.151971 × 10−7	2.35697 × 10−12	4.77049 × 10−17
	1.0	0.650199	0.63147	0.621632	0.613291	0.613291	1.810671 × 10−7	2.82376 × 10−10	9.99201 × 10−16
	2.0	0.902204	0.89594	0.892655	0.889867	0.889867	6.167394 × 10−8	5.74951 × 10−11	9.99201 × 10−16
	3.0	0.975328	0.97367	0.972799	0.972060	0.972060	1.165205 × 10−9	3.75726 × 10−11	1.11022 × 10−16
0.05	0.0	0.173216	0.11130	0.073843	0.035340	0.035340	1.306675 × 10−5	7.361971 × 10−9	3.72481 × 10−12
	1.0	0.715566	0.67742	0.65435	0.630632	0.630632	2.224480 × 10−5	1.736922 × 10−7	1.39521 × 10−11
	2.0	0.923783	0.91113	0.903476	0.895608	0.895608	7.794449 × 10−6	3.622408 × 10−8	1.85728 × 10−11
	3.0	0.981023	0.97768	0.975656	0.973577	0.973577	1.257660 × 10−7	2.328496 × 10−8	1.11566 × 10−12
0.08	0.0	0.229497	0.15456	0.10748	0.056508	0.056508	4.940148 × 10−5	7.713501 × 10−8	9.99082 × 10−11
	1.0	0.748476	0.70289	0.674246	0.643237	0.643237	8.990891 × 10−5	1.124520 × 10−6	2.17279 × 10−10
	2.0	0.934445	0.91940	0.909956	0.899727	0.899727	3.218897 × 10−5	2.387229 × 10−7	3.08464 × 10−10
	3.0	0.983821	0.97985	0.977359	0.974662	0.974662	4.548965 × 10−7	1.516340 × 10−7	1.89203 × 10−11
0.1	0.0	0.262219	0.18058	0.128409	0.070593	0.070593	9.109940 × 10−5	2.352262 × 10−7	4.76052 × 10−10
	1.0	0.767071	0.71788	0.686336	0.651452	0.651452	1.740220 × 10−4	2.722916 × 10−6	7.86751 × 10−10
	2.0	0.940392	0.92420	0.913853	0.902386	0.902386	6.321236 × 10−5	5.848640 × 10−7	1.16877 × 10−9
	3.0	0.985376	0.98111	0.978381	0.975358	0.975358	8.108096 × 10−7	3.686350 × 10−7	7.26751 × 10−11

and visualized in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6.

Problem 2. 

Consider the fractional Burgers–Poisson Equation [40]:

$$\begin{array}{l}\frac{{\partial }^{\alpha }\theta (\eta ,r)}{\partial r^{\alpha }}-\frac{{\partial }^{\alpha }}{\partial r^{\alpha }}\left(\frac{{\partial }^2\theta (\eta ,r)}{\partial {\eta }^2}\right)+\frac{\partial \theta (\eta ,r)}{\partial \eta }+\theta (\eta ,r)\frac{\partial \theta (\eta ,r)}{\partial \eta }-\left(3\frac{\partial \theta (\eta ,r)}{\partial \eta }\frac{{\partial }^2\theta (\eta ,r)}{\partial {\eta }^2}+\theta (\eta ,r)\frac{{\partial }^3\theta (\eta ,r)}{\partial {\eta }^3}\right)=0,\hspace{0.17em}\hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}r>0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0<\alpha \le 1\hspace{0.17em}.\end{array}$$

(30)

Then, it is subject to the initial condition:

$${\theta }_0\left(\eta ,0\right)=\eta$$

(31)

The exact solution of Equation (30) at $\alpha =1$ is

$$\theta \left(\eta ,r\right)=\frac{1+\eta }{1+r}-1$$

(32)

According to the non-linear operator, we choose $E_1$ and $E_2$ for problem 2.

$$\begin{array}{l}{E}_1=-\left(C_1+3C_2\left(r\right)+2C_3\left(r^2\right)+2C_4{\left(r\right)}^3+2C_5{\left(r\right)}^4\right)\\ E_2=0\end{array}$$

(33)

Using the same procedure of the OAFM method, we obtain the zero-order and first-order solution as

$$\hspace{0.17em}{\theta }_0\left(\eta ,r\right)=\eta ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}$$

(34)

$${\theta }_1\left(\eta ,r\right)=\frac{r^{\alpha }\left(C_1+r\left(3C_2+2r\left(C_3+r\left(C_4+C_{5}r\right)\right)\right)\right)\left(1+\eta \right)}{\alpha \hspace{0.17em}\mathsf{\Gamma }\left[\alpha \right]},$$

(35)

In combining Equations (34) and (35), we obtain the OAFM solution given by the following expression:

$$\begin{array}{l}\widehat{\theta }\hspace{0.17em}(\eta ,r)={\theta }_0\left(\eta ,r\right)+{\theta }_1(\eta ,r,C_1,C_2,C_3,C_4,C_5,C_6).\\ \widehat{\theta }(\eta ,r)=\eta +\frac{r^{\alpha }\left(C_1+r\left(3C_2+2r\left(C_3+r\left(C_4+C_{5}r\right)\right)\right)\right)\left(1+\eta \right)}{\alpha \hspace{0.17em}\mathsf{\Gamma }\left[\alpha \right]}.\end{array}$$

(36)

Below, the results of Problem 2 are presented in

Table 3. Numerical values of convergence control parameters obtained by the Galerkin method for different values of $\alpha$ for Equation (30).

$\mathit{\alpha }=1$		$\mathit{\alpha }=0.7$	$\mathit{\alpha }=0.6$
$C_1$	−0.9999870166919811	−1.1189030434611396	−1.208213789561331
$C_2$	0.33313615275841246	0.9544139953390676	1.341269491284959
$C_3$	−0.4945108518067398	−6.4337069205146395	−10.808604146484456
$C_4$	0.4503627958849215	20.772596091262024	37.507958650955004
$C_5$	−0.2727180163051805	−29.1932425180932	−54.58723568118725

and

Table 4. Comparison of the OAFM, with exact and HPM solutions, and comparison with abs. errors at $\alpha =1$.

$\mathit{\eta }$	$\mathit{r}$	HPM Sol. [40]	OAFM Solution			Exact Solution	Abs. HPM [40]	Abs. OAFM
			$\mathit{\alpha }=1$	$\mathit{\alpha }=0.7$	$\mathit{\alpha }=0.6$
0.9	0.2	0.5804	0.583333	0.342866	0.217758	0.583333	2.933 × 10−3	9.39489 × 10−8
1.2	0.2	0.8304	0.833333	0.554898	0.410035	0.833333	2.933 × 10−3	1.08783 × 10−7
1.5	0.2	1.08	1.08333	0.766929	0.602313	1.08333	3.33 × 10−3	1.23617 × 10−7
1.8	0.2	1.3296	1.33333	0.978961	0.79459	1.33333	3.73 × 10−3	1.38451 × 10−7
2	0.2	1.496	1.5	1.12032	0.922775	1.5	4.0 × 10−3	1.4834 × 10−7
1	0.3	0.526	0.538435	0.271101	0.124476	0.538462	2.462 × 10−3	2.64625 × 10−5
1	0.35	0.45925	0.481359	0.168277	−0.03724	0.481481	2.223 × 10−2	1.2237 × 10−4
1	0.4	0.392	0.428173	−0.00623	−0.35796	0.428571	3.657 × 10−2	3.98493 × 10−4
1	0.45	0.32275	0.378264	−0.31461	−0.96298	0.37931	5.656 × 10−2	1.04634 × 10−3
1	0.5	0.25	0.330963	−0.84290	−2.02297	0.333333	8.333 × 10−2	2.3706 × 10−3

and visualized in Figure 7, Figure 8, Figure 9 and Figure 10.

Problem 3. 

The BBM–Burger equation can be written as [41]

$$\frac{{\partial }^{\alpha }\theta (\eta ,r)}{\partial r^{\alpha }}-\frac{{\partial }^3\theta (\eta ,r)}{\partial {\eta }^2\partial r}+\frac{\partial \theta (\eta ,r)}{\partial \eta }+{\left(\frac{{\theta }^2(\eta ,r)}{2}\right)}_{\eta }=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}r>0,\hspace{0.17em}\eta \in \hspace{0.17em}\hspace{0.17em}I\hspace{0.17em}\hspace{0.17em}\ge R,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\alpha \in \left(0,\hspace{0.17em}1\right],$$

(37)

with the initial conditions of

$${\theta }_0\left(\eta ,0\right)=Sec{h}^2\left(\frac{\eta }{4}\right)\hspace{0.17em}\hspace{0.17em}.$$

(38)

The exact solution of Equation (37) at α = 1 is:

$$\theta (\eta ,r)=\sec h^2\left(\frac{\eta }{4}-\frac{r}{4}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}.$$

(39)

According to the non-linear operator, we choose $E_1$ and $E_2$ for problem number 3:

$$\begin{array}{l}{E}_1=-\left(C_1\left(\frac{1}{2}Sec{h}^2\left(\frac{\eta }{4}\right)\tanh \left(\frac{\eta }{4}\right)\right)\right)\\ E_2=0,\hspace{0.17em}\end{array}$$

(40)

Using the same procedure of the OAFM method, we obtain the zero-order and first-order solution as

$${\theta }_0(\eta ,r)=Sec{h}^2\left(\frac{\eta }{4}\right)\hspace{0.17em}.$$

(41)

$${\theta }_1\left(\eta ,r,C_1\right)=\frac{C_1\hspace{0.17em}{r}^{\alpha }Sec{h}^2\left(\frac{\eta }{4}\right)\tanh \left(\frac{\eta }{4}\right)}{2\hspace{0.17em}\alpha \hspace{0.17em}\mathsf{\Gamma }\left(\alpha \right)}$$

(42)

We obtain the first-order approximate solution by combining Equations (41) and (42).

$$\begin{array}{l}\widehat{\theta }(\eta ,r)={\theta }_0(\eta ,r)+{\theta }_1(\eta ,r,C_1),\\ \widehat{\theta }(\eta ,r)=Sec{h}^2\left(\frac{\eta }{4}\right)+\frac{C_1{r}^{\alpha }Sec{h}^2\left(\frac{\eta }{4}\right)\tanh \left(\frac{\eta }{4}\right)}{2\alpha \hspace{0.17em}\mathsf{\Gamma }\left(\alpha \right)}.\end{array}$$

(43)

Below, the results of Problem 3 are presented in

Table 5. Comparison of OAFM solutions with RPSM solutions and the exact solution for different values of $\alpha$ and comparison with the abs. error of RPSM at $\alpha =1$ for the BBM–Burger equation.

$\mathit{\eta }$	$\mathit{r}$	RPSM [41] Sol. at $\mathit{\alpha }=1$	OAFM Solution at			Exact Solution	Abs. RPSM [41]	Abs. OAFM
			$\mathit{\alpha }=1$	$\mathit{\alpha }=0.7$	$\mathit{\alpha }=0.6$
−15	0.001	0.0022087864	0.00220879	0.0021987	0.00218763	0.00220879	2.44 × 10−9	9.17574 × 10−10
	0.01	0.0021988584	0.00219885	0.0021540	0.00212126	0.00219888	2.41 × 10−8	3.39094 × 10−8
	0.1	0.0021019994	0.00209946	0.0019299	0.00185702	0.00210223	2.29 × 10−7	2.77232 × 10−6
−10	0.001	0.0265787633	0.0265791	0.0264598	0.0263276	0.02657911	349 × 10−7	1.08249 × 10−8
	0.01	0.0264579103	0.026461	0.0259285	0.0255388	0.02646140	3.45 × 10−6	3.94981 × 10−7
	0.1	0.0252801959	0.0252796	0.0232655	0.0223983	0.02531183	3.16 × 10−5	3.21803 × 10−5
−5	0.001	0.2802626306	0.280296	0.2792140	0.2780160	0.28029508	3.33 × 10−5	8.95206 × 10−8
	0.01	0.2788971245	0.279225	0.2743971	0.2708640	0.27922753	3.30 × 10−4	2.72217 × 10−6
	0.1	0.2656871698	0.268514	0.2502531	0.2423990	0.26872355	3.30 × 10−3	2.09181 × 10−4
0	0.001	0.9999997500	1	1	1	1	1.88 × 10−7	6.25000 × 10−8
	0.01	0.9999750000	1	1	1	0.99999475	1.88 × 10−6	6.24997 × 10−6
	0.1	0.9974997396	1	1	1	0.99937526	1.88 × 10−3	6.24740 × 10−4
5	0.001	0.2805672046	0.280534	0.2816160	0.282814	0.28053822	3.34 × 10−5	4.89041 × 10−8
	0.01	0.2819428891	0.281605	0.2864333	0.289966	0.28160625	3.37 × 10−4	1.33948 × 10−6
	0.1	0.2961686088	0.292315	0.3105766	0.318439	0.29251226	3.66 × 10−3	1.96911 × 10−4
10	0.001	0.0266056972	0.0266054	0.0267247	0.0268568	0.02660534	3.49 × 10−7	4.44199 × 10−9
	0.01	0.0267272511	0.0267235	0.0272560	0.0276457	0.02672372	3.52 × 10−6	2.43308 × 10−7
	0.1	0.0279748560	0.0279048	0.0299189	0.0307862	0.02793650	3.84 × 10−5	3.16596 × 10−5
15	0.001	0.0022109987	0.002211	0.0022210	0.0022321	0.00221109	2.44 × 10−9	3.66932 × 10−10
	0.01	0.0022209817	0.00222094	0.0022657	0.0022985	0.00222096	2.46 × 10−8	2.11548 × 10−8
	0.1	0.0023233232	0.00232033	0.0024898	0.0025628	0.00232306	2.59 × 10−7	2.73523 × 10−6

and

Table 6. Comparison between the abs. error of the OAFM with RPSM and error of FHATM [42] at $\alpha =1$ for the BBM–Burger equation.

$\mathit{\eta }$	$\mathit{r}$	Exact Solution	OAFM Solution	Error RPSM [42]	Error [42] FHATM	Error OAFM
10	0.01	2.672 × 10−2	2.672 × 10−2	3523 × 10−6	4529 × 10−5	2.4331 × 10−7
	0.001	2.661 × 10−2	2.661 × 10−2	3492 × 10−7	4501 × 10−6	4.4419 × 10−9
15	0.01	2.221 × 10−3	2.221 × 10−3	2464 × 10−8	3717 ×10−6	2.1155 × 10−8
	0.001	2.211 × 10−3	2.211 × 10−3	2441 × 10−9	3697 × 10−7	3.6690 × 10−10
20	0.01	1.825 × 10−4	1.825 × 10−4	1663 × 10−10	3034 × 10−7	1.7446 × 10−9
	0.001	1.817 × 10−8	1.817 × 10−4	1640 × 10−11	3018 × 10−8	3.0135 × 10−11

and visualized in Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17 and Figure 18.

5. Discussion

This section discusses the results of the optimal auxiliary function method (OAFM) for solving fractional-order equations. In this study, Mathematica 9 was employed for all of our computational work. The accuracy and the validity of the method were evaluated by comparing the results obtained with other analytical methods available in the literature. Table 1 show the numerical values of convergence control parameters obtained using the collocation technique. Table 2 represents the absolute error of the OAFM, comparing it with NIM and q-HAM using the exact solution $\left(\alpha =\mu =1\right)$. Table 3 shows the numerical values of convergence control parameters obtained using the Galerkin method (for the fractional Burgers-Poisson equation, we used the Galerkin method because, for this problem, the collocation method did not provide us with a rapid result).

Similarly, Table 4 compares OAFM results with exact and HPM solutions; the absolute error shows that our method is more accurate than the HPM. These tables demonstrate the accuracy of the OAFM by evaluating the absolute errors for different problems with their exact solutions. Table 5 and Table 6 present a comparative analysis of OAFM solutions with RPSM solutions and exact solutions for various $\alpha$ values, along with a comparison of absolute errors in RPSM at $\alpha =1$, respectively. These comparisons further validate the effectiveness of the OAFM in delivering precise solutions.

Visual representations are provided through 2D and 3D graphs to complement the numerical results. Figure 1 and Figure 2 display 2D graphs for problem 1 with varying values of $\alpha$, while Figure 3, Figure 4, Figure 5 and Figure 6 show 3D graphs for problem 1. Additionally, Figure 7 and Figure 8 show 2D graphs for different values of $\alpha$ at $r=0.2$ for problem 2, and Figure 9 and Figure 10 display 3D plots comparing the exact and OAFM solutions for problem 2. Furthermore, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16 show 3D graphs for problem 3 with different values of $\alpha$, while Figure 17 and Figure 18 show 2D plots for varying $\alpha$ at $r=0.01$ for problem 3.

Overall, the results presented in Table 2, Table 4, Table 5 and Table 6 indicate that, as the value of $\alpha$ approaches 1, the OAFM solution rapidly converges to the exact solution. Based on the presented results, we can confidently say that the OAFM approach delivers remarkably precise solutions.

6. Conclusions

The optimal auxiliary function method (OAFM) is a strong and reliable analytical tool, providing solutions that closely match the exact solutions for various fractional-order equations. The time-fractional Cahn–Hilliard equation, fractional Burgers–Poisson equation, and Benjamin–Bona–Mahony–Burger equations have all been successfully used for finding the approximate solutions by using the OAFM. We can conclude from the numerical results and the presented figures that the proposed method for fractional-order non-linear partial differential equations is very reliable and easy to use. In comparison to NIM, FHATM, q-HAM, RPSM, and HPM, the OAFM approach solutions converge quickly to the exact solution. Based on the mathematical findings, we determined that the proposed method is simple, quick, and effective.