Numerical Solution of High-Dimensional Shockwave Equations by Bivariate Multi-Quadric Quasi-Interpolation

: Radial basis function-based quasi-interpolation performs efﬁciently in high-dimensional approximation and its applications, which can attain the approximant and its derivatives directly without solving any large-scale linear system. In this paper, the bivariate multi-quadrics (MQ) quasi-interpolation is used to simulate two-dimensional (2-D) Burgers’ equation. Speciﬁcally, the spatial derivatives are approximated by using the quasi-interpolation, and the time derivatives are approximated by forward ﬁnite difference method. One advantage of the proposed scheme is its simplicity and easy implementation. More importantly, the proposed scheme opens the gate to meshless adaptive moving knots methods for the high-dimensional partial differential equations (PDEs) with shock or soliton waves. The scheme is also applicable to other non-linear high-dimensional PDEs. Two numerical examples of Burgers’ equation (shock wave equation) and one example of the Sine–Gordon equation (soliton wave equation) are presented to verify the high accuracy and efﬁciency of this method. meshless method


Introduction
Studying the numerical solution of high-dimensional non-linear Burgers' equation is of great significance.One reason is that it is a fundamental shockwave PDE from fluid mechanics which often occurs in applied and computational mathematics, such as modeling of dynamics, heat conduction, and acoustic wave [1,2].More importantly, it often plays the role of a testing equation to check the feasibility of the scheme, i.e., once the scheme is efficient for Burgers' equation, there is great possibility that it is applicable to other equations with shock or soliton waves, such as non-linear Schrodinger equation [3], Sine-Gordon equation [4], and Klein-Gordon equation [5].
In this paper, the proposed scheme is applicable to other non-linear high-dimensional PDEs.Burgers' equation has been studied by many researchers not only because its shockwave behavior when the coefficient R is large, but also because it is the simplest form of non-linear advection equation.Hence we take the 2-D coupled Burgers' equation as one example: The initial conditions are: u(x, y, 0) = f 1 (x, y), (x, y) ∈ D, v(x, y, 0) = f 2 (x, y), (x, y) ∈ D, and the boundary conditions are: u(x, y, t) = g 1 (x, y, t), (x, y) ∈ ∂D, v(x, y, t) = g 2 (x, y, t), (x, y) ∈ ∂D, (3) where t ∈ [0, T], D ⊂ R 2 is a bounded domain and ∂D is the boundary.u(x, y, t) and v(x, y, t) are the coupled solutions to be determined, f 1 (x, y), f 2 (x, y), g 1 (x, y, t) and g 2 (x, y, t) are given functions, and R is the Reynolds number.Fletcher [6] obtained its exact solutions, so that one could evaluate the quality of numerical solutions.
During the past decades, various numerical techniques have been developed for the solution of 2-D Burgers' equation.Bahadir [7] introduced one fully implicit finite difference scheme, then Zhu et al. [8] proposed the discrete Adomian decomposition method (ADM), Yu et al. [9] applied the bivariate quintic spline to solve numerically.
Compared with these mesh-based methods, which depend on a topology shape of the mesh, the meshless method performs satisfactorily for problems with complicated and irregular geometries [10].In recent years, the meshless collocation method has drawn considerable attention for solving PDEs [11][12][13].Mittal and Tripathi [14] proposed a collocation method based on Modified bi-cubic B-Spline functions.However, it requires to solve the inverse of a large-scale matrix on each time step, which costs computational time and might be ill-conditioned.
Considering the problems mentioned above, we try to apply the MQ quasi-interpolation method for the numerical solution of PDEs because of the following reasons:

•
Quasi-interpolation could give the approximant and its derivatives directly needing no solution of any large-scale system, hence it is easy to be implemented by the engineers in applications.

•
Quasi-interpolation could filter the noise of the data which often occurs in applications because of various reasons, e.g., measure error, data pollution, and so forth.

•
Ma [15] proved the superiority of MQ quasi-interpolation over finite difference and other interpolation methods in both stability and accuracy, especially when approximating scattered data.
The ambition of this paper is to construct a numerical scheme for solving 2-D Burgers' equation applying the bivariate MQ quasi-interpolation.Firstly, we develop a bivariate MQ quasi-interpolation and analyze its approximation order to a given function and its high order derivatives.Then the numerical scheme for solving the 2-D coupled Burgers' equation is proposed applying the bivariate MQ quasi-interpolation.In the numerical scheme, the spatial derivatives are approximated by using the derivatives of the quasi-interpolation, and the time derivatives of the dependent variables are approximated by the forward finite difference method.The truncation error and the total error of our scheme are estimated to show the accuracy of the proposed scheme.At last, two numerical examples of Burgers' equation (shockwave equation) and one example of the Sine-Gordon equation (soliton wave equation) are presented to verify the efficiency of the proposed method.
One advantage of the proposed scheme is its simplicity, easy implementation, but high accuracy and stability.More importantly, since the proposed method requires no triangulation of the region, adaptive refinement can be added in the scheme directly.In other words, it allows more flexible knots movement and one need not to worry the problems of mesh over gathering or condition numbers.In addition, the proposed scheme is applicable to other high-dimensional PDEs with shock or even soliton waves.
The rest of the paper is organized as follows.In Section 2, the bivariate MQ quasi-interpolant is introduced.In Section 3, we present the numerical scheme to solve 2-D Burgers' equation and give the error estimations.In Section 4, three numerical examples are proposed to verify our method.Section 5 concludes the whole paper.

Bivariate MQ Quasi-Interpolations
In this section, we firstly introduce the univariate MQ quasi-interpolation.On the basis, we develop the bivariate MQ quasi-interpolation and show its approximation order to the function and its derivatives.
Given the data {x j , f j }, f j = f (x j ), the univariate MQ quasi-interpolation is 1 is named MQ function, c 1 is shape parameter, which is often a small constant.As c 1 goes to zero, φ j (x) will converge to |x − x j |, consequently the univariate MQ interpolation could be taken as piecewise linear interpolation but possesses smoothness property.Denoting by h 1 = max j (x j+1 − x j ), the approximation error of Q f (x) to the original function f (x) and the k-th derivatives were given in [19]: holds, provided that c 1 = O(h Besides the high accuracy, the MQ quasi-interpolation has been proven to possess many other favorable properties, such as monotonicity preservation, shape preservation, and variation diminishing [17,18].Because of these positive properties, researchers began to extend the univariate quasi-interpolation to multi-dimension [24][25][26]. In this paper, we propose one kind of bivariate MQ quasi-interpolation using the tensor product technique.Specifically, given data (x i , y j , f ij ) and f ij = f (x i , y j ), the bivariate MQ quasi-interpolation is defined as: where ψ i (x) are described in ( 5) and ψ j (y) are represented as follows: is a small constant.Denoting h 2 = max j (y j+1 − y j ), the error estimation of the proposed bivariate MQ quasi-interpolation ( 7) is given in Theorem 2.
is bounded by a polynomial p 1 (x)p 2 (y), p 1 (x) and p 2 (y) are polynomials of degree α 1 + 2 − i and α 2 + 2 − j respectively (α 1 , α 2 , α 1 + 2 − i, and α 2 + 2 − j are all nonnegative integers), then the error of bivariate MQ quasi-interpolation (7) satisfies provided that c 1 = O(h Proof.By adding one middle term, the following inequalities are obtained: Since ∑ j f ij ψ j (y) and ∑ i f (x i , y)ψ i (x) are the approximation of f (x i , y) and f (x, y) separately, when c 1 = O(h 1 ) and c 2 = O(h 2 ), they are estimated as: ).
In addition, | ∑ i ψ i (x) |≤ 1, the result are as follows: Similarly, the derivative errors of (Q f ) (k) to f (k) could be estimated based on Theorem 1.This ends the proof.

Numerical Scheme Applying Bivariate MQ Quasi-Interpolant
In this section, we develop a numerical scheme for solving 2-D Burgers' Equations ( 1)-(3) based on bivariate MQ quasi-interpolation in Section 3.1.Then we provide the error estimations of Algorithm 1 in Section 3.2.

The Main Algorithm
We represent the approximations of u(x, y, t) and v(x, y, t) at the knots (x i , y j , t k ) = (ih 1 , jh 2 , kτ) by u k ij and v k ij respectively, where h 1 and h 2 are the mesh sizes in x and y direction separately, τ is the time step sizes, i.e., u k ij ≈ u(x i , y j , t k ), v k ij ≈ v(x i , y j , t k ).We use the similar expressions for the derivative approximations, e.g., (u x ) k ij ≈ u x (x i , y j , t k ), (u y ) k ij ≈ u y (x i , y j , t k ) and so forth.The numerical solution could be realized in the following Algorithm 1.

Error Estimations of Algorithm 1
In Section 3.2, the truncation error of Algorithm 1 is proposed in Lemma 1. Then based on Lemma 1, the total error of Algorithm 1 is provided in Theorem 3.
Algorithm 1 Numerical scheme of 2-D Burgers' equation applying the bivariate MQ quasi-interpolation.
Approximate the solutions' (u and v) first and second space derivatives by using the bivariate MQ quasi-interpolation.For example, the first derivatives of function u could be approximated as follows, the other derivatives can be approximated similarly based on the bivariate MQ quasi-interpolation.
2: Discrete the Burgers' equation in time and substitute the points (x i , y j ) into the derivatives, compute u k+1 ij and v k+1 ij : 3: Return to Step 1.
For the error estimation, we denote the following operators Proof.Firstly, based on the Taylor expansion, we have: Then based on Theorem 2, we get the following error estimations: 2 ).Finally, we get the truncation errors: This ends the proof.
Suppose u k (x, y) and v k (x, y) are the simulations of the exact solutions u(x, y, t k ) and v(x, y, t k ) using Algorithm 1, the total error of Algorithm 1 is provided in Theorem 3: Theorem 3. Defining the total errors of Algorithm 1 as e k u (x, y) = u k (x, y) − u(x, y, t k ), e k v (x, y) = v k (x, y) − v(x, y, t k ), if u(x, y, t) ∈ C 4 (R), v(x, y, t) ∈ C 4 (R), the error estimation of Algorithm 1 satisfies: where the Soblev norm for the function g(x, y) is defined as: Proof.Since u(x, y, t), v(x, y, t) are the exact solutions of (1), they satisfy: According to Lemma 1, u k (x, y), v k (x, y) satisfy Subtracting ( 10) from (11), one can get the following equation by adding some middle terms: Therefore, the Soblev norm of the errors satisfy B is a positive definite matrix, hence there exists a nonsingular matrix X satisfying Since e 0 u 2,∞ = 0, e 0 v 2,∞ = 0, we can get This completes the proof.
Remark 1.In this section, we apply the bivariate MQ quasi-interpolation for solving the 2-D coupled Burgers' equation.As is shown in Ma [15], the MQ method is more accurate, more stable, and simpler for simulating the high order derivatives than finite difference method, especially for scatted data or the data with noise.In addition, the proposed method is simple and easy to implement, since one need not to solve any large-scale linear system and could get the original function and derivatives' approximation directly.

2-D Burgers' Equation
In this section, we select two presentative test examples to illustrate the performance of the method described in the previous section.There are several ways to process the boundary, here the popular bivariate operator which appears in Ling [24] is employed.All these experiments are carried out on a computer with an Intel (R) Core (TM) i7-7500 CPU, 3.20 GHZ processor and 8.00 GB RAM.To check the validity of the scheme, the root mean square (RMS) and L ∞ error norms are applied to make comparisons with the previous methods [8].The error norms are defined as where u exact ij and u num ij are the exact and numerical results of u at the knot (x i , y j ).
Example 1.In this problem, we choose the computational domain as D = {(x, y) : 0 ≤ x ≤ 0.5, 0 ≤ y ≤ 0.5}, and the spatial knots distribute uniformly on the computational domain with a mesh width h 1 = h 2 = 0.025.For the Burgers' Equations (1) and (2), the initial conditions at t = 0 is Under the above conditions, the exact solutions could be given [6]: In this example, we observe the performance of our method with R = 100 for the convenience of comparison, and the numerical solutions are simulated using uniform knots, with a space width h 1 = h 2 = 0.025, the time step is τ = 10 −4 .The parameter is set as c = 2.65 × 10 −6 .
The study compares the errors of the presented method and Zhu et al.'s method [8].The numerical errors at time t = 0.1 and t = 0.4 are listed in Tables 1 and 2 respectively.It is observed that at time t = 0.1, our errors are similar(even a little larger) than Zhu et al.'s method [8].However, at time t = 0.4, our method is far smaller than Zhu et al.'s method [8].That is to say, our method could maintain a longer time simulation.[8] and the present method with R = 100, τ = 10 −4 , at t = 0.1.

Mesh Grid Error of u Error of v
Zhu et al. [8] The Present Zhu et al. [8] The Present In Table 3, the L ∞ errors and the CPU time of the present method at t = 0.1 and t = 0.4 are listed.One can conclude that our method is efficient in sense that it could get a satisfying error with rather little CPU time.Figure 1 describes the numerical solutions of u(x, y, t), v(x, y, t) by the present method, as expected, the present method could simulate Example 1 quickly and rather precisely.Example 2. When applying the Hopf-Cole transformation in [6], one could get the exact solutions of Burgers' Equations (1) and ( 2) In this example, we choose the computational domain as D = {(x, y) : 0 ≤ x ≤ 1, 0 ≤ y ≤ 1}, and we take the initial and boundary conditions from the above exact solution.The knots are distributed uniformly on the computational domain with a width h 1 = h 2 = 0.05, and R = 100.However the time step sizes of our method (τ = 10 −3 ) in this example is chosen to be ten times of the method in Zhu et al. [8], and Bahadir [7] (τ = 10 −4 ).The parameter c is chosen as c = 1.53 × 10 −5 .The numerical results will show that we could get the similar accuracy while with little CPU time.
In Table 4, the L ∞ numerical errors of the solutions at some typical mesh points at time t = 0.01, 0.05 and 0.2 with τ = 10 −3 .It is observed that our method could get satisfactory accuracy.In addition, the CPU time is rather little.The numerical results using the proposed algorithm are in good agreement with the exact solutions.In Table 5, we compare the numerical errors of the present method with time step τ = 10 −3 , and that of Zhu et al. [8], and Bahadir [7] with time step size τ = 10 −4 .It is observed that even using the time step that is ten times of the method in Zhu et al. [8], and Bahadir [7], our method could get similar accuracy.[8], Bahadir [7] with R = 100, τ = 10 −4 , at t = 0.01, and the present method with τ = 10 −3 .

Mesh Grid Error of u
Error of v Zhu et al. [8] Bahadir [7] The Present Zhu et al. [8] Bahadir [7] The Present The numerical solutions of u(x, y, t), v(x, y, t) by the present method are shown in Figure 2. In conclusion, Example 2 shows that the present method is easy to implement, accurate, and rather time saving.Example 3. In this example, we will show that the proposed scheme is also applicable to other high-dimensional PDEs with shock waves or even soliton waves.We take the typical soliton equation Sine-Gordon equation as the example: The initial conditions are: u(x, y, 0) = 4arctan(exp(x)) + 4arctan(exp(y)), − 6 ≤ x, y ≤ 6, v(x, y, 0) = 0, and the boundary conditions are: The example depicts the superposition of two orthogonal line solitons, which is also used by many researchers to verify their schemes [4,29].The numerical computations are simulated using equi-distributed knots, with a space width h 1 = h 2 = 0.2 and time steps τ = 0.01.In this example, the parameter c is chosen as c = 0.0075.Figure 3 shows that our methods could simulate the solitons rather exactly.That is to say, our scheme can solve a wild range of high-dimensional PDEs efficiently.

Conclusions
In this paper, we present a numerical scheme for 2-D Burgers' equation by bivariate MQ quasi-interpolation.The error estimations of the algorithm are also provided.Compared with the previous method (such as Zhu et al. [8]), the proposed method owns similar accuracy, while using far less computational effort.The present algorithm owns the following advantages: Firstly, it is very simple and easy to be implemented by the engineers in the applications.Secondly, it is of satisfactory accuracy and efficient for the time dependent PDEs with shock waves.Thirdly, it uses less computational effort, because it need not to solve any large-scale linear system.The method is also applicable to other high-dimensional time dependent PDEs [22,30], such as the Burgers-Fisher, the Klein-Gordon, and so forth.
Remark 2. Since the quasi-interpolation method have no requirement of the knots' topology structure and allows flexible knots moving, there are many works could be done to improve the efficiency and accuracy of the proposed algorithm with similar computational effort.One method is to introduce the moving knots strategy into the algorithm and let the nodes concentrated in the areas with shockwave, such as [31].Hence better approximation accuracy could be achieved with the same computational effort.Another method is using the symmetric quasi-interpolation to construct the algorithm that preserve the energy or some important property of the original PDE, e.g., [12,22].In this way, one could simulate the PDE for a longer time.

Table 1 .
Example 1: Comparisons of numerical errors by Zhu et al.

Table 2 .
Example 1: Comparison of the numerical errors by Zhu et al.

Table 4 .
Example 2:The L ∞ numerical errors and CPU time by the present method with R = 100, τ = 10 −3 at different time t.

Table 5 .
Example 2: The error of numerical results by Zhu et al.