A modified homotopy analysis method for solution of fractional wave equations

Abstract

The aim of this article is to introduce a modified analytical approach to obtain quick and accurate solution of wave-like fractional physical models. This modified analytical approach is an innovative adjustment in Laplace transform algorithm and homotopy analysis method for fractional partial differential equations. The proposed technique solves the problems using Adomian’s polynomials. The homotopy analysis transform method utilizes a simple and powerful method to adjust and control the convergence region of the infinite series solution using an auxiliary parameter. The numerical solutions obtained by this modified proposed method indicate that the approach is easy to implement, highly accurate, and computationally very attractive. A good agreement between the obtained solutions and some well-known results has been obtained.

Introduction

Fractional differential equations have drawn the interest of many researchers^1–4 due to their important applications in science and engineering. Fractional derivatives provide an excellent instrument for the description of memory and hereditary properties of various materials and processes.

Our concern in this work is to consider the numerical solution of the time-fractional wave equations. Wave-like model can describe many physical problems in different fields of science and engineering. These physical problems describe some nonlinear phenomena such as earthquake stresses ⁵ and non-homogeneous elastic waves in soils. ⁶ Recently, there has been a growing interest for obtaining the explicit solutions to heat-like and wave-like models by analytic techniques. The different types of heat-like and wave-like physical models in physics have been solved by Wazwaz and Gorguis, ⁷ Shou and He, ⁸ and Özis and Agırseven. ⁹

In this article, the homotopy analysis transform method (HATM) basically illustrates how the Laplace transform can be used to approximate analytical solutions of the wave-like fractional models by manipulating the homotopy analysis method (HAM). The proposed method is coupling of the HAM and Laplace transform. The main advantage of this proposed method is its capability of combining two powerful methods for obtaining rapid convergent series for fractional partial differential equations. HAM was first proposed and applied by Liao^10–13 based on homotopy, a fundamental concept in topology and differential geometry. The HAM has been successfully applied by many researchers for solving linear and nonlinear partial differential equations.^14–22 In recent years, many researchers have paid attention to obtaining solutions to linear and nonlinear differential and integral equations using various methods by combining the Laplace transform method. Among these, we may mention the following: Laplace decomposition methods^23,24 and homotopy perturbation transform method.^25–28 Recently, Khan et al. ²⁹ have applied to obtain the solutions of the Blasius flow equation on a semi-infinite domain by coupling homotopy analysis with the Laplace transform method.

The main aim of this article presents approximate analytical solutions of wave-like physical models with time-fractional derivative $\alpha$ ( $1<\alpha <2$ ) in the form of a rapidly convergent series with easily computable components using HATM.

Definition 1

The Laplace transform $L[u(x,t)]$ of the Riemann–Liouville fractional integral $I_t^{\alpha }$ is defined as ³

$$L[I_t^{\alpha }u(x,t)]=s^{-\alpha }L[u(x,t)]$$ (1)

Definition 2

The Laplace transform $L[u(x,t)]$ of Caputo’s fractional derivative $D_t^{\alpha }$ is defined as ³

$$\begin{array}{c}\multicolumn{1}{c}{L[D_t^{\alpha }u(x,t)]=s^{\alpha }L[u(x,t)]}\\ \multicolumn{1}{c}{-\sum _{k=0}^{n-1}{s}^{(\alpha -k-1)}{u}^{(k)}(x,0),n-1<\alpha \le n}\end{array}$$ (2)

Basic idea of HATM

To illustrate the basic idea of the HATM for the fractional partial differential equation, we consider the following fractional partial differential equation as

$$D_t^{n\alpha }u(r,t)+R[r]u(r,t)+N[r]u(r,t)=g(r,t),t>0,r\in R^3,n-1<n\alpha \le n$$ (3)

where $D_t^{n\alpha }=({\partial }^{n\alpha }/\partial t^{n\alpha }),$ $R[r]$ is the linear operator in $r\in R^3,$ $N[r]$ is the general nonlinear operator in $r\in R^3$ , and $g(r,t)$ are continuous functions. For simplicity, we ignore all initial and boundary conditions, which can be treated in a similar way. Now, the methodology consists of applying Laplace transform first on both sides of equation (3), and we get

$$L[D_t^{n\alpha }u(r,t)]+L[R[r]u(r,t)+N[r]u(r,t)]=L[g(r,t)]$$ (4)

Now, using the differentiation property of the Laplace transform, we have

$$\begin{array}{c}\multicolumn{1}{c}{L[u(r,t)]-\frac{1}{s^{n\alpha }}\sum _{k=0}^{n-1}{s}^{(n\alpha -k-1)}{u}^k(r,0)}\\ \multicolumn{1}{c}{+\frac{1}{s^{n\alpha }}L(R[r]u(r,t)+N[r]u(r,t)-g(r,t))=0}\end{array}$$ (5)

We define the nonlinear operator

$$\begin{array}{c}\multicolumn{1}{c}{N[\varphi (r,t;q)]=L[\varphi (r,t;q)]-\frac{1}{s^{n\alpha }}\sum _{k=0}^{n-1}{s}^{(n\alpha -k-1)}{u}^k(r,0)}\\ \multicolumn{1}{c}{+\frac{1}{s^{n\alpha }}L(R[r]\varphi (r,t)+N[r]\varphi (r,t)-g(r,t;q))}\end{array}$$ (6)

where $q\in [0,1]$ is an embedding parameter and $\varphi (r,t;q)$ is the real function of $r,t,\mathrm{and}q$ . By means of generalizing the traditional homotopy method, Liao^10–13 constructs the zero-order deformation equation

$$(1-q)\ell [\varphi (r,t;q)-u_0(r,t)]=\hslash \mathrm{qH}(r,t)N[\varphi (r,t;q)],$$ (7)

where $\hslash$ is a nonzero auxiliary parameter, $H(r,t)\ne 0$ is an auxiliary function, $u_0(r,t)$ is an initial guess of $u(r,t)$ , and $\varphi (r,t;q)$ is an unknown function. It is important that one has great freedom to choose auxiliary thing in HATM. Obviously, when $q=0$ and $q=1$ in equation (7), it holds

$$\varphi (r,t;0)=u_0(r,t),\varphi (r,t;1)=u(r,t)$$ (8)

Thus, as $q$ increases from $0$ and $1,$ the solution varies from the initial guess $u_0(x,t)$ to the solution $u(r,t).$ Expanding $\varphi (r,t;q)$ in Taylor’s series with respect to $q,$ we have

$$\varphi (r,t;q)=u_0(r,t)+\sum _{m=1}^{\infty }{q}^m{u}_m(r,t)$$ (9)

where

$$u_m(r,t)=\frac{1}{m!}{\frac{{\partial }^m\varphi (r,t;q)}{\partial q^m}|}_{q=0}$$ (10)

If the auxiliary linear operator $\ell$ , the initial guess $u_0(r,t)$ , the auxiliary parameter $\hslash ,$ and the auxiliary function are properly chosen, the series (9) converges at $q=1,$ and we have

$$u(r,t)=u_0(r,t)+\sum _{m=1}^{\infty }{u}_m(r,t)$$ (11)

which must be one of the solutions of original nonlinear equations.

Define the vectors

$${\overrightarrow{u}}_n=\{u_0(r,t),u_1(r,t),u_2(r,t),\dots ,u_n(r,t)\}$$ (12)

Differentiating equation (7) $m$ time with respect to embedding parameter $q$ and then setting $q=0$ and finally dividing them by $m!,$ we obtain the mth-order deformation equation

$$L[u_m(r,t)-{\chi }_m{u}_{m-1}(r,t)]=\hslash \mathrm{qH}(r,t)R_m({\overrightarrow{u}}_{m-1},x,t)$$ (13)

Operating the inverse Laplace transform on both sides, we get

$$u_m(r,t)={\chi }_m{u}_{m-1}(r,t)+\hslash q{L}^{-1}[H(r,t)R_m({\overrightarrow{u}}_{m-1},r,t)]$$ (14)

where

$$R_m({\overrightarrow{u}}_{m-1},r,t)=\frac{1}{(m-1)!}{\frac{{\partial }^{m-1}\varphi (r,t;q)}{\partial q^{m-1}}|}_{q=0}$$ (15)

and

$${\chi }_m=\{\begin{array}{c}\hfill 0,m\le 1\hfill \\ \multicolumn{1}{c}{1,m>1}\end{array}$$

In this way, it is easy to obtain $u_m(r,t)$ for $m\ge 1,$ and at Mth order, we have

$$u(r,t)=\sum _{m=0}^M{u}_m(r,t)$$ (16)

when $M\to \infty$ we get an accurate approximation of the original equation (3)

$$\begin{array}{c}\multicolumn{1}{c}{u_m(r,t)=({\chi }_m+\hslash )u_{,m-1}-\hslash (1-{\chi }_m)L^{-1}\left(\sum _{k=0}^{j-1}\frac{u^k(x,0)}{s^{k+1}}\right)}\\ \multicolumn{1}{c}{+\hslash L^{-1}\left(\frac{1}{s^{j\alpha }}L(R[t]u(r,t)+N[t]u(r,t)-g(r,t))\right)}\end{array}$$ (17)

In this way, it is easily to obtain $u_m(t)$ for $m\ge 1,$ and at Mth order, we have

$$u_n(r,t)=u_0(r,t)+\underset{k\to \infty }{lim}\sum _{m=1}^k{u}_n(r,t)$$ (18)

where $u_0(r,t)=u(r,0)+\mathrm{tu}\prime (r,0)+t^{2}u″(r,0)+t^{3}u″\prime (r,0)+\cdots +t^{j-1}{u}^{(j-1)}(r,0)=\sum _{i=0}^{j-1}{t}^i{u}^{(i-1)}(r,0).$

The nonlinear term $N_n[t]u_n(t)$ is expanded in terms of Adomian’s polynomials as

$$\begin{array}{c}\multicolumn{1}{c}{N[u(r,t)]=N(\sum _{m=0}^k{u}_m(t))}\\ \multicolumn{1}{c}{=\sum _{m=0}^k{A}_m(u_0,u_1,u_2,u_3,\dots ,u_m)}\end{array}$$ (19)

where $A_m$ are Adomian’s polynomials which are calculated by algorithm (20) constructed by Adomian ³⁰

$$=\frac{1}{n!}{\left[\frac{d^k}{d{t}^k}N\left(\sum _{p=0}^k{\lambda }^p{u}_p\right)\right]}_{\lambda =0},k\ge 0$$ (20)

Combining equations (17) and (20), we get

$$\begin{array}{c}\multicolumn{1}{c}{u_m(r,t)=({\chi }_m+\hslash )u_{m-1}-\hslash (1-{\chi }_m)\sum _{i=0}^{j-1}{t}^i{u}^{(i-1)}(r,0)}\\ \multicolumn{1}{c}{+\hslash L^{-1}\left(\frac{1}{s^{j\alpha }}L(R[t]u_{m-1}(t)+\sum _{k=0}^{m-1}{A}_k(u_0,u_1,u_2,u_3,\dots ,u_m)-g(r,t))\right)·m=1,2,3,\dots }\end{array}$$ (21)

From equation (21), we calculate the various $u_m(t)$ for $m\ge 1$ , and substituting these values in equation (18), we obtain the analytical approximate solution of equation (3). The novelty of our proposed algorithm is that a new correction functional (21) is constructed and the nonlinear term as a series of Adomian’s polynomials in equation (21) is expanded. We combine the new hybrid iterative algorithm for nonlinear fractional partial differential equations arising in science and engineering.

In order to select the optimal value of auxiliary parameter $(\hslash ),$ the averaged residual error is introduced as

$$E_{\alpha }=\frac{1}{(Q+1)(W+1)}\sum _{i=0}^Q\sum _{j=0}^W\left[\mathrm{Re}s(\sum _{z=0}^{\alpha }u(i\Delta x,j\Delta t))\right]$$ (22)

with $Q=W=1,\Delta x=0.05/Q,\mathrm{and}\Delta t=0.05/W$ . It should be stated that the selection of optimal value of the auxiliary parameter in the HATM and HAM has an approximately similar procedures (or more details, see Cheng et al., ³¹ Liao^32,33, and Rashidi et al. ³⁴ ). For the current approximation, the optimal value of $\hslash$ is given by the minimum value of $E_{\alpha }$ corresponding to nonlinear algebraic equation

$$\frac{d{E}_{\alpha }}{d\hslash }=0.$$ (23)

Numerical experiments

In this section, two examples on time-fractional wave equations are solved to demonstrate the performance and efficiency of the HATM.

Example 1

We consider the one-dimensional wave-like equation with variable coefficients^7–9 as

$$D_t^{2\alpha }u(x,t)=\frac{1}{2}{x}^2{u}_{\mathrm{xx}},0<x<1,t>0,0<\alpha \le 1$$ (24)

subject to the Neumann boundary conditions

$$u(0,t)=0,u(1,t)=1+\sin \mathrm{h}t$$ (25)

and the initial conditions

$$u(x,0)=x,u_t(x,0)=x^2$$ (26)

We choose linear operator as

$$\pounds [\varphi (x,t;q)]=L[\varphi (x,t;q)]$$ (27)

with property $\pounds [c]=0,$ where $c$ is constant.

We now define a nonlinear operator as

$$N[\varphi (x,t;q)]=L[\varphi (x,t;q)]-\frac{x}{s}-\frac{x^2}{s^2}-\frac{s^{-2\alpha }}{2}L\left[x^2{\varphi }_{\mathrm{xx}}\right]$$ (28)

Thus, we obtain the mth_-order deformation equation

$$\pounds [u_m(x,t)-{\chi }_m{u}_{m-1}(x,t)]=\hslash R_m({\overrightarrow{u}}_{m-1},x,t)$$ (29)

Operating the inverse Laplace transform on both sides in equation (29), we get

$$u_m(x,t)={\chi }_m{u}_{m-1}(x,t)+\hslash q{L}^{-1}[R_m({\overrightarrow{u}}_{m-1},x,t)]$$ (30)

where

$$\begin{array}{c}\multicolumn{1}{c}{R_m({\overrightarrow{u}}_{m-1},x,t)=L[u_{m-1}(x,t)]}\\ \multicolumn{1}{c}{-(1-{\chi }_m)(\frac{x}{s}+\frac{x^2}{s^2})-\frac{s^{-2\alpha }}{2}L\left[x^2{\left(u_{m-1}\right)}_{\mathrm{xx}}\right]}\end{array}$$ (31)

Now the solution of mth-order deformation equations

$$\begin{array}{c}\multicolumn{1}{c}{u_m(x,t)=({\chi }_m+\hslash )u_{m-1}-(1-{\chi }_m)(x+x^{2}t)}\\ \multicolumn{1}{c}{-\frac{\hslash }{2}{L}^{-1}(s^{-2\alpha }L[x^2{(u_{m-1})}_{\mathrm{xx}}),m\ge 1}\end{array}$$ (32)

Using the initial approximation $u_0(x,t)=u(x,0)+t{u}_t(x,0)=x+x^{2}t$ and the iterative scheme (32), we obtain the various iterates

$$u_1(x,t)=-\frac{x^2\hslash t^{2\alpha +1}}{\Gamma (2\alpha +2)}$$

$$u_2(x,t)=-\frac{\hslash (1+\hslash )x^2{t}^{2\alpha +1}}{\Gamma (2\alpha +2)}+\frac{{\hslash }^2{x}^2{t}^{4\alpha +1}}{\Gamma (4\alpha +2)}$$

$$\begin{array}{c}\multicolumn{1}{c}{u_3(x,t)=-\frac{\hslash {(1+\hslash )}^2{x}^2{t}^{2\alpha +1}}{\Gamma (2\alpha +2)}+\frac{{2\hslash }^2(1+\hslash )x^2{t}^{4\alpha +1}}{\Gamma (4\alpha +2)}}\\ \multicolumn{1}{c}{-\frac{{\hslash }^3{x}^2{t}^{6\alpha +1}}{\Gamma (6\alpha +2)}}\end{array}$$

$$\begin{array}{c}\multicolumn{1}{c}{u_4(x,t)=-\frac{\hslash {(1+\hslash )}^3{x}^2{t}^{2\alpha +1}}{\Gamma (2\alpha +2)}+\frac{{3\hslash }^2{(1+\hslash )}^2{x}^2{t}^{4\alpha +1}}{\Gamma (4\alpha +2)}}\\ \multicolumn{1}{c}{-\frac{{3\hslash }^3(1+\hslash )x^2{t}^{6\alpha +1}}{\Gamma (6\alpha +2)}+\frac{{\hslash }^4{x}^2{t}^{8\alpha +1}}{\Gamma (8\alpha +2)},\dots }\end{array}$$

Proceeding in this manner, the rest of the components $u_n(x,t)$ for $n\ge 5$ can be completely obtained, and the series solutions are thus entirely determined. Consequently, we obtain the solution of equation (24) in a series form

$$u(x,t)=u_0(x,t)+\sum _{m=0}^{\infty }{u}_m(x,t)$$ (33)

However, mostly, the results given by the Laplace decomposition method and homotopy perturbation transform method converge to the corresponding numerical solutions in a rather small region. But, different from those two methods, the HATM provides us with a simple way to adjust and control the convergence region of solution series by choosing a proper value for the auxiliary parameter $\hslash .$ If we select $\hslash =-1,$ then

$$\begin{array}{c}\multicolumn{1}{c}{u(x,t)=x+x^2}\\ \multicolumn{1}{c}{(t+\frac{t^{2\alpha +1}}{\Gamma (2\alpha +2)}+\frac{t^{4\alpha +1}}{\Gamma (4\alpha +2)}+\frac{t^{6\alpha +1}}{\Gamma (6\alpha +2)}+\cdots )}\\ \multicolumn{1}{c}{=x+x^2\sin \mathrm{h}(t^{\alpha },\alpha )}\end{array}$$ (34)

where

$$\sinh (t^{\alpha },\alpha )=\left(t+\frac{t^{2\alpha +1}}{\Gamma (2\alpha +2)}+\frac{t^{4\alpha +1}}{\Gamma (4\alpha +2)}+\frac{t^{6\alpha +1}}{\Gamma (6\alpha +2)}+\cdots \right)$$

When we choose $\alpha =1$ , then clearly we can conclude that the obtained solution $\sum _{m=0}^{\infty }{u}_m(x,t)$ converges to the exact solution $u(x,t)=x+x^2\sin \mathrm{ht}$ The above result is in complete agreement with Özis and Agırseven. ⁹

Figure 1 shows the behavior of obtained solution $u(x,t)$ by the proposed method for different fraction Brownian motions $\alpha =0.7,0.8,0.9$ and for standard motion, that is, at $\alpha =1.$ It is seen from Figure 1 that the solution obtained by new fractional homotopy analysis transform method (FHATM) increases very rapidly with increases in $t$ at the value of $x=0.5,y=0.5,\mathrm{and}\hslash =-1.$

Figure 1.

Plot of solutions at different values of $\alpha =0.7,0.8,0.9,\mathrm{and}1$ for Example 1.

The convergence and rate of approximation for the HAM solution strongly depend on the value of auxiliary parameter $\hslash .$ Even if the initial approximation $u_0(x,t)$ , the auxiliary linear operator $L,$ and the auxiliary function $H(x,t)$ are given, we still have great freedom to choose the value of the auxiliary parameter $\hslash .$ So, the auxiliary parameter $\hslash$ provides us with an additional way to conveniently adjust and control the convergence region and rate of solution series. By means of the so-called $\hslash -\mathrm{curves}$ it is easy to find out the so-called valid regions of $\hslash$ to gain a convergent solution series. When the valid region of $\hslash$ is a horizontal line segment, then the solution is converged.

Figure 2 shows the $\hslash -\mathrm{curve}$ obtained from the 20th-order FHATM approximation solution of time-fractional wave equation (24) at $x=0.5.$ In our study, it is obvious from Figure 2 that the acceptable range of auxiliary parameter $\hslash$ is $-2.0\le \hslash <0.$ We still have freedom to choose the auxiliary parameter according to $\hslash -\mathrm{curve}$ . From Figure 2, the valid regions of convergence correspond to the line segments nearly parallel to the horizontal axis.

Figure 2.

Plot of $\hslash -\mathrm{curve}$ for different values of $\alpha$ for Example 1.

Example 2

In this example, we consider the following two-dimensional wave-like equation with variable coefficients^7–9 as

$$D_t^{2\alpha }u(x,y,t)=\frac{1}{12}\left(x^2{u}_{xx}+y^2{u}_{yy}\right),0<x,y<1,t>0,0<\alpha \le 1$$ (35)

subject to the Neumann conditions

$$\{\begin{array}{c}\hfill u_x(0,y,t)=0,u_x(1,y,t)=4\cos \mathrm{h}t\hfill \\ \multicolumn{1}{c}{u_y(x,0,t)=0,u_y(x,1,t)=4\sin \mathrm{h}t}\end{array}$$ (36)

with initial condition

$$u(x,y,0)=x^4,u_t(x,y,0)=y^4$$ (37)

Applying the Laplace transform on both sides in equation (35) and after using the differentiation property of Laplace transform, we get

$$L[u(x,y,t)]-s^{-1}{x}^4-s^{-2}{y}^4-\frac{1}{12}{s}^{-2\alpha }L[x^2{u}_{\mathrm{xx}}+y^2{u}_{\mathrm{yy}}]=0$$ (38)

We choose the linear operator as

$$\pounds [\varphi (x,y,t;q)]=L[\varphi (x,y,t;q)]$$ (39)

with property $\pounds [c]=0,$ where $c$ is constant.

We now define a nonlinear operator as

$$\begin{array}{c}\multicolumn{1}{c}{N[\varphi (x,y,t;q)]=L[\varphi (x,y,t;q)]-(\frac{x^4}{s}+\frac{y^4}{s^2})}\\ \multicolumn{1}{c}{-\frac{s^{-2\alpha }}{12}L[x^2{\varphi }_{\mathrm{xx}}+y^2{\varphi }_{\mathrm{yy}}]}\end{array}$$ (40)

Using equation (7), with assumption $H(x,y,t)=1,$ we construct the zero-order deformation equation

$$(1-\mathrm{q})\pounds [\varphi (x,y,t;q)-u_0(x,y,t)]=q\hslash N[\varphi (x,y,t;q)]$$ (41)

Obviously, when $q=0$ and $q=1$

$$\varphi (x,y,t;0)=u_0(x,y,t),\varphi (x,y,t;1)=u(x,y,t)$$ (42)

Thus, we obtain the mth-order deformation equation

$$\pounds [u_m(x,y,t)-{\chi }_m{u}_{m-1}(x,y,t)]=\hslash R_m({\overrightarrow{u}}_{m-1},x,y,t)$$ (43)

Operating the inverse Laplace transform on both sides in equation (43), we get

$$u_m(x,y,t)={\chi }_m{u}_{m-1}(x,y,t)+\hslash q{L}^{-1}[R_m({\overrightarrow{u}}_{m-1},x,y,t)]$$ (44)

where

$$\begin{array}{c}\multicolumn{1}{c}{R_m({\overrightarrow{u}}_{m-1},x,y,t)=L[u_{m-1}(x,y,t)]-(1-{\chi }_m)}\\ \multicolumn{1}{c}{(\frac{x^4}{s}+\frac{y^4}{s^2})-\frac{s^{-2\alpha }}{12}L[x^2{\left(u_{m-1}\right)}_{\mathrm{xx}}+y^2{\left(u_{m-1}\right)}_{\mathrm{yy}}]}\end{array}$$ (45)

Now, the solution of mth_-order deformation equations

$$\begin{array}{c}\multicolumn{1}{c}{u_m(x,y,t)=({\chi }_m+\hslash )u_{m-1}-\hslash (1-{\chi }_m)(x^4+y^{4}t)}\\ \multicolumn{1}{c}{-\frac{\hslash }{12}{L}^{-1}\left(s^{-\alpha }L[x^2{\left(u_{m-1}\right)}_{\mathrm{xx}}+y^2{\left(u_{m-1}\right)}_{\mathrm{yy}}]\right),m\ge 1}\end{array}$$ (46)

Using the initial approximation $u_0(x,y,t)=u(x,y,0)+t{u}_t(x,y,0)=x^4+t{y}^4$ and the iterative scheme (46), we obtain the various iterates

$$u_1(x,y,t)=-\frac{\hslash x^4{t}^{2\alpha }}{\Gamma (2\alpha +1)}-\frac{\hslash y^4{t}^{2\alpha +1}}{\Gamma (2\alpha +2)}$$

$$\begin{array}{c}\multicolumn{1}{c}{u_2(x,y,t)=-\frac{\hslash (1+\hslash )x^4{t}^{2\alpha }}{\Gamma (2\alpha +1)}-\frac{\hslash (1+\hslash )y^4{t}^{2\alpha +1}}{\Gamma (2\alpha +2)}}\\ \multicolumn{1}{c}{+\frac{{\hslash }^2{x}^4{t}^{4\alpha }}{\Gamma (4\alpha +1)}+\frac{{\hslash }^2{y}^4{t}^{4\alpha +1}}{\Gamma (4\alpha +2)}}\end{array}$$

$$\begin{array}{c}\multicolumn{1}{c}{u_3(x,y,t)=-\frac{\hslash {(1+\hslash )}^2{x}^4{t}^{2\alpha }}{\Gamma (2\alpha +1)}-\frac{\hslash {(1+\hslash )}^2{y}^4{t}^{2\alpha +1}}{\Gamma (2\alpha +2)}}\\ \multicolumn{1}{c}{+\frac{{2\hslash }^2(1+\hslash )x^4{t}^{4\alpha }}{\Gamma (4\alpha +1)}+\frac{{2\hslash }^2(1+\hslash )y^4{t}^{4\alpha +1}}{\Gamma (4\alpha +2)}}\\ \multicolumn{1}{c}{-\frac{{\hslash }^3{x}^4{t}^{6\alpha }}{\Gamma (6\alpha +1)}-\frac{{\hslash }^3{y}^4{t}^{6\alpha +1}}{\Gamma (6\alpha +2)},\dots }\end{array}$$

Proceeding in this manner, the rest of the components $u_n(x,y,t)$ for $n\ge 4$ can be completely obtained, and the series solutions are thus entirely determined. Hence, the solution of the given problem is given as

$$u(x,y,t)=u_0(x,y,t)+\sum _{m=0}^{\infty }{u}_m(x,y,t)$$ (47)

For $\hslash =-1,$ we have the following solution

$$\begin{array}{c}\multicolumn{1}{c}{u(x,y,t)=x^4(1+\frac{t^{2\alpha }}{\Gamma (2\alpha +1)}+\frac{t^{4\alpha }}{\Gamma (4\alpha +1)}+\cdots )}\\ \multicolumn{1}{c}{+y^4(t+\frac{t^{2\alpha +1}}{\Gamma (2\alpha +2)}+\frac{t^{4\alpha +1}}{\Gamma (4\alpha +2)}+\frac{t^{6\alpha +1}}{\Gamma (6\alpha +2)}+\cdots )}\\ \multicolumn{1}{c}{=x^4\sum _{k=0}^{\infty }\frac{t^{2k\alpha }}{\Gamma (2k\alpha +1)}+y^4\sum _{k=0}^{\infty }\frac{t^{(2k+1)\alpha }}{\Gamma ((2k+1)\alpha +1)}}\\ \multicolumn{1}{c}{=x^4\cos \mathrm{h}(t^{\alpha },\alpha )+y^4\sin \mathrm{h}(t^{\alpha },\alpha )}\end{array}$$ (48)

where the functions $\sin \mathrm{h}(z,\alpha )$ and $\cos \mathrm{h}(z,\alpha )$ are defined as follows

$$\begin{array}{c}\multicolumn{1}{c}{\sinh (z,\alpha )=\frac{E_{\alpha }(z)-E_{\alpha }(-z)}{2}\mathrm{and}\cos \mathrm{h}(z,\alpha )}\\ \multicolumn{1}{c}{=\frac{E_{\alpha }(z)+E_{\alpha }(-z)}{2}}\end{array}$$ (49)

As $\alpha =1,$ this series has the closed form $u(x,y,t)=x^4\cos \mathrm{h}(t)+y^4\sin \mathrm{h}(t),$ which is an exact solution of the classical heat equation. The above result is in complete agreement with Özis and Agırseven. ⁹

Figure 3 shows the behavior of the approximate solution $u(x,y,t)$ for different fraction Brownian motions $\alpha =0.7,0.8,0.9$ and for standard motion, that is, at $\alpha =1$ for Example 2. It is seen from Figure 3 that the solution obtained by new FHATM increases very rapidly with increases in $t$ at the value of $x=0.5,y=0.5,\mathrm{and}\hslash =-1.$ Figure 4 shows the $\hslash -\mathrm{curve}$ obtained from the 20th-order FHATM approximation solution of two-dimensional time-fractional heat equation at $x=0.5$ and $y=0.5.$ In our study, it is obvious from Figure 4 that the acceptable range of auxiliary parameter $\hslash$ is $-2.0\le \hslash <0.$ We still have freedom to choose the auxiliary parameter according to $\hslash -\mathrm{curve}$ . From Figure 4, the valid regions of convergence correspond to the line segments nearly parallel to the horizontal axis.

Figure 3.

Plot of solutions at different values of α = 0.7, 0.8, 0.9, and 1 for Example 2.

Figure 4.

Plot of $\hslash$ curves for different values of $\alpha$ for Example 2.

Concluding remarks

In this article, we carefully proposed a reliable modification of the HAM which introduces a promising tool for solving wave equations of fractional order. Different from all other analytic methods, it provides us with a simple way to adjust and control the convergence region of solution series by introducing an auxiliary parameter. This work illustrates the validity and great potential of the HATM for fractional differential equations. An excellent agreement is achieved. The solution is very rapidly convergent by utilizing the modified HAM by modification of Laplace operator. It may be concluded that the modified HATM methodology is very powerful and efficient in finding approximate solutions as well as analytical solutions of many fractional physical models.