2. Algorithm of the FDTD Method

The FDTD method is a computational method for analyzing the space-time dependence of electromagnetic fields by discretizing space-time variables. This method utilizes a dual lattice called a Yee lattice [1].

Figure 1 shows the Yee lattice. In the figure, there are two cubes called cells. The component parallel to the edge of the electric field is at the center of each edge of a yellow cell. The component perpendicular to the surface of the magnetic field is at the center of each surface of the cell. The cyan cell is placed in such a manner that the component parallel to the edge of the magnetic field is at the center Electro-magnetic Simulation Based on the Integral Form of Maxwell's Equations DOI: http://dx.doi.org/10.5772/intechopen.81338

Figure 1. Yee lattice used in the FDTD method.

divergence of both sides of Eq. (2) and using Eq. (4), law of charge conservation

Maxwell's equations combined with law of Lorentz are the foundation of electronics, optics, and electric circuits used to understand the physical structure dependence of an electromagnetic field distribution, the interaction between the structure and field, and other relevant characteristics. However, situations having analytical solutions of them are rare. Thus, computational method of electromag-

For computational methods of electromagnetics, there are two major types, time domain and frequency domain. In a time-domain method, time is discretized. The field distribution of a particular time step is determined by Maxwell's equations and by the distribution of the previous time step. In a frequency-domain method, the time derivative is replaced by iω, where i is the imaginary unit and ω is the angular frequency. Thus, Maxwell's equations are solved. A user chooses a method by considering the analysis object, calculation accuracy, specifications of his or her

The finite-difference time-domain (FDTD) method is a time-domain method used to analyze high-frequency electromagnetic phenomena in optical devices, antennae, and similar devices [1]. Its algorithm is based on the laws of Faraday (1), Ampére-Maxwell (2), and charge conservation (5). In the FDTD method, Gauss's laws (3) and (4) are not considered except for the initial condition. The reason can be easily understood by taking divergence of both sides of Eq. (1) and (2), and

the time derivatives of Eq. (3) and (4), respectively. This means that Gauss's laws of electric and magnetic flux densities are always satisfied when they are

The FDTD method is a computational method for analyzing the space-time dependence of electromagnetic fields by discretizing space-time variables. This

Figure 1 shows the Yee lattice. In the figure, there are two cubes called cells. The component parallel to the edge of the electric field is at the center of each edge of a yellow cell. The component perpendicular to the surface of the magnetic field is at the center of each surface of the cell. The cyan cell is placed in such a manner that the component parallel to the edge of the magnetic field is at the center

In the next section, an algorithm of the original FDTD method is shown. Next, a corrected algorithm of the FDTD method based on the integral form of Maxwell's equations is shown [2, 3]. Then, a numerical result of the propagation of electromagnetic waves in a two-dimensional slab waveguide is shown. In the subsequent section, the accuracy of the original and corrected FDTD methods is compared by showing the differences between the computational and analytical methods. The analytical method is shown in the appendix. The last section is devoted to

<sup>∂</sup>t<sup>ρ</sup> <sup>þ</sup> <sup>∇</sup> � <sup>i</sup> <sup>¼</sup> <sup>0</sup>: (5)

<sup>∂</sup>t<sup>∇</sup> � <sup>B</sup> <sup>¼</sup> <sup>0</sup>, (6) <sup>∂</sup>tð Þ¼ <sup>∇</sup> � <sup>D</sup> � <sup>ρ</sup> <sup>0</sup>, (7)

is derived:

Recent Advances in Integral Equations

netics is important.

initially satisfied.

conclusions.

64

computer, and other relevant factors.

combining the charge conservation law (5) yields

2. Algorithm of the FDTD Method

method utilizes a dual lattice called a Yee lattice [1].

of its edge, and the component normal to the surface of the electric field is at the center of its surface.

The yellow cell is used to calculate the magnetic field at time t ¼ t<sup>0</sup> þ Δt using the magnetic field at <sup>t</sup> <sup>¼</sup> <sup>t</sup><sup>0</sup> and the electric field at <sup>t</sup> <sup>¼</sup> <sup>t</sup><sup>0</sup> <sup>þ</sup> <sup>1</sup> <sup>2</sup> Δt by applying Eq. (1), where t<sup>0</sup> is a particular time and Δt is the time step. Let us now consider the top surface of the cell. At the center of the surface whose coordinates are x0; y0; z<sup>0</sup> , Eq. (1) becomes

$$
\partial\_t B\_\mathbf{z}(t, \mathbf{x}\_0, y\_0, \mathbf{z}\_0) = -\partial\_\mathbf{x} E\_\mathbf{y}(t, \mathbf{x}\_0, y\_0, \mathbf{z}\_0) + \partial\_\mathbf{y} E\_\mathbf{x}(t, \mathbf{x}\_0, y\_0, \mathbf{z}\_0), \tag{8}
$$

where the variables x, y, and z of B and E represent the x-, y-, and z-components of the B and E fields, respectively, and ∂<sup>x</sup> and ∂<sup>y</sup> represent the partial derivatives in the x- and y-directions, respectively. Replacing the partial derivatives by the central differences yields

$$\partial\_t B\_x \left( t, \mathbf{x}\_0, \mathbf{y}\_0, \mathbf{z}\_0 \right) = \frac{B\_x \left( t\_0 + \Delta t, \mathbf{x}\_0, \mathbf{y}\_0, \mathbf{z}\_0 \right) - B\_x \left( t\_0, \mathbf{x}\_0, \mathbf{y}\_0, \mathbf{z}\_0 \right)}{\Delta t} + O(\Delta t)^2,\tag{9}$$

$$\partial\_x E\_\gamma \left( t, \mathbf{x}\_0, \mathbf{y}\_0, \mathbf{z}\_0 \right) = \frac{E\_\gamma \left( t\_0 + \frac{1}{2} \Delta t, \mathbf{x}\_0 + \frac{1}{2} \Delta \mathbf{x}, \mathbf{y}\_0, \mathbf{z}\_0 \right) - E\_\gamma \left( t\_0 + \frac{1}{2} \Delta t, \mathbf{x}\_0 - \frac{1}{2} \Delta \mathbf{x}, \mathbf{y}\_0, \mathbf{z}\_0 \right)}{\Delta t} + O(\Delta \mathbf{y})^2,\tag{10}$$

$$\partial\_{\mathbf{y}}E\_{\mathbf{x}}\left(\mathbf{t}, \mathbf{x}\_{0}, \mathbf{y}\_{0}, \mathbf{z}\_{0}\right) = \frac{E\_{\mathbf{x}}\left(\mathbf{t}\_{0} + \frac{1}{2}\Delta \mathbf{t}, \mathbf{x}\_{0}, \mathbf{y}\_{0} + \frac{1}{2}\Delta \mathbf{y}, \mathbf{z}\_{0}\right) - E\_{\mathbf{x}}\left(\mathbf{t}\_{0} + \frac{1}{2}\Delta \mathbf{t}, \mathbf{x}\_{0}, \mathbf{y}\_{0} - \frac{1}{2}\Delta \mathbf{y}, \mathbf{z}\_{0}\right)}{\Delta \mathbf{y}} + O(\Delta \mathbf{x})^{2},\tag{11}$$

where <sup>t</sup> <sup>¼</sup> <sup>t</sup><sup>0</sup> <sup>þ</sup> <sup>1</sup> <sup>2</sup> Δt. Then, Bz t<sup>0</sup> þ Δt; x0; y0; z<sup>0</sup> are derived from Eqs. (9), (10), and (11) as the following:

Recent Advances in Integral Equations

$$\begin{split} B\_{\pi}(t\_{0} + \Delta t, \mathbf{x}\_{0}, y\_{0}, z\_{0}) &= B\_{\pi}(t\_{0}, \mathbf{x}\_{0}, y\_{0}, z\_{0}) + \Delta t \left[ -\frac{E\_{\mathcal{V}}\left(t\_{0} + \frac{1}{2}\Delta t, \mathbf{x}\_{0} + \frac{1}{2}\Delta \mathbf{x}, y\_{0}, z\_{0}\right)}{\Delta \mathbf{x}} \right. \\\\ &\left. + \frac{E\_{\mathcal{V}}\left(t\_{0} + \frac{1}{2}\Delta t, \mathbf{x}\_{0} - \frac{1}{2}\Delta \mathbf{x}, y\_{0}, z\_{0}\right)}{\Delta \mathbf{x}} + \frac{E\_{\mathbf{x}}\left(t\_{0} + \frac{1}{2}\Delta t, \mathbf{x}\_{0}, y\_{0} + \frac{1}{2}\Delta \mathbf{y}, z\_{0}\right)}{\Delta \mathbf{y}} \right] \\\\ &\left. - \frac{E\_{\mathbf{x}}\left(t\_{0} + \frac{1}{2}\Delta t, \mathbf{x}\_{0}, y\_{0} - \frac{1}{2}\Delta \mathbf{y}, z\_{0}\right)}{\Delta \mathbf{y}} + O(\Delta \mathbf{x}, \Delta \mathbf{y})^{2} \right] + O(\Delta t)^{3}. \end{split} \tag{12}$$

where Oð Þ Δx; Δy 2 means <sup>∑</sup><sup>m</sup>þn≥2,m,n≥<sup>2</sup><sup>O</sup> ð Þ <sup>Δ</sup><sup>x</sup> <sup>m</sup>ð Þ <sup>Δ</sup><sup>y</sup> <sup>n</sup> ð Þ. Bx and By at <sup>t</sup> <sup>¼</sup> <sup>t</sup><sup>0</sup> <sup>þ</sup> <sup>Δ</sup><sup>t</sup> can also be derived similarly. Usually, the H field can be derived from the B field. For example, in vacuum, air, or a dielectric

$$H = \frac{1}{\mu\_0} \mathbf{B},\tag{13}$$

<sup>∂</sup>zHx <sup>t</sup>; <sup>x</sup>1; <sup>y</sup>1; <sup>z</sup><sup>1</sup>

Dy t<sup>0</sup> þ

�

material

67

1 2

Hx <sup>t</sup>0; <sup>x</sup>1; <sup>y</sup>1; <sup>z</sup><sup>1</sup> � <sup>1</sup>

�iy t0; x1; y1; z<sup>1</sup>

can be approximated by

called the index, is used.

� � <sup>¼</sup> Hx <sup>t</sup><sup>0</sup> <sup>þ</sup> <sup>1</sup>

where t ¼ t0. Then, Dy t0; x1; y1; z<sup>1</sup>

DOI: http://dx.doi.org/10.5772/intechopen.81338

Δt; x1; y1; z<sup>1</sup> � �

� �

� � <sup>þ</sup> <sup>O</sup>ðΔx; <sup>Δ</sup>zÞ<sup>2</sup>

2 Δz

Δz �

<sup>2</sup> <sup>Δ</sup>t; <sup>x</sup>1; <sup>y</sup>1; <sup>z</sup><sup>1</sup> <sup>þ</sup> <sup>1</sup>

Electro-magnetic Simulation Based on the Integral Form of Maxwell's Equations

<sup>¼</sup> Dy <sup>t</sup><sup>0</sup> � <sup>1</sup>

#

<sup>2</sup> <sup>Δ</sup><sup>z</sup> � � � Hx <sup>t</sup><sup>0</sup> <sup>þ</sup> <sup>1</sup>

2

Hz t0; x<sup>1</sup> þ

<sup>þ</sup> <sup>O</sup>ð Þ <sup>Δ</sup><sup>t</sup> <sup>3</sup>

be derived from the D flux density. For example, in vacuum, air, or magnetic

<sup>E</sup> <sup>¼</sup> <sup>1</sup> ε0

<sup>E</sup> <sup>¼</sup> <sup>1</sup> ε

<sup>ε</sup><sup>∗</sup> <sup>¼</sup> <sup>ε</sup> ε0

<sup>n</sup> <sup>¼</sup> ffiffiffiffi

t ¼ t<sup>0</sup> � Δt=2 and the H field at t ¼ t0. The H field at t ¼ t<sup>0</sup> þ Δt is calculated using Eqs. (12) and (14), given the H field at t ¼ t0, having determined the E field at t ¼ t<sup>0</sup> þ Δt=2. If the time t is less than tfin, then the time becomes t þ Δt and the flow

Figure 2 shows the algorithm of the FDTD method for the case in which Eqs. (14) and (22) are satisfied. Initially, distributions of the E and H fields are given which in turn satisfy Eqs. (3) and (4). When the E field distribution at t ¼ t<sup>0</sup> � Δt=2 and the H field distribution at t ¼ t<sup>0</sup> are known, the E field at t ¼ t<sup>0</sup> þ Δt=2 is calculated using Eqs. (20) and (22), given the E field at

where ε depends on the material. Often a value

called the relative permittivity and a value

repeats. If the time t exceeds tfin, the algorithm terminates.

� � is derived as follows:

; x1; y1; z<sup>1</sup> � �

> 1 2

Δx

:

Dx and Dz at t ¼ t<sup>0</sup> þ Δt can also be derived similarly. Typically, the E field can

where ε<sup>0</sup> is the vacuum permittivity, whose value is <sup>8</sup>:<sup>85418782</sup> � <sup>10</sup>�<sup>12</sup> ½ � As=Vm . In a dielectric material, the relationship between E and D often becomes nontrivial, but this is beyond the scope of this book. However, in small E and D regions, it

� �

Δx; y1; z<sup>1</sup>

þ Δt

þ

<sup>2</sup> <sup>Δ</sup>t; <sup>x</sup>1; <sup>y</sup>1; <sup>z</sup><sup>1</sup> � <sup>1</sup> <sup>2</sup> <sup>Δ</sup><sup>z</sup> � � <sup>Δ</sup><sup>z</sup> <sup>þ</sup> <sup>O</sup>ð Þ <sup>Δ</sup><sup>z</sup> <sup>2</sup>

Hx t0; x1; y1; z<sup>1</sup> þ

Hz <sup>t</sup>0; <sup>x</sup><sup>1</sup> � <sup>1</sup>

D, (21)

D, (22)

, (23)

<sup>ε</sup><sup>∗</sup> <sup>p</sup> , (24)

� �

Δz

2

Δx

� �

,

(19)

1 2 Δz

Δx; y1; z<sup>1</sup>

(20)

where <sup>μ</sup><sup>0</sup> is the vacuum permeability with the value 4<sup>π</sup> � <sup>10</sup>�<sup>7</sup> ½ � Vs=Am . In a magnetic material, the relationship between H and B often becomes nontrivial, but this is beyond the scope of this book. However, in small H and B regions, it can be approximated by the following equation:

$$H = \frac{1}{\mu} \mathbf{B},\tag{14}$$

where μ depends on the material. Often, a value

$$
\mu^\* = \frac{\mu}{\mu\_0},
\tag{15}
$$

called the relative permeability is used. In an optical wavelength region, μ<sup>0</sup> is 1.

The cyan cell is used to calculate the electric field at <sup>t</sup> <sup>¼</sup> <sup>t</sup><sup>0</sup> <sup>þ</sup> <sup>1</sup> <sup>2</sup> Δt using the electric field at <sup>t</sup> <sup>¼</sup> <sup>t</sup><sup>0</sup> � <sup>1</sup> <sup>2</sup> Δt and the magnetic field at t ¼ t � Δt applying Eq. (2) representing Ampére-Maxwell's law. Let us consider the right-hand surface of the cell. At the center of the surface whose coordinates are x1; y1; z<sup>1</sup> � �, Eq. (2) becomes

$$\partial\_t D\_\mathcal{\mathcal{Y}}(t, \mathbf{x}\_1, \mathbf{y}\_1, \mathbf{z}\_1) = \partial\_\mathbf{z} H\_\mathbf{x}(t, \mathbf{x}\_1, \mathbf{y}\_1, \mathbf{z}\_1) - \partial\_\mathbf{x} H\_\mathbf{z}(t, \mathbf{x}\_1, \mathbf{y}\_1, \mathbf{z}\_1) - i\_\mathbf{y}(t, \mathbf{x}\_1, \mathbf{y}\_1, \mathbf{z}\_1). \tag{16}$$

Replacing partial derivatives with central differences yields

$$\partial\_t D\_\mathbf{y} \left( t, \mathbf{x}\_1, \mathbf{y}\_1, \mathbf{z}\_1 \right) = \frac{D\_\mathbf{y} \left( t\_0 + \Delta t, \mathbf{x}\_1, \mathbf{y}\_1, \mathbf{z}\_1 \right) - D\_\mathbf{y} \left( t\_0, \mathbf{x}\_1, \mathbf{y}\_1, \mathbf{z}\_1 \right)}{\Delta t} + O \left( \Delta t \right)^2,\tag{17}$$

$$\partial\_{\mathbf{x}}H\_{\mathbf{z}}(\mathbf{t},\mathbf{x}\_{1},\mathbf{y}\_{1},\mathbf{z}\_{1}) = \frac{H\_{\mathbf{z}}\left(\mathbf{t}\_{0} + \frac{1}{2}\Delta t, \mathbf{x}\_{1} + \frac{1}{2}\Delta \mathbf{x}, \mathbf{y}\_{1},\mathbf{z}\_{1}\right) - H\_{\mathbf{z}}\left(\mathbf{t}\_{0} + \frac{1}{2}\Delta \mathbf{t}, \mathbf{x}\_{1} - \frac{1}{2}\Delta \mathbf{x}, \mathbf{y}\_{1},\mathbf{z}\_{1}\right)}{\Delta \mathbf{x}} + O\left(\Delta \mathbf{z}\right)^{2},\tag{18}$$

Electro-magnetic Simulation Based on the Integral Form of Maxwell's Equations DOI: http://dx.doi.org/10.5772/intechopen.81338

$$\partial\_{\mathbf{z}}H\_{\mathbf{x}}\left(\mathbf{t},\mathbf{x}\_{1},\mathbf{y}\_{1},\mathbf{z}\_{1}\right) = \frac{H\_{\mathbf{x}}\left(\mathbf{t}\_{0} + \frac{1}{2}\Delta t, \mathbf{x}\_{1}, \mathbf{y}\_{1}, \mathbf{z}\_{1} + \frac{1}{2}\Delta \mathbf{z}\right) - H\_{\mathbf{x}}\left(\mathbf{t}\_{0} + \frac{1}{2}\Delta t, \mathbf{x}\_{1}, \mathbf{y}\_{1}, \mathbf{z}\_{1} - \frac{1}{2}\Delta \mathbf{z}\right)}{\Delta \mathbf{z}} + O\left(\Delta \mathbf{z}\right)^{2},\tag{19}$$

where t ¼ t0. Then, Dy t0; x1; y1; z<sup>1</sup> � � is derived as follows:

$$D\_{\mathcal{V}}\left(t\_0 + \frac{1}{2}\Delta t, \mathbf{x}\_1, \mathbf{y}\_1, \mathbf{z}\_1\right) = D\_{\mathcal{V}}\left(t\_0 - \frac{1}{2}, \mathbf{x}\_1, \mathbf{y}\_1, \mathbf{z}\_1\right) + \Delta t \left[\frac{H\_{\mathbf{x}}\left(t\_0, \mathbf{x}\_1, \mathbf{y}\_1, \mathbf{z}\_1 + \frac{1}{2}\Delta \mathbf{z}\right)}{\Delta \mathbf{z}}\right]$$

$$-\frac{H\_x\left(t\_0, \mathbf{x}\_1, \mathbf{y}\_1, \mathbf{z}\_1 - \frac{1}{2}\Delta z\right)}{\Delta x} - \frac{H\_x\left(t\_0, \mathbf{x}\_1 + \frac{1}{2}\Delta \mathbf{x}, \mathbf{y}\_1, \mathbf{z}\_1\right)}{\Delta x} + \frac{H\_x\left(t\_0, \mathbf{x}\_1 - \frac{1}{2}\Delta \mathbf{x}, \mathbf{y}\_1, \mathbf{z}\_1\right)}{\Delta x}$$

$$\begin{bmatrix} -i\_{\eta} \left( t\_0, \mathbf{x}\_1, y\_1, \mathbf{z}\_1 \right) + O \left( \Delta \mathbf{x}, \Delta \mathbf{z} \right)^2 \end{bmatrix} + O \left( \Delta t \right)^3. \tag{20}$$

Dx and Dz at t ¼ t<sup>0</sup> þ Δt can also be derived similarly. Typically, the E field can be derived from the D flux density. For example, in vacuum, air, or magnetic material

$$\mathbf{E} = \frac{1}{\varepsilon\_0} \mathbf{D},\tag{21}$$

where ε<sup>0</sup> is the vacuum permittivity, whose value is <sup>8</sup>:<sup>85418782</sup> � <sup>10</sup>�<sup>12</sup> ½ � As=Vm . In a dielectric material, the relationship between E and D often becomes nontrivial, but this is beyond the scope of this book. However, in small E and D regions, it can be approximated by

$$\mathbf{E} = \frac{1}{\varepsilon} \mathbf{D},\tag{22}$$

where ε depends on the material. Often a value

$$
\varepsilon^\* = \frac{\varepsilon}{\varepsilon\_0},
\tag{23}
$$

called the relative permittivity and a value

$$
\mathfrak{n} = \sqrt{\mathfrak{e}^\*},
\tag{24}
$$

called the index, is used.

Figure 2 shows the algorithm of the FDTD method for the case in which Eqs. (14) and (22) are satisfied. Initially, distributions of the E and H fields are given which in turn satisfy Eqs. (3) and (4). When the E field distribution at t ¼ t<sup>0</sup> � Δt=2 and the H field distribution at t ¼ t<sup>0</sup> are known, the E field at t ¼ t<sup>0</sup> þ Δt=2 is calculated using Eqs. (20) and (22), given the E field at t ¼ t<sup>0</sup> � Δt=2 and the H field at t ¼ t0. The H field at t ¼ t<sup>0</sup> þ Δt is calculated using Eqs. (12) and (14), given the H field at t ¼ t0, having determined the E field at t ¼ t<sup>0</sup> þ Δt=2. If the time t is less than tfin, then the time becomes t þ Δt and the flow repeats. If the time t exceeds tfin, the algorithm terminates.

Bz t<sup>0</sup> þ Δt; x0; y0; z<sup>0</sup>

where Oð Þ Δx; Δy

electric field at <sup>t</sup> <sup>¼</sup> <sup>t</sup><sup>0</sup> � <sup>1</sup>

<sup>∂</sup>tDy <sup>t</sup>; <sup>x</sup>1; <sup>y</sup>1; <sup>z</sup><sup>1</sup>

<sup>∂</sup>tDy <sup>t</sup>; <sup>x</sup>1; <sup>y</sup>1; <sup>z</sup><sup>1</sup>

� � <sup>¼</sup> Hz <sup>t</sup><sup>0</sup> <sup>þ</sup> <sup>1</sup>

<sup>∂</sup>xHz <sup>t</sup>; <sup>x</sup>1; <sup>y</sup>1; <sup>z</sup><sup>1</sup>

66

becomes

� � <sup>¼</sup> Bz <sup>t</sup>0; <sup>x</sup>0; <sup>y</sup>0; <sup>z</sup><sup>0</sup>

Recent Advances in Integral Equations

þ

�

2

For example, in vacuum, air, or a dielectric

approximated by the following equation:

� � <sup>¼</sup> <sup>∂</sup>zHx <sup>t</sup>; <sup>x</sup>1; <sup>y</sup>1; <sup>z</sup><sup>1</sup>

� � <sup>¼</sup> Dy <sup>t</sup><sup>0</sup> <sup>þ</sup> <sup>Δ</sup>t; <sup>x</sup>1; <sup>y</sup>1; <sup>z</sup><sup>1</sup>

where μ depends on the material. Often, a value

Ey t<sup>0</sup> þ

Ex t<sup>0</sup> þ

� � <sup>þ</sup> <sup>Δ</sup><sup>t</sup> �

<sup>Δ</sup>t; <sup>x</sup><sup>0</sup> � <sup>1</sup> 2

� �

<sup>Δ</sup>t; <sup>x</sup>0; <sup>y</sup><sup>0</sup> � <sup>1</sup>

Δy

� �

1 2

1 2 Ey t<sup>0</sup> þ

Δx; y0; z<sup>0</sup>

Δx þ

2 Δy; z<sup>0</sup>

can also be derived similarly. Usually, the H field can be derived from the B field.

<sup>H</sup> <sup>¼</sup> <sup>1</sup> μ0

where <sup>μ</sup><sup>0</sup> is the vacuum permeability with the value 4<sup>π</sup> � <sup>10</sup>�<sup>7</sup> ½ � Vs=Am . In a magnetic material, the relationship between H and B often becomes nontrivial, but this is beyond the scope of this book. However, in small H and B regions, it can be

> <sup>H</sup> <sup>¼</sup> <sup>1</sup> μ

<sup>μ</sup><sup>∗</sup> <sup>¼</sup> <sup>μ</sup> μ0

representing Ampére-Maxwell's law. Let us consider the right-hand surface of

� � � <sup>∂</sup>xHz <sup>t</sup>; <sup>x</sup>1; <sup>y</sup>1; <sup>z</sup><sup>1</sup>

� � � Dy <sup>t</sup>0; <sup>x</sup>1; <sup>y</sup>1; <sup>z</sup><sup>1</sup>

Δt

<sup>2</sup> Δx; y1; z<sup>1</sup> � � � Hz <sup>t</sup><sup>0</sup> <sup>þ</sup> <sup>1</sup>

The cyan cell is used to calculate the electric field at <sup>t</sup> <sup>¼</sup> <sup>t</sup><sup>0</sup> <sup>þ</sup> <sup>1</sup>

the cell. At the center of the surface whose coordinates are x1; y1; z<sup>1</sup>

Replacing partial derivatives with central differences yields

<sup>2</sup> <sup>Δ</sup>t; <sup>x</sup><sup>1</sup> <sup>þ</sup> <sup>1</sup>

called the relative permeability is used. In an optical wavelength region, μ<sup>0</sup> is 1.

<sup>2</sup> Δt and the magnetic field at t ¼ t � Δt applying Eq. (2)

1 2 Δt; x<sup>0</sup> þ 1 2

1 2

B, (13)

B, (14)

, (15)

� � � iy <sup>t</sup>; <sup>x</sup>1; <sup>y</sup>1; <sup>z</sup><sup>1</sup>

<sup>2</sup> <sup>Δ</sup>t; <sup>x</sup><sup>1</sup> � <sup>1</sup>

� � <sup>Δ</sup><sup>x</sup> <sup>þ</sup> <sup>O</sup>ð Þ <sup>Δ</sup><sup>x</sup> <sup>2</sup>

� �

<sup>2</sup> Δt using the

� �: (16)

, (17)

,

(18)

� �, Eq. (2)

<sup>þ</sup> <sup>O</sup>ð Þ <sup>Δ</sup><sup>t</sup> <sup>2</sup>

<sup>2</sup> Δx; y1; z<sup>1</sup>

� �

Δx

Ex t<sup>0</sup> þ

<sup>þ</sup> <sup>O</sup>ð Þ <sup>Δ</sup>x; <sup>Δ</sup><sup>y</sup> <sup>2</sup>

means <sup>∑</sup><sup>m</sup>þn≥2,m,n≥<sup>2</sup><sup>O</sup> ð Þ <sup>Δ</sup><sup>x</sup> <sup>m</sup>ð Þ <sup>Δ</sup><sup>y</sup> <sup>n</sup> ð Þ. Bx and By at <sup>t</sup> <sup>¼</sup> <sup>t</sup><sup>0</sup> <sup>þ</sup> <sup>Δ</sup><sup>t</sup>

Δx; y0; z<sup>0</sup>

Δt; x0; y<sup>0</sup> þ

Δy

<sup>þ</sup> <sup>O</sup>ð Þ <sup>Δ</sup><sup>t</sup> <sup>3</sup> :

� �

1 2 Δy; z<sup>0</sup>

(12)

An algorithm in which differentials are not replaced with differences must be

Electro-magnetic Simulation Based on the Integral Form of Maxwell's Equations

<sup>d</sup><sup>a</sup> � <sup>B</sup> ¼ � <sup>Z</sup>

Z

∂S

da � Bð Þ� t0; x; y; z

where S is a compact and connected surface, ∂S is the boundary curve of the surface, da is a surface element normal to the surface, and ds is a line element parallel to the curve. Moreover, integrating both sides of Eq. (25) over t from t ¼ t<sup>0</sup> to t<sup>0</sup> þ Δt and those of Eq. (26) over t from t ¼ t<sup>0</sup> � Δt=2 to t<sup>0</sup> þ Δt=2 yields

As a result of Stokes' theorem, Faraday's and Ampére-Maxwell's laws in Eqs. (1)

∂S

Z

S

t Z <sup>0</sup>þΔt

t0

da � Dð Þ t<sup>0</sup> � Δt=2; x; y; z

dt <sup>Z</sup>

Note that no derivative is used in Eqs. (27) and (28), with the result that problems such as fermion doubling cannot occur. Our problem is how to approxi-

> Z ξ0þΔξ=2

ξ0�Δξ=2

and when g is analytical in a region including ξ<sup>0</sup> � Δξ=2 ≤ξ≤ξ<sup>0</sup> þ Δξ=2 and

Z<sup>η</sup>0þΔη=<sup>2</sup>

η0�Δη=2

The lowest-order approximations of Eqs. (29) and (30) are, respectively,

2 4

In general, when f is analytical in a region including ξ<sup>0</sup> � Δξ=2≤ ξ≤ξ<sup>0</sup> þ Δξ=2,

dξ ∑ ∞ n¼0

1 n! dn f ξ<sup>0</sup> ð Þ

dη ∑ ∞ m, <sup>n</sup>¼<sup>0</sup>

<sup>d</sup><sup>ξ</sup> <sup>f</sup>ð Þ¼ <sup>ξ</sup> <sup>f</sup> <sup>ξ</sup><sup>0</sup> ð ÞΔ<sup>ξ</sup> <sup>þ</sup> <sup>O</sup>ð Þ <sup>Δ</sup><sup>ξ</sup> <sup>3</sup>

1 m!n!

∂S

dt Z

ds � Hð Þ� t; x; y; z

∂S

Z

S

<sup>d</sup>ξ<sup>n</sup> <sup>ξ</sup> � <sup>ξ</sup><sup>0</sup> ð Þ<sup>n</sup>

<sup>∂</sup><sup>m</sup>þng <sup>ξ</sup>0; <sup>η</sup><sup>0</sup> ð Þ

ds � H �

ds � E, (25)

da � i, (26)

ds � Eð Þ t; x; y; z , (27)

da � ið Þ t; x; y; z

, (29)

<sup>∂</sup>ξ<sup>m</sup>∂η<sup>n</sup> <sup>ξ</sup> � <sup>ξ</sup><sup>0</sup> ð Þ<sup>m</sup> <sup>η</sup> � <sup>η</sup><sup>0</sup> ð Þ<sup>n</sup>

, (31)

:

(30)

3 5:

(28)

considered in order to avoid such problems.

DOI: http://dx.doi.org/10.5772/intechopen.81338

and (2) can be written in an integral form as

<sup>d</sup><sup>a</sup> � <sup>B</sup>ðt<sup>0</sup> <sup>þ</sup> <sup>Δ</sup>t; <sup>x</sup>; <sup>y</sup>; <sup>z</sup>Þ ¼ <sup>Z</sup>

<sup>d</sup><sup>a</sup> � <sup>D</sup>ðt<sup>0</sup> <sup>þ</sup> <sup>Δ</sup>t=2; <sup>x</sup>; <sup>y</sup>; <sup>z</sup>Þ ¼ <sup>Z</sup>

mate the integrals in Eqs. (27) and (28).

the following relationship is satisfied:

Z ξ0þΔξ=2

ξ0�Δξ=2

dηgð Þ¼ ξ; η

Z

S

Z

S

Z ξ0þΔξ=2

dξ

Z<sup>η</sup>0þΔη=<sup>2</sup>

η0�Δη=2

ξ0�Δξ=2

69

d dt Z S

d dt Z S

S

S

Z<sup>t</sup>0þΔt=<sup>2</sup>

t0�Δt=2

þ

dξ fð Þ¼ ξ

η<sup>0</sup> � Δη=2 ≤η≤η<sup>0</sup> þ Δη=2, the following relationship is satisfied:

dξ

Z ξ0þΔξ=2

ξ0�Δξ=2

Z ξ0þΔξ=2

ξ0�Δξ=2

da � D ¼

Figure 2. Flow of the FDTD algorithm.
