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.
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 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 denote the radial basis function, which can expressed as The corresponding first- and second-order derivative of the radial basis function can be derived as follows: Obviously, RBF is a compact support function with the interpolation property. is the support domain; and are constants.
A discrete point sequence of is defined as The subspaces of constructed from these discrete point sequence are defined as possesses the interpolation property, such as , so an interpolation operator can be defined as Without loss of generality, we set which has interpolation property as well; that is, where . Based on the interpolation theory, it is possible to define an interpolating transform mapping any continuous function into a sequence of its coefficients . Any function may be reconstructed from its transform by means of Under a minimal regularity assumption the coefficients and have the following meaning:
So, the multi-scale interpolating operator can be deduced as follows: where The operator maps a set of functional values at the level of resolution into a set of interpolation transform coefficients at level. For where the operator is the restriction operator defined as Since the restriction operator is known, it is easy to calculate the numerical solution of the interpolation operator from (9)~(11). The th derivative of is
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 and . is the set of the indices of collocation points whose interpolation transform coefficients are larger than a certain threshold at the level of resolution, and 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,
So, the adaptive interpolating operator and its th derivative can be written as where The th derivative of is
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: with the initial and boundary conditions where denotes time and Re denotes Reynolds number.
Following the classical collocation approach and the definition of the operator in (15), the approximating formulation of a function can be written as Substituting (20) into (18) leads to a system of nonlinear ordinary differential equations as follows: where . The corresponding vector expression is where According HAM, we can choose the auxiliary linear operator as follows: where is the embedded variable, is the function with respect to and , and is the known initial value . Based on (22), the nonlinear operator can be defined as follows: and 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: that is, where , . Let According to the Taylor series theory, can be expressed as the following power series: For simplification, (27) can be rewritten as which can be rewritten as In order to obtain the analytical solution of this nonlinear matrix differential equation, it is necessary to choose a suitable basis function to approximate 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 in (31) can be identified as Then Equation (31) can be rewritten as So, can be expressed as the Taylor series with respect to as follows: Substituting (35) into (31) and rearranging based on powers of -terms, we have Equation (36) is the system of homogeneous linear ODEs and its general solution is where . Both (37) and (38) are the system of inhomogeneous linear ODEs, and the general solutions are respectively, where , . Substituting (39) and (40) into (35), the numerical solution of (18) can be obtained. To improve the computational accuracy, the time interval can be divided evenly as Taking the solution at time instead of “” as the initial value in (40), the recurrence formula can be obtained. The matrix exponential function can be calculated accurately as follows: Let , where is a positive integer (usually take , and then ). Taking account of that is a small time interval, is a very small value, and then which is the Taylor series expansion of . In order to calculate the matrix more accurately, it is necessary to factorize the matrix as After times of factorization as above, a more accurate solution can be obtained.
The mesh parameter in the extrapolation method is usually taken as , where . It is easily noticed from (42) to (44) that the precise algorithm of exponent matrix is a -type algorithm. Let , and then Equation (45) can be described in program statement as At the end of the routine, it follows that Obviously, the sum of the resultant matrix at the end of th routine and the identity matrix is the matrix . And so, for calculating all , it needs merely to rewrite the above program statement as The amount of calculation is the same as the original program for . It should be pointed out that the computation accuracy is remained due to the decreasing of the time step although the cycle index in calculating the matrix is fewer than .
The calculation of the exponent matrix at different time steps is needed in solving nonlinear equations through iteration based on the precise integration method, and the algorithm of the matrix 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: where , .
4. Numerical Result and Discussion
The Burgers equation has analytical solution as follows: where is the th 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 . The figures of the analytic solution with at are shown in Figure 1, respectively. It is easy to find that the analytic solution evolves into a shock wave near , 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.
In our method, RBF was employed to construct the adaptive meshless interpolation operator and the regular width parameter , scale parameter , time step , and threshold parameter .
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. Figure 2 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.
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 and , respectively; the condition of stop iteration is . 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.
The 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 Figures 3(a) and 3(b), respectively. It should be pointed that the error obtained by HAM-based MAMIO occurs merely near 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 . 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.
(a) HAM
(b) Runge-Kutta method
The comparison of the classical adaptive precise integration method (APIM) and the HAM-based APIM is shown in Figure 5.
| (a) |
| (b) |
| (c) |
Figure 5(b) shows that APIM is failed when the time step is larger (). However, the HAM-based APIM has enough precision even when the time step increases to 0.2 and is more than 0.4 (Figure 5(c)).
As the time step becomes small enough (Figure 5, ), 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.