Subsections

• 2.4.1 Basic Concepts
• 2.4.2 Analysis of the Finite Difference Method
• 2.4.3 Further Analysis

2.4 Finite Differences

The finite difference discretization scheme is one of the simplest forms of discretization and does not easily include the topological nature of equations. A classical finite difference approach approximates the differential operators constituting the field equation locally. Therefore a structured grid is required to store local field quantities. For each of the points of the structured grid the differential operators appearing in the main problem specification are rendered in a discrete expression. The order of the differential operator of the original problem formulation directly dictates the number of nodes to be involved.

Here, the main drawback of finite differences can already be seen. The association of physical field values only to points cannot handle higher dimensional geometrical objects. Furthermore, the n-point discretization given by the order of the original formulation only includes the logically direct orthogonal neighbors while the other neighbors are neglected, depicted in Figure 2.8.

Figure 2.8: Finite difference requirements for the mesh.

The advantages of this method are that it is easy to understand and to implement, at least for simple material relations. The finite difference method optimizes the approximation for the differential operator in the central node of the considered patch. Enhancements related to the use of non-orthogonal grids and the low order of accuracy were developed but have not proven successful.

2.4.1 Basic Concepts

The derivatives of the partial differential equation are approximated by linear combinations of function values at the structured grid points. Arbitrary order approximations can be derived from a Taylor series expansion:

$$u(x) = \sum_{n=0}^{\infty} \frac{(x-x_i)^n}{n!} \left ( \frac{\partial^n u}{\partial x^n} \right )_i, \; u \in C^\infty$$ (2.47)

A geometric interpretation of the different equations is shown in Figure 2.9.

Figure 2.9: Different geometric interpretations of the first-order finite difference approximation related to forward, backward, and central difference approximation.

$$\color{red} \left ( \frac{\partial u}{ \partial x} \right ) \approx \frac{u_{i+1}-u_i}{\Delta x}$$ $$\color{green}\left ( \frac{\partial u}{ \partial x} \right ) \approx \frac{u-u_{i-1}}{\Delta x}$$ $$\color{blue} \left ( \frac{\partial u}{ \partial x} \right ) \approx \frac{u_{i+1}-u_{i-1}}{2 \Delta x}$$

For second-order derivatives the central difference scheme can be used:

$$T1: \; u_{i+1} = u_i + \Delta x \left ( \frac{\partial u}{\partia... ...i + \frac{(\Delta x)^3}{6} \left ( \frac{\partial^3 u}{\partial x^3} \right )_i$$ (2.48)

$$T2: \; u_{i\color{red}-\color{black}1} = u_i \color{red}-\color{b... ...ck} \frac{(\Delta x)^3}{6} \left ( \frac{\partial^3 u}{\partial x^3} \right )_i$$ (2.49)

$$T1 + T2 \Rightarrow \left ( \frac{\partial^2 u}{\partial x^2} \right )_i = \frac{u_{i+1} - 2 u_i + u_{i-1}}{(\Delta x)^2} + \mathcal{O}(\Delta x)^2$$ (2.50)

The accuracy of the finite difference approximations is given by:

• forward difference: truncation error: $\mathcal{O}(\Delta x)$
• backwards difference: truncation error: $\mathcal{O}(\Delta x)$
• central difference: truncation error: $\mathcal{O}(\Delta x^2)$

Mixed derivatives, illustrated in Figure 2.10 can be approximated, e.g., for two dimensions by:

$$\frac{\partial^2 u}{\partial x \partial y} = \frac{\partial}{\par... ... ) = \frac{\partial}{\partial y} \left ( \frac{\partial u}{\partial x} \right )$$ (2.51)

Figure 2.10: Approximation of two-dimensional mixed derivatives.

$$\left (\frac{\partial^2 u}{\partial x \partial y} \right )_{i,j} = \frac{ \left( \frac{\partial u}{\partial y} \right )_{i+1,j} - ... ...artial u}{\partial y} \right )_{i-1,j} }{ 2 \Delta x} + \mathcal{O}(\Delta x)^2$$ $$\left (\frac{\partial u}{\partial y} \right )_{i+1,j} = \frac {u_{i+1,j+1}-u_{i+1,j-1}}{2 \Delta y} + \mathcal{O}(\Delta y)^2$$ $$\left (\frac{\partial u}{\partial y} \right )_{i-1,j} = \frac {u_{i-1,j+1}-u_{i-1,j-1}}{2 \Delta y} + \mathcal{O}(\Delta y)^2$$

A second-order finite difference approximation for two dimensions results in:

$$\left ( \frac{\partial^2 u}{\partial x \partial y} \right )_{i,j}... ...1}}{4 \Delta x \Delta y} + \mathcal{O}\left( (\Delta x)^2, (\Delta y)^2\right )$$ (2.52)

Boundaries also have to be handled by the finite difference method, see Figure 2.11,

Figure 2.11: Boundary treatment for the finite difference scheme. The necessary grid points are not available for all different types of approximation schemes.

$$\left ( \frac{\partial u}{\partial x} \right )_{x_0} = \frac{u_1 - u_0}{\Delta x} + \mathcal{O}(\Delta x)$$ (2.53)

but at the boundary on the left side, the backward and central difference approximations would need $x_{-1}$ which is not available. Therefore treatment of boundaries is more complex for the finite difference method compared to the other discretization schemes. This particular problem grows in complexity if higher order discretization schemes are used, because more grid points are required to evaluate the corresponding approximation. Higher order approximations with finite differences are given by, e.g.:

$$\left ( \frac{\partial u}{\partial x} \right )_i = \frac{2 u_{i+1... ...i-2}}{6 \Delta x} + \mathcal{O}(\Delta x)^3 \; \textnormal{backward difference}$$ (2.54)

$$\left ( \frac{\partial u}{\partial x} \right )_i = \frac{- u_{i+2... ...{i-1}}{6 \Delta x} + \mathcal{O}(\Delta x)^3 \; \textnormal{forward difference}$$ (2.55)

$$\left ( \frac{\partial u}{\partial x} \right )_i = \frac{- u_{i+2... ...i-2}}{12 \Delta x} + \mathcal{O}(\Delta x)^4 \; \textnormal{central difference}$$ (2.56)

$$\left ( \frac{\partial^2 u}{\partial x^2} \right )_i = \frac{- u_... ...}{(12 \Delta x)^2} + \mathcal{O}(\Delta x)^4 \; \textnormal{central difference}$$ (2.57)

2.4.2 Analysis of the Finite Difference Method

One method of directly transfering the discretization concepts (Section 2.1) is the finite difference time domain method. It is analyzed here related to time-dependent Maxwell equations, as was first introduced by Yee [30]. It is one of the exceptional examples of engineering illustrating great insights into discretization processes.

With this method, the partial spatial and time derivatives are replaced by a finite difference approximation. This system is solved using an explicit time evaluation. One of the main advantages of this method is that no matrix operations or algebraic solution methods have to be used.

The spatial domain is discretized by two dual orthogonal regular Cartesian grids based on cubes with spatial subdivisions of $\Delta x, \Delta y, \Delta z$ , whereas the time domain is subdivided into intervals of $\Delta t$ . The original formulation was based on half-step staggered grids in space and time. The quantities from the second complex are denoted with $\Delta \tilde x, \Delta \tilde y, \Delta \tilde z, \Delta \tilde t$ . It is important to highlight that two different grids are necessary, due to the fact that different quantities with different orientation reside on these two distinct grids, even if in this method the secondary quantities coincide numerically with the primary ones.

Maxwell's equations, given in Section 1.2, if projected onto a regular Cartesian structured grid, yield the following six coupled scalar equations:

$$\partial_t H_x = \frac{1}{\mu}(\partial_z E_y - \partial_y E_z )$$ (2.58)

$$\partial_t H_y = \frac{1}{\mu}(\partial_x E_z - \partial_z E_x )$$ (2.59)

$$\partial_t H_z = \frac{1}{\mu}(\partial_y E_x - \partial_x E_y )$$ (2.60)

$$\partial_t E_x = \frac{1}{\varepsilon}(\partial_y H_z - \partial_z H_y - \gamma E_x)$$ (2.61)

$$\partial_t E_y = \frac{1}{\varepsilon}(\partial_z H_x - \partial_x H_z - \gamma E_y)$$ (2.62)

$$\partial_t E_z = \frac{1}{\varepsilon}(\partial_x H_y - \partial_y H_x - \gamma E_z)$$ (2.63)

These equations are the basic expressions for the finite difference time domain method (FDTD). The divergence relations are fulfilled by this method implicitly.

The components of the electric and magnetic field $\ensuremath{\mathbf{E}}$ and $\ensuremath{\mathbf{H}}$ with their corresponding projections to the coordinate axes are the variables used. These variables and the local values of material properties are attached to the midpoints of the grid edges. The variable indexing scheme is also used consistently with [30]

$$E_x(i \Delta x,j \pm 1/2 \Delta y ,k \Delta z, n \Delta t) = E_{x\vert i,j \pm 1/2,k}^{n}$$ (2.64)

With this notation the following expressions are obtained with central difference approximation:

$$\partial_x E_z (i\Delta x, j\Delta y, k \Delta z, n \Delta t) = \frac{E_{z\vert i+1/2, j,k}^{n} - E_{z\vert i-1/2, j,k}^{n} }{\Delta x} + \mathcal{O}(\Delta x)^2$$ (2.65)

$$\partial_t E_z (i\Delta x, j\Delta y, k \Delta z, n \Delta t) = \frac{E_{z\vert i,j,k}^{n+1/2} - E_{z\vert i,j,k}^{n-1/2}}{\Delta t} + \mathcal{O}(\Delta t)^2$$ (2.66)

The time-stepping formulas for $H_x$ are:

$$H_{x\vert i+1/2,j,k}^{n+1} = H_{x\vert i+1/2,j,k}^{n} +$$ (2.67)

$$\frac{\Delta t}{\Delta z  \mu_{\vert i+1/2,j,k}} \left ( E_{y\vert(i+1/2,j,k+1/2)}^{n+1/2} - E_{y\vert(i+1/2,j,k-1/2)}^{n+1/2} \right ) -$$

$$\frac{\Delta t}{\Delta y  \mu_{\vert i+1/2,j,k}} \left ( E_{z\vert(i+1/2,j+1/2,k)}^{n+1/2} - E_{z\vert(i+1/2,j-1/2,k)}^{n+1/2} \right )$$

and the time-stepping for $E_x$ :

$$E_{x\vert(i,j+1/2,k+1/2)}^{n+1/2} = E_{x\vert(i,j+1/2,k+1/2)}^{n-1/2} +$$ (2.68)

$$\frac{\Delta \tilde t}{\Delta \tilde z  \varepsilon_{\vert(i,j+1... .../2)}} \left ( H_{y\vert(i,j+1/2,k+1)}^{n} -H_{y\vert(i,j+1/2,k)}^{n} \right ) -$$

$$\frac{\Delta \tilde t}{\Delta \tilde y  \varepsilon_{\vert(i,j+1... ...+1/2)}} \left ( H_{z\vert(i,j+1,k+1/2)}^{n} -H_{z\vert(i,j,k+1/2)}^{n} \right )$$

2.4.3 Further Analysis

As already introduced, the FDTD method does not use global quantities. Instead only local nodal values of the corresponding vector values are used. Based on the initial problem formulation, it can be seen that these local values are the projections of averaged field components onto 2-cells, and therefore are local representatives of the global quantities. With local constitutive relations (only the $x$ part is given):

$$B_{x\vert(i+1/2, j,k)}^{n} = \mu_{\vert(i+1/2,j,k)} H_{x\vert(i+1/2,j,k)}^{n}$$ (2.69)

and

$$D_{x\vert(i,j+1/2,k+1/2)}^{n+1/2} = \varepsilon_{\vert(i,j+1/2,k+1/2)} E_{x\vert(i,j+1/2,k+1/2)}^{n+1/2}$$ (2.70)

the following expressions are obtained using global quantities:

$$\Delta y \Delta z \mu_{\vert(i+1/2,j,k)} H_{x\vert(i+1/2,j,k)}^{n+1} = \Phi_{B,x\vert(i+1/2,j,k)}^{n+1}$$ (2.71)

$$\Delta y \Delta t E_{y\vert(i+1/2,j,k\pm 1/2)}^{n+1/2} = \Phi_{E,y\vert(i+1/2,j,k \pm 1/2)}^{n+1/2}$$ (2.72)

$$\Delta z \Delta t E_{z\vert(i+1/2,j\pm 1/2,k)}^{n+1/2} = \Phi_{E,z\vert(i+1/2,j \pm 1/2,k)}^{n+1/2}$$ (2.73)

and

$$\Delta \tilde y \Delta \tilde z \varepsilon_{\vert(i,j+1/2,k+1/2)} E_{x\vert(i,j+1/2,k+1/2)}^{n+1/2} = \Psi_{D,x\vert(i,j+1/2,k+1/2)}^{n+1/2}$$ (2.74)

$$\Delta \tilde z \Delta \tilde t H_{z\vert(i,j+1,k+1)}^{n} = \Psi_{H,z\vert(i,j+1,k+1/2)}^{n}$$ (2.75)

Rewriting Equation 2.72 and Equation 2.76 yields:

$$H_{x\vert(i+1/2,j,k)}^{n} = \frac{1}{ \Delta y \Delta z \mu_{\vert(i+1/2,j,k)} } \Phi_{B,x\vert(i+1/2,j,k)}^{n}$$ (2.76)

$$H_{x\vert(i+1/2,j,k)}^{n} = \frac{1}{ \Delta \tilde y \Delta \tilde t } \Psi_{H,x\vert(i+1/2,j,k)}^{n}$$ (2.77)

which can be expressed for $\Phi$ and $\Psi$ . For $\Phi$ it reads:

$$\frac{ \Phi_{B,x\vert(i+1/2,j,k)}^{n} } { \Delta y \Delta z } = \... ...k)} \frac{ \Psi_{H,x\vert(i+1/2,j,k)}^{n} } { \Delta \tilde x \Delta \tilde t }$$ (2.78)

while the result for $\Psi$ is:

$$\frac{ \Psi_{D,x\vert(i,j+1/2,k+1/2)}^{n+1/2} } { \Delta \tilde y... ...2,k+1/2)} \frac{ \Phi_{E,x\vert(i,j+1/2,k+1/2)}^{n+1/2} } { \Delta x \Delta t }$$ (2.79)

Therefore this equation represents a discrete constitutive equation of the simplest type, obtained by extending the local constitutive equations $\ensuremath{\mathbf{B}}=\mu \ensuremath{\mathbf{H}}$ , and $\ensuremath{\mathbf{D}} = \varepsilon \ensuremath{\mathbf{E}}$ with the assumption of planarity, regularity, and orthogonality of the cells.

Examining the time-stepping formulae for $H_x$ , the time-stepping for $\Phi_{B,x}$ becomes:

$$\Phi_{B,x\vert(i+1/2,j,k)}^{n+1} = \Phi_{B,x\vert(i+1/2,j,k)}^{n} +$$ (2.80)

$$\Phi_{E,y\vert(i+1/2,j,k+1/2)}^{n+1/2} - \Phi_{E,y\vert(i+1/2,j,k-1/2)}^{n+1/2} -$$

$$\Phi_{E,z\vert(i+1/2,j+1/2,k)}^{n+1/2} + \Phi_{E,z\vert(i+1/2,j-1/2,k)}^{n+1/2}$$

and from the $E_x$ and the corresponding $\Psi_{D,x}$ is obtained:

$$\Psi_{D,x\vert(i,j+1/2,k+1/2)}^{n+1/2} = \Psi_{D,x\vert(i,j+1/2,k+1/2)}^{n-1/2} -$$ (2.81)

$$\Psi_{H,y\vert(i,j+1/2,k+1)}^{n} + \Psi_{H,y\vert(i,j+1/2,k)}^{n} +$$

$$\Psi_{H,z\vert(i,j+1,k+1/2)}^{n} + \Psi_{H,z\vert(i,j,k+1/2)}^{n}$$

The expressions can then be transformed into equations depending only on two global values:

$$\Phi_{B,x\vert(i+1/2,j,k)}^{n+1} = \Phi_{B,x\vert(i+1/2,j,k)}^{n} + \frac{ \Delta x \Delta t }{ \Delta \tilde y \Delta \tilde z }$$ (2.82)

$$\Biggl( \frac{1}{ \varepsilon_{\vert(i+1/2,j,k+1/2)} } \Psi_{D,x\vert(i+1/2,j,k+1/2)}^{n+1/2} - \frac{1}{ \varepsilon_{\vert(i+1/2,j,k-1/2)} } \Psi_{D,x\vert(i+1/2,j,k-1/2)}^{n+1/2} -$$

$$\frac{1}{ \varepsilon_{\vert(i+1/2,j+1/2,k)} } \Psi_{D,x\vert(i+1/2,j+1/2,k)}^{n+1/2} + \frac{1}{ \varepsilon_{\vert(i+1/2,j-1/2,k)} } \Psi_{D,x\vert(i+1/2,j-1/2,k)}^{n+1/2} \Biggr )$$

For the special case of a dominant magnetic system, TM-mode, the following expressions can be derived:

$$\partial_t H_x = \frac{1}{\mu}(\partial_z E_y - \partial_y E_z )$$ (2.83)

$$\partial_t H_y = \frac{1}{\mu}(\partial_x E_z - \partial_z E_x )$$ (2.84)

$$\partial_t E_z = \frac{1}{\varepsilon}(\partial_x H_y - \partial_y H_x - \gamma E_z)$$ (2.85)

With these expressions the TM-mode at $t=n \Delta t$ is discretized by:

$$\partial_t H_x = \frac{1}{\mu}(\partial_z E_y - \partial_y E_z )$$ (2.86)

$$\partial_t H_y = \frac{1}{\mu}(\partial_x E_z - \partial_z E_x )$$ (2.87)

$$H_{x\vert i,j,k}^{n+1/2} - H_{x\vert i,j,k}^{n-1/2} = \frac{\Delt... ... \Delta y} \left ( E_{z\vert i,j+1/2,k}^{n} - E_{z\vert i,j-1/2,k}^{n} \right )$$ (2.88)

$$H_{y\vert i,j,k}^{n+1/2} - H_{y\vert i,j,k}^{n-1/2} = \frac{\Delt... ... \Delta x} \left ( E_{z\vert i+1/2,j,k}^{n} - E_{z\vert i-1/2,j,k}^{n} \right )$$ (2.89)

For a dominant electric system, the TE-mode is given by:

$$\partial_t E_x = \frac{1}{\varepsilon}(\partial_y H_z - \partial_z H_y - \gamma E_x)$$ (2.90)

$$\partial_t E_y = \frac{1}{\varepsilon}(\partial_z H_x - \partial_x H_z - \gamma E_y)$$ (2.91)

$$\partial_t H_z = \frac{1}{\mu}(\partial_y E_x - \partial_x E_y )$$ (2.92)

With these expressions the TE-mode at $t=(n+1/2) \Delta t$ is discretized by:

$$\partial_t E_z = \frac{1}{\varepsilon}(\partial_x H_y - \partial_y H_x - \gamma E_z)$$ (2.93)

$$E_{z\vert i,j,k}^{n+1} - E_{z\vert i,j,k}^{n} = \frac{\Delta t}{\... ...\vert i,j-1/2,k}^{n+1/2} }{\Delta y} - \gamma E_{z\vert i,j,k}^{n+1/2} \right )$$ (2.94)

As can be seen, the $E_z$ values for the full time-steps and the $H_x$ , $H_y$ for half time-steps are already available, and the the update for $H_x$ , $H_y$ can be done without any further calculation. On the contrary $E_{z\vert i,j,k}^{n+1/2}$ is not available and has to be approximated by:

$$E_{z\vert i,j,k}^{n+1/2} = \frac{1}{2}\left ( E_{z\vert i,j,k}^{n+1} + E_{z\vert i,j,k}^{n} \right )$$ (2.95)

From this the following expression is finally obtained:

$$E_{z\vert i,j,k}^{n+1} = \frac{1-\frac{\gamma \Delta t}{2 \vareps... ...H_{x\vert i,j+1/2,k}^{n+1/2} -H_{x\vert i,j-1/2,k}^{n+1/2} }{\Delta y} \right )$$ (2.96)

As can be seen only the neighboring values $H_x$ , $H_y$ are used to evaluate the spatial derivative of $E_z$ . Also, only the neighboring elements of $H_x$ , $H_y$ and $E_z$ are used to calculate the spatial derivatives of $H_x$ , $H_y$ . Therefore this method calculates both fields, $\ensuremath{\mathbf{E}}$ and $\ensuremath{\mathbf{H}}$ , based on the $\mathrm{curl} \left( \right)$ expressions with special requirements on the given field. Note this discretization can also be derived by the global integral discretization [33], which eases the evaluation of boundary conditions.