Finite difference method based on polynomial interpolation for solving Helmholtz equation

The generalized finite difference formula plays an important role in the meshless method for solving differential equations. The main purpose of this paper is to find the numerical solution of Helmholtz equation, an elliptic partial differential equation describing electromagnetic waves in physics, by using the generalized finite difference formula. Firstly, this paper introduces a simple and practical nodal distribution, which not only guarantees the uniqueness of multivariate polynomial interpolation, but also makes the matrix triangular, so that the constructed basic polynomials can be transformed into Lagrange basis polynomials. Secondly, the generalized finite difference formula is created by polynomial interpolation. Finally, a numerical example of Helmholtz equation under general boundary conditions is given.


Introduction
In solving the problem of partial differential equation in physical mechanics, if the equation has no analytical solution, we often hope to get the numerical solution of the equation. Finite element method (FEA) is a popular method to solve partial differential equations in mechanics. It is a modern calculation method developed rapidly with the development of computer. However, because the finite element method is based on the numerical algorithm with mesh, when solving problems involving large deformation such as metal forming and stamping process, as well as problems such as dynamic crack propagation, explosion impact, high-speed collision, the method is faced with the problems of mesh re division complexity and mesh distortion, which seriously affect the efficiency and accuracy of calculation. If the mesh remodeling is not good enough in the process of calculation, such as when the mesh is discontinuous or does not meet the requirements of pre-processing, it is difficult to carry on the whole process of mesh based solution and iteration. On this basis, in order to make up for the problems of traditional finite element calculation, the meshless method [1][2][3] is proposed in the academic circle.
Mesh less method is a kind of solution technology which does not consider the mesh division and only uses the node distribution to realize the modeling and Simulation of physical problems. This method not only gets rid of the constraints of mesh initialization and mesh reconstruction, but also ensures the accuracy of the solution [4]. At present, there are more than ten mature meshless methods, including smoothed particle hydrodynamics (SPH), meshless Galerkin method (EFGM), meshless local Galerkin method (mlpgm), diffusion element method (DEM), HP clouds meshless method. Because this method does not need to mesh and is flexible, it is widely used in computational mechanics, large deformation problems and emerging nanoscale multi-scale problems, cell infiltration, blood flow, biological microelectronic systems and other fields.
Helmholtz equation is an elliptic partial differential equation describing electromagnetic waves, named after Helmholtz, a German physicist. Because of its relationship with wave equation, Helmholtz 2 equation appears in the fields of electromagnetic radiation, seismology and acoustics in physics. In this paper, the generalized finite difference method, a meshless method proposed by Huang Xiaowei and Wu Chuansheng [5], is used to solve the numerical solution of Helmholtz equation under the general boundary condition, that is, the Dirichlet problem. This method is based on local polynomial interpolation method, and requires local nodes to have a special distribution proposed by Gasca and maeztu, so as to maintain the existence and uniqueness of Multivariable Interpolation [6][7][8]. On this basis, the Newton interpolation basis function is constructed directly, and the generalized finite difference formula with simple and stable algorithm is established. This method overcomes the singularity problem that may exist in the generalized finite difference method, and the determination of the finite difference formula only needs to solve the linear equations whose coefficient matrix is a triangular matrix, which reduces the amount of calculation to produce the generalized finite difference formula.
In this paper, we first introduce a simple method to determine the node distribution and interpolation node set, then derive the approximation function at the point ( ∈ , is a given bounded region of the problem under study), and finally give a numerical example of solving Helmholtz equation with the generalized finite difference method under general boundary conditions.

Interpolation node and interpolation node set
The distributed generation of nodes is very important to the numerical solution of partial differential equations by meshless method. We hope that the node distribution algorithm is simple enough and can satisfy the uniqueness and Solvability of multivariate polynomial interpolation. Then, according to the selected nodes, the generalized finite difference formula is constructed by polynomial interpolation. Finally, the finite difference formula is introduced into the partial differential equation and the numerical solution of the problem is obtained.
According to the theory of polynomial interpolation, if the generated nodes are on a series of parallel planes and the points on the plane are on a series of parallel lines, the interpolation nodes of n-ary polynomials can be easily selected. Based on this principle, we propose a two-dimensional linear node distribution model.

Distribution of linear nodes in two-dimensional region
For the two-dimensional bounded connected region ⊂ , we select the linear nodes parallel to the y-axis. In order to obtain a group of parallel linear nodes, we use to represent the discretization level of the planar region. The generation process of node distribution is as follows: First, the projection region , of region on the x-axis is defined as follows: Then, the number of columns m of the linear node is determined by the discretization level of the region.
Then a group of equidistant dividing points and a series of lines parallel to the y-axis are generated on the projection area , of the x-axis.
We subdivide the intersection ∩ of line and region into several straight line segments, and we make the following formula true. The column number of linear node is determined by discretization level of region.
Similarly, equidistant nodes are generated on these straight line segments.
In MATLAB, we choose parameter to get the node distribution of circular area as shown in Figure 1. The rotation of x-axis and y-axis in the algorithm can produce linear nodes parallel to y-axis as shown in Figure 2. In this paper, the two-dimensional discrete nodes are composed of linear nodes parallel to x-axis.

Interpolation node set
After obtaining the linear node set in the region parallel to the y-axis by the above method, this section will use the radial basis function interpolation method to select the interpolation node set near the node , so that the polynomial interpolation problem on ∏ has a unique feasible solution with respect to the interpolation node set . The generalized finite difference formula approximates the derivative of the calculation point by the linear combination of the function values of the adjacent nodes near the calculation point, so it is necessary to select the set of adjacent nodes near the calculation point to establish the generalized finite difference formula. The radial basis function interpolation method used in this section first finds the node closest to the calculation point , and then selects the interpolation node set near the node according to the given number of nodes n. Firstly, in the generated node set ⋃ , | 1, . . . , ; 1, . . . , , the 1-points nearest to node are determined to form the point set J, and then the node set G is transformed into a node set , , . . . , ,and we make the following formula true.
, | 1, . . . , ; 1, . . . , , Then, we apply the greedy strategy. Firstly, we define index set 1, . . . , , and the point set closest to point is selected as in , , . . . , . the point set , , . . . , is composed of n points which are different from and closest to point . And so on, until all nodes are selected on 1 different grid lines, the interpolation node set ⋃ is constructed, that is, through the following rule iteration: For any ∈ , select 1 points in which are closest to point to form point set .
(11) Finally, the interpolation node set is obtained as follows: (12) We implement it in MATLAB. Figure 2 shows the interpolation nodes of a given point in the circle.

Generalized finite difference formula
Based on the plane line type node distribution and interpolation node set selection method introduced in the above section, the Newton polynomial interpolation basis function and Lagrange polynomial interpolation [9][10][11] basis function are derived, and the generalized finite difference formula is finally derived. In this paper, we only consider the two-dimensional case, i.e. given region ∈ .
First, according to the method of the upper section, the plane linear node distribution ∈ is constructed. For any given point ∈ , the set of two-dimensional interpolation nodes is obtained： According to the conclusion of Gasca et al. [4,5], for the given set of interpolation nodes, the Newton polynomial interpolation basis functions of two-dimensional polynomial space are constructed as follows： ϕ x, y μ x λ y , α i, j ∈ I .
In this formula: It is easy to verify that satisfies the following equation: In this formula, ， ∈ . Moreover, the two-dimensional Newton polynomial interpolation basis function derived has the same properties as the one-dimensional Newton polynomial interpolation basis function, and the two-dimensional polynomial interpolated on is unique for . Similarly, for the interpolation node set , the Lagrange polynomial interpolation basis function is constructed so that the function is 1 at a given interpolation node and 0 at other interpolation nodes. The expression is as follows: ∈ and ∈ form column vectors and respectively according to the dictionary order of subscripts. The values of column vector at each interpolation node are arranged into a square matrix ∈ according to the dictionary order of node labels [12,13]. From (2)(3)(4)(5) formula, it is easy to prove that the matrix is a triangular matrix of unit, and the determinant value of the matrix is 1, that is, the matrix is nonsingular. According to the uniqueness of interpolation polynomial, Newton polynomial interpolation basis function can be expressed by Lagrange polynomial interpolation basis function in the following matrix form: Then the Lagrange polynomial interpolation basis function can be obtained According to the uniqueness of the polynomial interpolation problem, we can prove that the Lagrange polynomial interpolation basis function ∈ of our given interpolation node set ∈ satisfies the following conditions for ∀ ， ∈ , | | .
q p ∈ D l p γ!, γ γ 0, γ γ In the formula is the -th derivative of . Finally, the generalized finite difference formula is derived, that is, given : → , the Lagrange interpolation polynomial of given interpolation node set is: f q ∈ l p , q ∈ R .