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

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:

where Ω=[a,b]×[c,d] is the problem region, u(x,y) is an unknown function of x and y.k(x,y)is the wave number, which could be either a constant or a function of x and y.f (x,y) is a source term.The Dirichlet boundary conditions are given as:

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.[2-6].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 [8-11].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

Letxi,i=1,2,...,m, be a given interpolation node, and the corresponding function value isfi, then Lagrange polynomial interpolation formula can be expressed as

whereLi(x)is Lagrange interpolation basis function in the form of

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

Defines weight of the barycentric Lagrange interpolation function as

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

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

Figure 1:The algorithm flow diagram of ωi

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

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

By using Eq.(2.6) to interpolate the constant 1 with(xi,1),i=1,2,...,m, then we have

After simplification, it can be obtained

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

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 O(m).

2.2 Bivariate Barycenter Lagrange Interpolation Formula

Let theu(x,y)be a function about the variablexandy, and define the interval as(x,y)∈[a,b]×[c,d].When the interval [a,b] and [c,d] is discretized intomandnnodes, respectively, then there arem×ntensor product interpolation nodes in the whole domain1,2,...,n.Then the barycentric Lagrange interpolation formula [23] of theu(x,y)is

whereuij=u(xi,yj).Li(x)andMj(y)are the interpolation basis functions in thexandydirection,respectively.

where the weight is

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:

Chebyshev nodes of the second kind:

What needs to be explained is that for two types of Chebyshev node formulas, the defined interval is [-1,1].As for the general interval [a,b], the coordinate transformation formula ¯x=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.

Whenmnodes are inserted on the interval [a,b], there isa=x1＜x2＜···＜xm=b, to set the function value of the unknown functionu(x)at the nodexiisui=u(xi),i=1,2,...m, then the barycentric Lagrange interpolation formula of the functionu(x)is

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

Using Eq.(2.14), thekderivative of the functionu(x)at the nodexj(j=1,2,...,m)can be expressed as

In the form of a matrix as

whereandu=[u1,u2,...,um]Tare the vector.The matrixD(k)is as thekorder differential matrix of the functionu(x), andIn addition, we can calculateD(k)(see [25]) with

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 matrixand matrixisA⊗B, recorded as following

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 functionu(x,y)inxandydirections of the twodimensional Helmholtz equation aremandnthen by using Eq.(2.17), the expressions of the differential matrices corresponding to the partial derivative of the differential equation about the variablesxandycan be briefly expressed as

whereImandInrepresent the unit matrices ofmorder andnorder, respectively, and theD(2)andC(2)are barycentric interpolates of the second-order matrices at the nodex1,x2,...,xmandy1,y2,...,yn, 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

If Eq.(3.1) holds at the nodes of(xI,yJ),I=1,2,...,m;J=1,2,...,n, then Eq.(3.1) is replaced by the following

Let

Then Eq.(3.3) can be written into

That is

In which

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)

whererepresent the first andmlines of the first order differential matrix of themorder unit matrix,andrepresent the first andnlines of the first order differential matrix of thenorder unit matrix, andg1=[g1(y1),g1(y2),...,g1(yn)]T,g2=[g2(y1),g2(y2),...,g2(yn)]T,g3=[g3(x1),g3(x2),...,g3(xm)]T,g4=[g4(x1),g4(x2),...,g4(xm)]T.

Combining Eqs.(3.5) and (3.6), we can obtain the value of the functionu(x,y)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 (E∞) and relative error (Er) are defined as follows:

whereuijandueijrepresent numerical and exact solutions, respectively.NtandMtare the numbers of test nodes in thexandydirections, respectively.

Example 1[26]:

Consider the following non-homogeneous Helmholtz equation

uxx+uyy+k2u=sin(πx)sin(kπy),(x,y)∈[0,1]×[0,1]

with the boundary conditions:

The exact solution is

Tab.1 gives the calculation results ofE∞for the four kinds of interpolation nodes withNt=Mt=N,k=5.In which,Nis the number of interpolation nodes in the domain,NtandMtare 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 E∞with four kinds of interpolation nodes when Nt=Mt=N,k=5 for Example 1

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

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

Then theE∞with two kinds of Chebyshev interpolation nodes and three kinds of test nodes whenNt=Mt=100,k=5 andN=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 asNtandMtchange, which means that the number of test nodes has little influence on our numerical results.So we can take the value ofNtandMtwith 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 ofE∞with the 4th-order method in [26] and with the present scheme are compared whenNt=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 reachesand 12th-order,respectively when the number of collocation points isN=128, while that in the literature just reachesand 4th-order, respectively.It shows that the calculation accuracy of the present scheme is much better than that in the literature.

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

The node distribution for the four kinds of interpolation nodes withNt=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.

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 withNt=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.

Finally, the convergence order withNt=Mt=100,k=5 andNt=130,Mt=140,k=30 is shown in Fig.4.It can be seen that with the wave number ofkincreases, 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.

Example 2[27]:

Consider the following non-homogeneous Helmholtz equation

with the boundary conditions:

The exact solution is

u(x,y)=sinh(αx)sinh[α(π-x)]+sinh(βy)sinh[β(π-y)]

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

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

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

Tab.5 gives theE∞computed by the present scheme and compares with EB-6 method in [27] whenα=0.5,β=0.7,Nt=Mt=120,k=6.4.The results show that the error order of the proposed scheme reacheswhen the number of collocation points isN=16, while that in the literature, the error order only reacheswhenN=16 andwhenN=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 nodesN, 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 E∞when α=0.5,β=0.7,Nt=Mt=120,k=6.4 for Example 2

TheE∞,Erand CPU time withα=0.6,β=0.8,Nt=Mt=50,k=10 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 E∞,Er and CPU time when α=0.6,β=0.8,Nt=Mt=50,k=10 for Example 2

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α=0.6,β=0.8,Nt=Mt=50,k=10,N=16.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

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

Example 3[28]:

Consider the following non-homogeneous high wavenumber Helmholtz equation

uxx+uyy+k2u=-a2sin(ax)sin(ky),(x,y)∈[0,1]×[0,1]

with the boundary conditions:

u(0,y)=0,u(1,y)=sin(a)sin(ky)

u(x,0)=0,u(x,1)=sin(ax)sin(k)

The exact solution is

u(x,y)=sin(ax)sin(ky)

The calculation results with the 6th-order method in [28] and the present scheme are compared under different wave numbers withk=24 andk=48 whena=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 theE∞,Erand CPU time witha=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.

Fig.6.gives the exact (a) and numerical (b) solutions and error (c) with different test nodes:random (I), uniform (II) and Chebyshev (III) witha= 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.

Example 4[29]:

Consider the following non-homogeneous variable wavenumber Helmholtz equation

with the boundary conditions:

u(0,y)=sin(3y),u(1,y)=sin(2+3y)

u(x,0)=sin(2x),u(x,1)=sin(2x+3)

The exact solution is

u(x,y)=sin(2x+3y)

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

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

Tab.9 shows theE∞whenNt= 150,Mt= 130 with variable wavenumberk(x,y)=Compared with the MHADI method in [29], the error order of the present scheme iswithN=16, while that in the literature is onlyIn 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 theE∞,Erat different collocation points and CPU time withNt=Mt=50,k(x,y)=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.

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

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

Table 9:The E∞when Nt=150,Mt=130,k(x,y)=for Example 4

Table 9:The E∞when Nt=150,Mt=130,k(x,y)=for Example 4

N MHADI [29] Present scheme Random Uniform Chebyshev 16 1.70(-2) 7.064(-12) 8.911(-12) 9.255(-12)32 1.30(-3) 1.155(-10) 1.990(-10) 1.968(-10)64 7.10(-5) 3.513(-9) 4.694(-9) 6.344(-9)128 5.10(-6) 9.268(-8) 1.300(-7) 1.467(-7)

Table 10:The E∞,Er and CPU time when Nt=Mt=50,k(x,y)= for Example 4

Table 10:The E∞,Er and CPU time when Nt=Mt=50,k(x,y)= for Example 4

N E∞ Er CPU time Random Uniform Chebyshev Random Uniform Chebyshev 8 8.942(-7) 8.948(-7) 8.930(-7) 5.651(-7) 5.478(-8) 4.514(-7) 0.319 12 1.584(-11) 1.594(-11) 1.603(-11) 1.139(-11) 9.730(-12) 8.066(-12) 0.304 16 4.761(-12) 8.564(-12) 8.488(-12) 5.709(-13) 7.908(-13) 1.405(-12) 0.427 24 2.562(-11) 8.154(-11) 8.213(-11) 3.927(-12) 5.838(-12) 9.301(-12) 0.639

a) and numerical (b) solutions and error (c) with different test nodes:ran-dom (I), uniform (II) and Chebyshev (III) when Nt = Mt = 50,k(x,y) =,N = 16for Example 4

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

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

Example 5[30]:

Consider the following homogeneous Helmholtz equation

uxx+uyy+k2u=0,(x,y)∈Ω

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

ρ(θ)=exp(sin(θ))sin2(2θ)+exp(cos(θ))cos2(2θ),

x=ρ(θ)cos(θ),y=ρ(θ)sin(θ),0 ≤θ≤2π

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

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

Table 14:The E∞,Er and CPU time when Nt=Mt=50 for Example 6

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

The exact solution is

Tab.11 and Tab.12 give theE∞,Erand CPU time with an ameba-shape domain whenNt=Mt=50,k=5 andNt=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.

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

Then the Fig.8.gives the computational domain with different test nodes:random (a),uniform (b) and Chebyshev (c) whenNt=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)withNt=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.

Example 6[30]:

Consider the following non-homogeneous modified Helmholtz equation

with the boundary conditions:

The exact solution is

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 theE∞,Erand CPU time with different collocation points underNt=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) whenNt=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.

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 nodesNtandMt, 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.