Abstract
This paper presents an efficient weighted Laguerre polynomials based meshless finite-difference time domain (WLP-MFDTD). By decomposing the coefficients of the system matrix and adding a perturbation term, a factorization-splitting scheme is introduced. The huge sparse matrix is transformed into two matrices with 9 unknown elements in each row regardless of the duplicated ones. Consequently, compared with the conventional implementation, the CPU time and memory requirement can be saved greatly. The perfectly matched layer absorbing boundary condition is also extended to this approach. A numerical example demonstrates the capability and efficiency of the proposed method.
1. Introduction
Meshless methods first presented for vibration analysis have been applied to electromagnetic numerical computation due to their characteristics of adaptive and accurate description of nonuniform geometries and multiscale solutions [1]. Instead of a mesh to discretize the problem domain such as FDTD, FEM, and MOM, a set of random scattered points are constructed to approximate the field component. Consequently, no reconstruction of mesh lines is required when the problem domain is modified partially.
Among the existing electromagnetic meshless approaches, for example, smooth particle hydrodynamics method, least-squares-based methods, and radial basis functions, the radial point interpolation method (RPIM) emerged as an alternative one [1, 2]. This method is used for discretization of local space and is based on Gaussian and multiquadric radial functions [3, 4]. It has smoother approximation for derivative and unique approximation to the field component. The shape function possesses delta function property and is capable of modeling arbitrary boundaries.
However, the conventional time-domain RPIM suffers from late-time instability mainly because of the numerical truncation errors of the Gaussian basis function at the boundary. Optimal shape parameter, global basis function, and alternative direction implicit (ADI) FDTD are recommended to take a balance between minimum numerical dispersion errors and high accuracy [5, 6]. Besides, the weighted decaying Laguerre polynomial (WLP) is developed to incorporate with RPIM called Laguerre-RPIM [7, 8]. The WLP based time-domain analysis was originally proposed by Chung et al. [9] and extended by Yi et al. [10] and Zhao and Shen [11]. The unconditional stability, numerical accuracy, and efficiency have been verified.
Although the system matrix mentioned above is sparse and banded, the lower-upper (LU) factorization is very expensive [12, 13]. Generally, Voronoi decomposition or Delaunay tessellation is utilized to allocate the field nodes with irregular distribution [4]. For the former 2D case implemented directly by Matlab function Voronoi, each node is surrounded by four nodes in the interior of the problem domain. As a result, the final matrix (, where is the number of total nodes) has 18 unknown elements in each row despite duplicated nodes. The scheme is hardly usable for practical large-scale problems. More specifically, for complex or fine structure where additional local nodes are required to describe the dimension accurately, the unknown elements will increase dramatically. Therefore, it is desirable to introduce an efficient algorithm for implementing the Laguerre-RPIM.
In this paper, a factorization-splitting scheme is applied by decomposing the coefficients of the system matrix and adding a perturbation term. The proposed method solves two matrices with 9 unknown elements in each row regardless of the duplicated ones. Consequently, compared with the conventional implementation, the CPU time and memory requirement can be saved greatly. Besides, a new uniform nodal placement is proposed: each node is surrounded by four nodes and each node is surrounded by four nodes. The singularity of the polynomial basis function in the support domain is eliminated [14]. The implementation of PML absorbing boundary condition is also discussed.
2. Mathematical Formulation
2.1. Laguerre-RPIM
The Laguerre polynomials of order are defined as
Using the orthogonality of the above polynomials with respect to the weighted function , an orthogonal set of basis functions can be formed:where is a time-scale factor. The basis functions are absolutely convergent to zero as . Then, the time-domain field component is expanded as
The first derivative with respect to of field component can be found as
Taking as an example, assuming a linear, isotropic, nondispersive, and lossless medium, the differential Maxwell equation can be written as
By inserting (4) into (5) and using a temporal Galerkin testing procedure to eliminate the time variable, one can getThe other field qualities can be obtained in a similar way.
We assume that the function is approximated by a set of local random scattered points , where is the number of above points which are in the neighborhood of . RPIM is based on radial basis function and polynomial basis function :where and are the associated expansion coefficients. Gaussian form is adopted as radial basis function in the support domain and is built by utilizing the four-term () Pascal triangle:where shape parameter controls the flatness of RPIM and denotes the maximum distance between interpolation point and nodes involved in the local support domain. The vector form is defined asConstraints are usually imposed to ensure the unique solution:
Equations (7) and (10) can be expressed as matrix form:Coefficient matrices and depend only upon the positions of corresponding points:
Once the unknown and based on (11) are solved, the final interpolation is obtained:where is the shape function related to the basis functions which are only determined by the coordinates of interpolation points. It also can be observed that has no relation to field value . Detailed discussion of the shape function can be found in [1]. By substituting (13) into (6) and expanding magnetic component in local support domain with RPIM, one obtains
The final expression of can be obtained once is expanded in a similar manner. After some mathematical manipulations, a set of matrices are denoted aswhere , , , , , , and
By inserting (16) into (15), the final huge matrix of nodes is
2.2. High Efficient Laguerre-RPIM
As discussed in the previous section, the coefficient matrix is determined by medium parameter (, ) and the position of the interpolation points (shape functions). It is independent on the order . Equation (18) is usually computed in a recursive manner. Although the LU factorization can decompose the banded matrix at the beginning of the simulation, it is computationally expensive even with the unsymmetrical multifunctional sparse LU factorization package which accounts for a large proportion of the whole computational time.
As mentioned previously in the Introduction, the node and node would have the same scale of support domain nodes () for the new uniform nodal placement. The number of nodes reduces to about half of the conventional method [8], whereas the number of nodes remains the same.
The huge matrix has unknown elements which will increase rapidity for large-scale problems or complex structures. Consequently, the conventional Laguerre-RPIM is hardly suitable for practical situation. In this section, a high efficient Laguerre-RPIM is presented to implement the conventional method.
In order to solve the huge system matrix efficiently, we decompose the derivative term () of the left side into two diagonal matrices and , and the equation is rewritten asHere, is a lower diagonal matrix and is an upper diagonal matrix. () is difference operators for the derivatives. It should be noticed that the former derivative is for the shape function , while the latter one is for shape function . In the meantime, a perturbation term is added and the factorized form isWith the splitting scheme, the above equation can be computed into two substeps:here, is a nonphysical intermediate value. It should be pointed out that the first step contains 2 sparse banded matrices, each node is related to 9 nodes after accumulating the duplicated ones, and the total unknown element is . The proposed method greatly reduces the memory of the method and increases the computation efficiency while the accuracy is maintained. However, the second step has two explicit equations which can be computed in a direct way. Due to the unknown values at the right hand of (21), the expansion coefficients of components must be updated following the serial numbers of the above equations. Then, the components are updated according to (16). All the expansion coefficients are solved matching in order . The final difference equations of (21) can be expanded as
2.3. PML Absorbing Boundary Condition for the High Efficient Laguerre-RPIM
The PML absorbing boundary condition while maintaining the high efficient form of the proposed method is implemented in this section. According to the discussion in Sections 2.1 and 2.2, (15)-(16) can be rewritten aswhere
Following the deduced procedure in Section 2.3, (23) can be solved efficiently by factorizing it into two substeps as
3. Numerical Results
In order to validate the formulations of the efficient Laguerre-RPIM FDTD, a numerical example is implemented with the comparison of the conventional RPIM FDTD and Laguerre-RPIM. The performance is studied by simulating the radiation of a point source as illustrated in Figure 1. Only one-fourth of the total nodes are plotted for a clear view. The total numbers of nodes and nodes are 8100 and 7920, respectively. The exciting point and observing point are located at (2.5 mm, 10 mm) and (11.93 mm, 6.09 mm), respectively. The excitation wave is sinusoidal modulated Gaussian pulsewhere GHz, , and . The value of shape function parameter is chosen as 0.05 [2]. Time-scale factor and the order of the WLPs are chosen as and , respectively [6, 13].
The time-domain electric field of the observation point is shown in Figure 2 together with the comparison of other methods. It can be observed that the time-domain results well comply with each other. The associated computational time is 43.06 s for RPIM, 14.26 s for Laguerre-FDTD, and 2.84 s for the proposed method. All the results have demonstrated the accuracy of the proposed method.
4. Conclusion
A new algorithm referred to as efficient Laguerre-FDTD is proposed in this letter. By decomposing the coefficients of the system matrix and adding a perturbation term, the huge sparse matrix is transformed into two matrices with 9 unknown elements in each row regardless of the duplicated ones, instead of one matrix with 18 unknown elements. Compared with the conventional implementation, the CPU time and memory requirement can be saved greatly. A numerical example has proved the accuracy and efficiency of the presented method.
Competing Interests
The authors declare that there are no competing interests regarding the publication of this paper.
Acknowledgments
This work was supported by the Chinese National Science Foundation under Grants 51477182 and 51407199.