# On the Optimal Auxiliary Linear Operator for the Spectral Homotopy Analysis Method Solution of Nonlinear Ordinary Differential Equations

Authors: S. S. Motsa
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The purpose of this study is to identify the auxiliary linear operator that gives the best convergence and accuracy in the implementation of the spectral homotopy analysis method (SHAM) in the solution of nonlinear ordinary differential equations. The auxiliary linear operator is an essential element of the homotopy analysis method (HAM) algorithm that strongly influences the convergence of the method. In this work we introduce new procedures of defining the auxiliary linear operators and compare solutions generated using the new linear operators with solutions obtained using well-known linear operators. The applicability and validity of the proposed linear operators is tested on four highly nonlinear ordinary differential equations with fluid mechanics applications that have recently been reported in the literature. The results from the study reveal that the new linear operators give better results than the previously used linear operators. The identification of the optimal linear operator will direct future research on further applications of HAM-based methods in solving complicated nonlinear differential equations.


---

## Body

## 1. Introduction

The homotopy analysis method (HAM) is an analytic method that has been widely used for solving highly nonlinear equations with applications in computational and applied mathematics, economics and finance, engineering, and many other areas of fundamental science. A systematic exposition of features of the HAM and its application can be found in recent books by Liao [1–3] and Vajravelu and Van Gorder [4]. A distinctive characteristic of the HAM that sets it apart from all other analytical methods is the presence of a convergence-controlling parameter and the flexibility to select auxiliary functions and linear operators in order to guarantee convergence and improve accuracy of the approximate solutions. This makes the HAM suitable for solving many highly nonlinear problems including those that do not contain small or large embedded parameters.A discrete numerical approach that uses the principles of the traditional HAM and combines them with the Chebyshev spectral collocation method was introduced by Motsa et al. [5, 6] and called spectral homotopy analysis method (SHAM). In the traditional HAM approach the options of initial approximations and auxiliary linear operators and functions are limited by the requirement that the obtained approximate solution must be a continuous analytical series solution. In the SHAM application, this restriction is removed and, consequently, the SHAM can accommodate infinite options of initial approximations, auxiliary linear operators, and convergence-controlling auxiliary functions. This makes the SHAM more robust and capable of solving a wider range of complicated nonlinear equations than its analytical counterpart.The effect of different auxiliary linear operators on the accuracy and convergence of the HAM has been studied by different authors in the past. Van Gorder and Vajravelu [7] presented different methods for selecting auxiliary linear operators and gave the necessary and sufficient conditions for the convergence of the HAM series solutions. The aim of the illustrative study of [7] was to demonstrate how different linear operators can be obtained. A comparative analysis of the accuracy and convergence of the approximate solutions obtained using their proposed linear operators was not considered. Shan and Chaolu [8] proposed a linear operator constructed by forming a differential equation using the guessed initial solution that satisfies the boundary conditions. The method proposed in [8] was demonstrated using the Thomas Fermi equation and a nonlinear heat transfer equation. Van Gorder [9] discussed methods for selecting auxiliary linear operators in solving the Painlevé equations, Lane-Emden equation, and a third order equation that models a nonlinear stretching sheet in a fluid flow problem. It was noted that, for some auxiliary linear operators, the rate of convergence may be slow and for others the rate of convergence was rapid. All the studies reviewed so far, however, suffer from the fact that the linear operators have to be carefully selected in such a way that the solution method gives analytical solutions made up of simple functions such as polynomials, decaying exponential, sines, cosines, rational functions, and other elementary functions or products of such functions which are fairly easy to integrate. A particular method of selecting a linear operator such as the method of complete differential matching used in [7, 9] can only be used in very few nonlinear equations if the output solution is expected to be a simple series function. Despite being suggested in [7] as a choice of linear operator that, potentially, can give better convergence and accuracy, the method of complete differential matching has not been widely used in the HAM literature because it often leads to higher order decomposed equations that cannot be integrated analytically.In seeking to identify auxiliary parameters and operators that give the best accuracy and convergence of the HAM algorithm, a more suitable study would seek to identify the optimal parameters from the infinitely many possibilities and not to search from the small set of options that guarantee exact solutions of the decomposed equations. The aim of this paper is to determine the auxiliary linear operator and initial approximation that give the optimal accuracy and convergence in the HAM solution of highly nonlinear boundary value problems defined in bounded domains. The study gives a systematic way of selecting initial approximations, auxiliary linear operators, and convergence controlling parameter. New methods of defining linear operators in the context of the discrete spectral method based version of the HAM (SHAM) are proposed. The accuracy and convergence of the SHAM algorithms derived using the proposed new linear operators are validated on the well-known nonlinear differential equations that have been previously solved using the HAM/SHAM algorithms with the standard linear operators. In particular, we consider the Darcy-Brinkman-Forchheimer model [5, 10], Jeffery-Hamel equation [6, 11–13], the laminar viscous flow in a semiporous channel subject to a transverse magnetic field [14–16], and the model of two-dimensional viscous flow in a rectangular domain bounded by two moving porous walls [17–20].The remainder of the paper is organized as follows. In Section2, we give a brief description of the SHAM for general nonlinear ordinary differential equations. Section 3 introduces the auxiliary linear operators that are used in the SHAM algorithm. Section 4 describes the application of the SHAM with the proposed linear operators in the problems that are selected for numerical experimentation. The numerical simulations and results are presented in Section 5. Finally, we conclude and describe future work in Section 6.


## 2. Description of the Spectral Homotopy Analysis Method in a General Setting

In this section we give the basic description of the homotopy analysis method for the solution of a general nonlinear differential equation of orderp (where p is a positive integer) given by
(1)∑j=0p‍αj(x)y(j)(x)+N(y(x),y′(x),y′′(x),…,y(p)(x))=Φ(x),
where αj(x) and Φ(x) are known functions of the independent variable x and y(x) is an unknown function. The primes in (1) denote differentiation with respect to x and y(0)≡y(x). The function N represents the nonlinear component of the governing equation. For illustrative purposes, we assume that (1) is to be solved in the domain x∈[a,b] subject to the separated boundary conditions
(2)Ba(y(a))=0,Bb(y(b))=0,
where Ba and Bb are linear operators.In the framework of the homotopy analysis method (HAM) [1–4], we define the following zeroth order deformation equations:
(3)(1-q)L[Y(x;q)-y0(x)]=qħ{N[Y(x;q)]-Φ(x)},
where q∈[0,1] denotes an embedding parameter, Y(x;q) is a kind of continuous mapping function of y(x), and ħ is the convergence-controlling parameter. The nonlinear operator N is defined from the governing equation (1) as
(4)N[Y(x;q)]=L[Y(x;q)]+N[Y(x;q)],
with
(5)L[Y(x;q)]=∑j=0p‍αj(x)y(p)(x).By differentiating the zeroth order equations (3) m times with respect to q, setting q=0, and finally dividing the resulting equations by m!, we obtain the following mth order deformation equations:
(6)L[ym(x)-(χm+ħ)ym-1(x)]=ħRm-1[y0,y1,…,ym-1],
where
(7)Rm-1[y0,y1,…,ym-1]=1(m-1)!∂m-1{N[Y(x;q)]-Φ(x)}∂qm-1|q=0,χm={0,m⩽1,1,m>1.From the solutions of (6), the approximate solution for y(x) is determined as the series solution
(8)y(x)=∑k=0+∞‍yk(x).A HAM solution is said to be of orderK if the above series is truncated at k=K, that is, if
(9)y(x)=∑m=0K‍ym(x).For nontrivial linear operators, the higher order equations (6) cannot be integrated using analytical means. Consequently, numerical approaches are employed. When the Chebyshev spectral collocation method is used to solve (6), the method is called spectral homotopy analysis method [3, 5, 6].In using the SHAM, the initial guess is obtained simply as a solution of the linear part of the governing equation (1) subject to the underlying boundary conditions (2). That is, we solve
(10)∑j=0p‍αj(x)y(p)(x)=Φ(x).In essence, collocation approximates the solutiony(x) using an interpolating polynomial of degree M which satisfies the boundary conditions and the differential equations at all points, called the collocation points, xj, where j=0,1,…,M. In the Chebyshev spectral collocation method, the collocation points are chosen to be the extrema of Chebyshev polynomials of degree TM on the interval -1≤η≤1 defined as
(11)ηj=cos⁡(πjM),j=0,1,…,M.
We use the transformation x=(b-a)(η+1)/2 to map the interval [a,b] to [-1,1]. The so-called differentiation matrix D is used to approximate the derivatives of the unknown variables y(x) at the collocation points as the matrix vector product
(12)dydx|x=xj=∑k=0M‍Djky(ηk)=DY,j=0,1,…,M,
where D=2D/(b-a) and
(13)Y=[y(η0),y(η1),…,y(ηM)]T
is the vector function at the collocation points. Higher order derivatives are obtained as powers of D; that is,
(14)y(p)(xj)=DpY,j=0,1,…,M.The matrixD is of size (M+1)×(M+1) and its entries are defined in [21, 22]. In the SHAM algorithm the continuous derivatives of the higher order deformation equations are replaced by the discrete Chebyshev differentiation matrices. As a result, the higher order deformation equations reduce to matrix equations and are solved using standard techniques for solving linear systems of equations.


## 3. Definition and Selection of Linear Auxiliary Operators

In this section we describe the different approaches used to define the auxiliary linear operators in the SHAM algorithm. In order to highlight the subtle differences between the variety of linear operators, it is convenient to express the nonlinear operator of the governing equation (1) using a sum formula of the derivative products as follows:
(15)N(y,y′,…,y(p))=yβ0,0N0,0(y)+y′∑j=01‍β1,jN1,j(y,…,y(j))+y′′∑j=02‍β2,jN2,j(y,y′,…,y(j))+⋯+y(p-1)∑j=0p-1‍βp-1,jNp-1,j(y,y′,…,y(j))+y(p)∑j=0p‍βp,jNp,j(y,y′,…,y(j))=∑r=0p‍y(r)∑j=0r‍βr,jNr,j(y,y′,y′′,…,y(j)),
where βr,j(x), (j=0,1,…,p) are coefficients of the nonlinear terms containing yr as the highest derivative. In this way, (1) becomes
(16)∑j=0p‍αj(x)y(j)(x)+∑r=0p‍y(r)(x)∑j=0r‍βr,jNr,j(y(x),y′(x),…,y(j)(x))=Φ(x).The selection of auxiliary linear operators in the application of the homotopy analysis method was described by Van-Gorder and Vajravelu [7] in a fairly general setting. In particular, three methods, namely, the method ofhighest order differential matching,linear partition matching, andcomplete differential matching, were defined in [7]. Below, we give a definition of the linear operator selection of [7] and present other options that can be implemented under the SHAM algorithm for any nonlinear differential equation that can be represented by (16).(i) Method of the Highest Order Differential Matching. This approach defines the auxiliary linear operator using only the highest order derivative. This approach is the most widely used approach in the HAM solution of nonlinear differential equations defined in finite domains. With reference to (16), the linear operator is selected as
(17)L1[y]=dpydxp.(ii) Method of Linear Partition Matching. This approach defines the auxiliary linear operator to be the collection of all linear terms in the governing equation. With reference to (16), we set
(18)L2[y]=∑j=0p‍αj(x)djydxj.(iii) Method of Complete Differential Matching. In this approach, the collection of all linear operators and some nonlinear factors of the governing nonlinear equations are used to define the linear operator. The aim is to ensure that all the terms of the governing nonlinear differential equation contribute to the auxiliary linear operator selected. The following rules were defined in [7]:(a)
in the case where we have a term which is the product of derivatives, we take the higher order derivative in the term;(b)
if the term has a product of derivatives with functions of the unknown function, we again take the highest order derivative in the term;(c)
in the case where we have a nonlinear expression in just the unknown function, we take the function itself.Thus, applying the above rules to (16), we set
(19)L3[y]=∑j=0p‍αj(x)djydxj+∑r=0p‍drydxr∑j=0r‍βr,j.(iv) Linear Partition Mapping after Transformation. This approach uses linear partition mapping after the transformation y(x)=z(x)+y0(x), where the function y0(x) is carefully chosen to satisfy the underlying boundary conditions. In earlier versions of the SHAM [5, 6], it was suggested that y0(x) be defined as a solution of (10) which is solved subject to the given boundary conditions. Substituting in the governing equation (1) to express the equation in terms of z(x) gives
(20)∑j=0p‍αj(x)z(j)(x)+∂N∂z(y0,y0′,…,y0(p))z+∂N∂z′(y0,y0′,…,y0(p))z′+⋯+∂N∂z(p)(y0,y0′,…,y0(p))z(p)+NonlinearTerms=Φ(x)-∑j=0p‍αj(x)y0(j)(x)-N(y0,y0′,…,y0(p)).Thus, the linear operator is set to be(21)L4[z]=∑j=0p‍αj(x)z(j)(x)+∂N∂z(y0,y0′,…,y0(p))z+⋯+∂N∂z(p)(y0,y0′,…,y0(p))z(p).(v) Modified Complete Differential Matching Method. This is a new approach that is proposed as a modification of the method of complete differential matching of [7]. In the original approach [7], from the group of nonlinear terms of the governing equation, only the highest derivative is selected. In the proposed approach we want to ensure that all the terms of the nonlinear groups that make up the terms of the differential equation contribute to the linear operator. We define the following rules which are a modification of the rules set in [7].(a)
In the case where we have a term which is the product of derivatives, we take the higher order derivative in the term and approximate each function in the remaining derivative product byy0(x). For example, if the nonlinear product is y′(x)y′′(x)y′′′(x), we set y0′(x)y0′′(x)y′′′(x), where y0(x) is the solution of (10).(b)
If the term has a product of derivatives with functions of the unknown function, we again take the highest order derivative in the term and approximate each function in the remaining derivative product byy0(x). For example, if the nonlinear product is y2(x)y′′′(x), we set y02(x)y′′′(x), where y0(x) is the solution of (10).(c)
In the case where we have a nonlinear expression in just the unknown function, we take the function itself and approximate the rest of the functions byy0(x). For example, if the nonlinear product is y2(x), we set y0(x)y(x), where y0(x) is the solution of (10).


## 4. Numerical Experiments

In this section we illustrate how using a different linear operator can significantly improve convergence and accuracy of the SHAM. For illustration purposes we consider examples of nonlinear differential equations that have previously been solved in the published literature using the HAM/SHAM with standard linear operators. It is worth mentioning that in the previous HAM-based investigations the linear operators were chosen in such a way that(i)
the obtained approximate solution is a continuous analytical expression;(ii)
the HAM series solution conforms to a predefined rule of solution expression.The spectral method based approach of [5, 6] was aimed at removing the above restrictions by considering linear operators which are defined from part of the governing differential equation. In this case, it was found that the linearised deformation equations could not be solved exactly and hence the use of the spectral method. Below, we give a systematic procedure for defining the linear operators for various nondifferential equations selected for numerical experimentation.


### 4.1. Example 1: Darcy-Brinkman-Forchheimer Equation

We consider the following Darcy-Brinkman-Forchheimer equation that models the steady state pressure driven fully developed parallel flow through a horizontal channel that is filled with porous media:(22)y′′(x)-s2y(x)-Fsy(x)2+1M=0,y(-1)=0,y(1)=0,
where F is the dimensionless Forchheimer number and s is the porous media shape parameter. This problem was previously solved using the SHAM with the auxiliary linear operator selected using the method of linear partition matching L2 in [5]. The equivalent problem cast in cylindrical coordinates was subsequently solved in [10] where the linear operator was selected as the linear part of the governing equation excluding the linear terms with variable coefficients. For this reason, we remark that the linear operator used in [10] is not part of the class of linear operators considered in this study. Using the systematic choices described above, we set
(23)L1[y]=y′′,L2[y]=y′′-s2y,L3[y]=y′′-s2y-Fsy,L4[z]=z′′-s2z-2Fsy~0z,L5[y]=y′′-s2y-Fsy~0y.
The linear operator L4 is obtained by using the linear partition mapping after using the transformation y(x)=z(x)+y~0(x). The resulting equation is
(24)z′′-s2z-2Fsy~0z-Fsz2=-(y~0′′-s2y~0-Fsy~02+1M).
The function y~0(x) is chosen as a function that satisfies the boundary conditions and the linear part of the (22). That is, y~0(x) must be a solution of
(25)y~0′′-s2y~0=-1M.
Thus, the initial approximation is chosen to be
(26)y~0(x)=1s2M(1-cosh⁡(sx)cosh⁡(s)).
The initial approximation used in the SHAM algorithms that use the linear operators L1, L2, and L3 is obtained by solving the linear equations:
(27)Li(y0)=-1M,y0(-1)=y0(1)=0,fori=1,2,3.Thus, the initial approximations corresponding to the SHAM withL1, L2, and L3 are given, respectively, by
(28)y0(x)=12M(1-x2),y0(x)=1s2M(1-cosh⁡(sx)cosh⁡(s)),y0(x)=1s12M(1-cosh⁡⁡(s1x)cosh⁡⁡(s1)),
where s1=s2+Fs.The initial approximationz0(x) to use in conjunction with L4 is obtained as a solution of the differential equation formed from the linear part of the (24). That is, we solve
(29)z0′′-s2z0-2Fsy~0z0=-(y~0′′-s2y~0-Fsy~02+1M),z0(-1)=z0(1)=0.
This equation is a linear equation with variable coefficients. Thus, it is solved using the spectral method as described in the earlier section above. Similarly, the initial approximation to use with L5 is obtained as a solution of the differential equation
(30)y0′′-s2y0-Fsy~0y0=-1M,y0(-1)=y0(1)=0,
with y~0 given by (26).With the nonlinear operatorNi defined as the homogeneous part of the governing nonlinear equation in each case, we obtain the following higher order deformation equations:
(31)Li[ym(x)-χmym-1(x)]=ħRi,m-1[y0,y1,…,ym-1],i=1,2,3,5,
where
(32)Ri,m-1=ym-1′′-s2ym-1+1M(1-χm)-Fs∑n=0m-1‍ynym-1-n.The nonlinear operator and higher order deformation equation corresponding toi=4 is defined as
(33)N4(z)=L4(z)-Fsz2,L4[zm(x)-χmzm-1(x)]=ħR4,m-1[z0,z1,…,zm-1],
where
(34)R4,m-1=L4(zm-1)-(1-χm)ϕ1(x)-Fs∑n=0m-1‍znzm-1-n
and ϕ1(x) is the right hand side of (29).


### 4.2.  Example 2: Jeffery-Hamel Equation

Here, we consider the Jeffery-Hamel equation that models the steady two-dimensional flow of an incompressible conducting viscous fluid between two rigid plane walls that meet at an angle2α. The rigid walls are considered to be divergent if α>0 and convergent if α<0. The governing equation is defined as
(35)y′′′+2αReyy′+4α2y′=0,
subject to the boundary conditions
(36)y(0)=1,y′(0)=0,y(1)=0,
where Re is the Reynolds number. This problem was investigated using the standard HAM with the highest order differential matching linear operator L1(y) in [11, 12]. A related method, called optimal homotopy asymptotic method (OHAM), was used to solve the same problem in [13], again with L1(y). In Motsa et al. [6], the SHAM was used with the auxiliary linear operator defined using the transformed linear partition mapping method L4(z). In this example, we explore the auxiliary linear operators defined as
(37)L1[y]=y′′′,L2[y]=y′′′+4α2y′,L3[y]=y′′′+4α2y′+2αRey′,L4[z]=z′′′+4α2z′+2αRey~0z′+2αRey~0′z,L5[y]=y′′′+4α2y′+2αRey~0y′.The functiony~0(x) is chosen as a function that satisfies the boundary conditions and the linear part of (35). Thus, y~0(x) is
(38)y~0(x)=cos⁡(2α)-cos⁡(2αx)cos⁡(2α)-1.The linear operatorL4 is obtained by using the linear partition mapping after using the transformation y(x)=z(x)+y~0(x). The resulting equation is
(39)z′′′+4α2z′+2αRey~0z′+2αRey~0′z+2αRezz′=ϕ2(x)
subject to
(40)z(0)=0,z′(0)=0,z(1)=0,
where
(41)ϕ2(x)=-(y0′′′+2αRey0y0′+4α2y0′).
The initial approximation used in the SHAM algorithms that employ the linear operators L1, L2, and L3 is obtained by solving the linear equations:
(42)Li(y0)=0,y0(0)=1,y0′(0)=0y0(1)=0,fori=1,2,3.
Thus, the initial approximations corresponding to the SHAM with L1, L2, and L3 are given, respectively, by
(43)y0(x)=1-x2,y0(x)=cos⁡(2α)-cos⁡(2αx)cos⁡(2α)-1,y0(x)=cos⁡⁡(2α1)-cos⁡⁡(2α1x)cos⁡⁡(2α1)-1,
where s1=4α2+2αRe. The initial approximation z0(x) to use with L4 is obtained as a solution of the differential equation formed from the linear part of (39). That is, we solve
(44)z0′′′+4α2z0′+2αRey~0z0′+2αRey~0′z0=ϕ2(x).
Since the above equation has variable coefficients, it is solved using the spectral method. Similarly, the initial approximation to use with L5 is obtained as a solution of the differential equation
(45)y0′′′+4α2y0′+2αRey~0y0′=0,y0(0)=1,y0′(0)=0,y0(1)=0
with y~0 given by (38).With the nonlinear operatorNi defined as the homogeneous part of the governing nonlinear equation in each case, we obtain the following higher order deformation equations
(46)Li[ym(x)-χmym-1(x)]=ħRi,m-1[y0,y1,…,ym-1],i=1,2,3,5,
where
(47)Ri,m-1=ym-1′′′+4α2ym-1′+2αRe∑n=0m-1‍ynym-1-n′.The nonlinear operator and higher order deformation equation corresponding toi=4 is defined as
(48)N4(z)=L4(z)+2αRezz′,L4[zm(x)-χmzm-1(x)]=ħR4,m-1[z0,z1,…,zm-1],
where
(49)R4,m-1=L4(zm-1)-(1-χm)ϕ2(x)+2αRe∑n=0m-1‍znzm-1-n′.


### 4.3. Example 3: Laminar Viscous Flow in a Semiporous Channel Subject to a Transverse Magnetic Field

In this section, we consider the problem of laminar viscous flow in a semiporous channel subject to a transverse magnetic field given by(50)y′′′′-Ha2y′′+Re(yy′′′-y′y′′)=0,
subject to the boundary conditions
(51)y=0,y′=0,x=0,y=1,y′=0,x=1.This example was previously investigated using the standard HAM with the highest order differential matching linear operatorL1(y) in [14, 15]. In [16] a blend between the SHAM and a spectral method based quasilinearisation method was used to solve the same problem. In this investigation, we consider the following auxiliary linear operators:
(52)L1[y]=y′′′′,L2[y]=y′′′′-Ha2y′′,L3[y]=y′′′′-Ha2y′′+Re(y′′′-y′′),L4[z]=z′′′′-Ha2z′′+Re(y~0z′′′+y~0′′′z-y~0′z′′-y~0′′z′),L5[y]=y′′′′-Ha2y′′+Re(y~0y′′′-y~0′y′′)=0.
The initial approximation used in the SHAM algorithms that use the linear operators L1, L2, and L3 is obtained by solving the linear equations:
(53)Li(y0)=0,y0(0)=0,y0′(0)=0,y0(1)=1,y0′(1)=0,fori=1,2,3.
Thus, the initial approximations corresponding to the SHAM with L1 and L2 are given by
(54)y0(x)=3x2-2x2,y0(x)=Hax-eHax+eHa-Hax+eHa(Hax-1)+1eHa(Ha-2)+Ha+2.
We note that in this example y~0 is the second solution in (54). The initial approximations corresponding to L3 and L5 are obtained by solving (53) numerically using the spectral method. Similarly, z0(x) is obtained by numerically solving
(55)z0′′′′-Ha2z0′′+Re(y~0z0′′′+y~0′′′z0-y~0′z0′′-y~0′′z0′)=ϕ3(x)
subject to the boundary conditions
(56)z0(0)=0,z0′(0)=0,z0(1)=0,z0′(1)=0,
where
(57)ϕ3(x)=-(y~0′′′′-Ha2y~0′′+Re(y~0y~0′′′-y~0′y~0′′)).
The higher order deformation equations are defined as
(58)Li[ym(x)-χmym-1(x)]=ħRi,m-1[y0,y1,…,ym-1],i=1,2,3,5,
where
(59)Ri,m-1=ym-1′′′′-Ha2ym-1′′+Re∑n=0m-1‍(ynym-1-n′′′-yn′ym-1-n′′).The nonlinear operator and higher order deformation equation corresponding toi=4 is defined as
(60)N4(z)=L4(z)+Re(zz′′′-z′z′′),L4[zm(x)-χmzm-1(x)]=ħR4,m-1[z0,z1,…,zm-1],
where
(61)R4,m-1=L4(zm-1)-(1-χm)ϕ3(x)+Re∑n=0m-1‍(znzm-1-n′′′-zn′zm-1-n′′).


### 4.4. Example 4: Two-Dimensional Viscous Flow in a Rectangular Domain Bounded by Two Moving Porous Walls

In this section we consider the problem of two-dimensional viscous flow in a rectangular domain bounded by two moving porous walls. The governing equations are given in [17, 23] as
(62)y′′′′+α(xy′′′+3y′′)+Re(yy′′′-y′y′′)=0,
subject to the boundary conditions
(63)y=0,y′′=0,atx=0,y=1,y′=0,atx=1,
where α is the nondimensional wall dilation rate defined to be positive for expansion and negative for contraction and Re is the permeation Reynolds number defined positive for injection and negative for suction through the walls. The nonlinear equation (62) was investigated using the standard HAM with the highest order differential matching linear operator L1(y) in [17, 18]. Rashidi et al. [19] used the optimal homotopy asymptotic method (OHAM) to solve the same problem, again using L1(y) as linear operator. The SHAM with transformed linear partition method L4(z) was used in [20].In this investigation, we consider the following auxiliary linear operators:(64)L1[y]=y′′′′,L2[y]=y′′′′+α(xy′′′+3y′′),L3[y]=y′′′′+α(xy′′′+3y′′)+Re(y′′′-y′′),L4[z]=z′′′′+α(xz′′′+3z′′)+Re(y~0z′′′+y~0′′′z-y~0′z′′-y~0′′z′),L5[y]=y′′′′+α(xy′′′+3y′′)+Re(y~0y′′′-y~0′z′′).The initial approximations used in the SHAM algorithms that use the linear operatorsL1, L2, L3, and L5 are obtained by solving the linear equations:
(65)Li(y0)=0,y0(0)=0,y0′(0)=0,y0(1)=1,y0′(1)=0,fori=1,2,3,5.Thus, the initial approximationy0(x) corresponding to L1 is as follows:
(66)y0(x)=12(3x-x3).
The initial approximations corresponding to L2, L3, and L5 are obtained by numerically solving (65). Similarly, z0(x) is obtained by numerically solving
(67)z0′′′′+α(xz0′′′+3z0′′)+Re(y0z0′′′+y0′′′z0-y0′z0′′-y0′′z0′)=ϕ4(x),
where
(68)ϕ4(x)=-[y0′′′′+α(xy0′′′+3y0′′)+Re(y0y0′′′-y0′y0′′)].The higher order deformation equations are defined as(69)Li[ym(x)-χmym-1(x)]=ħRi,m-1[y0,y1,…,ym-1],i=1,2,3,5,
where
(70)Ri,m-1=ym-1′′′′+α(xym-1′′′+3ym-1′′)+Re∑n=0m-1‍(ynym-1-n′′′-yn′ym-1-n′′).The nonlinear operator and higher order deformation equation corresponding toi=4 is defined as
(71)N4(z)=L4(z)+Re(zz′′′-z′z′′),L4[zm(x)-χmzm-1(x)]=ħR4,m-1[z0,z1,…,zm-1],
where
(72)R4,m-1=L4(zm-1)-(1-χm)ϕ4(x)+Re∑n=0m-1‍(znzm-1-n′′′-zn′zm-1-n′′).


## 4.1. Example 1: Darcy-Brinkman-Forchheimer Equation

We consider the following Darcy-Brinkman-Forchheimer equation that models the steady state pressure driven fully developed parallel flow through a horizontal channel that is filled with porous media:(22)y′′(x)-s2y(x)-Fsy(x)2+1M=0,y(-1)=0,y(1)=0,
where F is the dimensionless Forchheimer number and s is the porous media shape parameter. This problem was previously solved using the SHAM with the auxiliary linear operator selected using the method of linear partition matching L2 in [5]. The equivalent problem cast in cylindrical coordinates was subsequently solved in [10] where the linear operator was selected as the linear part of the governing equation excluding the linear terms with variable coefficients. For this reason, we remark that the linear operator used in [10] is not part of the class of linear operators considered in this study. Using the systematic choices described above, we set
(23)L1[y]=y′′,L2[y]=y′′-s2y,L3[y]=y′′-s2y-Fsy,L4[z]=z′′-s2z-2Fsy~0z,L5[y]=y′′-s2y-Fsy~0y.
The linear operator L4 is obtained by using the linear partition mapping after using the transformation y(x)=z(x)+y~0(x). The resulting equation is
(24)z′′-s2z-2Fsy~0z-Fsz2=-(y~0′′-s2y~0-Fsy~02+1M).
The function y~0(x) is chosen as a function that satisfies the boundary conditions and the linear part of the (22). That is, y~0(x) must be a solution of
(25)y~0′′-s2y~0=-1M.
Thus, the initial approximation is chosen to be
(26)y~0(x)=1s2M(1-cosh⁡(sx)cosh⁡(s)).
The initial approximation used in the SHAM algorithms that use the linear operators L1, L2, and L3 is obtained by solving the linear equations:
(27)Li(y0)=-1M,y0(-1)=y0(1)=0,fori=1,2,3.Thus, the initial approximations corresponding to the SHAM withL1, L2, and L3 are given, respectively, by
(28)y0(x)=12M(1-x2),y0(x)=1s2M(1-cosh⁡(sx)cosh⁡(s)),y0(x)=1s12M(1-cosh⁡⁡(s1x)cosh⁡⁡(s1)),
where s1=s2+Fs.The initial approximationz0(x) to use in conjunction with L4 is obtained as a solution of the differential equation formed from the linear part of the (24). That is, we solve
(29)z0′′-s2z0-2Fsy~0z0=-(y~0′′-s2y~0-Fsy~02+1M),z0(-1)=z0(1)=0.
This equation is a linear equation with variable coefficients. Thus, it is solved using the spectral method as described in the earlier section above. Similarly, the initial approximation to use with L5 is obtained as a solution of the differential equation
(30)y0′′-s2y0-Fsy~0y0=-1M,y0(-1)=y0(1)=0,
with y~0 given by (26).With the nonlinear operatorNi defined as the homogeneous part of the governing nonlinear equation in each case, we obtain the following higher order deformation equations:
(31)Li[ym(x)-χmym-1(x)]=ħRi,m-1[y0,y1,…,ym-1],i=1,2,3,5,
where
(32)Ri,m-1=ym-1′′-s2ym-1+1M(1-χm)-Fs∑n=0m-1‍ynym-1-n.The nonlinear operator and higher order deformation equation corresponding toi=4 is defined as
(33)N4(z)=L4(z)-Fsz2,L4[zm(x)-χmzm-1(x)]=ħR4,m-1[z0,z1,…,zm-1],
where
(34)R4,m-1=L4(zm-1)-(1-χm)ϕ1(x)-Fs∑n=0m-1‍znzm-1-n
and ϕ1(x) is the right hand side of (29).


## 4.2.  Example 2: Jeffery-Hamel Equation

Here, we consider the Jeffery-Hamel equation that models the steady two-dimensional flow of an incompressible conducting viscous fluid between two rigid plane walls that meet at an angle2α. The rigid walls are considered to be divergent if α>0 and convergent if α<0. The governing equation is defined as
(35)y′′′+2αReyy′+4α2y′=0,
subject to the boundary conditions
(36)y(0)=1,y′(0)=0,y(1)=0,
where Re is the Reynolds number. This problem was investigated using the standard HAM with the highest order differential matching linear operator L1(y) in [11, 12]. A related method, called optimal homotopy asymptotic method (OHAM), was used to solve the same problem in [13], again with L1(y). In Motsa et al. [6], the SHAM was used with the auxiliary linear operator defined using the transformed linear partition mapping method L4(z). In this example, we explore the auxiliary linear operators defined as
(37)L1[y]=y′′′,L2[y]=y′′′+4α2y′,L3[y]=y′′′+4α2y′+2αRey′,L4[z]=z′′′+4α2z′+2αRey~0z′+2αRey~0′z,L5[y]=y′′′+4α2y′+2αRey~0y′.The functiony~0(x) is chosen as a function that satisfies the boundary conditions and the linear part of (35). Thus, y~0(x) is
(38)y~0(x)=cos⁡(2α)-cos⁡(2αx)cos⁡(2α)-1.The linear operatorL4 is obtained by using the linear partition mapping after using the transformation y(x)=z(x)+y~0(x). The resulting equation is
(39)z′′′+4α2z′+2αRey~0z′+2αRey~0′z+2αRezz′=ϕ2(x)
subject to
(40)z(0)=0,z′(0)=0,z(1)=0,
where
(41)ϕ2(x)=-(y0′′′+2αRey0y0′+4α2y0′).
The initial approximation used in the SHAM algorithms that employ the linear operators L1, L2, and L3 is obtained by solving the linear equations:
(42)Li(y0)=0,y0(0)=1,y0′(0)=0y0(1)=0,fori=1,2,3.
Thus, the initial approximations corresponding to the SHAM with L1, L2, and L3 are given, respectively, by
(43)y0(x)=1-x2,y0(x)=cos⁡(2α)-cos⁡(2αx)cos⁡(2α)-1,y0(x)=cos⁡⁡(2α1)-cos⁡⁡(2α1x)cos⁡⁡(2α1)-1,
where s1=4α2+2αRe. The initial approximation z0(x) to use with L4 is obtained as a solution of the differential equation formed from the linear part of (39). That is, we solve
(44)z0′′′+4α2z0′+2αRey~0z0′+2αRey~0′z0=ϕ2(x).
Since the above equation has variable coefficients, it is solved using the spectral method. Similarly, the initial approximation to use with L5 is obtained as a solution of the differential equation
(45)y0′′′+4α2y0′+2αRey~0y0′=0,y0(0)=1,y0′(0)=0,y0(1)=0
with y~0 given by (38).With the nonlinear operatorNi defined as the homogeneous part of the governing nonlinear equation in each case, we obtain the following higher order deformation equations
(46)Li[ym(x)-χmym-1(x)]=ħRi,m-1[y0,y1,…,ym-1],i=1,2,3,5,
where
(47)Ri,m-1=ym-1′′′+4α2ym-1′+2αRe∑n=0m-1‍ynym-1-n′.The nonlinear operator and higher order deformation equation corresponding toi=4 is defined as
(48)N4(z)=L4(z)+2αRezz′,L4[zm(x)-χmzm-1(x)]=ħR4,m-1[z0,z1,…,zm-1],
where
(49)R4,m-1=L4(zm-1)-(1-χm)ϕ2(x)+2αRe∑n=0m-1‍znzm-1-n′.


## 4.3. Example 3: Laminar Viscous Flow in a Semiporous Channel Subject to a Transverse Magnetic Field

In this section, we consider the problem of laminar viscous flow in a semiporous channel subject to a transverse magnetic field given by(50)y′′′′-Ha2y′′+Re(yy′′′-y′y′′)=0,
subject to the boundary conditions
(51)y=0,y′=0,x=0,y=1,y′=0,x=1.This example was previously investigated using the standard HAM with the highest order differential matching linear operatorL1(y) in [14, 15]. In [16] a blend between the SHAM and a spectral method based quasilinearisation method was used to solve the same problem. In this investigation, we consider the following auxiliary linear operators:
(52)L1[y]=y′′′′,L2[y]=y′′′′-Ha2y′′,L3[y]=y′′′′-Ha2y′′+Re(y′′′-y′′),L4[z]=z′′′′-Ha2z′′+Re(y~0z′′′+y~0′′′z-y~0′z′′-y~0′′z′),L5[y]=y′′′′-Ha2y′′+Re(y~0y′′′-y~0′y′′)=0.
The initial approximation used in the SHAM algorithms that use the linear operators L1, L2, and L3 is obtained by solving the linear equations:
(53)Li(y0)=0,y0(0)=0,y0′(0)=0,y0(1)=1,y0′(1)=0,fori=1,2,3.
Thus, the initial approximations corresponding to the SHAM with L1 and L2 are given by
(54)y0(x)=3x2-2x2,y0(x)=Hax-eHax+eHa-Hax+eHa(Hax-1)+1eHa(Ha-2)+Ha+2.
We note that in this example y~0 is the second solution in (54). The initial approximations corresponding to L3 and L5 are obtained by solving (53) numerically using the spectral method. Similarly, z0(x) is obtained by numerically solving
(55)z0′′′′-Ha2z0′′+Re(y~0z0′′′+y~0′′′z0-y~0′z0′′-y~0′′z0′)=ϕ3(x)
subject to the boundary conditions
(56)z0(0)=0,z0′(0)=0,z0(1)=0,z0′(1)=0,
where
(57)ϕ3(x)=-(y~0′′′′-Ha2y~0′′+Re(y~0y~0′′′-y~0′y~0′′)).
The higher order deformation equations are defined as
(58)Li[ym(x)-χmym-1(x)]=ħRi,m-1[y0,y1,…,ym-1],i=1,2,3,5,
where
(59)Ri,m-1=ym-1′′′′-Ha2ym-1′′+Re∑n=0m-1‍(ynym-1-n′′′-yn′ym-1-n′′).The nonlinear operator and higher order deformation equation corresponding toi=4 is defined as
(60)N4(z)=L4(z)+Re(zz′′′-z′z′′),L4[zm(x)-χmzm-1(x)]=ħR4,m-1[z0,z1,…,zm-1],
where
(61)R4,m-1=L4(zm-1)-(1-χm)ϕ3(x)+Re∑n=0m-1‍(znzm-1-n′′′-zn′zm-1-n′′).


## 4.4. Example 4: Two-Dimensional Viscous Flow in a Rectangular Domain Bounded by Two Moving Porous Walls

In this section we consider the problem of two-dimensional viscous flow in a rectangular domain bounded by two moving porous walls. The governing equations are given in [17, 23] as
(62)y′′′′+α(xy′′′+3y′′)+Re(yy′′′-y′y′′)=0,
subject to the boundary conditions
(63)y=0,y′′=0,atx=0,y=1,y′=0,atx=1,
where α is the nondimensional wall dilation rate defined to be positive for expansion and negative for contraction and Re is the permeation Reynolds number defined positive for injection and negative for suction through the walls. The nonlinear equation (62) was investigated using the standard HAM with the highest order differential matching linear operator L1(y) in [17, 18]. Rashidi et al. [19] used the optimal homotopy asymptotic method (OHAM) to solve the same problem, again using L1(y) as linear operator. The SHAM with transformed linear partition method L4(z) was used in [20].In this investigation, we consider the following auxiliary linear operators:(64)L1[y]=y′′′′,L2[y]=y′′′′+α(xy′′′+3y′′),L3[y]=y′′′′+α(xy′′′+3y′′)+Re(y′′′-y′′),L4[z]=z′′′′+α(xz′′′+3z′′)+Re(y~0z′′′+y~0′′′z-y~0′z′′-y~0′′z′),L5[y]=y′′′′+α(xy′′′+3y′′)+Re(y~0y′′′-y~0′z′′).The initial approximations used in the SHAM algorithms that use the linear operatorsL1, L2, L3, and L5 are obtained by solving the linear equations:
(65)Li(y0)=0,y0(0)=0,y0′(0)=0,y0(1)=1,y0′(1)=0,fori=1,2,3,5.Thus, the initial approximationy0(x) corresponding to L1 is as follows:
(66)y0(x)=12(3x-x3).
The initial approximations corresponding to L2, L3, and L5 are obtained by numerically solving (65). Similarly, z0(x) is obtained by numerically solving
(67)z0′′′′+α(xz0′′′+3z0′′)+Re(y0z0′′′+y0′′′z0-y0′z0′′-y0′′z0′)=ϕ4(x),
where
(68)ϕ4(x)=-[y0′′′′+α(xy0′′′+3y0′′)+Re(y0y0′′′-y0′y0′′)].The higher order deformation equations are defined as(69)Li[ym(x)-χmym-1(x)]=ħRi,m-1[y0,y1,…,ym-1],i=1,2,3,5,
where
(70)Ri,m-1=ym-1′′′′+α(xym-1′′′+3ym-1′′)+Re∑n=0m-1‍(ynym-1-n′′′-yn′ym-1-n′′).The nonlinear operator and higher order deformation equation corresponding toi=4 is defined as
(71)N4(z)=L4(z)+Re(zz′′′-z′z′′),L4[zm(x)-χmzm-1(x)]=ħR4,m-1[z0,z1,…,zm-1],
where
(72)R4,m-1=L4(zm-1)-(1-χm)ϕ4(x)+Re∑n=0m-1‍(znzm-1-n′′′-zn′zm-1-n′′).


## 5. Results and Discussion

The SHAM algorithm using the proposed auxiliary linear operators was applied to the several test problems presented in the last section. In this section we present the numerical results that highlight the difference in the accuracy of the approximate solutions and convergence of the different SHAM algorithms. All the simulations were conducted usingM=20 collocation points. The accuracy and convergence of the SHAM approximate series solution depend on careful selection of the auxiliary parameter ħ. The residual error was used to determine the optimal values of the auxiliary parameter. Using finite terms of the SHAM series we define Kth order approximation at the collocation points xj for j=0,1,2,…,M as
(73)Y=∑k=0K‍yk(xj)=∑k=0K‍Yk,
where Yk=[yk(x0)yk(x1)⋯yk(xM)]t. Given that Y is the SHAM approximate solution at the collocation points, the residual error at the collocation points is defined as
(74)Res(Y)=N[(Y;ħ)],
where N is the nonlinear operator defining the governing nonlinear equation (see (4)). We define the infinity norm of the residual error as
(75)Er(ħ)=∥N[(Y;ħ)]∥∞.
The optimal ħ that gives the best accuracy was identified to be the value of ħ located at the minimum of the graph of the variation of the infinity norm against ħ.Figure1 gives the typical residual error curves that are obtained using the 20th order SHAM approximate solution of the Darcy-Brinkman-Forchheimer equation (Example 1). The minimum of the curve gives the best possible residual error for each auxiliary linear operator used in the SHAM algorithm. A comparison of the minima of the residual error curves gives an indication of which auxiliary linear operator is likely to give better accuracy when used with the optimal auxiliary linear operator in the SHAM algorithm. It can be seen from Figure 1 that the best accuracy, at the 20th order SHAM approximation, is obtained when using the linear partition mapping after transformation (L4). The highest order differential mapping (L1) gives the least accuracy. Figure 2 shows the variation of the norm of the residual error at the optimal value of ħ as a function of the SHAM order. This graph gives an indication of how the SHAM series converges with an increase in the number of iterations (terms used in the series). Rapid convergence is determined by attaining the lowest possible residual error after few iterations. In this example it can be seen from Figure 2 that the highest order differential mapping linear operator L1 is the slowest in terms of convergence, followed by the complete differential mapping L3. The operator that yields the fastest convergence in this example is L4 which gives full convergence after only 10 iterations. It can also be seen from Figure 2 that L2,L3, and L5 all converge to the same level which is slightly lower than that of L1, but slightly higher than L4. This result suggests that the SHAM algorithm that applies L4 gives the best accuracy and when using L1 the accuracy is the lowest, in comparison with the other 4 choices of linear operators in this example.Figure 1
Residualħ-curve for Example 1 using the 20th order SHAM with M=1, s=1, and F=2.Figure 2
Variation of residual error with SHAM order for Example 1 whenM=1, s=1, and F=2.Table1 gives the minimum number of terms of the SHAM series that are required to attain full convergence of the SHAM algorithm using the optimal values of the convergence controlling parameter ħ. The corresponding maximum residual errors are also tabulated. It can be seen from Table 1 that the least error is obtained when using L1 and the other linear operators give more or less comparable residual error values. The SHAM with L4 requires the fewest number terms of the SHAM series to attain maximum accuracy. The SHAM with L5 comes second, followed by L2, L3, and L1 in terms of the minimum SHAM order required to give best accuracy. This observation is in accord with the trends displayed in Figure 2. The table also indicates that the number of terms of the SHAM series required for the SHAM approximation to attain maximum accuracy increases with an increase in the value of the Forchheimer number parameter F. It is worth noting from the governing equation (22) that the parameter F multiplies the nonlinear component of the equation. Table 2 shows the effect of increasing the porous media shape parameter s on the SHAM order required to attain the best accuracy. It can be seen from the table that the effect of s on convergence and accuracy is similar to that of F. In particular, it can be noted that as s increases more SHAM series terms are required to attain the maximum accuracy in all linear operators.Table 1
Example  1 withs=1 and M=1.
F
L
1
L
2
L
3
L
4
L
5
SHAM order
1
35
25
25
17
11
5
81
45
59
18
24
10
130
72
90
38
39
15
174
101
132
33
46
20
254
122
166
37
61
Maximum Residual errorEr
1
1.599
e
-
14
4.108
e
-
15
5.329
e
-
15
4.441
e
-
16
8.993
e
-
15
5
3.375
e
-
14
6.439
e
-
15
3.442
e
-
15
1.443
e
-
15
4.663
e
-
15
10
3.952
e
-
14
9.770
e
-
15
3.442
e
-
15
1.776
e
-
15
5.218
e
-
15
15
5.529
e
-
14
3.397
e
-
14
2.887
e
-
15
3.553
e
-
15
3.331
e
-
15
20
1.612
e
-
13
2.620
e
-
14
4.885
e
-
15
3.664
e
-
15
3.331
e
-
15Table 2
Example  1 withF=1 and M=1.
s
L
1
L
2
L
3
L
4
L
5
SHAM order
1
35
25
25
17
11
5
81
45
59
18
24
10
130
72
90
38
39
15
174
101
132
33
46
20
254
122
166
37
61
Maximum residual errorEr
1
1.599
e
-
14
4.108
e
-
15
5.329
e
-
15
4.441
e
-
16
8.993
e
-
15
5
3.375
e
-
14
6.439
e
-
15
3.442
e
-
15
1.443
e
-
15
4.663
e
-
15
10
3.952
e
-
14
9.770
e
-
15
3.442
e
-
15
1.776
e
-
15
5.218
e
-
15
15
5.529
e
-
14
3.397
e
-
14
2.887
e
-
15
3.553
e
-
15
3.331
e
-
15
20
1.612
e
-
13
2.620
e
-
14
4.885
e
-
15
3.664
e
-
15
3.331
e
-
15Figure3 shows the ħ-curves for the Jeffery-Hamel problem (Example 2) corresponding to the different SHAM schemes generated using all the linear operators using 20 terms of the SHAM approximation series. It can be seen from Figure 3 that, after 20 iterations, L4 would give the best accuracy and L1 gives the least accuracy. This is in accord with the observation made earlier in Example 1. It can also be observed that the behaviour of the SHAM with L1 is approximately the same as that of the SHAM with L2. This can be noted from the observation that the ħ-curve for L1 lies on top of that for L2. In Figure 4 the variation of the residual error is shown as a function of the SHAM order. It can be observed that, in all test linear operators, the residual error decreases with the number of terms used in the SHAM series. This demonstrates that the SHAM approximation converges in all cases considered in this example. The convergence is fastest for L4 followed by L5 and is slowest when L1 is used in the SHAM algorithm. It can also be observed that the results for L1, L2, L3, and L5 converge to approximately the same level while L4 converges to a level lower than that of the other linear operators. Thus, in addition to converging much faster, L4 gives results that are slightly more accurate than those obtained using the other linear operators.Figure 3
Residualħ-curve for Example 2 using the 20th order SHAM with Re=100 and α=-5°.Figure 4
Variation of residual error with SHAM order for Example 2 whenRe=100 and α=-5°.Table3 presents the minimum SHAM order required to give maximum accuracy at the the optimal values of the convergence controlling parameter ħ and the corresponding residual errors. The results are generated for increasing values of the Reynolds number (Re). It can be seen from Table 3 that the minimum order of the SHAM required gradually increases with an increase in Re for both L1 and L2. There is no significant increase in the minimum SHAM order required in the case of L3, L4, and L5. The results from Table 3 also reveal that the required minimum SHAM order of L1 and L2 is comparable particularly for larger values of Re. This confirms the observation made in Figures 3 and 4. Thus, it can be inferred that the behaviour of the SHAM with higher order differential mapping (L1) is comparable with that of SHAM with linear partition mapping (L2) in the solution of the Jeffery-Hamel equation.Table 3
Example  2 withα=-5°.
Re
L
1
L
2
L
3
L
4
L
5
SHAM order
1
10
6
6
4
4
50
26
20
15
10
14
100
34
33
18
19
14
150
43
41
17
28
17
200
53
53
17
17
20
Maximum residual errorEr
1
8.451
e
-
12
8.987
e
-
12
1.074
e
-
11
4.886
e
-
14
1.256
e
-
11
50
8.636
e
-
12
8.887
e
-
12
1.148
e
-
11
1.357
e
-
12
7.468
e
-
12
100
7.930
e
-
12
8.759
e
-
12
1.303
e
-
11
1.886
e
-
12
1.633
e
-
11
150
7.815
e
-
12
8.092
e
-
12
2.546
e
-
11
2.973
e
-
12
1.808
e
-
11
200
9.128
e
-
12
8.421
e
-
12
3.209
e
-
11
4.435
e
-
12
1.176
e
-
11Figure5 shows the ħ-curves corresponding to the example of laminar viscous flow in a semiporous channel subject to a transverse magnetic field (Example 3). Again, it can be seen from Figure 5 that the residual error is lowest when the SHAM is used with L4 and largest when used with L1 which, as was observed in the case of Jeffery-Hamel flow, is comparable with L2. The same observation can be made in Figure 6 where the residual error is plotted against the SHAM order. From Figure 6, it can also be seen that L4 gives the best convergence followed by L5 and the slowest convergence is obtained when using L1 which is comparable to L2.Figure 5
Residualħ-curve for Example 3 using the 20th order SHAM with Re=10 and Ha=1.Figure 6
Variation of residual error with SHAM order for Example 3 whenRe=10 and Ha=1.Table4 shows effect of increasing the Reynolds number (Re) on the minimum SHAM orders and the corresponding residual errors. It can be observed, from Table 4, that the minimum SHAM orders increase with an increase in Re when the SHAM is used with L1, L2, and L3. It can also be seen that there is no significant increase in the minimum SHAM orders when using L4 and L5. Furthermore, the residual error obtained using L4 is lower than that of the other linear operators. This is in line with similar observations made earlier.Table 4
Example  3 with Ha=1.
Re
L
1
L
2
L
3
L
4
L
5
SHAM order
1
11
14
10
5
7
5
26
24
21
11
12
10
49
46
30
12
13
15
81
76
38
13
16
20
118
112
44
15
18
25
160
145
49
16
20
Maximum residual errorEr
1
1.744
e
-
09
9.366
e
-
09
2.376
e
-
09
2.130
e
-
12
2.081
e
-
09
5
1.329
e
-
09
5.971
e
-
09
1.947
e
-
09
1.197
e
-
11
3.024
e
-
09
10
1.285
e
-
09
5.408
e
-
09
1.746
e
-
09
1.636
e
-
11
5.125
e
-
09
15
1.880
e
-
09
6.254
e
-
09
1.977
e
-
09
6.425
e
-
11
1.686
e
-
09
20
1.666
e
-
09
5.110
e
-
09
1.150
e
-
09
3.407
e
-
11
4.478
e
-
09
25
1.857
e
-
09
6.505
e
-
09
2.491
e
-
09
5.090
e
-
11
3.529
e
-
09Table5 shows the effect of increasing the Hartmann number (Ha) on the minimum SHAM orders and the corresponding residual errors in Example 3. In this case, it can be observed that the minimum SHAM orders increase with an increase in Ha only in the case of L1. In the case of the other linear operators the minimum SHAM orders do not increase with an increase in Ha. This means that convergence of the SHAM series is independent of the value of Ha when the SHAM is used with L2, L3, L4, and L5. It can also be noted from Table 5 that the residual error is smallest when L4 is used. If can be inferred that the SHAM is most accurate when used with L4.Table 5
Example  3 with Re=1.
Ha
L
1
L
2
L
3
L
4
L
5
SHAM order
2
10
11
10
5
7
4
16
9
9
5
6
6
26
10
11
5
5
8
32
10
8
6
11
10
45
18
8
5
5
15
102
14
8
5
5
Maximum Residual errorEr
2
1.484
e
-
08
1.323
e
-
08
9.815
e
-
09
3.391
e
-
12
2.954
e
-
09
4
8.618
e
-
09
3.334
e
-
09
3.888
e
-
09
5.698
e
-
12
2.855
e
-
09
6
2.258
e
-
09
6.625
e
-
09
1.281
e
-
08
1.461
e
-
11
6.068
e
-
09
8
1.172
e
-
09
2.022
e
-
09
1.012
e
-
08
1.690
e
-
11
7.852
e
-
09
10
8.887
e
-
09
5.535
e
-
09
1.261
e
-
08
5.424
e
-
11
2.253
e
-
09
15
3.325
e
-
09
1.870
e
-
09
5.521
e
-
09
2.964
e
-
10
5.871
e
-
09Figure7 shows the ħ-curves corresponding to the example of two-dimensional viscous flow in a rectangular domain bounded by two moving porous walls (Example 4). It can be noted from Figure 7 that the residual error is the lowest when the SHAM is used with L4 and the largest when used with L1 after 20 iterations. Figure 8 illustrates the variation of the residual error against the SHAM order. It can be seen from Figure 8 that L4 gives the best convergence followed by L5 and the slowest convergence is obtained when using L1. The linear operators L1, L2, L3, and L5 converge to the same level of residual error which is higher than that of L4. Again, it can be inferred that L4 gives superior accuracy compared to the other linear operators.Figure 7
Residualħ-curve for Example 4 using the 20th order SHAM with Re=10 and α=1.Figure 8
Variation of residual error with SHAM order for Example 4 whenRe=10 and α=1.Table6 shows the effect of increasing the Reynolds number (Re) on the minimum SHAM orders and corresponding residual errors in Example 4. It can be observed that the minimum SHAM orders increase with an increase in Re when the SHAM is used with L1, L2, and L3. In the case of L4 and L5 it seems possible that an increase in Re has no strong influence in the convergence speed of the SHAM. It can also be noted from Table 6 that the residual error is the smallest when L4 is used. This is in accord with the trends observed in all the other examples considered in this work. Table 7 shows the effect of the nondimensional wall dilation rate α on the minimum SHAM order and corresponding residual errors. It can be observed from Table 7 that an increase in α results in a corresponding significant increase in the minimum SHAM order for L1. The minimum SHAM orders are not significantly influenced by α when the SHAM is used with L2, L3, L4, and L5. Furthermore, the best accuracy is obtained when using L4.Table 6
Example  4 withα=1.
Re
L
1
L
2
L
3
L
4
L
5
SHAM order
1
16
8
12
5
5
5
29
16
16
5
6
10
40
21
24
7
8
20
57
37
35
6
8
50
117
76
65
7
10
Maximum Residual errorEr
1
2.103
e
-
09
4.065
e
-
09
1.354
e
-
08
4.946
e
-
13
6.441
e
-
09
5
2.107
e
-
09
3.635
e
-
09
1.016
e
-
08
1.837
e
-
12
9.894
e
-
09
10
2.219
e
-
09
3.233
e
-
09
5.365
e
-
09
1.990
e
-
12
2.995
e
-
09
20
2.408
e
-
09
4.137
e
-
09
5.088
e
-
09
2.459
e
-
12
1.235
e
-
08
50
2.160
e
-
09
3.528
e
-
09
1.098
e
-
08
8.652
e
-
12
5.574
e
-
09Table 7
Example  4 with Re=1.
α
L
1
L
2
L
3
L
4
L
5
SHAM order
1
16
8
12
5
5
2
21
10
13
6
7
4
32
11
8
8
6
6
52
21
10
9
8
8
99
22
11
14
11
Maximum residual errorEr
1
2.103
e
-
09
4.065
e
-
09
1.354
e
-
08
4.946
e
-
13
6.441
e
-
09
2
1.198
e
-
09
1.109
e
-
08
7.087
e
-
10
7.009
e
-
12
9.531
e
-
09
4
1.954
e
-
09
1.017
e
-
08
1.604
e
-
08
3.667
e
-
11
3.063
e
-
09
6
4.419
e
-
09
5.935
e
-
09
1.772
e
-
08
7.262
e
-
11
1.388
e
-
08
8
1.930
e
-
09
2.864
e
-
09
1.449
e
-
08
9.716
e
-
11
1.221
e
-
09


## 6. Conclusion

This study gives a systematic way of choosing initial approximations and linear operators that can be employed in the spectral homotopy analysis method (SHAM) algorithm. The main aim of the study was to determine the linear operator definition that gives approximate solutions with optimal convergence and accuracy. A comparison between the accuracy and convergence of SHAM algorithms using commonly used linear operators and newly proposed linear operators was conducted. Simulations were conducted on selected nonlinear boundary value problems defined on bounded domains. One of the more significant findings to emerge from this study is that the linear operator defined by the highest order differential matching, which is the most commonly used linear operator for solving boundary value problems defined on bounded domains, results in the slowest converging and least accurate SHAM algorithm. The SHAM algorithm that uses a linear operator defined by the method of linear partition matching was found to perform slightly better than the SHAM with the highest order differential matching in some cases. However, in other cases the methods were found to be comparable. It was also shown that the SHAM with linear operator defined by the method of complete differential matching performed much better than the SHAM with linear partition mapping or the highest order differential matching. The study proposed a modification of the complete differential method approach which, unlike the original approach, does not ignore the nonlinear factors of the derivative terms in the definition of the linear operator. This new definition of linear operator was found to perform better than the original complete differential method approach. The most optimal linear operator was found to be one that uses the linear partition mapping after a transformation that expresses the unknown function as a sum of a new unknown plus some initial guess that satisfies the underlying boundary conditions.The current findings add substantially to our understanding of how homotopy analysis based methods can be modified to yield approximate solutions with optimum accuracy and convergence. Further work needs to be done to investigate the application of the linear operators proposed in this study in solving nonlinear partial differential equations. It would also be very interesting to extend the current study to systems of nonlinear differential equations.



---

*Source: 697845-2014-08-14.xml*