# HAM-Based Adaptive Multiscale Meshless Method for Burgers Equation

Authors: Shu-Li Mei
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

Based on the multilevel interpolation theory, we constructed a meshless adaptive multiscale interpolation operator (MAMIO) with the radial basis function. Using this operator, any nonlinear partial differential equations such as Burgers equation can be discretized adaptively in physical spaces as a nonlinear matrix ordinary differential equation. In order to obtain the analytical solution of the system of ODEs, the homotopy analysis method (HAM) proposed by Shijun Liao was developed to solve the system of ODEs by combining the precise integration method (PIM) which can be employed to get the analytical solution of linear system of ODEs. The numerical experiences show that HAM is not sensitive to the time step, and so the arithmetic error is mainly derived from the discrete in physical space.


---

## Body

## 1. Introduction

Burgers equation is a typical nonlinear partial differential equation, which was constructed to describe a kind of hydromechanical phenomenon. As there is a sharp wave in the solution of Burgers whose grads can vary with the time parameter, it is often employed to test the merits of the numerical algorithm. In recent years, there are many analytical methods for nonlinear problems that have been proposed in recent years. In these methods, the homotopy analysis method (HAM) proposed by Liao [1–6] is constantly being developed and applied to solve various nonlinear problems [7–13]. In order to further improve the properties of HAM, Liao [14] and some other researchers [15, 16] studied the choice rules about the auxiliary parameter ℏ and the auxiliary function H(t) on different nonlinear problems, which requires the users to have high-level skill about it. However, an effective way to develop HAM to solve the matrix differential equations is decoupling which is complicated obviously. In our research, we try to develop HAM to solve the matrix ODEs without decoupling by combining the precise integration method (PIM) and meshless method, which can be used to solve the nonlinear PDEs.Meshless methods eliminate some or all of the traditional mesh-based view of the computational domain and rely on a particle (either Lagrangian or Eulerian) view of the field problem. Compared with the traditional numerical method for PDEs such as the finite element method (FEM) and finite difference methods (FDM), meshless method is better suited to the problems associated with extremely large deformation and problems associated with frequency remeshing. During the past several decades, there are about 10 different meshless methods that have been developed, such as the Smooth Particle Hydrodynamics (SPH) [17], the Element-Free Galerkin (EFG) method [18], the Reproducing Kernel Particle (RKP) method [19], the Finite Point (FP) method [20], the hp clouds method [21], Meshless Local Petrov-Galerkin (MLPG) [22], Local Boundary Integral Equation (LBIE) [23, 24], and several others.Collocation method and Galerkin method are the two discretization methods which have been dominant in existing meshless methods, in which the radial basis function is often employed to construct the interpolation operator. Radial basis functions are functions which have no preferred direction but only depend on norms in space. As a powerful tool in both of pure and applied mathematics, multi-scale analysis theory has been developed rapidly in recent years. Multilevel systems are by now used widely in many fields of science and technology such as signal analysis, data compression, pattern recognition, and solutions of the partial differential equations [25, 26]. The major advantage of the multiscale analysis method is that it can improve the algorithm efficiency effectively. In recent years, many effective numerical methods for PDEs based on wavelets have been proposed, such as the wavelet Galerkin method (WGM) [27, 28], the wavelet finite element method [29], and the wavelet collocation method [30, 31]. So, construction of multi-scale interpolation operator with RBF is meaningful for meshless method.After the discretization to the nonlinear PDEs, a nonlinear system of ODEs can be obtained. HAM has been developed to solve the system of ODEs by Sami, Noorani, and Hashim, which can be employed to solve PDEs combining the various numerical discretization methods mentioned above.The purpose of this study is to construct a HAM-based multi-scale meshless method for Burgers equation. First, we construct a multi-scale interpolation operator with RBF, which can discretize PDEs adaptively into a system of nonlinear ODEs. Then, we develop HAM to solve the system of nonlinear ODEs obtained in the first step, in which the adaptability was introduced to the HAM. So, a new adaptive numerical-analytical method for nonlinear PDEs was obtained. At last, we test the algorithm by comparing with other algorithms. The numerical experiments show that this is an efficient second-order time-marching solver for time-dependent problems as long as a factorization of the differential operator is available.


## 2. Construction of Meshless Adaptive Multiscale Interpolation Operator

Letϕ(x) denote the radial basis function, which can expressed as
(1)ϕ(x)=exp(-r2x2σ2-x2).
The corresponding first- and second-order derivative of the radial basis function can be derived as follows:
(2)ddxϕ(x)=[-2r2xσ2-x2-2r2x3(σ2-x2)2]ϕ(x),d2dx2={-2r2σ2,x=00,x=±σ-2r2σ2(x2-σ2)4×(-σ4-2σ2x2+3x4+2r2x2σ2)×ϕ(x),otherwise.
Obviously, RBF is a compact support function with the interpolation property. [-σ,σ] is the support domain; r and σ are constants.A discrete point sequence ofϕ(x) is defined as
(3)mmmmmmmmmmlllllφj,k(x)=φ(2jx-k),x∈[xmin,xmax],j∈Z,k=0,1,2,…,2j.
The subspaces of L2(0,1) constructed from these discrete point sequence are defined as
(4)Vj:=span{ϕj,k∣k=0,1,2,…,2j}⊂L2(0,1),j∈Z,ϕj,k(x) possesses the interpolation property, such as ϕjk(n2-j)=δn,k, so an interpolation operator Ij:C0(0,1)→Vj can be defined as
(5)Ijf=∑k=02jf(xj,k)φj,k,xj,k=xmin+kxmax-xmin2j.
Without loss of generality, we set ψj,k=ϕj+1,2k+1 which has interpolation property as well; that is,
(6)ψjk(yjn)=δkn,ψjk(yj′n)=0,∀j′<j,
where yjk=xj+1,2k+1. Based on the interpolation theory, it is possible to define an interpolating transform mapping any continuous function f into a sequence of its coefficients {{βj0,k},{αj0,k},{αj0+1,k},…}. Any function f∈C0(0,1) may be reconstructed from its transform by means of
(7)f≈fj=∑k=02j0βj0kϕj0k+∑j≥j0∑k=02jαjkψjk.
Under a minimal regularity assumption the coefficients βjk and αjk have the following meaning:
(8)βj0k=f(xj0k),αjk=f(yjk)-Ijf(yjk).So, the multi-scale interpolating operator can be deduced as follows:(9)Ii(x)=∑k0=02j0Rk0,ij0,Jϕj0k(x)+∑j=j0J-1∑k∈ZjCk,ij,Jψj,k(x),
where
(10)Zj={0,1,2,…,2j},Ck,nj,J=R2k+1,nj+1,J-∑k0=02j0Rk0,nj0,Jϕj0,k0(xj+1,2k+1)-∑j1=j0j-1∑k1=02j1-1Ck1,nj1,Jψj1,k1(xj+1,2k+1)j=j0+1,j0+2,…,J-1,k=0,1,2,…,2j-1.
The operator Ck,nj,J maps a set of functional values at the J level of resolution into a set of interpolation transform coefficients at j level. For j=j0(11)Ck,ij0,J=R2k+1,ij+1,J-∑k0=02j0Rk0,ij0,Jθj0,k0(xj0+1,2k+1),
where the operator Ri,ml,j is the restriction operator defined as
(12)Ri,ml,j={1,xil=xmj,0,otherwise.
Since the restriction operator is known, it is easy to calculate the numerical solution of the interpolation operator from (9)~(11). The mth derivative of Ii(x) is
(13)Di(m)(x)=∑k0=02j0Rk0,ij0,Jϕj0k0(m)(x)+∑j=j0J-1∑k∈ZjCk,ij,Jψj,k(m)(x).For collocation method, it is necessary to setup a thresholdε so that the collocation points corresponding to interpolation transform coefficients less than ε can be neglected. For convenience, we define two integer subsets of integers Z≥j and Zc. Z≥j⊂Zj is the set of the indices of collocation points whose interpolation transform coefficients are larger than a certain threshold ε at the j level of resolution, and Zc⊂ZJ is the set of the indices of the ultimate set of collocation points, which are larger than a certain threshold used in the interpolation, that is,
(14)Zc=Z≥1∪Z≥2∪⋯∪Z≥J.So, the adaptive interpolating operator and itsmth derivative can be written as
(15)Ii(x)=∑k0=02j0Rk0,ij0,Jφj0k(x)+∑j=j0J-1∑(k∈Z≥j)∩(k∈[N1,N2])Ck,ij,Jψj,k(x),i∈Zc,
where
(16)Ck,nj,J=R2k+1,nj+1,J-∑k0=02j0Rk0,nj0,Jφj0,k0(xj+1,2k+1)-∑j1=j0j-1∑(k1∈Z≥j1)∩(k1∈[N1,N2])Ck1,nj1,Jψj1,k1(xj+1,2k+1),j=j0+1,j0+2,…,J-1,k∈Z≥j,n∈Zc.
The mth derivative of Ii(x) is
(17)Di(m)(x)=∑k0=02j0Rk0,ij0,Jφj0k0(m)(x)+∑j=j0J-1∑(k∈Z≥j)∩(k∈[N1,N2])Ck,ij,Jψj,k(m)(x),i∈Zc.


## 3. Coupling Technique of HAM and Adaptive Multilevel Interpolator on Nonlinear PDEs

HAM is an analytic approximation method for highly nonlinear equations in science, finance, and engineering. It was first proposed by Dr. Shijun LIAO in 1992 in his Ph.D, dissertation and modified and developed by Dr. Liao with his team and researchers in many other countries. As a result, the HAM overcomes the restrictions of all other analytic approximation methods mentioned above and is valid for highly nonlinear problems.Consider the Burgers equation:(18)∂u∂t+u∂u∂x=1Re∂2u∂2x,x∈[0,2],t≥0,
with the initial and boundary conditions
(19)u(x,0)=sin(πx),u(0,t)=u(2,t)=0,
where t denotes time and Re denotes Reynolds number.Following the classical collocation approach and the definition of the operatorIi(x) in (15), the approximating formulation uJ(x) of a function u(x) can be written as
(20)uJ(x)=∑i∈ZcIi(x)uJi.
Substituting (20) into (18) leads to a system of nonlinear ordinary differential equations as follows:
(21)∑n∈ZcuJ(xn,t)[1ReDn′′(xk)-uj(xk,t)Dn′(xk)]=∂uj(xk,t)∂t,
where k∈Zc. The corresponding vector expression is
(22)∂∂tVJ=M0VJ+M1(VJ)VJ,
where
(23)VJ=(uJ(x0,t),uJ(x1,t),…,uJ(x2J,t))T,M0(k,n)=mk,n0=1ReDn′′(xk),k,n∈Zc,M1(k,n)=mk,n1=-uJ(xk,t)Dn′(xk),k,n∈Zc.
According HAM, we can choose the auxiliary linear operator as follows:
(24)ℒ[Φ(t;q)]=∂Φ(t;q)∂t-(M0+M1(A))Φ(t;q),
where q∈[0,1] is the embedded variable, Φ(t;q) is the function with respect to t and q, and A is the known initial value Φ(0;q). Based on (22), the nonlinear operator can be defined as follows:
(25)𝒩[Φ(t;q)]=∂Φ(t;q)∂t-(M0+M1(Φ(t;q)))Φ(t;q).ℏ≠1 and H(t)≠0 denote the auxiliary parameter and the auxiliary function, respectively.According to the homotopy analysis theory, we can construct the 0th order deformed equation as follows:(26)(1-q)ℒ[Φ(t;q)-Φ(t;0)]=ℏqH(t)𝒩[Φ(t;q)]
that is,
(27)(1-q)[∂Φ(t;q)∂t-(M0+M1(A))Φ(t;q)]=ℏqH(t)[∂Φ(t;q)∂t-(M0+M1(Φ(t;q)))Φ(t;q)],
where Φ(t;0)=Vj0(t), Φ(t;1)=Vj(t). Let
(28)Vjm(t)=1m!∂mΦ(t;q)∂qm|q=0.
According to the Taylor series theory, Φ(t;q) can be expressed as the following power series:
(29)Φ(t;q)=Φ(t;0)+∑m=1∞1m!∂mΦ(t;q)∂qm|q=0qm=Vj0(t)+∑m=1∞Vjm(t)qm.
For simplification, (27) can be rewritten as
(30)[∂Φ(t;q)∂t-(M0+M1(A))Φ(t;q)][1-q-ℏqH(t)]+ℏqH(t)[M1(Φ(t;q))-M1(A)]Φ(t;q)=0
which can be rewritten as
(31)∂Φ(t;q)∂t=(M0+M1(A))Φ(t;q)-ℏqH(t)1-q-ℏqH(t)×[M1(Φ(t;q))-M1(A)]Φ(t;q).
In order to obtain the analytical solution of this nonlinear matrix differential equation, it is necessary to choose a suitable basis function to approximate Φ(t;q) based on the homotopy analysis method. It is well known that the matrix differential equation is quite different from the common differential equation, as the matrix one needs decoupling. To avoid the complicated decoupling, without loss of generality, the term ℏqH(t)/(1-q-ℏqH(t)) in (31) can be identified as
(32)p=τt=ℏqH(t)1-q-ℏqH(t),τ∈[0,t],ε∈[0,1].
Then
(33)Φ(t;q)≡Vj(t).
Equation (31) can be rewritten as
(34)dVjdt=M0Vj+M1(A)Vj-ε[M1(A)Vj-M1(Vj)Vj].
So, Vj(t) can be expressed as the Taylor series with respect to t as follows:
(35)Vj(t)=∑m=0∞1m!dmVj(t)dt.
Substituting (35) into (31) and rearranging based on powers of p-terms, we have
(36)p0:dVj0dt=(M0+M1(A))Vj0,(37)p1:dVj1dt=[M0+M1(A)]Vj1+[M1(Vj0)-M1(A)]Vj0,(38)p2:dVj2dt=[M0+M1(A)]Vj2+[M1(Vj0)-M1(A)]Vj1+M1(Vj1)Vj0.
Equation (36) is the system of homogeneous linear ODEs and its general solution is
(39)Vj0(t)=eHtA,
where H=M0+M1(A). Both (37) and (38) are the system of inhomogeneous linear ODEs, and the general solutions are
(40)Vj1(t)=eHt(H-1r0)-H-1r0,Vj2(t)=eHt(H-1r1)-H-1r1
respectively, where r0=(M1(Vj0)-M1(A))Vj0, r1=(M1(Vj0)-M1(A))Vj1+M1(Vj1)Vj0. Substituting (39) and (40) into (35), the numerical solution of (18) can be obtained. To improve the computational accuracy, the time interval [0,t] can be divided evenly as
(41)t0=0,t1=τ,…,tk=k·τ,….
Taking the solution at time tk instead of  “A” as the initial value in (40), the recurrence formula can be obtained. The matrix exponential function eHt can be calculated accurately as follows:
(42)exp(H·τ)=[exp(H·τ2N)]2N.
Let Δt=τ/2N, where N is a positive integer (usually take N=20, and then Δt=τ/1048576). Taking account of that τ is a small time interval, Δt is a very small value, and then
(43)exp(H·Δt)=I+Ta=I+HΔt+(HΔt)2[I+(HΔt)/3+(HΔt)2/12]2
which is the Taylor series expansion of exp(HΔt). In order to calculate the matrix exp(H·τ) more accurately, it is necessary to factorize the matrix as
(44)[exp(H·Δt)]2N=(I+Ta)2(N-1)(I+Ta)2(N-1).
After N times of factorization as above, a more accurate solution can be obtained.The mesh parameter in the extrapolation method is usually taken asτi=2-iτ, where τ=tk+1-tk. It is easily noticed from (42) to (44) that the precise algorithm of exponent matrix is a 2N-type algorithm. Let T(τ)=[exp(H·Δt)]2N, and then
(45)T(τi)=[exp(H·Δt)]2N-i.
Equation (45) can be described in program statement as
(46)for(iter=0;iter<N;iter++)Ta=2*Ta+Ta*Ta.
At the end of the routine, it follows that
(47)T=I+Ta.
Obviously, the sum of the resultant matrix Ta at the end of nth routine and the identity matrix I is the matrix T(τN-n). And so, for calculating all T(τi), it needs merely to rewrite the above program statement as
(48)for(iter=0;iter<N;iter++){Ta=2Ta+Ta×Ta;T(N-iter)=Ta+I;}
The amount of calculation is the same as the original program for T(τi). It should be pointed out that the computation accuracy is remained due to the decreasing of the time step τi although the cycle index in calculating the matrix T(τi) is fewer than N.The calculation of the exponent matrixT(ih) at different time steps is needed in solving nonlinear equations through iteration based on the precise integration method, and the algorithm of the matrix T(ih) presented here can obtain all the matrices at different time steps for once.It should be noticed that the nonlinear term in (22) can be expressed in Taylor series directly, and then another iterative format can be obtained. For comparison, we give this iterative format as follows directly:
(49)Vjk+1=T[Vjk+M0-1(s0+M0-1s1)]-M0-1(s0+M0-1s1+s1τ),
where s0=M1Vj, s1=ds0/dt.


## 4. Numerical Result 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)kkkkkkk×In(Re2π)cos(nπx)+I0(Re2π)(2∑n=1∞nexp(-n2π2tRe)))-1,
where In(x) is the nth order modified Besssel function of the first kind.In this section, the adaptive method which is proposed in this paper is used to calculate Burgers equation (18) with Reynolds number Re=1000. The figures of the analytic solution with Re=100 at t=0.0,0.4,0.6 are shown in Figure 1, respectively. It is easy to find that the analytic solution evolves into a shock wave near t=1, and the gradient becomes larger and larger with the increasing of Re. So the Burgers equation is often used to test the validity of numerical methods.Figure 1
Analytical solutions of the Burgers equation at different times.In our method, RBF was employed to construct the adaptive meshless interpolation operator and the regular width parameterσ=3Δ, scale parameter j=15, time step τ=0.02, and threshold parameter e=0.005.The evolution of the solution of Burgers equation from the uniformly smooth distribution to the shock structure causes the growth of the interpolation coefficients corresponding to the smaller scales,which in turn results in the refinement of the grids. Figure2 shows the progressive refinement of the irregular grid with the decreasing of the shock thickness. This illustrates that the meshless adaptive multi-level interpolation operator (MAMIO) constructed with RBF is effective, which is helpful to improve the efficiency of the algorithm.Figure 2
Evolution of the solution and collocation points for the solution of the Burgers equation with HAM-based MAMIO.Next, we will discuss the efficiency of the HAM-based adaptive precise integration method on nonlinear dynamics equations. Letε denote the maximum difference between the numerical solutions derived from the time steps Δt and Δt/2, respectively; the condition of stop iteration is |ε|<1×10-6. Table 1 shows the iteration times of the method. It is obvious that iteration times are adaptive with the evolvement of Burgers equation. That is, we do not have to set the smallest time step to satisfy the precision requirement, which is helpful to decrease the computation costs effectively.Table 1
Iteration times at different evolution time domain of the Burgers equation.
Time domain
0–0.13
0.13–0.22
0.22–0.31
0.31–0.39
0.39–0.4
Iteration times
2
3
3
4
3The numerical method for nonlinear PDEs constructed by combining RBF collocation method and Runge-Kutta method is a more common method. To test the validity of the multi-level adaptive PIM, these two methods are being compared in the following.The computational error curves of them are shown in Figures3(a) and 3(b), respectively. It should be pointed that the error obtained by HAM-based MAMIO occurs merely near x=1 and that at other points is almost equal to zero; however, the error obtained by Runge-Kutta method-based MAMIO is evident at all points, especially near the boundary and x=1. Furthermore, it should be noted that the error that occurs near the boundary is also evident even if the multistep methods, such as the Adams-Bashforth-Moulton method, are used (see Figure 4). The largest errors of these methods at different times are listed in Table 2. It can be seen that the computational precision of the HAM-based MAMIO is higher than that of the other two methods.Table 2
Largest computational errors of the Burgers equation.
t (s)
HAM
Runge-Kutta method
Adams-Bashforth-Moulton method
0.4
8.0949 × 10−5
9.1795 × 10−5
6.8626 × 10−5
0.8
1.7884 × 10−4
1.8425 × 10−4
1.4890 × 10−4
1.2
8.4061 × 10−6
1.1917 × 10−5
3.9004 × 10−5
1.6
1.6539 × 10−6
6.2157 × 10−6
3.4486 × 10−5
2.0
3.2393 × 10−7
9.6786 × 10−6
9.5912 × 10−6The computational errors of HAM the Runge-Kutta method (t=0.4).
(a)
HAM
(b)
Runge-Kutta methodFigure 4
The computational error of the Adams-Bashforth-Moulton method (t=0.4).The comparison of the classical adaptive precise integration method (APIM) and the HAM-based APIM is shown in Figure5.Numerical results of the classical APIM and the HPM-based APIM.
(a)
t
=
0.0
1
,
t
=
0
.
6
(b)
τ
=
0.02
,
t
=
0.08
(c)
τ
=
0.2
,
t
=
0
.
6Figure5(b) shows that APIM is failed when the time step τ is larger (τ=0.02). However, the HAM-based APIM has enough precision even when the time step τ increases to 0.2 and t is more than 0.4 (Figure 5(c)).As the time stepτ becomes small enough (Figure 5, τ=0.01), both of the APIM and the HAM-based APIM have higher precision however, the latter is still better evidently than the APIM.The above comparisons show that the HAM-based APIM is not sensitive to the time step compared to the APIM, which can improve the efficiency of the algorithm to some extent.


## 5. Conclusions

Coupling technique of HAM and the meshless multi-level adaptive interpolating operator constructed with RBF was developed to solve partial equations in a finite domain. In this method, the adaptive interpolating operator must be reconstructed at different moment because the amount and the position of collocation points are different at different time. The constructing efficiency of the interpolating operator is crucial to the efficiency of the adaptive method. Compared with multi-level wavelet operator, the dynamic operator proposed in this paper is more simple and efficient, as RBF possesses almost all the excellent numerical properties such as the compact support, analytical expression, and interpolation. This is helpful to improve calculation efficiency. Furthermore, the method becomes more efficient and accurate by means of HAM which possesses uniform convergence.



---

*Source: 248246-2013-09-05.xml*