# Wavelet-Based Homotopy Analysis Method for Nonlinear Matrix System and Its Application in Burgers Equation

Authors: Xing Ruyi
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

To generalize the homotopy analysis method (HAM) to multidegree-of-freedom nonlinear system, the adaptive precise integration method (APIM) is introduced into the HAM, with which the almost exact value of the exponential matrix can be obtained. Combining the interval interpolation wavelet collocation method, HAM-based APIM can be employed to solve the nonlinear PDEs. As an example, Burgers equation is spatially discretized by the interval quasi-Shannon wavelet collocation method and solved by the proposed method to illustrate the effectiveness and great potential of the homotopy analysis method in nonlinear problems.


---

## Body

## 1. Introduction

The homotopy analysis method (HAM) is developed in 1992 by Liao in [1–5]. This method has been successfully applied in various nonlinear problems in science and engineering in [6–11] and references therein. Based on homotopy of topology, the validity of the HAM is independent of whether or not there exist small parameters in the considered equation. Therefore, the HAM can overcome the foregoing restrictions and limitations of perturbation techniques [12, 13]. In fact, He [14, 15] proposed another asymptotic numerical method for solving nonlinear problems named homotopy perturbation method (HPM). The better improvement is adding an auxiliary parameter into the homotopy equation, which is helpful to eliminate the secular term in the perturbation solution. This can improve the rate of convergence greatly. Unlike analytical perturbation methods, HPM does not depend on small parameter which is difficult to find. Most researchers think of HPM as a special case of HAM, which is much simpler than classical HAM and has been widely used in solving various nonlinear problems [16]. However, the HAM applies mainly in single-degree-of-freedom systems; as for multidegree-of-freedom nonlinear differential equations, due to the correlation of the equations, the solving process is very tedious, and strong calculating technique is needed.In an attempt to overcome the previous difficulties, the precise integration method (PIM) proposed by Zhong [17] in 1994 is employed, because it is effective in computation precision and stability for the solutions of linear differential equations [18, 19]. However, like other numerical methods, PIM has its own limitation. When it is applied to the problems of time-variant or nonlinear time integration, PIM usually expands the nonlinear term in Taylor series firstly, which is then solved by using an iteration method. In order to obtain high precision, a very small time step should be taken for PIM, which results in huge computational complexity and greatly restricts the application of PIM in engineering.Burgers’ model of turbulence is a crucial fluid dynamic model, and the study of this model and theory of shock waves has been considered by many authors both for conceptual understanding of a class of physical flows and for testing various numerical methods. In this work, the Burgers equation is firstly transformed into a set of nonlinear matrix differential equations using wavelet collocation method [20–22]. Since the definition domain of wavelet transformation is an infinite interval, which can deduce the boundary effect, this can be eliminated by introducing an interval wavelet [23] instead of the common wavelet function.Taking full advantage of the HAM and PIM, a novel approximate method for multidegree-of-freedom nonlinear matrix differential equation based on interval wavelet is presented in this paper.


## 2. The HAM Based on the PIM

### 2.1. Basic Idea of the Homotopy Analysis Method

Taking as an example, we consider the following nonlinear algebraic equation:(1)f(x)=0.
First of all, we construct such a homotopy
(2)H(x;q)=(1-q)[f(x)-f(x0)]+q·f(x),
where x0 is an initial guess of x and q∈[0,1] is called embedding parameter. Obviously, at q=0 and q=1, one has
(3)H(x;0)=f(x)-f(x0),H(x;1)=f(x),
respectively. Thus, as q increases from 0 to 1, H(x;q) varies continuously from f(x)-f(x0) to f(x). Such kind of continuous variation is called deformation in topology [20]. Now, enforcing H(x;q)=0, that is,
(4)(1-q)[f(x)-f(x0)]+q·f(x)=0,
we have now a family of algebraic equations. Obviously, the solution of the previous family of algebraic equations is dependent upon the embedding parameter q. So, the family of equations can be rewritten as
(5)(1-q){f[ϕ(x;q)]-f(x0)}+q·f[ϕ(x;q)]=0.As the embedding parameterq increases from 0 to 1, ϕ(x;q) varies from the initial guess x0 to the solution x of f(x)=0. We call the family of equations like (5) the zeroth-order deformation equationBecauseϕ(x;q) is now a function of the embedding parameter q, we can expand it into Maclaurin series
(6)ϕ(x;q)=x0+∑m=1+∞xm·qm,
where ϕ(x;0)=x0 is employed, and
(7)xm=1m!·∂mϕ(x;q)∂qm|q=0.
Here, the series (6) is called homotopy series, and xm is called the kth-order homotopy derivative of ϕ(x;q). If the convergence radius of the previous Maclaurin series is no less than 1, using the relationship ϕ(x;1)=x, one has
(8)x=x0+∑m=1+∞xm.According to the fundamental theorem of calculus about Taylor series, the coefficientxm of the homotopy series (6) is unique.Substituting (6) into (5) and rearranging based on powers of q-terms, one has
(9)x1f′(x0)+f(x0)=0,x2f′(x0)+12x12f′′(x0)=0.
Whose solutions are
(10)x1=-f(x0)f′(x0),(11)x2=-x12f′′(x0)2f′(x0)=-f2(x0)f′′(x0)2[f′(x0)]3.In this way, one obtainsxm one by one in the order m=1,2,3,…. Then one has 1st-order homotopy-series approximation
(12)x≈x0+x1=x0-f(x0)f′(x0),
2nd-order homotopy-series approximation
(13)x≈x0+x1+x2=x0-f(x0)f′(x0)-f2(x0)f′′(x0)2[f′(x0)]3.


### 2.2. The Precise Integration Method for Nonlinear Matrix Differential Equations

Taking as an example, we consider the following matrix differential equations:(14)dvdt=Hv+f(v,t),
where v is an n-dimensional unknown vector, H and f(v,t) are n×n matrix, and t is a time variable. The general solution of (14) is
(15)v(t)=eH(t-tk)v(tk)+eHt∫tkte-Ht·f(v,t)dt.
To solve the integral part of (15), the linear interpolation approximation is employed
(16)f(v,t)=r0+r1(t-tk),
where r0=f(v(tk),tk), r1=r1=(∂f/∂t)|t=tk+∑i=12n[(∂f/∂vi)(Hv(tk)+f(v(tk),tk))i], and vi is the ith element of vector v. Based on the previous derivations, the recurrence formula of (14) can be obtained as
(17)v(tk+1)=eH(tk+1-tk)[v(tk)+H-1(r0+H-1r1)]-H-1[r0+H-1r1+(t-tk)r1].Supposing the time stepτ=tk+1-tk is constant, the matrix exponential function can be calculated accurately by using PIM. The calculation of the exponential function is based on the identity
(18)eH·τ=[eH·(τ/m)]m,
where m is an integer. It is suggested to select m=2N, such as N=20, m=1048576.Becauseτ is a small time interval, Δt=τ/m is an extremely small time interval. Hence for the interval, the truncated Taylor expansion is applied with high precision
(19)T=eH·Δt=I+Ta,
where Ta=H·Δt+(H·Δt)2/2+(H·Δt)3/3!+(H·Δt)4/4!BecauseΔt is very small, the first five-term series expansion should be enough.For computing the matrixT, (17) could be factored as
(20)T=[I+Ta]2N=[I+2·Ta+Ta·Ta]2N-1.
Such factorization should be iterated N times. Noted that it is equivalent to do the following instruction:
(21)for(iter=0;iter!=N;++iter)Ta=2·Ta+Ta·Ta,
after the execution of the above instruction, T=I+Ta.Because the nonlinear termf(v,t) is expanded in Taylor series and the PIM uses an iterative method, the time step τ must be very small.


## 2.1. Basic Idea of the Homotopy Analysis Method

Taking as an example, we consider the following nonlinear algebraic equation:(1)f(x)=0.
First of all, we construct such a homotopy
(2)H(x;q)=(1-q)[f(x)-f(x0)]+q·f(x),
where x0 is an initial guess of x and q∈[0,1] is called embedding parameter. Obviously, at q=0 and q=1, one has
(3)H(x;0)=f(x)-f(x0),H(x;1)=f(x),
respectively. Thus, as q increases from 0 to 1, H(x;q) varies continuously from f(x)-f(x0) to f(x). Such kind of continuous variation is called deformation in topology [20]. Now, enforcing H(x;q)=0, that is,
(4)(1-q)[f(x)-f(x0)]+q·f(x)=0,
we have now a family of algebraic equations. Obviously, the solution of the previous family of algebraic equations is dependent upon the embedding parameter q. So, the family of equations can be rewritten as
(5)(1-q){f[ϕ(x;q)]-f(x0)}+q·f[ϕ(x;q)]=0.As the embedding parameterq increases from 0 to 1, ϕ(x;q) varies from the initial guess x0 to the solution x of f(x)=0. We call the family of equations like (5) the zeroth-order deformation equationBecauseϕ(x;q) is now a function of the embedding parameter q, we can expand it into Maclaurin series
(6)ϕ(x;q)=x0+∑m=1+∞xm·qm,
where ϕ(x;0)=x0 is employed, and
(7)xm=1m!·∂mϕ(x;q)∂qm|q=0.
Here, the series (6) is called homotopy series, and xm is called the kth-order homotopy derivative of ϕ(x;q). If the convergence radius of the previous Maclaurin series is no less than 1, using the relationship ϕ(x;1)=x, one has
(8)x=x0+∑m=1+∞xm.According to the fundamental theorem of calculus about Taylor series, the coefficientxm of the homotopy series (6) is unique.Substituting (6) into (5) and rearranging based on powers of q-terms, one has
(9)x1f′(x0)+f(x0)=0,x2f′(x0)+12x12f′′(x0)=0.
Whose solutions are
(10)x1=-f(x0)f′(x0),(11)x2=-x12f′′(x0)2f′(x0)=-f2(x0)f′′(x0)2[f′(x0)]3.In this way, one obtainsxm one by one in the order m=1,2,3,…. Then one has 1st-order homotopy-series approximation
(12)x≈x0+x1=x0-f(x0)f′(x0),
2nd-order homotopy-series approximation
(13)x≈x0+x1+x2=x0-f(x0)f′(x0)-f2(x0)f′′(x0)2[f′(x0)]3.


## 2.2. The Precise Integration Method for Nonlinear Matrix Differential Equations

Taking as an example, we consider the following matrix differential equations:(14)dvdt=Hv+f(v,t),
where v is an n-dimensional unknown vector, H and f(v,t) are n×n matrix, and t is a time variable. The general solution of (14) is
(15)v(t)=eH(t-tk)v(tk)+eHt∫tkte-Ht·f(v,t)dt.
To solve the integral part of (15), the linear interpolation approximation is employed
(16)f(v,t)=r0+r1(t-tk),
where r0=f(v(tk),tk), r1=r1=(∂f/∂t)|t=tk+∑i=12n[(∂f/∂vi)(Hv(tk)+f(v(tk),tk))i], and vi is the ith element of vector v. Based on the previous derivations, the recurrence formula of (14) can be obtained as
(17)v(tk+1)=eH(tk+1-tk)[v(tk)+H-1(r0+H-1r1)]-H-1[r0+H-1r1+(t-tk)r1].Supposing the time stepτ=tk+1-tk is constant, the matrix exponential function can be calculated accurately by using PIM. The calculation of the exponential function is based on the identity
(18)eH·τ=[eH·(τ/m)]m,
where m is an integer. It is suggested to select m=2N, such as N=20, m=1048576.Becauseτ is a small time interval, Δt=τ/m is an extremely small time interval. Hence for the interval, the truncated Taylor expansion is applied with high precision
(19)T=eH·Δt=I+Ta,
where Ta=H·Δt+(H·Δt)2/2+(H·Δt)3/3!+(H·Δt)4/4!BecauseΔt is very small, the first five-term series expansion should be enough.For computing the matrixT, (17) could be factored as
(20)T=[I+Ta]2N=[I+2·Ta+Ta·Ta]2N-1.
Such factorization should be iterated N times. Noted that it is equivalent to do the following instruction:
(21)for(iter=0;iter!=N;++iter)Ta=2·Ta+Ta·Ta,
after the execution of the above instruction, T=I+Ta.Because the nonlinear termf(v,t) is expanded in Taylor series and the PIM uses an iterative method, the time step τ must be very small.


## 3. Coupling Technique of the HAM and the PIM for Burgers Equation Based on Interval Wavelet

### 3.1. Spatial Discretization of the Burgers Equation

Consider that the one-dimensional Burgers equation has the form [21–23]
(22)∂u∂t+u∂u∂x=1Re∂2u∂2x,x∈[0,2],t≥0,
subject to the initial and boundary conditions
(23)u(x,0)=sin(πx),u(0,t)=u(2,t)=0,
where t represents the time, while Re denotes the Reynolds number.According to the basic idea of wavelet collocation method, the solutionu(x,t) could be approximately represented as
(24)u(x,t)≈uj(x,t)=∑n=02ju(xn,t)φj(x-xn),
where xn(n=0,1,2,…,2j) are uniform discrete points and φ(x) is the scaling function. Here we take the quasi-Shannon wavelet [24, 25] as the scaling function
(25)φ(x)=sin(πx)πxexp(-x22σ2),σ>0,
where σ is the window size parameter.In order to take the quasi-Shannon scaling function as the basis function,u(x,t) is discretized in its definition domain
(26)xn=21-j·n,j∈Z,n=0,1,2,…,2j.
And thereby the basis function is obtained
(27)φj(x-xn)=sin[2j-1π(x-xn)]2j-1π(x-xn)exp[-22j-3(x-xn)2r2],
where r=2j-1·σ.It is easy to see thatφj(x-xn) has the interpolation property
(28)φj(xk-xn)=δk,n={1k=n0k≠n.Using its interpolation property, the first derivative and second derivative ofφj(x-xn) can be obtained as follows:
(29)φj′(x-xn)={0k=n2j-1cos[π(k-n)]exp[-(k-n)2/2r2](k-n)k≠n,φj′′(x-xn)={-(3+π2r2)·22j-23r2k=n-22j-2·[1(k-n)2+1r2]·cos[π(k-n)]·exp[-(k-n)22r2]k≠n.Substituting (4) into (1), a series of nonlinear ODEs can be obtained
(30)du(xk,t)dt=∑n=02ju(xn,t)[1Reφj′′(xk-xn)-u(xk,t)φj′(xk-xn)].For simplicity, we can write (10) as matrix differential equation
(31)dV(t)dt=M0·V(t)+M1(V(t))·V(t),
where(32)V(t)=[u(x0,t),u(x1,t),…,u(x2j,t)]T,M0=1Re(φj′′(x0-x0)φj′′(x0-x1)⋯φj′′(x0-x2j)φj′′(x1-x0)φj′′(x1-x1)⋯φj′′(x1-x2j)⋮⋮⋯⋮φj′′(x2j-x0)φj′′(x2j-x1)⋯φj′′(x2j-x2j)),M1=-diag(V(t))·(φj′(x0-x0)φj′(x0-x1)⋯φj′(x0-x2j)φj′(x1-x0)φj′(x1-x1)⋯φj′(x1-x2j)⋮⋮⋮φj′(x2j-x0)φj′(x2j-x1)⋯φj′(x2j-x2j)).


### 3.2. Construction of Quasi-Shannon Interval Wavelet

The Lagrange polynomialslj,k1, lj,k2 are introduced
(33)lj,k1=∏i=0i≠kLx-xj,ixj,k-xj,i,lj,k2=∏i=2j-Li≠k2jx-xj,ixj,k-xj,i.
The following two weight coefficients are defined
(34)ank=lj,k1(xjn),bnk=lj,k2(xjn).
Then the interpolation basic function can be expressed as
(35)wj,k(x)=wj,k(x)+∑n=-N+1-1ank·wj,n(x),k=0,…,L,(36)wj,k(x)=wj,k(x),k=L+1,…,2j-L-1,(37)wj,k(x)=wj,k(x)+∑n=2j+12j+N-1bnk·wj,n(x),k=2j-L,…,2j,
where L is the number of he external points and N is the support domain of the wavelet function.It is easy to know that the quasi-Shannon interval wavelet is a linear combination of the quasi-Shannon scaling functionwj,k(x). Therefore, the quasi-Shannon interval wavelet function possesses all the properties of the quasi-Shannon scaling function.


### 3.3. The PIM for Matrix Differential Equations

We choose the time step asτ, assuming that the nonlinear term in (31) is linear in the relatively short time segment (tk,tk+1), the nonlinear term M1(V(t))·V(t) of (31) can be expressed in Taylor series, and the first two terms are taken. Therefore, (31) could be represented as
(38)dV(t)dt=M0V(t)+r0+r1(t-tk),
where
(39)r0=M1(tk)V(tk),r1=[M2(tk)+M1(tk)M0+M1(tk)M1(tk)]V(tk),M2(t)=dM1(t)dt.Based on the previous derivation, the iterative solution of (38) can be obtained
(40)V(tk+1)=eM0·τ[V(tk)+M0-1(r0+M0-1r1)]-M0-1(r0+M0-1r+r1τ),
where the matrix exponential function eM0·τ can be calculated accurately by using PIM.


### 3.4. The Homotopy Analysis Method Based on the Precise Integration Method

There are various ways to construct a homotopy function. For (11), a linear homotopy can be constructed as
(41)q[dϕ(x,t;q)dt-M0-M1(ϕ(x,t;q))ϕ(x,t;q)]+(1-q)×[dϕ(x,t;q)dt-M0ϕ(x,t;q)-M1(A)ϕ(x,t;q)]=0,
where A=V(0) is the initial approximation of V(t) and q∈[0,1] is the embedding parameter. As q increases from 0 to 1, the solution ϕ(x,t;q) varies from the initial guess A to the solution V(t), and it can be expanded in Taylor series about the embedding parameter
(42)ϕ(x,t;q)=A+∑m=1+∞Vm(t)·qm,
where
(43)Vm(t)=1m!∂mϕ(x,t;q)∂qmthat the approximate solutions are generally not valid for larget. In order to improve the computation accuracy, we divide the time interval [0-t] into a number of sections uniformly by setting the time step as τ, and use the iterative method. In such a way, a series of time nodes can be expressed as
(44)t0=0,t1=τ,…,tk=k·τ,….Since the only requirement forq is that from q=0 to q=1, ϕ(x,t;q) should change continuously from the initial guess A to V(t). Therefore, it is feasible to improve the computation accuracy by assuming that
(45)q=Δtτ,Δt∈[0,τ].
Then ϕ(x,t;q) can be expanded in the form
(46)ϕ(x,t;q)=V0(t)+V1(t)·(Δtτ)+V2(t)·(Δtτ)2+V3(t)·(Δtτ)3+⋯.Substituting (46) into (41) and rearranging based on powers of q-terms, the following are obtained:
(47)V0(t)=A,dV0(t)dt=H·V0(t)=H·A,V1(t)=τ·dV0(t)dt=τ·H·A,dV1(t)dt=H·V1(t)0=τ·H2·A,V2(t)=τ2·dV1(t)dt=(τH)22·A,dV2(t)dt=H·V2(t)+s0=τ2H32·A+s0,V3(t)=τ3·dV2(t)dt=(τH)33!·A+τ3·s0,dV3(t)dt=H·V3(t)+s1=τ3H43!·A+τH3·s0+s1,V4(t)=τ4·dV4(t)dt=(τH)44!·A+τ2H3·4·s0+τ4·s1,
where s0=M1(V1(t))·A, s1=M1(V1(t))·V1(t)+M1(V2(t))·V0(t) and H=M0+M1(A).Substituting (47) into (46) and assuming q=1, one has
(48)V(t)=V0(t)+V1(t)+V2(t)+V3(t)+⋯=[I+τH+(τH)22!+⋯]A+[τ3+τ2H3·4+τ3H23·4·5+⋯]s0+[τ4+τ2H4·5+τ3H24·5·6+⋯]s1+⋯=T·A+2·(H-1)3τ2[(τH)33!+(τH)44!+⋯]s0+3!·(H-1)4τ3[(τH)44!+(τH)55!+⋯]s1+⋯=T·A+2·(H-1)3τ2·R·s0+3!·(H-1)4τ3·Q·s1+⋯,
where R=[T-I-τH-(τH)2/2], Q=R-(τH)3/6, T=eτH. The matrix exponential function T can be calculated by using PIM.For simplicity, only the first two terms are taken as the approximate solution. Accordingly, the approximate solution is obtained(49)V(t)≈V0(t)+V1(t)=T·A+2·(H-1)3τ2·R·s0.


## 3.1. Spatial Discretization of the Burgers Equation

Consider that the one-dimensional Burgers equation has the form [21–23]
(22)∂u∂t+u∂u∂x=1Re∂2u∂2x,x∈[0,2],t≥0,
subject to the initial and boundary conditions
(23)u(x,0)=sin(πx),u(0,t)=u(2,t)=0,
where t represents the time, while Re denotes the Reynolds number.According to the basic idea of wavelet collocation method, the solutionu(x,t) could be approximately represented as
(24)u(x,t)≈uj(x,t)=∑n=02ju(xn,t)φj(x-xn),
where xn(n=0,1,2,…,2j) are uniform discrete points and φ(x) is the scaling function. Here we take the quasi-Shannon wavelet [24, 25] as the scaling function
(25)φ(x)=sin(πx)πxexp(-x22σ2),σ>0,
where σ is the window size parameter.In order to take the quasi-Shannon scaling function as the basis function,u(x,t) is discretized in its definition domain
(26)xn=21-j·n,j∈Z,n=0,1,2,…,2j.
And thereby the basis function is obtained
(27)φj(x-xn)=sin[2j-1π(x-xn)]2j-1π(x-xn)exp[-22j-3(x-xn)2r2],
where r=2j-1·σ.It is easy to see thatφj(x-xn) has the interpolation property
(28)φj(xk-xn)=δk,n={1k=n0k≠n.Using its interpolation property, the first derivative and second derivative ofφj(x-xn) can be obtained as follows:
(29)φj′(x-xn)={0k=n2j-1cos[π(k-n)]exp[-(k-n)2/2r2](k-n)k≠n,φj′′(x-xn)={-(3+π2r2)·22j-23r2k=n-22j-2·[1(k-n)2+1r2]·cos[π(k-n)]·exp[-(k-n)22r2]k≠n.Substituting (4) into (1), a series of nonlinear ODEs can be obtained
(30)du(xk,t)dt=∑n=02ju(xn,t)[1Reφj′′(xk-xn)-u(xk,t)φj′(xk-xn)].For simplicity, we can write (10) as matrix differential equation
(31)dV(t)dt=M0·V(t)+M1(V(t))·V(t),
where(32)V(t)=[u(x0,t),u(x1,t),…,u(x2j,t)]T,M0=1Re(φj′′(x0-x0)φj′′(x0-x1)⋯φj′′(x0-x2j)φj′′(x1-x0)φj′′(x1-x1)⋯φj′′(x1-x2j)⋮⋮⋯⋮φj′′(x2j-x0)φj′′(x2j-x1)⋯φj′′(x2j-x2j)),M1=-diag(V(t))·(φj′(x0-x0)φj′(x0-x1)⋯φj′(x0-x2j)φj′(x1-x0)φj′(x1-x1)⋯φj′(x1-x2j)⋮⋮⋮φj′(x2j-x0)φj′(x2j-x1)⋯φj′(x2j-x2j)).


## 3.2. Construction of Quasi-Shannon Interval Wavelet

The Lagrange polynomialslj,k1, lj,k2 are introduced
(33)lj,k1=∏i=0i≠kLx-xj,ixj,k-xj,i,lj,k2=∏i=2j-Li≠k2jx-xj,ixj,k-xj,i.
The following two weight coefficients are defined
(34)ank=lj,k1(xjn),bnk=lj,k2(xjn).
Then the interpolation basic function can be expressed as
(35)wj,k(x)=wj,k(x)+∑n=-N+1-1ank·wj,n(x),k=0,…,L,(36)wj,k(x)=wj,k(x),k=L+1,…,2j-L-1,(37)wj,k(x)=wj,k(x)+∑n=2j+12j+N-1bnk·wj,n(x),k=2j-L,…,2j,
where L is the number of he external points and N is the support domain of the wavelet function.It is easy to know that the quasi-Shannon interval wavelet is a linear combination of the quasi-Shannon scaling functionwj,k(x). Therefore, the quasi-Shannon interval wavelet function possesses all the properties of the quasi-Shannon scaling function.


## 3.3. The PIM for Matrix Differential Equations

We choose the time step asτ, assuming that the nonlinear term in (31) is linear in the relatively short time segment (tk,tk+1), the nonlinear term M1(V(t))·V(t) of (31) can be expressed in Taylor series, and the first two terms are taken. Therefore, (31) could be represented as
(38)dV(t)dt=M0V(t)+r0+r1(t-tk),
where
(39)r0=M1(tk)V(tk),r1=[M2(tk)+M1(tk)M0+M1(tk)M1(tk)]V(tk),M2(t)=dM1(t)dt.Based on the previous derivation, the iterative solution of (38) can be obtained
(40)V(tk+1)=eM0·τ[V(tk)+M0-1(r0+M0-1r1)]-M0-1(r0+M0-1r+r1τ),
where the matrix exponential function eM0·τ can be calculated accurately by using PIM.


## 3.4. The Homotopy Analysis Method Based on the Precise Integration Method

There are various ways to construct a homotopy function. For (11), a linear homotopy can be constructed as
(41)q[dϕ(x,t;q)dt-M0-M1(ϕ(x,t;q))ϕ(x,t;q)]+(1-q)×[dϕ(x,t;q)dt-M0ϕ(x,t;q)-M1(A)ϕ(x,t;q)]=0,
where A=V(0) is the initial approximation of V(t) and q∈[0,1] is the embedding parameter. As q increases from 0 to 1, the solution ϕ(x,t;q) varies from the initial guess A to the solution V(t), and it can be expanded in Taylor series about the embedding parameter
(42)ϕ(x,t;q)=A+∑m=1+∞Vm(t)·qm,
where
(43)Vm(t)=1m!∂mϕ(x,t;q)∂qmthat the approximate solutions are generally not valid for larget. In order to improve the computation accuracy, we divide the time interval [0-t] into a number of sections uniformly by setting the time step as τ, and use the iterative method. In such a way, a series of time nodes can be expressed as
(44)t0=0,t1=τ,…,tk=k·τ,….Since the only requirement forq is that from q=0 to q=1, ϕ(x,t;q) should change continuously from the initial guess A to V(t). Therefore, it is feasible to improve the computation accuracy by assuming that
(45)q=Δtτ,Δt∈[0,τ].
Then ϕ(x,t;q) can be expanded in the form
(46)ϕ(x,t;q)=V0(t)+V1(t)·(Δtτ)+V2(t)·(Δtτ)2+V3(t)·(Δtτ)3+⋯.Substituting (46) into (41) and rearranging based on powers of q-terms, the following are obtained:
(47)V0(t)=A,dV0(t)dt=H·V0(t)=H·A,V1(t)=τ·dV0(t)dt=τ·H·A,dV1(t)dt=H·V1(t)0=τ·H2·A,V2(t)=τ2·dV1(t)dt=(τH)22·A,dV2(t)dt=H·V2(t)+s0=τ2H32·A+s0,V3(t)=τ3·dV2(t)dt=(τH)33!·A+τ3·s0,dV3(t)dt=H·V3(t)+s1=τ3H43!·A+τH3·s0+s1,V4(t)=τ4·dV4(t)dt=(τH)44!·A+τ2H3·4·s0+τ4·s1,
where s0=M1(V1(t))·A, s1=M1(V1(t))·V1(t)+M1(V2(t))·V0(t) and H=M0+M1(A).Substituting (47) into (46) and assuming q=1, one has
(48)V(t)=V0(t)+V1(t)+V2(t)+V3(t)+⋯=[I+τH+(τH)22!+⋯]A+[τ3+τ2H3·4+τ3H23·4·5+⋯]s0+[τ4+τ2H4·5+τ3H24·5·6+⋯]s1+⋯=T·A+2·(H-1)3τ2[(τH)33!+(τH)44!+⋯]s0+3!·(H-1)4τ3[(τH)44!+(τH)55!+⋯]s1+⋯=T·A+2·(H-1)3τ2·R·s0+3!·(H-1)4τ3·Q·s1+⋯,
where R=[T-I-τH-(τH)2/2], Q=R-(τH)3/6, T=eτH. The matrix exponential function T can be calculated by using PIM.For simplicity, only the first two terms are taken as the approximate solution. Accordingly, the approximate solution is obtained(49)V(t)≈V0(t)+V1(t)=T·A+2·(H-1)3τ2·R·s0.


## 4. Numerical Results and Discussion

The Burgers equation has analytical solution as follows:(50)u(x,t)=4πRe·∑n=1∞nexp(-n2π2tRe)In(Re2π)sin(nπx)×(2∑n=1∞nexp(-n2π2tRe)In(Re2π)cos(nπx)+I0(Re2π))-1,where In(x) is the nth order modified Bessel function of the first kind. Figure 1 shows the analytical solutions with Re = 100 at t=0, 0.4, 0.6, respectively. It is easy to see that the analytical solution evolves into a shock wave near x=1, and the gradient becomes larger and larger with the increasing of Re.Figure 1
Analytical solutions of the Burgers equation at different times (t=0,0.4,0.6).In the analysis of the numerical result, we choose Reynolds numberRe=100, the scaling parameter j=8, the artificial parameter r=3, and the number of the external points L=4.The comparison of the classical PIM and the HAM is shown in Figures2 and 3 and Tables 1 and 2.Table 1
Difference between the classical PIM and the HAM whenτ=0.001 and τ=0.005.
t
τ
=
0.001
τ
=
0.005
Maximum of absolute error
Maximum of absolute error
Δ
=
|
Exact
-
HAM
|
Δ
=
|
Exact
-
PIM
|
Δ
=
|
Exact
-
HAM
|
Δ
=
|
Exact
-
PIM
|
0.1
0.001825
0.001862
0.001542
0.004271
0.2
0.000725
0.002339
0.001797
0.010446
0.3
0.000979
0.003806
0.002857
0.019874
0.4
0.000873
0.001949
0.003236
0.010663
0.5
0.000818
0.002270
0.001147
0.010949
0.6
0.001011
0.001955
0.001128
0.008018Table 2
Difference between the classical PIM and the HAM whenτ=0.05 and τ=0.1.
t
τ
=
0.05
τ
=
0.1
Maximum of absolute error
Maximum of absolute error
Δ
=
|
Exact
-
HAM
|
Δ
=
|
Exact
-
PIM
|
Δ
=
|
Exact
-
HAM
|
Δ
=
|
Exact
-
PIM
|
0.1
0.006067
0.029001
0.011791
0.049254
0.2
0.015968
0.084366
0.033981
0.137331
0.3
0.036130
0.186369
0.077718
0.283644
0.4
0.031070
0.272813
0.086969
0.436117
0.5
0.007711
0.295470
0.077946
0.463685
0.6
0.008237
0.544970
0.100736
0.710333Numerical results of the classical PIM and the HAM.
(a)
τ
=
0.05,  t=0.6
(b)
τ
=
0.1,  t=0.6Errors of the classical PIM and the HAM.
(a)
τ
=
0.005,  t=0.6
(b)
τ
=
0.001,  t=0.6We can see clearly that when the time stepτ is relatively large (τ=0.05 or τ=0.1), the classical PIM fails actually it fails when τ is 0.04. However, the HAM has enough precision even when the time step τ increases to 0.1 and t is more than 0.6.As the time stepτ becomes smaller, both the classical PIM and the HAM have higher precision; however, we can see clearly in Figure 3 that the HAM is still better than the classical PIM. Furthermore, from the comparison of (49) and (40), it can be seen that the computation cost of the proposed method is less than the classical PIM.


## 5. Conclusions

In this paper, it is shown that the proposed method can solve the Burgers equation successfully. The numerical comparisons presented demonstrate that the solutions obtained by the classical PIM are not valid for large time step; the method proposed in this paper can obtain higher accuracy than the classical PIM even if the time step is small. Although the method proposed aims at Burgers equation, it is a promising tool for other nonlinear problems.



---

*Source: 982810-2013-06-25.xml*