2 D Efficient Unconditionally Stable Meshless FDTD Algorithm

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 N × N matrices with 9 unknown elements in each row regardless of the duplicated ones. Consequently, comparedwith the conventional implementation, the CPU time andmemory 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.


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 (2×2, 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 2 International Journal of Antennas and Propagation 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.

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 { 0 ,  1 ,  2 , . ..} can be formed: where   > 0 is a time-scale factor.The basis functions are absolutely convergent to zero as  → ∞.Then, the timedomain field component is expanded as The first derivative with respect to  of field component can be found as Taking   (r, ) 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 get The other field qualities can be obtained in a similar way.We assume that the function (r) is approximated by a set of local random scattered points   (  ) ( = 1, 2, . . ., ), where  is the number of above points which are in the neighborhood of r.RPIM is based on radial basis function   (r) and polynomial basis function   (r): where   and   are the associated expansion coefficients.Gaussian form is adopted as radial basis function in the support domain and   (r) is built by utilizing the four-term ( = 4) Pascal triangle: where shape parameter  controls the flatness of RPIM and  max denotes the maximum distance between interpolation point and nodes involved in the local support domain.The vector form is defined as Constraints are usually imposed to ensure the unique solution: Equations ( 7) and ( 10) can be expressed as matrix form: Coefficient matrices  0 and  0 depend only upon the positions of corresponding points: Once the unknown   and   based on (11) are solved, the final interpolation is obtained: International Journal of Antennas and Propagation 3 where   (r) 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   (r) 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 as where By inserting (16) into (15), the final huge matrix of  nodes is (18)

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 (  =   = 4) 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 2 × 2 × 18 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 as Here,  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    (r), while the latter one   is for  shape function    (r  ).In the meantime, a perturbation term ( ) is added and the factorized form is With the splitting scheme, the above equation can be computed into two substeps: 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 2×××9.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 where Following the deduced procedure in Section 2. (25)

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 onefourth 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 where  0 = 8 GHz,  = 5 0 , and  0 = 0.9.The value of shape function parameter  is chosen as 0.05 [2].Time-scale factor and the order of the WLPs are chosen as   = 1 × 10 11 and   = 120, 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.

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 2 × 2 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.

Figure 1 :Figure 2 :
Figure 1: Node distribution for the simulated structure.