A Modified AH-FDTD Unconditionally Stable Method Based on High-Order Algorithm

The unconditionally stable method, Associated-Hermite FDTD, has attracted more and more attentions in computational electromagnetic for its time-frequency compact property. Because of the fewer orders of AH basis needed in signal reconstruction, the computational efficiency can be improved further. In order to further improve the accuracy of the traditional AH-FDTD, a highorder algorithm is introduced. Using this method, the dispersion error induced by the space grid can be reduced, which makes it possible to set coarser grid. The simulation results show that, on the condition of coarse grid, the waveforms obtained from the proposed method are matched well with the analytic result, and the accuracy of the proposed method is higher than the traditional AH-FDTD. And the efficiency of the proposed method is higher than the traditional FDTD method in analysing 2D waveguide problems with fine-structure.


Introduction
The unconditionally stable method is an important research field in computational electromagnetics.It has tremendous superiority, because it can break through the restraint of Courant-Friedrichs-Levy (CFL) stable condition and improve the efficiency of FDTD method greatly.In the past thirty years, a series of unconditionally stable methods spring up in computational electromagnetic field, such as Crank-Nicolson FDTD (CN-FDTD) method [1], Alternating-Direction Implicit FDTD (ADI-FDTD) method [2], and Weighted Laguerre polynomials FDTD (WLP-FDTD) method [3].
The CN-FDTD method adopted the average way and took the average value of  + 1th time step and the th time step as the value of the  + 1/2th time step.In general, the former iteration calculation problem can be converted into the solution of implicit equations, which breaks the limitation of the stable condition.But this method needs to solve the large matrix equations at each time step and consume huge memory, which makes the application of this method limited.ADI-FDTD can be seen as an improved algorithm based on CN-FDTD.The main idea of ADI-FDTD is to decouple the whole matrix equations into different dimension and convert the solution of the large matrix equations in the whole area into the solution of equations in each dimension, which reduces the dimensionality of matrix greatly and improves the efficiency.However, ADI-FDTD also brings one problem: the reduction of matrix dimension is obtained by abandoning the fixed term and the result is approximate.When the time step or the space step is big, the simulation results will suffer severe errors.
Chung et al. [3] put forward a new method based on weighted Laguerre polynomials orthogonal basis, which firstly expands the time-domain Maxwell equations using weighted Laguerre polynomials orthogonal basis and then eliminates every orthogonal term by Galerkin method so as to obtain the coefficients matrix equation.Finally, each point waveform in the spatial domain can be obtained by solving the coefficients matrix equation.The accuracy can be guaranteed for the signal reconstruction is based on the orthogonal basis.Later, WLP-FDTD method achieved comparatively large improvements in efficiency and boundary conditions [4][5][6].In 2014, Huang et al. [7] proposed an unconditionally stable method based on AH basis according to the most compactly supported characteristics of Associated-Hermite basis, which can reconstruct signals by fewer orders so as to increase the efficiency.In 2015, Huang et al. [8] conducted the eigenvalue decomposition on the former AH matrix and decoupled each order basis.Finally, the matrix solution can be converted into the solution of order-based submatrix equations, which increases the efficiency further and decreases the CPU memory.
Aiming at the problems of huge memory consumption and severe numerical dispersion in electrically large object simulation with conventional FDTD method, Fang firstly proposed a high-order FDTD method [9], which adopts high-order discrete scheme in space and normal FDTD scheme in time.This solution can maintain better dispersion characteristics and calculation accuracy in coarse mesh grid.Afterwards, this method is deeply researched on numerical dispersion, stability, and high-order PML absorbing boundary conditions and is further applied in analysis of waveguide field, resonance problem and electromagnetic compatibility, and so on.
This article applies the high-order algorithm to the unconditionally stable AH-FDTD method and conducts high-order discrete format in space.The accuracy can be improved and the mesh grid in non-fine-structure region can be dissected to be relatively coarse.At last, two 2D waveguide simulation models were established, which proved the validity of the proposed method.

The Associated-Hermite Orthogonal Basis
For the most compactly supported characteristics of AH orthogonal basis, it has a wide application in signal processing, image analysis, biological engineering, and so on.The AH orthogonal basis function can be expressed as where ) indicates th order Hermite polynomial and t = ( −   )/ is the time variable.  and  are time-shift factor and time-scale factor separately.The first four orders of AH orthogonal basis are shown in Figure 1.
The differential expression of the AH orthogonal basis can be written as For a signal() = ∑ ∞ =0     () expanded with AH orthogonal basis, the first-order time derivative of () can be expressed as Amplitude

Maxwell Equations Expanded by AH Basis.
Taking into account the lossless media, the time-domain Maxwell equations of the two-dimensional TEz model are where  is capacitivity and  is permeability.Substituting the time derivative term in ( 4) with (3), each order of the AH orthogonal basis can be eliminated using Galerkin method.By simplifying, it can be expressed as Discretizing ( 5) in space domain, the coefficient matrix equations of the discrete point (, ) can be written as [] [  ] , = where  6)-( 7) into ( 8) and gathering all points in space domain, the matrix equation on magnetic field can be obtained eventually.
where [] is the five-diagonal band nested matrix and {[]} is the source term.By solving the matrix equation, the field value of each point can be obtained in the space domain.Figure 2 is the diagram of the five-diagonal band matrix.
When  = 2, the fourth-order accuracy can be achieved in the space domain.Here,  1 = 9/8 and  2 = −1/24, at the same time substituting the result obtained from (12) and (13) in (5).The fourth-order space accuracy AH-FDTD format can be obtained: where the variables and parameters are in accordance with Section 3.1.Assembling the equations in the whole space, the matrix equation of magnetic field can eventually be expressed as where [] is the updated banded matrix and {[]} is the source term.By solving the equation, the field value of each point in the space domain can be obtained.Figure 3 shows the updated matrix diagram.

Case 1.
In order to verify the effectiveness of the proposed method, a two-dimensional parallel plate waveguide model is set up.In consideration of the convenience of comparing the numerical results with analytical results, the parallel waveguide is set to be filled with the homogeneous medium, and its size is 1.2 m × 0.08 m.Its permeability is taken as 4 × 10 −7 H/m and the capacitance is taken as 8.854187817 × 10 −12 F/m, which is the same as the free space.
In simulation, the whole region is divided into uniform grid of 100 × 8, and  1 and  2 are two observation points.A sinusoidal-modulated Gaussian pulse is chosen as  direction excitation source, which can be expressed as where   = 1/(2  ),   = 3  , and   = 0.8 GHz.The total time duration is   = 10.81 ns. Figure 4 is the sketch map of 2D parallel plate waveguide.

Result Analysis.
The parallel waveguide model is calculated, respectively, with AH-FDTD method and highorder AH-FDTD method.From the simulation, the results of AH-FDTD method and its modified method are both matched well with the analytical solution at points  1 and  2 .Figure 5 shows the waveform of  1 , where we can see that the numerical dispersion is smaller using the high-order AH-FDTD method, and its waveform is more close to the analytical results.Figure 6 shows the waveform of  2 , and   the numerical dispersion is smaller using the high-order AH-FDTD method.
Figure 7 shows the relative error of the results of AH-FDTD and high-order AH-FDTD to the analytic solution.Among them, Figure 7 where   () is  direction electric field obtained by AH-FDTD method or the modified AH-FDTD method, and    () is  direction electric field obtained by the analytical solution.From Figure 7(a), we can see that their error is small, and the relative error is under −40 dB.At the same time, the relative error of the proposed method in this paper is smaller.From Figure 7(b), we can see that the numerical dispersion is more serious than  1 because of the greater distance between point  2 and the excitation source, and the relative error is below −33 dB.The relative error of the proposed method in this paper is smaller.

Case 2.
In this case, a two-dimensional parallel plate waveguide with thin dielectric layer is investigated.The parallel waveguide size is 1.2 m × 0.08 m and a thin dielectric layer is embedded in the middle region shown in Figure 8.In order to clearly show the feature of the proposed algorithm, the relative permittivity of the thin layer is 2 and the width is set to three types: 2.4 × 10 −2 m, 2.4 × 10 −3 m, and 2.4 × 10 −5 m.The minimum cell size is 0.012 m × 0.01 m, 0.0012 m × 0.01 m, and 0.000012 m × 0.01 m, respectively.In the simulation, the nonuniform grid is used, and the whole calculation area is divided into 100 × 8, 108 × 8 and 110 × 8 grid. 1 and  2 are the two observation points, and other settings are shown in Section 4.1.

Result Analysis.
The simulation model is calculated, respectively, by traditional FDTD method and modified AH-FDTD method.Considering the similarity of these three   modified AH-FDTD method.As Figures 9 and 10 show, the results of traditional FDTD method and the modified AH-FDTD method are matched well with each other.Figure 9 shows the waveform of   at  1 , where the reflected wave from the dielectric interface can be clearly seen.And Figure 10 is the waveform of   at  2 , which displays mainly the transmission wave.
The relative error between traditional FDTD and highorder AH-FDTD is shown in Figure 11.The relative error size is defined by where   () is  direction electric field obtained by the modified AH-FDTD method and    () is  direction electric field obtained by the traditional FDTD method.From Figure 11, we can see that the relative errors of  1 and  2 are both under −50 dB, and the relative error of  1 is smaller than the result of  2 .
Table 1 shows the CPU resources cost under different width types.As Table 1 shows, the time step of the proposed method is not limited by the CLF stability condition and can be set as the same time step at these three situations.When the width of thin dielectric layer is  = 2.4 × 10 −5 m, the time step of the proposed method is 402 times that of the traditional FDTD method, and the total CPU time for the proposed method can be reduced to about 0.63% of the traditional FDTD method, with the accuracy still being guaranteed.And this advantage will be further enhanced as the width of the dielectric layer decreased.Considering the complexity of solving a banded nested matrix, the proposed method consumes more storage space than the traditional FDTD method, which is about 291 times the traditional FDTD method.

Conclusion
The AH-FDTD unconditionally stable method combined with high-order algorithm is proposed in this paper.Using this modified AH-FDTD method, the dispersion error induced by the space grid can be reduced, which makes it possible to set coarser grid and generally reduces the total number of grids.In simulation, two 2D parallel waveguide models are established and the simulation results show that the accuracy of the proposed method is higher than the traditional AH-FDTD when compared with the analytical solution and the efficiency of the proposed method is higher than the traditional FDTD method in analysing 2D waveguide problems with fine-structure.

Figure 1 :
Figure 1: The first four orders of AH orthogonal basis.

Figure 5 :
Figure 5: The waveform of   on point  1 .

Figure 6 :
Figure 6: The waveform of   on point  2 .

Figure 9 :
Figure 9: The waveform of   on point  1 .

Figure 11 :
Figure 11: The relative error between traditional FDTD and modified AH-FDTD.

Table 1 :
Comparison of the CPU resources for 2D waveguide.Figure 10: The waveform of   on point  2 .