[BACK]
images Computer Modeling in Engineering & Sciences images

DOI: 10.32604/cmes.2021.012575

ARTICLE

A Meshless Collocation Method with Barycentric Lagrange Interpolation for Solving the Helmholtz Equation

Miaomiao Yang, Wentao Ma and Yongbin Ge*

Institute of Applied Mathematics and Mechanics, Ningxia University, Yinchuan, 750021, China
*Corresponding Author: Yongbin Ge. Email: gyb@nxu.edu.cn
Received: 05 July 2020; Accepted: 21 September 2020

Abstract: In this paper, Chebyshev interpolation nodes and barycentric Lagrange interpolation basis function are used to deduce the scheme for solving the Helmholtz equation. First of all, the interpolation basis function is applied to treat the spatial variables and their partial derivatives, and the collocation method for solving the second order differential equations is established. Secondly, the differential matrix is used to simplify the given differential equations on a given test node. Finally, based on three kinds of test nodes, numerical experiments show that the present scheme can not only calculate the high wave numbers problems, but also calculate the variable wave numbers problems. In addition, the algorithm has the advantages of high calculation accuracy, good numerical stability and less time consuming.

Keywords: Helmholtz equation; Chebyshev interpolation nodes; Barycentric Lagrange interpolation; meshless collocation method; high wave number; variable wave number

1  Introduction

In the present work, we consider the following two-dimensional Helmholtz equation:

images

where images is the problem region, images is an unknown function of x and y. images is the wave number, which could be either a constant or a function of x and y. images is a source term. The Dirichlet boundary conditions are given as:

images

Helmholtz equation is an elliptical partial differential equation, which represents the solution of the wave equation independent of time and it arises from time-harmonic wave propagation [1]. This equation can simulate a variety of physical phenomena, including vibration analysis, water wave propagation, electromagnetic scattering, acoustic scattering, and radar scattering, etc. [26]. There is a long history about the development of wave propagation [7]. When the wavenumber is large, the solution of Helmholtz equation is highly oscillatory, so there are some difficulties in obtaining the numerical solution. Meanwhile, because of the complexity of variable wave number Helmholtz equation, it also brings great difficulties to numerical calculation, which still need to be further studied. Therefore, the research on numerical solution method of Helmholtz equation has important theoretical value and practical significance. A variety of numerical methods have been developed to solve the Helmholtz equation, such as the finite element method, boundary element method, finite difference method, meshless method and so on [811]. Among them, the finite element method is the most widely used method, but the calculation accuracy of the finite element method will sharply reduce with the increase of the wave number of the equation. At the same time, the meshless method is favored by many researchers because of its many advantages, such as simple discretization, easy derivation of format, easy realization of calculation program, high calculation accuracy and good numerical stability.

A new boundary method with method of external sources for eigenproblems was presented by Reutskiy [12] to solve Helmholtz equation in simply and multiply connected domains. However, this method depends on the grid, so when it applied to eigenproblem in irregular region, meshless technique is needed. Chen et al. [13] put forward a first-order system least square method for solving high wave number problems, and then a nontrivial decomposition is used to solve the first-order system least square problem, which was a new method different from the standard finite element method. Dogan et al. [14] studied the dispersion error of the Helmholtz equation by using the local meshless boundary integral equation method and the radial basis integral equation method. Both methods used the second-order polynomial and the frequency dependent radial basis function to interpolate the potential field, which can solve the high wave number problems. But the cost is to solve a larger linear system, which is computationally intensive. Based on the meshless and boundless analysis of Helmholtz equation, Chen et al. [15] proposed a boundless Burton Miller method by using Burton Miller formula. This method can find the unique solution for all wave numbers, and can also deal with Helmholtz equation under the high wave numbers. However, the cost of calculation is expensive and the derivation process is complicated. Qu et al. [16] developed a local basic solution method which was different from the traditional basic solution method. According to the nodes, the calculation domain was divided into several overlapping subdomains. Finally, the sparse strip matrix was used to solve large-scale problems. Whereas, the numerical implementation is complicated. By using an asymptotic expansion that removes the singularities up to several leading orders, Britt et al. [17] put forward a method that can keep high order convergence even if there is singular solution for Helmholtz equation. Whereas, it did not give out the results in the case of large wave number. A circular boundary knot method was studied by Lin et al. [18] by placing the collocation nodes on the boundary surface in a circular form to solve the high wave number Helmholtz equation, the BKM matrix obtained has the block structure of cyclic matrix and can be solved effectively. Which needs to solve a lot of small matrices, the calculation storage capacity is large. Wang et al. [19] considered the Helmholtz equation in the outer domain, and proposed a finite difference scheme in polar or spherical coordinates, which can solve high wave number problems. Whereas, the “pollution effect” in the domain of the origin of coordinates is inevitable. In addition, the common disadvantage of the above literature is that there is few mention of the case of variable wave number problem.

In this paper, based on the Chebyshev nodes, the scheme for solving the Helmholtz equation is derived by using barycentric Lagrange interpolation basis function. This method is simple in theory, needs few interpolation nodes, and the final discrete matrix is easy to be dealt with. It can not only calculate the high wave number problems, but also calculate the variable wave number problems. It has the advantages including short calculation time, high accuracy and good stability. In this meshless method, the point collocation technique is used for forming the final system of linear algebraic equations. Thus, the present method is a variant of the well-known finite point method [20]. And another analogous method is the localized method of fundamental solutions [21]. The remainder of this paper is arranged as follows: in Section 2, the calculation formula of the two-dimensional barycentric Lagrange interpolation and some theoretical knowledge of the simplified matrix are introduced. In Section 3, the interpolation collocation method and its discrete matrix for calculating the two-dimensional Helmholtz Eq. (1.1) are derived. In Section 4, the numerical results for some test problems are given and compared with the results in the literature. Finally, conclusions are included in Section 5.

2  Barycentric Lagrange Interpolation

2.1 Barycentric Lagrange Interpolation Formula

Let images be a given interpolation node, and the corresponding function value is fi, then Lagrange polynomial interpolation formula can be expressed as

images

where images is Lagrange interpolation basis function in the form of

images

According to Eq. (2.1), it can be concluded that when new interpolation nodes are added, the interpolation basis function needs to be recalculated, which increases the calculation amount. However, if we start from the basis function and extract the common part, its expression will be changed and the calculation amount will be reduced consequently.

Set

images

Defines weight of the barycentric Lagrange interpolation function as

images

As usual, we can get the following algorithm flow diagram (see Fig. 1) for computing images

images

Figure 1: The algorithm flow diagram of images

Then, the interpolation basis function (2.2) can be changed to

images

Substitute Eq. (2.5) into Eq. (2.1) to get

images

The Eq. (2.6) is called the improved Lagrange interpolation formula [22]. The part of main improvement is to construct Eq. (2.3), thus rewriting the interpolation basis function Eq. (2.2) into the form of Eq. (2.5). The purpose is with number of Lagrange interpolation nodes increases, the computation amount decreases a lot, from images to images.

By using Eq. (2.6) to interpolate the constant 1 with images then we have

images

After simplification, it can be obtained

images

Substituting Eq. (2.7) into Eq. (2.6), we have

images

The Eq. (2.8) is called the barycentric Lagrange interpolation formula [22]. Compared with Eq. (2.6), Eq. (2.8) overcomes the Runge phenomenon effectively while maintaining the same computational cost images.

2.2 Bivariate Barycenter Lagrange Interpolation Formula

Let the images be a function about the variable x and y, and define the interval as images. When the interval images and images is discretized into m and n nodes, respectively, then there are images tensor product interpolation nodes in the whole domain images. Then the barycentric Lagrange interpolation formula [23] of the images is

images

where images and Mj(y) are the interpolation basis functions in the x and y direction, respectively.

images

where the weight is

images

2.3 Interpolation Knot

We can see from the expression of the barycentric Lagrange interpolation weight function Eqs. (2.4) or (2.11) that the selection of the interpolation weight is only related to the distribution of the nodes. Therefore we consider the three kinds of nodes as interpolation nodes and test nodes. They are random nodes, uniform nodes and Chebyshev nodes. Because random nodes and uniforms nodes can be directly implemented in programming language, this part only gives the Chebyshev nodes formula:

Chebyshev nodes of the first kind:

images

Chebyshev nodes of the second kind:

images

What needs to be explained is that for two types of Chebyshev node formulas, the defined interval is images As for the general interval images the coordinate transformation formula images of the interval can be used.

2.4 Differential Matrix of Barycenter Lagrange Interpolation

The differential matrix was first found in the study of Chebyshev quasi-spectral method [24]. The barycentric Lagrange interpolated differential matrix can directly obtain the derivative function of the unknown function in the Helmholtz equation at the computing node. Therefore, the differential matrix is a very important part in the solution of the barycentric interpolation Lagrange collocation method.

When m nodes are inserted on the interval images, there is images, to set the function value of the unknown function images at the node xi is images, then the barycentric Lagrange interpolation formula of the function images is

images

In which, images is the barycentric interpolation basis function, and images is the weight of gravity interpolation.

Using Eq. (2.14), the k derivative of the function images at the node images can be expressed as

images

In the form of a matrix as

images

where images and images are the vector. The matrix images is as the k order differential matrix of the function images, and images In addition, we can calculate images (see [25]) with

images

When the solution of some partial differential equation is obtained by numerical method, it is also necessary to use the differential matrix to calculate the derivative values of the unknown function on calculating nodes.

2.5 Kronecker Product of a Matrix

Kronecker product is an operation form between two matrices of any size and a special form of tensor product. The Kronecker product of any matrix images and matrix images is images recorded as following

images

where images

Using the Kronecker product of the barycentric Lagrange interpolation differential matrix, the discrete equation of the Helmholtz equation can be expressed as a simple matrix form.

2.6 Relationship between Partial Derivative of Differential Equation and Differential Matrix

The Kronecker product with the barycentric Lagrange interpolated differential matrix can also be used. The unknown function in the Helmholtz equation has the following corresponding relation between the partial derivative of the variable and the Kronecker product of the corresponding differential matrix.

The number of nodes for the unknown function images in x and y directions of the two-dimensional Helmholtz equation are m and n then by using Eq. (2.17), the expressions of the differential matrices corresponding to the partial derivative of the differential equation about the variables x and y can be briefly expressed as

images

where Im and In represent the unit matrices of m order and n order, respectively, and the images and images are barycentric interpolates of the second-order matrices at the node images and images respectively.

Using the above marks, the discrete forms of partial differential equations and boundary conditions can be written directly into matrix form, which makes programming much easier.

3  Discrete Equations of the Two-Dimensional Helmholtz Equation

By substituting the barycenter Lagrange interpolation Eq. (2.9) into Eq. (1.1), we can obtain

images

If Eq. (3.1) holds at the nodes of images, then Eq. (3.1) is replaced by the following

images

Set images. From the properties of the Lagrange interpolation basis, images and using Eqs. (2.15), (2.16) and (2.18), then Eq. (3.2) can also be expressed as

images

Let

images

Then Eq. (3.3) can be written into

images

That is

images

In which images.

The above two-dimensional Helmholtz equation is discretized by the barycentric Lagrange interpolation method, and the final scheme is Eq. (3.5). In addition, the boundary conditions need to be discretized, and the following discrete formulas for boundary conditions are elaborated.

By employing Eqs. (2.17) and (2.18), we obtain the matrix form of the boundary condition in Eq. (1.2)

images

where images and images represent the first and m lines of the first order differential matrix of the m order unit matrix, images and images represent the first and n lines of the first order differential matrix of the n order unit matrix, and images

images

Combining Eqs. (3.5) and (3.6), we can obtain the value of the function images on each node.

4  Numerical Examples

In this section, we solve some test problems to demonstrate the effectiveness and accuracy of the barycentric Lagrange interpolation method. Six numerical examples of Helmholtz equation including high wavenumber and variable wavenumber problems are given. In which, the maximum absolute error (images) and relative error ( Er) are defined as follows:

images

where uij and ueij represent numerical and exact solutions, respectively. Nt and Mt are the numbers of test nodes in the x and y directions, respectively.

Example 1 [26]:

Consider the following non-homogeneous Helmholtz equation

images

with the boundary conditions:

images

The exact solution is

images

Tab. 1 gives the calculation results of images for the four kinds of interpolation nodes with Nt = Mt = N, k = 5. In which, N is the number of interpolation nodes in the domain, Nt and Mt are the numbers of test nodes, and the test nodes type are the second kind of Chebyshev nodes (that is Chebyshev II) rather than the first kind of Chebyshev nodes (that is Chebyshev I). The results show that the random interpolation nodes are bad for our proposed method, the uniform interpolation nodes are not good and the two kinds of Chebyshev interpolation nodes are very good. What’s more, the numerical results show that the accuracy of the Chebyshev node is highest, the stability is best, and the calculation error is smallest. Therefore, we use the two kinds of Chebyshev nodes rather than the random or uniform nodes as the interpolation nodes to calculate the Tabs. 2 and 3 for Example 1. In addition, the numerical results of the Chebyshev II is a little better than that of the Chebyshev I.

Table 1: The images with four kinds of interpolation nodes when Nt = Mt = N, k = 5 for Example 1

images

Table 2: The images with two kinds of interpolation nodes and three kinds of test nodes when Nt = Mt = 100, k = 5 for Example 1

images

Table 3: The images with two kinds of interpolation nodes and three kinds of test nodes when N = 24, k = 5 for Example 1

images

Then the images with two kinds of Chebyshev interpolation nodes and three kinds of test nodes when Nt = Mt = 100, k = 5 and N = 24, k = 5 is given in Tabs. 2 and 3, respectively. What should we know is that the two kinds of interpolation nodes here are Chebyshev I and the Chebyshev II. Three kinds of test nodes are random test nodes, uniform test nodes and Chebyshev II test nodes all the time. According to the calculation results, firstly, we can know that the results of Chebyshev II is better than Chebyshev I and there is almost no difference between three kinds of test nodes. In addition, the error, error order and the convergence order is almost constant as Nt and Mt change, which means that the number of test nodes has little influence on our numerical results. So we can take the value of Nt and Mt with any reasonable number in the left numerical examples. What’s more, the numerical results of the Chebyshev II is a little bit better than the Chebyshev I and the convergence order with the two kinds of interpolation nodes is almost identical, therefore, we use the Chebyshev II interpolation nods and three kinds of test nodes to calculate the left numerical examples in our work.

Finally, the calculation results of images with the 4th-order method in [26] and with the present scheme are compared when Nt = 130, Mt = 140, k = 30 in Tab. 4. The results show that the error order and the convergence order of the scheme in this paper reaches images and 12th-order, respectively when the number of collocation points is N = 128 , while that in the literature just reaches images and 4th-order, respectively. It shows that the calculation accuracy of the present scheme is much better than that in the literature.

Table 4: The images when Nt = 130, Mt = 140, k = 30 for Example 1

images

The node distribution for the four kinds of interpolation nodes with Nt = Mt = N, k = 5 is shown in Fig. 2. We can observe that the random nodes are distributed randomly within the interval, uniform nodes are distributed equably in the interval and the two kinds of Chebyshev nodes are concentrated at the boundaries while the middle is sparse relatively.

images

Figure 2: Random (a), uniform (b), Chebyshev I (c) and Chebyshev II (d) interpolation nodes when Nt = Mt = N, k = 5 for Example 1

Then the exact solution with Chebyshev test nodes is plotted in Fig. 3a-III, numerical solutions with random test nodes are plotted in Fig. 3b-I, uniform test nodes in Fig. 3b-II and Chebyshev test nodes in Fig. 3b-III; errors with random test nodes in Fig. 3c-I, uniform test nodes in Fig. 3c-II and Chebyshev test nodes in Fig. 3c-III with Nt = 130, Mt = 140, k = 30, N = 12. Among them, because the exact solutions with three kinds of test nodes are almost the same, we only give the figure with Chebyshev test nodes. In addition, it can be observed that all numerical solutions with the three kinds of test nodes agree well with the exact solution. What’s more, because of the addition of test nodes on the boundaries to calculate the numerical solution with proposed method, these nodes’ values on the boundaries are no longer given by the exact solution and resulting in some error’s appearance on the boundaries. It can also be observed from Fig. 3c that the errors are equivalent with that in the interior region.

imagesimages

Figure 3: Exact (a) and numerical (b) solutions and error (c) with different test nodes: random (I), uniform (II) and Chebyshev (III) when Nt = 130, Mt = 140, k = 30, N = 128 for Example 1

Finally, the convergence order with Nt = Mt = 100, k = 5 and Nt = 130, Mt = 140, k = 30 is shown in Fig. 4. It can be seen that with the wave number of k increases, the order of convergence decreases. Meanwhile, the convergence order of the two kinds of interpolation nodes Chebyshev I and Chebyshev II is almost identical. And then the convergence order of the proposed method is much higher than that in the literature.

images

Figure 4: The convergence order with Nt = Mt = 100, k = 5 and Nt = 130, Mt = 140, k = 30 for Example 1

Example 2 [27]:

Consider the following non-homogeneous Helmholtz equation

images

with the boundary conditions:

images

The exact solution is

images

Tab. 5 gives the images computed by the present scheme and compares with EB-6 method in [27] when images. The results show that the error order of the proposed scheme reaches images when the number of collocation points is N = 16 , while that in the literature, the error order only reaches images when N = 16 and images when N = 64. Therefore, the number of calculation points needed in the present scheme is much less than that in the literature and we can save more computational amount. In addition, we notice the error order of the proposed scheme decreases gradually and the error increase gradually with the increase of the number of nodes N , but the overall calculation result is better than that in the literature. What’s more, convergence order of the proposed method is the 21th-order, while that in the literature is the 6th-order, so the precision of present method is also better than the method in the literature. Compared with the finite difference method, the barycentric Lagrange interpolation collocation method does not mean that the more nodes there are, the better the calculation results will be. Sometimes, we can obtain the better results with comparably a few nodes.

Table 5: The images when images for Example 2

images

The images and CPU time with images are given in Tab. 6. It shows that the proposed scheme can approximate the exact solution with a few nodes. It is necessary to note that for the determination of CPU time in our work, we select the longest time taken by the three test nodes method. Even so, we still find that for the two-dimensional Helmholtz equation, the CPU time with the present scheme is very small when the number of collocation points is less, which shows that this method is less time-consuming and can save calculation time effectively.

Table 6: The images and CPU time when images for Example 2

images

Fig. 5. gives the exact (a) and numerical (b) solutions and error (c) with different test nodes: random (I), uniform (II) and Chebyshev (III) with images It shows that the numerical solution with three kinds of test nodes have the similar shape and are in good agreement with the exact solution as well. What’s more, Figs. 3c-I, 3c-II and 3c-III show that the error order reaches images.

imagesimages

Figure 5: Exact (a) and numerical (b) solutions and error (c) with different test nodes: random (I), uniform (II) and Chebyshev (III) when images for Example 2

Example 3 [28]:

Consider the following non-homogeneous high wavenumber Helmholtz equation

images

with the boundary conditions:

images

The exact solution is

images

The calculation results with the 6th-order method in [28] and the present scheme are compared under different wave numbers with k = 24 and k = 48 when a = 1, Nt = 130, Mt = 150 in Tab. 7. The results show that the calculation error and error order of the present scheme are much better than that in the literature. Therefore, it is more suitable for solving high wavenumber problems as well. Then the images and CPU time with a = 0.5, Nt = Mt = 60, k = 12 are given in Tab. 8. It shows the advantages of the present scheme, such as less number of collocation points, short calculation time, high precision and small error.

Table 7: The images with different wave numbers when a = 1, Nt = 130, Mt = 150 for Example 3

images

Table 8: The images and CPU time when a = 0.5, Nt = Mt = 60, k = 12 for Example 3

images

Fig. 6. gives the exact (a) and numerical (b) solutions and error (c) with different test nodes: random (I), uniform (II) and Chebyshev (III) with a = 1, Nt = 130, Mt = 150, k = 48, N = 128, respectively. It shows that the numerical solutions agree very well with the exact solution, and the error of three kinds of test nodes has little difference.

imagesimages

Figure 6: Exact (a) and numerical (b) solutions and error (c) with different test nodes: random (I), uniform (II) and Chebyshev (III) when a = 1, Nt = 130, Mt = 150, k = 48, N = 128 for Example 3

Example 4 [29]:

Consider the following non-homogeneous variable wavenumber Helmholtz equation

images

with the boundary conditions:

images

The exact solution is

images

Tab. 9 shows the images when Nt = 150, Mt = 130 with variable wavenumber images. Compared with the MHADI method in [29], the error order of the present scheme is images with N = 16 , while that in the literature is only images. In addition, the accuracy in the literature is the 4th-order, while that in this paper is up to the 9th-order. So the present scheme has higher calculation accuracy compared with that in the literature. We notice that although the error order of the proposed scheme decreases with the increase of the number of nodes, the overall calculation results are still better than that in the literature. Then the images at different collocation points and CPU time with Nt = Mt = 50, images are given in Tab. 10. The numerical results show that compared with the finite difference, the barycentric Lagrange interpolation method is more suitable for solving the problem of variable wavenumber. At the same time, we can obtain the better results with relatively a few nodes.

Table 9: The images when images for Example 4

images

Table 10: The images and CPU time when images for Example 4

images

Fig. 7. gives the exact (a) and numerical (b) solutions and error (c) with different test nodes: random (I), uniform (II) and Chebyshev (III) with images respectively. It shows that although k is the variable wavenumber, the numerical solutions are also in good agreement with the exact solution, and the error is steady as a whole.

imagesimages

Figure 7: Exact (a) and numerical (b) solutions and error (c) with different test nodes: random (I), uniform (II) and Chebyshev (III) when images for Example 4

Table 12: The images and CPU time when Nt = Mt = 50, k = 10 for Example 5

images

Example 5 [30]:

Consider the following homogeneous Helmholtz equation

images

with the Dirichlet boundary condition. The computational domain is an ameba-shape domain (see Fig. 7) with the boundary defined by the parametric equation:

images

The exact solution is

images

Tab. 11 and Tab. 12 give the images and CPU time with an ameba-shape domain when Nt = Mt = 50, k = 5 and Nt = Mt = 50, k = 10, respectively. What should be noted is that the [30] only computed a square domain, so we do not make comparisons with it and just give the results of our proposed method. The results show that even in an arbitrary region, the accuracy and convergence order of the present method are still very high. It is obvious that the present method is an effective method with high precision.

Table 11: The images and CPU time when Nt = Mt = 50, k = 5 for Example 5

images

Then the Fig. 8. gives the computational domain with different test nodes: random (a), uniform (b) and Chebyshev (c) when Nt = Mt = 50, k = 10, N = 48. Similar to Fig. 1, we can observe the distribution of the three test nodes. Fig. 9. gives the exact (a) and numerical (b) solutions and error (c) with different test nodes: random (I), uniform (II) and Chebyshev (III) with Nt = Mt = 50, k = 10, N = 48. It shows that the even though under the general domain, the numerical solutions agree very well with the exact solution, and the errors of three kinds of test nodes have little difference.

images

Figure 8: The ameba-shape domain with different test nodes: random (a), uniform (b) and Chebyshev (c) when Nt = Mt = 50, k = 10, N = 48 for Example 5

imagesimages

Figure 9: Exact (a) and numerical (b) solutions and error (c) with different test nodes: random (I), uniform (II) and Chebyshev (III) when Nt = Mt = 50, k = 10, N = 48 for Example 5

Example 6 [30]:

Consider the following non-homogeneous modified Helmholtz equation

images

with the boundary conditions:

images

The exact solution is

images

Lastly, we consider the non-homogeneous modified Helmholtz equation for Example 6. Tab. 13 gives calculation results with the present scheme and compared with [30]. Then the images and CPU time with different collocation points under Nt = Mt = 50 is given in Tab. 14. We find that the error and convergence order of our proposed method are still much better than that in the literature. Fig. 10 gives the exact (a), numerical (b) solutions and error (c) with different test nodes: random (I), uniform (II) and Chebyshev (III) when Nt = Mt = 50, N = 16 for Example 6. It is interesting that the proposed scheme is also suitable for non-homogeneous modified Helmholtz equation, and the numerical solution and the exact solution are in good agreement as well, at the same time the error is stable as a whole.

Table 13: The images when Nt = Mt = 50 for Example 6

images

Table 14: The images and CPU time when Nt = Mt = 50 for Example 6

images

imagesimages

Figure 10: Exact (a) and numerical (b) solutions and error (c) with different test nodes: random (I), uniform (II) and Chebyshev (III) when Nt = Mt = 50, N = 16 for Example 6

5  Conclusions

Based on the Chebyshev interpolation nodes, we have developed a meshless collocation method for solving the two-dimensional Helmholtz equations with the Dirichlet boundary condition, by using the barycentric Lagrange interpolation basis function. Firstly, this meshless collocation method applies the interpolation basis function to treat the spatial variables and their partial derivatives, and establishes the collocation method for solving the differential equations. Secondly, at a given node, the multivariate basis function is substituted into the equation and its boundaries and the discrete algebraic equations are obtained. Then, the differential matrix is used to simplify it. To investigate the accuracy and efficiency of the proposed method, we conduct some numerical experiments by using three kinds of test nodes. Based on our research, we can draw some conclusions as follows:

1.    As for the choice of interpolation nodes, the numerical results show that Chebyshev II is better than Chebyshev I for the present scheme. As for test nodes, there is almost no difference for the numerical results with three kinds of test nodes. In addition, with the change of the number of test nodes Nt and Mt, the error is almost unchanged. That is say that the number of test nodes has little influence for the numerical results.

2.    By using comparably fewer nodes, the present scheme can obtain high-precision numerical results and keep good stability compared with those methods in the literature. It shows the advantages of the present scheme, such as smaller number of collocation points needed, short calculation time, high precision, small error and high efficiency. What’s more, the numerical solutions agree well with the exact solutions.

3.    Compared with the methods in the literature, the present method can not only calculate the high wavenumber problems, but also calculate the variable wavenumber problems, and the numerical results are better than those in the literature. For modified Helmholtz equation, the results computed by the present scheme are far better than those in the literature as well.

4.    The present method can be extended to solve the 3D Helmholtz equations. We will report it in the near future.

Funding Statement: The authors would like to thank very mush to the anonymous reviewers whose constructive comments are helpful for this paper’s revision. This work is partially supported by National Natural Science Foundation of China (11772165, 11961054, 11902170), Key Research and Development Program of Ningxia (2018BEE03007), National Natural Science Foundation of Ningxia (2018AAC02003, 2020AAC03059), and Major Innovation Projects for Building First-class Universities in China’s Western Region (Grant No. ZKZD2017009).

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

 1.  Assari, P., Fatemeh, A., Cuomo, M. S. (2019). A numerical scheme for solving a class of logarithmic integral equations arisen from two-dimensional Helmholtz equations using local thin plate splines. Applied Mathematics and Computation, 356, 157–172. DOI 10.1016/j.amc.2019.03.042. [Google Scholar] [CrossRef]

 2.  Wang, K. Y., Wang, Q. S. (2013). Study on vibration of membranes with Taylor polynomial method and error analysis for Helmholtz equation. Advanced Materials Research, 763, 234–237. DOI 10.4028/www.scientific.net/AMR.763.234. [Google Scholar] [CrossRef]

 3.  Jo, C. H., Shin, C., Suh, J. H. (1996). An optimal 9-point, finite-difference, frequency-space, 2-D scalar wave extrapolato. Geophysics, 61(2), 529–537. DOI 10.1190/1.1443979.

 4.  Du, K. (2011). A composite preconditioner for the electromagnetic scattering from a large cavity. Journal of Computational Physics, 230(22), 8089–8108. DOI 10.1016/j.jcp.2011.07.011.

 5.  Wang, S., Hoop, M. V. D., Xia, J. L. (2010). Acoustic inverse scattering via Helmholtz operator factorization and optimization. Journal of Computational Physics, 229(22), 8445–8462. DOI 10.1016/j.jcp.2010.07.027.

 6.  Liu, G., Jayathilake, P. G., Khoo, B. C., Han, F., Liu, D. K. (2012). Conformal mapping for the Helmholtz equation: Acoustic wave scattering by a two dimensional inclusion with irregular shape in an ideal fluid. Journal of the Acoustical Society of America, 131(2), 1055–1065. DOI 10.1121/1.3675947. [Google Scholar] [CrossRef]

 7.  Linton, C. M. (2010). Lattice sums for the Helmholtz equation. SIAM Review, 52(4), 630–674. DOI 10.1137/09075130X. [Google Scholar] [CrossRef]

 8.  Swager, M. R., Gray, L. J., Fata, S. N. (2010). Galerkin boundary integral analysis for the 3D Helmholtz equation. Computer Modeling in Engineering & Sciences, 58(3), 297–234. [Google Scholar]

 9.  Tsai, C. C., Young, D. L. (2013). Using the method of fundamental solutions for obtaining exponentially convergent Helmholtz eigensolutions. Computer Modeling in Engineering & Sciences, 94(1), 175–205. [Google Scholar]

10. Wu, T. T. (2017). A dispersion minimizing compact finite difference scheme for the 2D Helmholtz equation. Journal of Computational and Applied Mathematics, 311, 497–512. DOI 10.1016/j.cam.2016.08.018. [Google Scholar] [CrossRef]

11. Liu, G. R., Gu, Y. T. (2005). An introduction to meshfree method and their programming, Berlin, New York: Springer. [Google Scholar]

12. Reutskiy, S. Y. (2006). The method of external sources (MES) for eigenvalue problems with Helmholtz equation. Computer Modeling in Engineering & Sciences, 12(1), 27–39. [Google Scholar]

13. Chen, H. X., Qiu, W. F. (2017). A first order system least squares method for the Helmholtz equation. Journal of Computational and Applied Mathematics, 309, 145–162. DOI 10.1016/j.cam.2016.06.019. [Google Scholar] [CrossRef]

14. Dogan, H., Popov, V., Ooi, E. H. (2015). Dispersion analysis of the meshless local boundary integral equation and radial basis integral equation methods for the Helmholtz equation. Engineering Analysis with Boundary Elements, 50, 360–371. DOI 10.1016/j.enganabound.2014.09.009. [Google Scholar] [CrossRef]

15. Chen, L. C., Li, X. L. (2020). A Burton–Miller boundary element-free method for Helmholtz problems. Applied Mathematical Modelling, 83, 386–399. DOI 10.1016/j.apm.2020.03.009. [Google Scholar] [CrossRef]

16. Qu, W. Z., Fan, C. M., Gu, Y. (2019). Localized method of fundamental solutions for interior Helmholtz problems with high wave number. Engineering Analysis with Boundary Elements, 107, 25–32. DOI 10.1016/j.enganabound.2019.06.018. [Google Scholar] [CrossRef]

17. Britt, S., Petropavlovsky, S., Tsynkov, S., Turkel, E. (2015). Computation of singular solutions to the Helmholtz equation with high order accuracy. Applied Numerical Mathematics, 93, 215–241. DOI 10.1016/j.apnum.2014.10.006. [Google Scholar] [CrossRef]

18. Lin, J., Chen, W., Chen, C. S., Jiang, X. R. (2020). Fast bundary knot method for solving axisymmetric Helmholtz problem with high wave number. Computer Modeling in Engineering & Sciences, 94(6), 485–505. [Google Scholar]

19. Wang, K., Wong, Y. S., Huang, J. Z. (2017). Solving Helmholtz equation at high wave numbers in exterior domains. Applied Mathematics and Computation, 298, 221–235. DOI 10.1016/j.amc.2016.11.015. [Google Scholar] [CrossRef]

20. Li, X. L., Dong, H. Y. (2020). Error analysis of the meshless finite point method. Applied Mathematics and Computation, 382, 125326. DOI 10.1016/j.amc.2020.125326. [Google Scholar] [CrossRef]

21. Qu, W. Z., Fan, C. M., Li, X. L. (2020). Analysis of an augmented moving least squares approximation and the associated localized method of fundamental solutions. Computers and Mathematics with Applications, 80(1), 13–30. DOI 10.1016/j.camwa.2020.02.015. [Google Scholar] [CrossRef]

22. Berrut, J. P., Trefethen, L. N. (2004). Baryccentric lagrange interpolation. SIAM Review, 46(3), 501–517. DOI 10.1137/S0036144502417715. [Google Scholar] [CrossRef]

23. Ma, W. T., Zhang, B. W., Ma, H. L. (2016). A meshless collocation approach with barycentric rational interpolation for two-dimensional hyperbolic telegraph equation. Applied Mathematics and Computation, 279, 236–248. DOI 10.1016/j.amc.2016.01.022. [Google Scholar] [CrossRef]

24. Kumar, S. (2014). A new analytical modeling for fractional telegraph equation via Laplace transform. Applied Mathematical Modelling, 38(13), 3154–3163. DOI 10.1016/j.apm.2013.11.035. [Google Scholar] [CrossRef]

25. Berrut, J. P., Baltensperger, R., Mittelmann, H. D. (2005). Recent developments in barycentric rational interpolation. International Series of Numerical Mathematics, 151, 27–51. [Google Scholar]

26. Fu, Y. P. (2008). Compact fourth-order finite difference schemes for Helmholtz equation with high wave numbers. Journal of Computational Mathematics, 26, 98–111. [Google Scholar]

27. Singer, I., Turkel, E. (2006). Sixth-order accurate finite difference schemes for the Helmholtz equation. Journal of Computational Acoustics, 214(03), 339–351. DOI 10.1142/S0218396X06003050. [Google Scholar] [CrossRef]

28. Wu, T. T., Xu, R. M. (2018). An optimal compact sixth-order finite difference scheme for the Helmholtz equation. Computers and Mathematics with Applications, 75(7), 2520–2537. DOI 10.1016/j.camwa.2017.12.023. [Google Scholar] [CrossRef]

29. Zhuang, Y., Sun, X. H. (2000). A high order ADI method for separable generalized Helmholtz equations. Advances in Engineering Software, 31(8–9), 585–591. DOI 10.1016/S0965-9978(00)00026-0. [Google Scholar] [CrossRef]

30. Liu, C. S., Kuo, C. L. (2016). A multiple-scale Pascal polynomial triangle solving elliptic equations and inverse Cauchy problems. Engineering Analysis with Boundary Elements, 62, 35–43. DOI 10.1016/j.enganabound.2015.09.003. [Google Scholar] [CrossRef]

images This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.