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

Authors: Zheng Pan; Zhikai Fu; Lihua Shi; Zhengyu Huang; Zheng Sun
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

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 high-order 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.


---

## Body

## 1. 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 ofn+1th time step and the nth time step as the value of the n+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–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.


## 2. 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(1)ϕ¯nt~=2nn!σπ1/2-1/2e-t~2/2Hnt~,where Hnt~=-1net~2(dn/dtn)e-t~2 indicates nth order Hermite polynomial and t~=t-Tf/σ is the time variable. Tf and σ are time-shift factor and time-scale factor separately. The first four orders of AH orthogonal basis are shown in Figure 1.Figure 1
The first four orders of AH orthogonal basis.The differential expression of the AH orthogonal basis can be written as(2)ddtϕ¯nt~=-1σ12ϕ¯1t~n=01σn2ϕ¯n-1t~-1σn+12ϕ¯n+1t~n≥1.For a signalf(t)=∑n=0∞αnϕ¯n(t) expanded with AH orthogonal basis, the first-order time derivative of f(t) can be expressed as(3)∂∂tft=1σ∑n=0∞αn+1n+12-αn-1n2ϕ¯nt~.


## 3. AH-FDTD Method Based on High-Order Algorithm

### 3.1. Maxwell Equations Expanded by AH Basis

Taking into account the lossless media, the time-domain Maxwell equations of the two-dimensional TEz model are(4)∂∂tExr,t=1εr∂∂yHzr,t∂∂tEyr,t=-1εr∂∂xHzr,t∂∂tHzr,t=1μr∂∂yExr,t-1μr∂∂xEyr,t,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(5)1σExq+1rq+12-Exq-1rq2=1εr∂∂yHzqr1σEyq+1rq+12-Eyq-1rq2=-1εr∂∂xHzqr1σHzq+1rq+12-Hzq-1rq2=1μr∂∂yExqr-1μr∂∂xEyqr.Discretizing (5) in space domain, the coefficient matrix equations of the discrete point (i,j) can be written as(6)αExi,j=C¯yEi,jβHzi,j-Hzi,j-1(7)αEyi,j=-C¯xEi,jβHzi,j-Hzi-1,j(8)αHzi,j=C¯yHi,jβExi,j+1-Exi,j-C¯xHi,jβEyi+1,j-Eyi,j+Hzt0i,j,where(9)α=221-12-2⋱⋱N-1-N-1N×Nand βN×N is the unit matrix.(10)C¯xEi,j=σεi,jΔx¯i,C¯yEi,j=σεi,jΔy¯jC¯xHi,j=σμi,jΔxi,C¯yHi,j=σμi,jΔyj.Exi,j,Eyi,j, and Hzi,j, respectively, represent the column vector of each field composed by AH basis expansion coefficients. Substituting (6)-(7) into (8) and gathering all points in space domain, the matrix equation on magnetic field can be obtained eventually.(11)AHz=J,where [A] is the five-diagonal band nested matrix and {[J]} 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.Figure 2
The five-diagonal banded coefficient matrix[A].


### 3.2. The High-Order Algorithm

In rectangular coordinate system, we setf(x,y,z,t) as a vector of electric field or magnetic field. The difference format with 2n order accuracy in space domain can be expressed as(12)∂f∂x=1Δx∑n=1Nanfx+Δx22n-1-fx-Δx22n-1+oΔx2N.Expandingf[x+(Δx/2)(2n-1)] and f[x-(Δx/2)(2n-1)] by Taylor series, (13) can be obtained using the method of undetermined coefficients.(13)13⋯2N-11333⋯2N-13⋮⋮⋱⋮1N3N⋯2N-1Na1a2⋮aN=10⋮0.WhenN=2, the fourth-order accuracy can be achieved in the space domain. Here, a1=9/8anda2=-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:(14)αExi,j=C¯yEi,jβ98Hzi,j-Hzi,j-1-124Hzi,j+1-Hzi,j-2αEyi,j=-C¯xEi,jβ98Hzi,j-Hzi-1,j-124Hzi+1,j-Hzi-2,jαHzi,j=C¯yHi,jβExi,j+1-Exi,j-C¯xHi,jβEyi+1,j-Eyi,j,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(15)A-Hz=J,where [A-] is the updated banded matrix and [J] 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.Figure 3
The banded coefficient matrix[A-].


## 3.1. Maxwell Equations Expanded by AH Basis

Taking into account the lossless media, the time-domain Maxwell equations of the two-dimensional TEz model are(4)∂∂tExr,t=1εr∂∂yHzr,t∂∂tEyr,t=-1εr∂∂xHzr,t∂∂tHzr,t=1μr∂∂yExr,t-1μr∂∂xEyr,t,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(5)1σExq+1rq+12-Exq-1rq2=1εr∂∂yHzqr1σEyq+1rq+12-Eyq-1rq2=-1εr∂∂xHzqr1σHzq+1rq+12-Hzq-1rq2=1μr∂∂yExqr-1μr∂∂xEyqr.Discretizing (5) in space domain, the coefficient matrix equations of the discrete point (i,j) can be written as(6)αExi,j=C¯yEi,jβHzi,j-Hzi,j-1(7)αEyi,j=-C¯xEi,jβHzi,j-Hzi-1,j(8)αHzi,j=C¯yHi,jβExi,j+1-Exi,j-C¯xHi,jβEyi+1,j-Eyi,j+Hzt0i,j,where(9)α=221-12-2⋱⋱N-1-N-1N×Nand βN×N is the unit matrix.(10)C¯xEi,j=σεi,jΔx¯i,C¯yEi,j=σεi,jΔy¯jC¯xHi,j=σμi,jΔxi,C¯yHi,j=σμi,jΔyj.Exi,j,Eyi,j, and Hzi,j, respectively, represent the column vector of each field composed by AH basis expansion coefficients. Substituting (6)-(7) into (8) and gathering all points in space domain, the matrix equation on magnetic field can be obtained eventually.(11)AHz=J,where [A] is the five-diagonal band nested matrix and {[J]} 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.Figure 2
The five-diagonal banded coefficient matrix[A].


## 3.2. The High-Order Algorithm

In rectangular coordinate system, we setf(x,y,z,t) as a vector of electric field or magnetic field. The difference format with 2n order accuracy in space domain can be expressed as(12)∂f∂x=1Δx∑n=1Nanfx+Δx22n-1-fx-Δx22n-1+oΔx2N.Expandingf[x+(Δx/2)(2n-1)] and f[x-(Δx/2)(2n-1)] by Taylor series, (13) can be obtained using the method of undetermined coefficients.(13)13⋯2N-11333⋯2N-13⋮⋮⋱⋮1N3N⋯2N-1Na1a2⋮aN=10⋮0.WhenN=2, the fourth-order accuracy can be achieved in the space domain. Here, a1=9/8anda2=-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:(14)αExi,j=C¯yEi,jβ98Hzi,j-Hzi,j-1-124Hzi,j+1-Hzi,j-2αEyi,j=-C¯xEi,jβ98Hzi,j-Hzi-1,j-124Hzi+1,j-Hzi-2,jαHzi,j=C¯yHi,jβExi,j+1-Exi,j-C¯xHi,jβEyi+1,j-Eyi,j,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(15)A-Hz=J,where [A-] is the updated banded matrix and [J] 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.Figure 3
The banded coefficient matrix[A-].


## 4. Simulation

### 4.1. 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 as4π×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 P1 and P2 are two observation points. A sinusoidal-modulated Gaussian pulse is chosen as y direction excitation source, which can be expressed as (16)Jyt=e-t-tc2/td2sin⁡2πfct-tc,where td=1/2fc, tc=3td, and fc = 0.8 GHz. The total time duration is Ts = 10.81 ns. Figure 4 is the sketch map of 2D parallel plate waveguide.Figure 4
The 2D parallel plate waveguide.


#### 4.1.1. Result Analysis

The parallel waveguide model is calculated, respectively, with AH-FDTD method and high-order 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 pointsP1 and P2. Figure 5 shows the waveform of P1, 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 P2, and the numerical dispersion is smaller using the high-order AH-FDTD method.Figure 5
The waveform ofEy on point P1.Figure 6
The waveform ofEy on point P2.Figure7 shows the relative error of the results of AH-FDTD and high-order AH-FDTD to the analytic solution. Among them, Figure 7(a) gives the relative error results of P1, and Figure 7(b) shows the result of point P2. The relative error is defined by R=20log10⁡Eyat-Eyt/max⁡Eyat, where Eyt is Y direction electric field obtained by AH-FDTD method or the modified AH-FDTD method, and Eyat is Y 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.Figure 7
The relative error between AH-FDTD and high-order AH-FDTD.
(a)
P 1 (b)
P 2From Figure7(b), we can see that the numerical dispersion is more serious than P1 because of the greater distance between point P2 and the excitation source, and the relative error is below −33 dB. The relative error of the proposed method in this paper is smaller.


### 4.2. 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 Figure8. 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−2m, 2.4 × 10−3m, and 2.4 × 10−5m. 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. P1 and P2 are the two observation points, and other settings are shown in Section 4.1.Figure 8
The 2D waveguide with dielectric material embedded in the middle region.


#### 4.2.1. Result Analysis

The simulation model is calculated, respectively, by traditional FDTD method and modified AH-FDTD method. Considering the similarity of these three results, we just show the simulation result under the widthd=2.4×10-2 m and Figures 9, 10, and 11 show the comparison result between traditional FDTD method and 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 Ey at P1, where the reflected wave from the dielectric interface can be clearly seen. And Figure 10 is the waveform of Ey at P2, which displays mainly the transmission wave.Figure 9
The waveform ofEy on point P1.Figure 10
The waveform ofEy on point P2.Figure 11
The relative error between traditional FDTD and modified AH-FDTD.The relative error between traditional FDTD and high-order AH-FDTD is shown in Figure11. The relative error size is defined by R=20log10⁡Eyat-Eyt/max⁡Eyat, where Eyt is y direction electric field obtained by the modified AH-FDTD method and Eyat is y direction electric field obtained by the traditional FDTD method. From Figure 11, we can see that the relative errors of P1 and P2 are both under −50 dB, and the relative error of P1 is smaller than the result of P2.Table1 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 d=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.Table 1
Comparison of the CPU resources for 2D waveguide.
d = 2.4  ×  10 - 2 m d = 2.4 × 10−4 m d = 2.4 × 10−5 m Traditional FDTD Modified AH-FDTD Traditional FDTD Modified AH-FDTD Traditional FDTD Modified AH-FDTD Δ t 1.6 × 10−11 s 1.6 × 10−11 s 3.98 × 10−13 s 1.6 × 10−11 s 3.98 × 10−14 s 1.6 × 10−11 s CPU time 0.0963 s 2.5423 s 5.0036 s 2.8457 s 312.7724 s 2.9816 s Memory 0.79 M 216.72 M 0.82 M 232.69 M 0.83 M 236.81 M


## 4.1. 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 as4π×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 P1 and P2 are two observation points. A sinusoidal-modulated Gaussian pulse is chosen as y direction excitation source, which can be expressed as (16)Jyt=e-t-tc2/td2sin⁡2πfct-tc,where td=1/2fc, tc=3td, and fc = 0.8 GHz. The total time duration is Ts = 10.81 ns. Figure 4 is the sketch map of 2D parallel plate waveguide.Figure 4
The 2D parallel plate waveguide.


### 4.1.1. Result Analysis

The parallel waveguide model is calculated, respectively, with AH-FDTD method and high-order 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 pointsP1 and P2. Figure 5 shows the waveform of P1, 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 P2, and the numerical dispersion is smaller using the high-order AH-FDTD method.Figure 5
The waveform ofEy on point P1.Figure 6
The waveform ofEy on point P2.Figure7 shows the relative error of the results of AH-FDTD and high-order AH-FDTD to the analytic solution. Among them, Figure 7(a) gives the relative error results of P1, and Figure 7(b) shows the result of point P2. The relative error is defined by R=20log10⁡Eyat-Eyt/max⁡Eyat, where Eyt is Y direction electric field obtained by AH-FDTD method or the modified AH-FDTD method, and Eyat is Y 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.Figure 7
The relative error between AH-FDTD and high-order AH-FDTD.
(a)
P 1 (b)
P 2From Figure7(b), we can see that the numerical dispersion is more serious than P1 because of the greater distance between point P2 and the excitation source, and the relative error is below −33 dB. The relative error of the proposed method in this paper is smaller.


## 4.1.1. Result Analysis

The parallel waveguide model is calculated, respectively, with AH-FDTD method and high-order 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 pointsP1 and P2. Figure 5 shows the waveform of P1, 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 P2, and the numerical dispersion is smaller using the high-order AH-FDTD method.Figure 5
The waveform ofEy on point P1.Figure 6
The waveform ofEy on point P2.Figure7 shows the relative error of the results of AH-FDTD and high-order AH-FDTD to the analytic solution. Among them, Figure 7(a) gives the relative error results of P1, and Figure 7(b) shows the result of point P2. The relative error is defined by R=20log10⁡Eyat-Eyt/max⁡Eyat, where Eyt is Y direction electric field obtained by AH-FDTD method or the modified AH-FDTD method, and Eyat is Y 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.Figure 7
The relative error between AH-FDTD and high-order AH-FDTD.
(a)
P 1 (b)
P 2From Figure7(b), we can see that the numerical dispersion is more serious than P1 because of the greater distance between point P2 and the excitation source, and the relative error is below −33 dB. The relative error of the proposed method in this paper is smaller.


## 4.2. 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 Figure8. 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−2m, 2.4 × 10−3m, and 2.4 × 10−5m. 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. P1 and P2 are the two observation points, and other settings are shown in Section 4.1.Figure 8
The 2D waveguide with dielectric material embedded in the middle region.


### 4.2.1. Result Analysis

The simulation model is calculated, respectively, by traditional FDTD method and modified AH-FDTD method. Considering the similarity of these three results, we just show the simulation result under the widthd=2.4×10-2 m and Figures 9, 10, and 11 show the comparison result between traditional FDTD method and 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 Ey at P1, where the reflected wave from the dielectric interface can be clearly seen. And Figure 10 is the waveform of Ey at P2, which displays mainly the transmission wave.Figure 9
The waveform ofEy on point P1.Figure 10
The waveform ofEy on point P2.Figure 11
The relative error between traditional FDTD and modified AH-FDTD.The relative error between traditional FDTD and high-order AH-FDTD is shown in Figure11. The relative error size is defined by R=20log10⁡Eyat-Eyt/max⁡Eyat, where Eyt is y direction electric field obtained by the modified AH-FDTD method and Eyat is y direction electric field obtained by the traditional FDTD method. From Figure 11, we can see that the relative errors of P1 and P2 are both under −50 dB, and the relative error of P1 is smaller than the result of P2.Table1 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 d=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.Table 1
Comparison of the CPU resources for 2D waveguide.
d = 2.4  ×  10 - 2 m d = 2.4 × 10−4 m d = 2.4 × 10−5 m Traditional FDTD Modified AH-FDTD Traditional FDTD Modified AH-FDTD Traditional FDTD Modified AH-FDTD Δ t 1.6 × 10−11 s 1.6 × 10−11 s 3.98 × 10−13 s 1.6 × 10−11 s 3.98 × 10−14 s 1.6 × 10−11 s CPU time 0.0963 s 2.5423 s 5.0036 s 2.8457 s 312.7724 s 2.9816 s Memory 0.79 M 216.72 M 0.82 M 232.69 M 0.83 M 236.81 M


## 4.2.1. Result Analysis

The simulation model is calculated, respectively, by traditional FDTD method and modified AH-FDTD method. Considering the similarity of these three results, we just show the simulation result under the widthd=2.4×10-2 m and Figures 9, 10, and 11 show the comparison result between traditional FDTD method and 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 Ey at P1, where the reflected wave from the dielectric interface can be clearly seen. And Figure 10 is the waveform of Ey at P2, which displays mainly the transmission wave.Figure 9
The waveform ofEy on point P1.Figure 10
The waveform ofEy on point P2.Figure 11
The relative error between traditional FDTD and modified AH-FDTD.The relative error between traditional FDTD and high-order AH-FDTD is shown in Figure11. The relative error size is defined by R=20log10⁡Eyat-Eyt/max⁡Eyat, where Eyt is y direction electric field obtained by the modified AH-FDTD method and Eyat is y direction electric field obtained by the traditional FDTD method. From Figure 11, we can see that the relative errors of P1 and P2 are both under −50 dB, and the relative error of P1 is smaller than the result of P2.Table1 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 d=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.Table 1
Comparison of the CPU resources for 2D waveguide.
d = 2.4  ×  10 - 2 m d = 2.4 × 10−4 m d = 2.4 × 10−5 m Traditional FDTD Modified AH-FDTD Traditional FDTD Modified AH-FDTD Traditional FDTD Modified AH-FDTD Δ t 1.6 × 10−11 s 1.6 × 10−11 s 3.98 × 10−13 s 1.6 × 10−11 s 3.98 × 10−14 s 1.6 × 10−11 s CPU time 0.0963 s 2.5423 s 5.0036 s 2.8457 s 312.7724 s 2.9816 s Memory 0.79 M 216.72 M 0.82 M 232.69 M 0.83 M 236.81 M


## 5. 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.



---

*Source: 3595947-2017-02-13.xml*