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.
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.
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:
A geometric interpretation of the different equations is shown in Figure 2.9.
For second-order derivatives the central difference scheme can be used:
The accuracy of the finite difference approximations is given by:
Mixed derivatives, illustrated in Figure 2.10 can be approximated, e.g., for two dimensions by:
A second-order finite difference approximation for two dimensions results in:
Boundaries also have to be handled by the finite difference method, see Figure 2.11,
but at the boundary on the left side, the backward and central difference approximations would need 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.:
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 . 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 , whereas the time domain is subdivided into intervals of . The original formulation was based on half-step staggered grids in space and time. The quantities from the second complex are denoted with . 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:
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 and 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 
With this notation the following expressions are obtained with central difference approximation:
The time-stepping formulas for are:
and the time-stepping for :
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 part is given):
Rewriting Equation 2.72 and Equation 2.76 yields:
which can be expressed for and . For it reads:
while the result for is:
Therefore this equation represents a discrete constitutive equation of the simplest type, obtained by extending the local constitutive equations , and with the assumption of planarity, regularity, and orthogonality of the cells.
Examining the time-stepping formulae for , the time-stepping for becomes:
and from the and the corresponding is obtained:
The expressions can then be transformed into equations depending only on two global values:
For the special case of a dominant magnetic system, TM-mode, the following expressions can be derived:
With these expressions the TM-mode at is discretized by:
For a dominant electric system, the TE-mode is given by:
With these expressions the TE-mode at is discretized by:
As can be seen, the values for the full time-steps and the , for half time-steps are already available, and the the update for , can be done without any further calculation. On the contrary is not available and has to be approximated by:
From this the following expression is finally obtained:
As can be seen only the neighboring values , are used to evaluate the spatial derivative of . Also, only the neighboring elements of , and are used to calculate the spatial derivatives of , . Therefore this method calculates both fields, and , based on the expressions with special requirements on the given field. Note this discretization can also be derived by the global integral discretization , which eases the evaluation of boundary conditions.