# 2D Efficient Unconditionally Stable Meshless FDTD Algorithm

Authors: Kang Luo; Yun Yi; Yantao Duan; Boao Xu; Bin Chen
Journal: International Journal of Antennas and Propagation
 (2016)

---

## 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 twoN×N 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.


---

## Body

## 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 E/H node is surrounded by four H/E nodes in the interior of the problem domain. As a result, the final matrix (2N×2N, where N is the number of total E 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 twoN×N 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 E node is surrounded by four H nodes and each H node is surrounded by four E nodes. The singularity of the polynomial basis function in the E 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 orderp are defined as(1)Lpt=etp!dpdtptpe-t,forp≥0,t≥0.Using the orthogonality of the above polynomials with respect to the weighted functione-t, an orthogonal set of basis functions f0,f1,f2,… can be formed:(2)∫0∞e-tLptLqt=δpq,fpt,s=e-sc·t/2Lpsc·t,where sc>0 is a time-scale factor. The basis functions are absolutely convergent to zero as t→∞. Then, the time-domain field component is expanded as(3)Ur,t=∑p=0∞Uprfpt.The first derivative with respect tot of field component can be found as(4)∂Ur,t∂t=s∑p=0∞0.5Upr+∑k=0p-1Ukrfpsct.TakingEx(r,t) as an example, assuming a linear, isotropic, nondispersive, and lossless medium, the differential Maxwell equation can be written as(5)ε∂Exr,t∂t=∂Hzr,t∂y-Jxr,t.By inserting (4) into (5) and using a temporal Galerkin testing procedure to eliminate the time variable, one can get(6)Exqr=2scε∂∂yHzqr-2∑k=0q-1Exkr-2εJxqr.The other field qualities can be obtained in a similar way.We assume that the functionu(r) is approximated by a set of local random scattered points ui(ri)(i=1,2,…,n), where n is the number of above points which are in the neighborhood of r. RPIM is based on radial basis function ri(r) and polynomial basis function pj(r):(7)ur=∑i=1nrirai+∑j=1mpjrbj,where ai and bj are the associated expansion coefficients. Gaussian form is adopted as radial basis function in the support domain and pj(r) is built by utilizing the four-term (m=4) Pascal triangle:(8)rir=exp⁡-cx-xi2+y-yi2+z-zi2rmax2,pjr=1,xj,yj,zj,where shape parameter c controls the flatness of RPIM and rmax denotes the maximum distance between interpolation point and nodes involved in the local support domain. The vector form is defined as(9)aT=a1,a2,…,an,bT=b1,b2,…,bm,RTr=r1r,r2r,…,rnr,pTr=p1r,p2r,…,pmr.Constraints are usually imposed to ensure the unique solution:(10)∑i=1npjrai=0,j=1,2,…,m.Equations (7) and (10) can be expressed as matrix form:(11)R0P0P0T0ab=u0.Coefficient matrices R0 and P0 depend only upon the positions of corresponding points:(12)R0=RTr1,RTr2,…,RTrnT,P0=pTr1,pTr2,…,pTrmT,u=u1,u2,…,unT.Once the unknownai and bj based on (11) are solved, the final interpolation is obtained:(13)ur=RTrPTrabT=RTrPTrR0P0P0T0-1u→0=∑i=1nϕirui,where ϕi(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 ϕi(r) has no relation to field value ui. 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(14)Es,xq=2scε∑i=1NH∂ϕs,iH∂yHs,i,zq-2∑k=0q-1Es,xk-2scεJs,xq.The final expression ofEs,xq can be obtained once Hs,zq is expanded in a similar manner. After some mathematical manipulations, a set of matrices are denoted as(15)WEq=br∑i=1NHDϕiHHi,zq+γq-1+JEq,(16)Hzq=cr∑j=1NEDTφjEWEq+βq-1,where WEq=Es,xqEs,yqT, D=∂y-∂xT, ar=4/sc2με, br=2/scε, cr=2/scμ, JEq=br-Js,xq-Js,yqT, and (17)
γ
q
-
1
=
-
2
∑
k
=
0
q
-
1
E
s
,
x
k
-
2
∑
k
=
0
q
-
1
E
s
,
y
k
T
,
β
q
-
1
=
-
2
∑
k
=
0
q
-
1
H
s
,
z
k
.By inserting (16) into (15), the final huge matrix of E nodes is (18)WEq-ar∑i=1NH∑j=1NEDϕiHDTφi,jEWEq=br∑i=1NHDϕiHβq-1+γq-1+JEq.


### 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 q. 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, theE node and H node would have the same scale of support domain nodes (NE=NH=4) for the new uniform nodal placement. The number of H nodes reduces to about half of the conventional method [8], whereas the number of E nodes remains the same.The huge matrix has2N×2N×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 (DDT) of the left side into two diagonal matrices A and B, and the equation is rewritten as(19)A=ar∑i=1NH∑j=1NEϕiHφi,jE∂y,y0-∂x,y∂x,x,B=ar∑i=1NH∑j=1NEϕiHφi,jE0-∂y,x00,WEq-A+BWEq=br∑i=1NHDϕiHβq-1+γq-1+JEq.Here, A is a lower diagonal matrix and B is an upper diagonal matrix. Dκ,ζ (κ/ζ=x,y) is difference operators for the derivatives. It should be noticed that the former derivative Dκ is for the H shape function ϕiHr, while the latter one Dζ is for E shape function φjEri. In the meantime, a perturbation term ABWEq-2∑k=0q-1WEk is added and the factorized form is(20)I-AI-BWEq=I+ABγq-1+br∑i=1NHDϕiHβq-1+JEq.With the splitting scheme, the above equation can be computed into two substeps:(21)I-AWE∗=I+Bγq-1+br∑i=1NHDϕiHβq-1+JEq,I-BWEq=WE∗-Bγq-1;here, WE∗=Ex∗Ey∗T 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×N×N×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 E components must be updated following the serial numbers of the above equations. Then, the H components are updated according to (16). All the expansion coefficients are solved matching in order q. The final difference equations of (21) can be expanded as(22)Es,x∗-ar∑i=1NH∑j=1NE∂yϕs,iH∂yφs,i,jEEs,j,x∗=-αr∑i=1NH∑j=1NE∂yϕs,iH∂xφs,i,jEγs,j,yq-1+γs,xq-1+br∑i=1NH∂yϕs,iHβs,zq-1+Js,xq,Es,y∗-ar∑i=1NH∑j=1NE∂xϕs,iH∂xφs,i,jEEs,j,y∗=-ar∑i=1NH∑j=1NE∂xϕs,iH∂yφs,i,jEEs,j,x∗+γs,yq-1-br∑i=1NH∂xϕs,iHβs,zq-1+Js,yq,Ey=Ey∗,Es,xq=Es,x∗+ar∑i=1NH∑j=1NE∂yϕs,iH∂xφs,i,jEγs,j,yk-Es,j,yq,Hs,zq=cr∑j=1NE∂yφs,jEEs,j,xq-∂xφs,jEEs,j,yq-2∑k=0q-1Hs,zk.


### 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 Sections2.1 and 2.2, (15)-(16) can be rewritten as(23)Ex,sq=σy,sE2sε0∑i=1NH∂ϕs,iH∂yHz,s,iq-2σy,sE∑k=0q-1Ex,sk,Ey,sq=-σx,sE2sε0∑i=1NH∂ϕs,iH∂xHz,s,iq-2σx,sE∑k=0q-1Ey,sk,Hzx,sq=-σx,sH2sμ0∑i=1NE∂φs,iE∂xEy,s,iq-2σx,sH∑k=0q-1Hzx,sk,Hzy,sq=σy,sH2sμ0∑i=1NE∂φs,iE∂yEx,s,iq-2σy,sH∑k=0q-1Hzy,sk,Hz,sq=Hzx,sq+Hzy,sq,where(24)σsH=sε0sε0+2ρs,σsE=sε0sε0+2σs,y.Following the deduced procedure in Section2.3, (23) can be solved efficiently by factorizing it into two substeps as(25)Ex,s∗-αr∑i=1NH∑j=1NEσy,sEσy,s,iH∂ϕs,iH∂y∂φs,i,jE∂yEx,s,i,j∗=-2σy,sE∑k=0q-1Ex,sk+2αr∑i=1NH∑j=1NEσy,sEσx,s,iH∂ϕs,iH∂y∂φs,i,jE∂xσx,sE∑k=0q-1Ey,sk-2brσy,sE∑i=1NH∂ϕs,iH∂yσx,sH∑k=0q-1Hzx,sk+σy,sH∑k=0q-1Hzy,sk,Ey,s∗-αr∑i=1NH∑j=1NEσx,sEσx,sH∂ϕs,iH∂x∂φs,i,jE∂xEy,s,i,j∗=-2∑k=0q-1σx,sEEy,sk-αr∑i=1NH∑j=1NEσx,sEσy,sH∂ϕs,iH∂x∂φs,i,jE∂yEx,s,i,j∗+2br∑i=1NHσx,sE∂ϕs,iH∂xσx,sH∑k=0q-1Hzx,sk+σy,sH∑k=0q-1Hzy,sk,Ey=Ey∗,Ex,sq=Ex,s∗-αr∑i=1NH∑j=1NEσy,sEσx,s,iH∂ϕs,iH∂y∂φs,i,jE∂xEy,s,i,jq+2σx,sE∑k=0q-1Ey,sk.


## 2.1. Laguerre-RPIM

The Laguerre polynomials of orderp are defined as(1)Lpt=etp!dpdtptpe-t,forp≥0,t≥0.Using the orthogonality of the above polynomials with respect to the weighted functione-t, an orthogonal set of basis functions f0,f1,f2,… can be formed:(2)∫0∞e-tLptLqt=δpq,fpt,s=e-sc·t/2Lpsc·t,where sc>0 is a time-scale factor. The basis functions are absolutely convergent to zero as t→∞. Then, the time-domain field component is expanded as(3)Ur,t=∑p=0∞Uprfpt.The first derivative with respect tot of field component can be found as(4)∂Ur,t∂t=s∑p=0∞0.5Upr+∑k=0p-1Ukrfpsct.TakingEx(r,t) as an example, assuming a linear, isotropic, nondispersive, and lossless medium, the differential Maxwell equation can be written as(5)ε∂Exr,t∂t=∂Hzr,t∂y-Jxr,t.By inserting (4) into (5) and using a temporal Galerkin testing procedure to eliminate the time variable, one can get(6)Exqr=2scε∂∂yHzqr-2∑k=0q-1Exkr-2εJxqr.The other field qualities can be obtained in a similar way.We assume that the functionu(r) is approximated by a set of local random scattered points ui(ri)(i=1,2,…,n), where n is the number of above points which are in the neighborhood of r. RPIM is based on radial basis function ri(r) and polynomial basis function pj(r):(7)ur=∑i=1nrirai+∑j=1mpjrbj,where ai and bj are the associated expansion coefficients. Gaussian form is adopted as radial basis function in the support domain and pj(r) is built by utilizing the four-term (m=4) Pascal triangle:(8)rir=exp⁡-cx-xi2+y-yi2+z-zi2rmax2,pjr=1,xj,yj,zj,where shape parameter c controls the flatness of RPIM and rmax denotes the maximum distance between interpolation point and nodes involved in the local support domain. The vector form is defined as(9)aT=a1,a2,…,an,bT=b1,b2,…,bm,RTr=r1r,r2r,…,rnr,pTr=p1r,p2r,…,pmr.Constraints are usually imposed to ensure the unique solution:(10)∑i=1npjrai=0,j=1,2,…,m.Equations (7) and (10) can be expressed as matrix form:(11)R0P0P0T0ab=u0.Coefficient matrices R0 and P0 depend only upon the positions of corresponding points:(12)R0=RTr1,RTr2,…,RTrnT,P0=pTr1,pTr2,…,pTrmT,u=u1,u2,…,unT.Once the unknownai and bj based on (11) are solved, the final interpolation is obtained:(13)ur=RTrPTrabT=RTrPTrR0P0P0T0-1u→0=∑i=1nϕirui,where ϕi(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 ϕi(r) has no relation to field value ui. 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(14)Es,xq=2scε∑i=1NH∂ϕs,iH∂yHs,i,zq-2∑k=0q-1Es,xk-2scεJs,xq.The final expression ofEs,xq can be obtained once Hs,zq is expanded in a similar manner. After some mathematical manipulations, a set of matrices are denoted as(15)WEq=br∑i=1NHDϕiHHi,zq+γq-1+JEq,(16)Hzq=cr∑j=1NEDTφjEWEq+βq-1,where WEq=Es,xqEs,yqT, D=∂y-∂xT, ar=4/sc2με, br=2/scε, cr=2/scμ, JEq=br-Js,xq-Js,yqT, and (17)
γ
q
-
1
=
-
2
∑
k
=
0
q
-
1
E
s
,
x
k
-
2
∑
k
=
0
q
-
1
E
s
,
y
k
T
,
β
q
-
1
=
-
2
∑
k
=
0
q
-
1
H
s
,
z
k
.By inserting (16) into (15), the final huge matrix of E nodes is (18)WEq-ar∑i=1NH∑j=1NEDϕiHDTφi,jEWEq=br∑i=1NHDϕiHβq-1+γq-1+JEq.


## 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 q. 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, theE node and H node would have the same scale of support domain nodes (NE=NH=4) for the new uniform nodal placement. The number of H nodes reduces to about half of the conventional method [8], whereas the number of E nodes remains the same.The huge matrix has2N×2N×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 (DDT) of the left side into two diagonal matrices A and B, and the equation is rewritten as(19)A=ar∑i=1NH∑j=1NEϕiHφi,jE∂y,y0-∂x,y∂x,x,B=ar∑i=1NH∑j=1NEϕiHφi,jE0-∂y,x00,WEq-A+BWEq=br∑i=1NHDϕiHβq-1+γq-1+JEq.Here, A is a lower diagonal matrix and B is an upper diagonal matrix. Dκ,ζ (κ/ζ=x,y) is difference operators for the derivatives. It should be noticed that the former derivative Dκ is for the H shape function ϕiHr, while the latter one Dζ is for E shape function φjEri. In the meantime, a perturbation term ABWEq-2∑k=0q-1WEk is added and the factorized form is(20)I-AI-BWEq=I+ABγq-1+br∑i=1NHDϕiHβq-1+JEq.With the splitting scheme, the above equation can be computed into two substeps:(21)I-AWE∗=I+Bγq-1+br∑i=1NHDϕiHβq-1+JEq,I-BWEq=WE∗-Bγq-1;here, WE∗=Ex∗Ey∗T 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×N×N×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 E components must be updated following the serial numbers of the above equations. Then, the H components are updated according to (16). All the expansion coefficients are solved matching in order q. The final difference equations of (21) can be expanded as(22)Es,x∗-ar∑i=1NH∑j=1NE∂yϕs,iH∂yφs,i,jEEs,j,x∗=-αr∑i=1NH∑j=1NE∂yϕs,iH∂xφs,i,jEγs,j,yq-1+γs,xq-1+br∑i=1NH∂yϕs,iHβs,zq-1+Js,xq,Es,y∗-ar∑i=1NH∑j=1NE∂xϕs,iH∂xφs,i,jEEs,j,y∗=-ar∑i=1NH∑j=1NE∂xϕs,iH∂yφs,i,jEEs,j,x∗+γs,yq-1-br∑i=1NH∂xϕs,iHβs,zq-1+Js,yq,Ey=Ey∗,Es,xq=Es,x∗+ar∑i=1NH∑j=1NE∂yϕs,iH∂xφs,i,jEγs,j,yk-Es,j,yq,Hs,zq=cr∑j=1NE∂yφs,jEEs,j,xq-∂xφs,jEEs,j,yq-2∑k=0q-1Hs,zk.


## 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 Sections2.1 and 2.2, (15)-(16) can be rewritten as(23)Ex,sq=σy,sE2sε0∑i=1NH∂ϕs,iH∂yHz,s,iq-2σy,sE∑k=0q-1Ex,sk,Ey,sq=-σx,sE2sε0∑i=1NH∂ϕs,iH∂xHz,s,iq-2σx,sE∑k=0q-1Ey,sk,Hzx,sq=-σx,sH2sμ0∑i=1NE∂φs,iE∂xEy,s,iq-2σx,sH∑k=0q-1Hzx,sk,Hzy,sq=σy,sH2sμ0∑i=1NE∂φs,iE∂yEx,s,iq-2σy,sH∑k=0q-1Hzy,sk,Hz,sq=Hzx,sq+Hzy,sq,where(24)σsH=sε0sε0+2ρs,σsE=sε0sε0+2σs,y.Following the deduced procedure in Section2.3, (23) can be solved efficiently by factorizing it into two substeps as(25)Ex,s∗-αr∑i=1NH∑j=1NEσy,sEσy,s,iH∂ϕs,iH∂y∂φs,i,jE∂yEx,s,i,j∗=-2σy,sE∑k=0q-1Ex,sk+2αr∑i=1NH∑j=1NEσy,sEσx,s,iH∂ϕs,iH∂y∂φs,i,jE∂xσx,sE∑k=0q-1Ey,sk-2brσy,sE∑i=1NH∂ϕs,iH∂yσx,sH∑k=0q-1Hzx,sk+σy,sH∑k=0q-1Hzy,sk,Ey,s∗-αr∑i=1NH∑j=1NEσx,sEσx,sH∂ϕs,iH∂x∂φs,i,jE∂xEy,s,i,j∗=-2∑k=0q-1σx,sEEy,sk-αr∑i=1NH∑j=1NEσx,sEσy,sH∂ϕs,iH∂x∂φs,i,jE∂yEx,s,i,j∗+2br∑i=1NHσx,sE∂ϕs,iH∂xσx,sH∑k=0q-1Hzx,sk+σy,sH∑k=0q-1Hzy,sk,Ey=Ey∗,Ex,sq=Ex,s∗-αr∑i=1NH∑j=1NEσy,sEσx,s,iH∂ϕs,iH∂y∂φs,i,jE∂xEy,s,i,jq+2σx,sE∑k=0q-1Ey,sk.


## 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 Figure1. Only one-fourth of the total nodes are plotted for a clear view. The total numbers of E nodes and H 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 pulse(26)Ex,s=sin⁡2πf0t-t0exp⁡-t-t02τ2,where f0=8 GHz, τ=5f0, and t0=0.9τ. The value of shape function parameter c is chosen as 0.05 [2]. Time-scale factor and the order of the WLPs are chosen as sc=1×1011 and NL=120, respectively [6, 13].Figure 1
Node distribution for the simulated structure.The time-domain electric field of the observation point is shown in Figure2 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.Figure 2
The electric field recorded at observation point.


## 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 twoN×N matrices with 9 unknown elements in each row regardless of the duplicated ones, instead of one 2N×2N 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.



---

*Source: 8749071-2016-11-27.xml*