The Improved Element-Free Galerkin Method for 3D Helmholtz Equations

: The improved element-free Galerkin (IEFG) method is proposed in this paper for solving 3D Helmholtz equations. The improved moving least-squares (IMLS) approximation is used to establish the trial function, and the penalty technique is used to enforce the essential boundary conditions. Thus, the ﬁnal discretized equations of the IEFG method for 3D Helmholtz equations can be derived by using the corresponding Galerkin weak form. The inﬂuences of the node distribution, the weight functions, the scale parameters of the inﬂuence domain, and the penalty factors on the computational accuracy of the solutions are analyzed, and the numerical results of three examples show that the proposed method in this paper can not only enhance the computational speed of the element-free Galerkin (EFG) method but also eliminate the phenomenon of the singular matrix.


Introduction
As an important elliptic differential equation, the Helmholtz equation has been widely applied in many different fields, such as mechanics, acoustics, physics, electromagnetics, engineering, and so on. It is well known that how to achieve the numerical solutions of Helmholtz equations effectively and accurately is one of the important directions in the scientific research.
As an important meshless method, the EFG method [16] was studied by Belytschko et al. In this method, a trial function is established by using the moving least-squares (MLS) approximation. Cheng et al. analyzed the error estimates of EFG method for potential problems [17]. Because the MLS approximation is based on the least-squares method [18][19][20][21][22], the disadvantages of the least-squares method also exist in the MLS approximation, in which sometimes ill-conditional or singular matrices occur.
In order to eliminate the singular matrices, the improved moving least-squares (IMLS) approximation [23] was proposed by Cheng et al., in which the orthogonal function system with a weight function is used as basis function, and thus can make up for the deficiency of the MLS approximation and has greater computational efficiency, using the IMLS approximation to establish the trial function. Thus, the improved element-free Galerkin (IEFG) method was applied for potential [24], transient heat conduction problems [25], the wave equation [26], the Schrödinger equation [27], advection-diffusion [28], elastodynamics [29], In Section 4, the influences of the node distribution, the weight functions, the scale parameters, and the penalty factors on the computational accuracy of the solutions are analyzed by giving examples. It is shown that the IEFG method for Helmholtz equations is convergent. Compared with the EFG method, the IEFG method has greater computational speed. Moreover, the singular matrix can be eliminated.

The IMLS Approximation
The approximation of a function u(x) is where p T (x) is the basis function vector, m is the basis function number, and a T (x) = (a 1 (x), a 2 (x), · · · , a m (x)) (2) is the coefficient vector of p T (x). In general, The local approximation is Define where w(x − x I ) is a weighting function, and x I (I = 1, 2, · · · , n) are the nodes with influence domains covering point x. Equation (6) can be written as where u T = (u 1 , u 2 , · · · , u n ), and From we have where Equation (12) sometimes forms a singular or ill-conditional matrix. In order to make up for this deficiency, for basis functions using the Gram-Schmidt process, we can obtain and Then, from Equation (12), a(x) can be obtained as where Substituting Equation (18) into Equation (5), we have where is the shape function. This is the IMLS approximation [23], in which the shape function can be obtained more easily than the MLS approximation. Moreover, the IMLS approximation can also avoid the singular matrix. Thus, it can enhance the computational efficiency of the MLS approximation.

The IEFG Method for 3D Helmholtz Equations
The governing equation is and the boundary conditions are where k 2 is the wave number, f (x) is the given function, u and q are the given values, and is the unit outward normal to the boundary Γ in direction x i . For 3D Helmholtz equations, the equivalent functional is By introducing the penalty technique to apply the boundary conditions, we can obtain the modified functional where α is the penalty factor. Let We can obtain the following equivalent integral weak form where In the cubic domain Ω, we employ M nodes x I (I = 1, 2, · · · , M). Thus, we have From the IMLS approximation, we can obtain where u = (u 1 , u 2 , · · · , u n ) T .
From Equations (29) and (31), we have where Substituting Equations (31) and (33) into Equation (28), we have In Equation (36), the form of u is the same as Equation (32), and n = M. All integral terms in Equation (36) are analyzed as follows: Let Substituting Equations (37)- (42) into Equation (36), we can obtain The δu T is arbitrary; thus we can obtain where This is the IEFG method for 3D Helmholtz equations.

Numerical Examples
The formula of the relative error is where In order to illustrate the advantages of the IFFG method, we chose three examples from other literature. The nodes distributed in the problem domains of these numerical examples were regular, the linear basis function was selected, and 3 × 3 × 3 Gaussian points were selected in each integral cell. The IEFG and the EFG methods are used to solve these examples.
The following equation is considered in the first example: The boundary conditions are The problem domain is is the analytical solution.
In order to study the convergence of the EFG and the IEFG methods for Helmholtz equations, all parameters of both methods were kept the same. The cubic spline weight functions were used, d max = 1.35, α = 2.0 × 10 4 . Table 1 shows the relationship between relative errors and node distribution. It is shown that, with the increase in nodes, the precision of numerical solutions improves as well, but the computational efficiency is reduced gradually. Therefore, the two methods in this paper are convergent. Both the computational accuracy and efficiency are considered, and 15 × 15 × 15 regularly distributed nodes are selected. The effects of the weight function, the scale parameter of the influence domain, and the penalty factor on solution of the IEFG method will be discussed.
(1) Weight function When the cubic spline function is used, 15 × 15 × 15 regularly distributed nodes and 14 × 14 × 14 background integral cells are selected, α = 2.0 × 10 4 , d max = 1.35. Thus, the smaller relative error is 0.6296%. When the quartic spline function is used, and the same regularly distributed nodes and background integral grids are used, α = 2.2 × 10 4 , d max = 1.28, the smaller relative error is 0.6274%. It is shown that the similar relative errors can be obtained when using two weight functions.
In addition, the singular matrix can be avoided in the IEFG method when using the cubic spline function. If d max = 1.0, the quartic spline function is selected. Unfortunately, the singular matrix occurs and the final result cannot be obtained. When the cubic spline function is used, the relative error is 0.6451%.
Thus, the cubic spline function is selected.
(2) Scale parameter The same node distribution and background integral grids are selected, α = 2.0 × 10 4 , and the cubic spline function is used. Figure 1 shows the relationship between d max and relative errors. Because of the error of computer itself, the relative error become larger when d max = 1.2. It is shown that when d max = 1.35, the relative error is smaller.
Mathematics 2021, 9, x FOR PEER REVIEW 9 of 22 1.28, the smaller relative error is 0.6274%. It is shown that the similar relative errors can be obtained when using two weight functions. In addition, the singular matrix can be avoided in the IEFG method when using the cubic spline function. If dmax = 1.0, the quartic spline function is selected. Unfortunately, the singular matrix occurs and the final result cannot be obtained. When the cubic spline function is used, the relative error is 0.6451%.
Thus, the cubic spline function is selected.
(2) Scale parameter The same node distribution and background integral grids are selected, α = 2.0 × 10 , and the cubic spline function is used. Figure 1 shows the relationship between dmax and relative errors. Because of the error of computer itself, the relative error become larger when d max = 1.2. It is shown that when dmax = 1.35, the relative error is smaller. (3) Penalty factor The same node distribution, background integral grids, and weight function are selected, dmax = 1.35. Figure 2 shows the relationship between α and relative errors. It is shown that when α = 2.0 × 10 , the relative error is smaller.  (3) Penalty factor The same node distribution, background integral grids, and weight function are selected, d max = 1.35. Figure 2 shows the relationship between α and relative errors. It is shown that when α = 2.0 × 10 4 , the relative error is smaller.
1.28, the smaller relative error is 0.6274%. It is shown that the similar relative errors can be obtained when using two weight functions.
In addition, the singular matrix can be avoided in the IEFG method when using th cubic spline function. If dmax = 1.0, the quartic spline function is selected. Unfortunately the singular matrix occurs and the final result cannot be obtained. When the cubic splin function is used, the relative error is 0.6451%.
Thus, the cubic spline function is selected.
(2) Scale parameter The same node distribution and background integral grids are selected, α = 2.0 × 10 , and the cubic spline function is used. Figure 1 shows the relationship between dmax and relative errors. Because of the error of computer itself, the relative error becom larger when d max = 1.2. It is shown that when dmax = 1.35, the relative error is smaller. (3) Penalty factor The same node distribution, background integral grids, and weight function are se lected, dmax = 1.35. Figure 2 shows the relationship between α and relative errors. It i shown that when α = 2.0 × 10 , the relative error is smaller.  The IEFG method is selected to solve it, 15 × 15 × 15 regularly distributed nodes and 14 × 14 × 14 background integral cells are selected, and the cubic spline function is used, α = 2.0 × 10 4 , d max = 1.35. When using the EFG method to solve it, the same parameters are selected, and thus the relative errors of two methods are equal to 0.6296%.  show the comparison between numerical solutions and analytical ones, and the CPU times of the IEFG method and the EFG method are 92.1 s and 98.0 s, respectively. Obviously, higher computational efficiency can be obtained when using the IEFG method.
Mathematics 2021, 9, x FOR PEER REVIEW 10 of 2 The IEFG method is selected to solve it, 15 × 15 × 15 regularly distributed nodes and 14 × 14 × 14 background integral cells are selected, and the cubic spline function is used α = 2.0 × 10 , dmax = 1.35. When using the EFG method to solve it, the same parameters are selected, and thus the relative errors of two methods are equal to 0.6296%.  show the comparison between numerical solutions and analytical ones and the CPU times of the IEFG method and the EFG method are 92.1 s and 98.0 s, respec tively. Obviously, higher computational efficiency can be obtained when using the IEFG method.  Additionally, the singular matrix can be avoided when constructing the shape func tions when the IEFG is used. If dmax = 1.0 and other parameters are the same, two methods are used to solve it, and two different results are obtained. When the EFG method is used the singular matrix occurs and the final result cannot be obtained. However, using the The IEFG method is selected to solve it, 15 × 15 × 15 regularly distributed nodes and 14 × 14 × 14 background integral cells are selected, and the cubic spline function is used α = 2.0 × 10 , dmax = 1.35. When using the EFG method to solve it, the same parameters are selected, and thus the relative errors of two methods are equal to 0.6296%. Figures 3-5 show the comparison between numerical solutions and analytical ones and the CPU times of the IEFG method and the EFG method are 92.1 s and 98.0 s, respectively. Obviously, higher computational efficiency can be obtained when using the IEFG method. Additionally, the singular matrix can be avoided when constructing the shape functions when the IEFG is used. If dmax = 1.0 and other parameters are the same, two methods are used to solve it, and two different results are obtained. When the EFG method is used the singular matrix occurs and the final result cannot be obtained. However, using the IEFG method to solve it, the relative error of the numerical solutions is 0.6451%. The nu merical and analytical results are compared in Figure 6; it is shown that the numerica results are in good agreement with the analytical ones. The second example [67] is The boundary conditions are Additionally, the singular matrix can be avoided when constructing the shape functions when the IEFG is used. If d max = 1.0 and other parameters are the same, two methods are used to solve it, and two different results are obtained. When the EFG method is used, the singular matrix occurs and the final result cannot be obtained. However, using the IEFG method to solve it, the relative error of the numerical solutions is 0.6451%. The numerical and analytical results are compared in Figure 6; it is shown that the numerical results are in good agreement with the analytical ones. IEFG method to solve it, the relative error of the numerical solutions is 0.6451%. The numerical and analytical results are compared in Figure 6; it is shown that the numerical results are in good agreement with the analytical ones. The second example [67] is The boundary conditions are ) ( The second example [67] is ∆u − k 2 u = 0.
The boundary conditions are u(x 1 , 1, The problem domain is is the analytical solution.
We set k = 2, ξ 1 = 1, and ξ 2 = 0.5. The IEFG method is used to solve it, α = 1.7 × 10 3 , d max = 1.21. The 15 × 15 × 15 regularly distributed nodes and 14 × 14 × 14 background integral grids are used. When using the EFG method to solve it, the same parameters are selected, and thus the same computational accuracy can be obtained. The relative errors of both methods are equal to 0.0844%. Figures 7-9 show the comparison of the numerical solutions of the two methods and the analytical ones. The CPU times of the IEFG method and the EFG method are 92.1 s and 98.0 s, respectively. We can see that the computational results of both methods are in very good agreement with the analytical ones.
Mathematics 2021, 9, x FOR PEER REVIEW 12 of 2 The problem domain is is the analytical solution.
We set k = 2, 1 ξ = 1, and 2 ξ = 0.5. The IEFG method is used to solve it, α = 1.7 × 10 dmax = 1.21. The 15×15×15 regularly distributed nodes and 14×14×14 background integra grids are used. When using the EFG method to solve it, the same parameters are selected and thus the same computational accuracy can be obtained. The relative errors of both methods are equal to 0.0844%. Figures 7-9 show the comparison of the numerical solu tions of the two methods and the analytical ones. The CPU times of the IEFG method and the EFG method are 92.1 s and 98.0 s, respectively. We can see that the computationa results of both methods are in very good agreement with the analytical ones.  When different parameters are selected, k = 5, 1 ξ = 3, and 2 ξ = 2.7. Using two meth ods to solve it, the same parameters are used. Thus, the relative errors of both method are equal to 0.5295%. Figures 10-12 show the comparison of the numerical solutions of th two methods and the analytical ones, and the CPU times of the IEFG method and the EFG method are 92.1 s and 98.0 s, respectively. We can see that the computational results o both methods are in good agreement with the analytical ones. When different parameters are selected, k = 5, 1 ξ = 3, and 2 ξ = 2.7. Using two meth ods to solve it, the same parameters are used. Thus, the relative errors of both methods are equal to 0.5295%. Figures 10-12 show the comparison of the numerical solutions of the two methods and the analytical ones, and the CPU times of the IEFG method and the EFG method are 92.1 s and 98.0 s, respectively. We can see that the computational results o both methods are in good agreement with the analytical ones. When different parameters are selected, k = 5, ξ 1 = 3, and ξ 2 = 2.7. Using two methods to solve it, the same parameters are used. Thus, the relative errors of both methods are equal to 0.5295%. Figures 10-12 show the comparison of the numerical solutions of the two methods and the analytical ones, and the CPU times of the IEFG method and the EFG method are 92.1 s and 98.0 s, respectively. We can see that the computational results of both methods are in good agreement with the analytical ones.
We can select k = 10, ξ 1 = 5.8, and ξ 2 = 6.2. Using the two methods to solve it, the same parameters are used, and the relative errors of both methods are equal to 2.3884%. Figures 13-15 show the comparison of the numerical solutions of the two methods and the analytical ones. The computational results of both methods are in good agreement with the analytical ones.   We can select k = 10, 1 ξ = 5.8, and 2 ξ = 6.2. Using the two methods to solve it, the same parameters are used, and the relative errors of both methods are equal to 2.3884%. Figures 13-15 show the comparison of the numerical solutions of the two methods and the analytical ones. The computational results of both methods are in good agreement with the analytical ones.

x3-axis.
We can select k = 10, 1 ξ = 5.8, and 2 ξ = 6.2. Using the two methods to solve it, the same parameters are used, and the relative errors of both methods are equal to 2.3884%. Figures 13-15 show the comparison of the numerical solutions of the two methods and the analytical ones. The computational results of both methods are in good agreement with the analytical ones.  From this example, we can draw two conclusions: On the one hand, the IEFG method has greater computational efficiency; on the other hand, the bigger the wave numbers are, the lower the computational accuracy.
Similarly, if dmax = 1.0, we select k = 10, 1 ξ = 5.8, and 2 ξ = 6.2. When the EFG method From this example, we can draw two conclusions: On the one hand, the IEFG method has greater computational efficiency; on the other hand, the bigger the wave numbers are, the lower the computational accuracy.
Similarly, if dmax = 1.0, we select k = 10, 1 ξ = 5.8, and 2 ξ = 6.2. When the EFG method From this example, we can draw two conclusions: On the one hand, the IEFG method has greater computational efficiency; on the other hand, the bigger the wave numbers are, the lower the computational accuracy.
Similarly, if d max = 1.0, we select k = 10, ξ 1 = 5.8, and ξ 2 = 6.2. When the EFG method is used, unfortunately, the singular matrix occurs. When the IEFG method is used, the relative error is 2.4229%. The numerical solutions and analytical ones are compared in Figure 16. It is shown that the numerical results are in good agreement with the analytical ones. The third example [68] is ) π sin( ) π sin( ) π cos( ) π 3 ( The boundary conditions are The problem domain is is the analytical solution. The IEFG method is used to solved it. The wave number is selected as 100, and 19 × 19 × 19 regularly distributed nodes and 18 × 18 × 18 background integral cells are used, α = 1.9 × 10 , dmax = 1.1. When using the EFG method to solve it, the same parameters are selected, and the relative errors of both methods are equal to 0.8646%. Figures 17-19 show the comparison of the numerical solutions and the analytical ones. We can see that numerical solutions are in good agreement with the analytical ones. The CPU times of the IEFG method and the EFG method are 200.6 s and 208.1 s, respectively. The third example [68] is ∆u + k 2 u = (k 2 − 3π 2 ) cos(πx 1 ) sin(πx 2 ) sin(πx 3 ).
The boundary conditions are u(x 1 , 0, The problem domain is 1], and is the analytical solution.
The IEFG method is used to solved it. The wave number is selected as 100, and 19 × 19 × 19 regularly distributed nodes and 18 × 18 × 18 background integral cells are used, α = 1.9 × 10 7 , d max = 1.1. When using the EFG method to solve it, the same parameters are selected, and the relative errors of both methods are equal to 0.8646%. Figures 17-19 show the comparison of the numerical solutions and the analytical ones. We can see that numerical solutions are in good agreement with the analytical ones. The CPU times of the IEFG method and the EFG method are 200.6 s and 208.1 s, respectively.   A similar computational accuracy can be obtained when using the two methods, but the higher computational speed can be obtained when using the IEFG method.
Similarly, if dmax = 1.0, when the EFG method is used, the singular matrix occurs and the final result cannot be obtained. However, when the IEFG method is selected, the relative error is 0.8648%. The numerical solutions and the analytical one are compared in Figure 20, where it is shown that the numerical results are in good agreement with the ana- A similar computational accuracy can be obtained when using the two methods, but the higher computational speed can be obtained when using the IEFG method.
Similarly, if d max = 1.0, when the EFG method is used, the singular matrix occurs and the final result cannot be obtained. However, when the IEFG method is selected, the relative error is 0.8648%. The numerical solutions and the analytical one are compared in Figure 20, where it is shown that the numerical results are in good agreement with the analytical ones. A similar computational accuracy can be obtained when using the two methods, but the higher computational speed can be obtained when using the IEFG method.
Similarly, if dmax = 1.0, when the EFG method is used, the singular matrix occurs and the final result cannot be obtained. However, when the IEFG method is selected, the relative error is 0.8648%. The numerical solutions and the analytical one are compared in Figure 20, where it is shown that the numerical results are in good agreement with the analytical ones.

Conclusions
In order to solve 3D Helmholtz equations efficiently, the IEFG method is proposed in this paper.
Some numerical examples are given in Section 4, and the convergence of the IEFG method is proven numerically. From these examples, we can see that the IEFG method in this paper can not only enhance the computational speed of the traditional EFG method, but also eliminate the phenomenon of the singular matrix.

Conclusions
In order to solve 3D Helmholtz equations efficiently, the IEFG method is proposed in this paper.
Some numerical examples are given in Section 4, and the convergence of the IEFG method is proven numerically. From these examples, we can see that the IEFG method in this paper can not only enhance the computational speed of the traditional EFG method, but also eliminate the phenomenon of the singular matrix.