Printer Friendly



Numerical solutions of partial differential equations have grown parallel to the popularization of digital computers. Solutions of Helmholtz equation by the classical second-order accurate finite-difference (FD) method had been investigated [2, 3]. Classical high-accuracy FD approximations to Helmholtz equation can be realized by non-compact stencils at the expense of increased computational costs. Having a compact stencil which connects a point of interest only to its neighboring points is key to the success in the frequency-domain FD method. Compact stencils produce block tri-diagonal matrices which require much less computing resources compared to non-compact stencils.

In last two decades, many high-order compact formulae for discretizing 2-D and 3-D homogeneous Helmholtz equation have been proposed. Based on the method of steepest descent, Jo et al. increased the order of accuracy with a new compact 9-point FD-like formulation for the 2D Helmholtz operator in 1996 [4]. In 1998, Nehrbass et al. improved the classical 3D FD2-7 (second-order accurate, 7-point) formula and derived the new numerical scheme called RD-FD [5]. Although it remains second-order accurate, RD-FD is able to half the dispersion error. It is identical to our LFE-FC-7 stencil in Section 4.1 derived from using a spherical Fourier-Bessel series (SFBS). By replacing the Laplace operator with -[k.sup.2], Singer and Turkel derived 4th-order accurate, compact 2D FD-like stencil for the Helmholtz equation in 1998 [6]. Using similar techniques, they published 6th-order accurate, compact 2D formula in 2006 [7]. Later in 2007, Sutmann extended Singer and Turkel's results and published the 3D compact Helmholtz stencil [8]. The works of Singer-Turkel and Sutmann are based on Taylor series expansion. Both Hadley [9, 10] and Chang-Mu [11, 12] used Fourier-Bessel Series as the basis function for the local field in deriving a 2D compact FD-like stencil which we refer to as LFE2D-9. Along these lines, Tsukerman [13] proposed other choices of basis functions for 2D Helmholtz equation including plane waves and harmonic polynomials. He called this new class of difference schemes FLAME (Flexible Local Approximation Methods). In 2010 Fernandes and Loula derived the sixth-order accurate scheme for 3-D Helmoltz equation using a quasi optimal finite difference method (QOFD) [14]. Their 3D stencil was given (see Equation (6) below) with truncated Taylor series with numerical coefficients. Recently, we derived and published a sixth-order accurate, 3D compact Helmholtz stencil [1]. Like our 2D formula, it is given in closed-form expression as weighted products of spherical Bessel functions (see Equation (10) below).

The method of connected local fields (CLF), which we recently proposed, is a new approach for obtaining semi-analytical solutions of Helmholtz equation. CLF is a mixture of the integral-equation (IE) based and PDE-based methods. In CLF, the entire solution consists of overlapping local fields. The fundamental building block (patch) of the 3-D CLF is a cube consisting of a central point and twenty six points on the cube's surface. Each local field (patch) is represented by a truncated SFBS satisfying Helmholtz equation. CLF connects these patches and forms a system of linear equations just like the FD-FD method.

Using these newly developed methods with the increased computing power, we can systematically apply them to study various EM wave problems even though there are still many difficult issues to be solved. For example, applying a high-order FD-FD method to analyze devices with high dielectric contrasts, such as SOI-based waveguide devices [15], one needs various FD stencils on uniform regions, regions with interfaces, regions with corners, or regions located near computational boundaries (i.e., absorbing boundary condition [16, 17]). The simulation error comes from the imperfections of all the various approximations. Since most waveguide devices are made of blocks of homogeneous materials, the propagation error in homogenous regions must be minimized. In this paper we shall focus on LFE-27 error analysis in free space. This includes Fourier-Bessel Series truncation error analysis of the local field inside a given patch, and the global dispersion characteristics of a given plane wave in the infinite region.

It is well known that the infinite homogeneous region supports isotropic non-dispersive plane wave solutions that satisfy Helmholtz equation. Plane wave solutions also exist for the discretized Helmholtz equation in infinite uniform media. While propagating numerical plane waves like analytic solutions, do not suffer from amplitude attenuation, they are dispersive and also anisotropic. Numerical dispersion means that wave fields of different frequencies propagate at different velocities. The direction of propagation with respect to the grid layout will affect the propagation speed. The latter phenomenon is called numerical anisotropy.

We classify numerical errors from FD-like methods into two types: local and global errors. Global errors are errors of a numerical plane wave propagating in infinite grids. The global error is usually referred as numerical dispersion which has also been called spurious dispersion by Harari and Turkel [18] because the dispersion phenomenon arises completely from the discretization process instead of actual physical phenomena. Numerical errors generated from all possible excitations including inhomogeneous (evanescent) plane waves are called local errors because their propagation is locally contained. The local error includes deviations from the exact analytic solution in both amplitude and phase.

Numerical dispersion analyses of compact formulae on discretizing Helmholtz equation in uniform regions can be found in references [12, 18-23]. In particular, thorough dispersion research on classical FD formulae for 1-D Helmholtz equation were given in Harari and Turkel's work in 1995 [18]. We are the first to derive the dispersion equations of compact formulae for homogeneous 3-D Helmholtz equation. Dispersion characteristics of the classical 7/27-point, second-order accurate formulae FD2-7, FD2-27 [24] and six-order accurate FD6-27 formula [8] will be briefly reviewed. Our focus will be on the in-depth investigation of the numerical phase and group velocity dispersion characteristics of the proposed SFBS-based formulae including three second-order accurate LFE-FC-7, LFE-EC-13, LFE-CR-9 equations, and the main sixth-order accurate LFE-27 formula [1, 25]. Here FC, EC and CR respectively denote face-centered, edge-centered and corner points on the twenty-seven point CLF basic cube/patch.


The general form of a compact 27-point FD-like formula for 3-D homogeneous Helmholtz equation can be expressed as [1, Equation (2a)]:

[mathematical expression not reproducible], (1)

The reduction of Equation (1) has taken into consideration the isotropy of four symmetry groups, the central node [u.sub.0], the six face-centered nodes, the twelve edge-centered nodes and the eight remain nodes on the corners. The summations of these fields on the same group are denoted by [u.sup.f.sub.[SIGMA]], [u.sup.e.sub.[SIGMA]] and [u.sup.c.sub.[SIGMA]] respectively. The three subscripts for the field [] specify its location relative to the central point [u.sub.0]. The plus, minus or zero sign indicates a unit's forward/backward or null displacement ([DELTA]) along the x, y or z direction.

Generally, the FD-like coefficients are functions of the normalized analytic wavenumber V, which is defined as the product of the analytic free-space wavenumber k and the grid spacing [DELTA]. [N.sub.[lambda]] is the field sampling density, number of sampling points per wavelength, given by:

[N.sub.[lambda]] [??] [lambda]/[DELTA] = 2[pi]/V, (2)

where [lambda] = 2[pi]/k is the free space wavelength.

2.1. Second-order Accurate 3D FD Formulae

2.1.1. Classical 7-point Formula: FD2-7

The classical means for discretizing a partial differential equation is to replace the partial derivative terms with a finite difference approximation. With Laplace operation in 3D Helmholtz equation replaced by the classical second-order accurate seven-point formula [2, 3], Equation (1) is discretized as:

[W.sub.f] = 1/[[DELTA].sup.2], [W.sub.e] = 0, [W.sub.c] = 0, [W.sub.0] = [6/[[DELTA].sup.2]] - [k.sup.2]. (FD2-7) (3)

The first number 2 in the above equation signifies that it is second-order accurate and the second number 7 implies that there are a total of 7 nonzero coefficients out of the total 27 3D compact stencil coefficients. This is further illustrated in Fig. 1. Solving problems with the classical FD2-7 needs more than thirteen sampling grids per wavelength to reduce the unwanted numerical dispersion (less than 1% phase velocity error) [9].

2.1.2. 27-point Formula: FD2-27

Using all 27 points, the discrete Laplace operator can be derived similar to how the seven-point formula was derived [24] Replacing Laplace operation in 3D Helmholtz equation by the 27-point FD approximation we have:

[mathematical expression not reproducible] (4)

While FD2-27 equation is of the same order of accuracy as the FD2-7 is, it however, possesses better anisotropy characteristics than the simple 7-point stencil. Note that by setting the spatial frequency k to zero this very FD2-27 formula ascends to a sixth-order accurate discretized Laplace equation [2,248-250].

2.2. High Accuracy FD Formulae

2.2.1. Sutmann's Sixth-order Formula: FD6-27

Using the information from Helmholtz equation itself, Singer and Turkel succeeded in deriving the six-order accurate scheme for 2-D Helmoltz equation [7]. Sutmann extends the work and derived the sixth-order accurate FD-like formula for the 3-D [8] case given by:

[mathematical expression not reproducible]. (FD6-27) (5)

It is worth noting that, for their specific applications, well-behaved source functions were considered in both Sutmann's and SingerTurkel's works. In EM and optical applications, source terms in 2D and 3D Helmholtz equations are usually replaced by incidental fields, allowing us to further improve accuracy of our LFE-based compact stencils.

2.2.2. Fernandes and Loula's Formula: QOFD-6-27

The sixth-order accurate QOFD-6-27 (quasi optimal finite difference) coefficients, given below, were numerically obtained by Fernandes and Loula [14] in 2010 by minimizing a least-squares functional of the local truncation error.

[mathematical expression not reproducible], (QOFD-6-27) (6)

Here [A.sub.1,6], [A.sub.1,8], [A.sub.2,6], [A.sub.2,8], [A.sub.3,6], [A.sub.3,8] are complex expressions made of integer numerators and denominators. This QOFD-6-27 stencil is remarkably accurate when compared to the Taylor series expansion of our LFE-27 analytic formula. They both agree up to fifteen significant digits in the first two powers of the normalized wavenumber V. They begin to differ after the 6th significant digit in all three [W.sub.x]- coefficients of [V.sup.4], [V.sup.6] and after the 7th significant digit in [V.sup.8].

2.3. CLF Formulae

The work on 3D CLF formulae has been detailed in our prior accompanying paper [1] except for the study of their global plane wave dispersion characteristics and local truncation error, which are the main focuses of this paper. We list below the four closely related 3D CLF stencils derived from using a local SFBS.

2.3.1. Face-centered LFE Formula: LFE-FC-7

In 1998, Nehrbass et al. derived a new numerical scheme called RDFD [5]. Their 3D RD-FD is identical to our face-centered SFBS LFE-FC-7 stencil which halves the dispersion error of the classical FD2-7 formula. However, RD-FD remains second-order accurate. Their 3D RD-FD formula is given by:

[W.sub.f] = 1, [W.sub.e] = 0, [W.sub.c] = 0,

[W.sub.0] = [6.sub.j0] (V). (RD-FD or LFE-FC-7) (7)

2.3.2. Edge-centered LFE Formula: LFE-EC-13

The edge-centered based SFBS formula is called LFE-EC-13 and is given by [1, Equation (23a)]:

[W.sub.f] = 0, [W.sub.e] = 1, [W.sub.c] = 0, [W.sub.0] = [12.sub.j0]([square root of 2]V). (LFE-EC-13) (8)

The twelve edge points are shown in Fig. 2.

2.3.3. Corner-point Formula: LFE-CR-9

The LFE-CR-9 stencil containing all corner points illustrated in Fig. 3 is given by [1, Equation (29a)]:

[W.sub.f] = 0, [W.sub.e] = 0, [W.sub.c] = 1, [W.sub.0] = [8.sub.j0]([square root of 3]V). (LFE-CR-9) (9)

2.3.4. Complete 27-point Formula: LFE-27

The coefficients of LFE-27 are obtained from a weighted superposition of LFE-FC-7, LFE-EC-13 and LFE-CR-9, which is derived in the accompanying paper (Ref. [1,Equations (43a)-(43b)], illustrated in Fig. 1). It is given below to be used in subsequent analysis.

[mathematical expression not reproducible]. (LFE-27) (10)


3.1. Derivation of Numerical Dispersion Equation

To analyze the numerical dispersion characteristics of Equation (1) as we have done in previous 2D LFE-9 stencil [12], we consider the plane wave solution defined as:

[mathematical expression not reproducible]. (11)

In Equation (11), we denote the numerical wavenumber by [kappa], instead of k. The zenith angle and azimuthal angle of the plane wave are denoted by [[theta].sub.k] and [[phi].sub.k] respectively. Substitute Equation (11) into Equation (1) we have:

[u.sup.f.sub.[SIGMA],n] = 2 (cos [B.sub.x] + cos [B.sub.y] + cos [B.sub.z]), (12a)

[u.sup.e.sub.[SIGMA],n] = 4 (cos [B.sub.x] cos [B.sub.y] + cos [B.sub.y] cos [B.sub.z] + cos [B.sub.z] cos [B.sub.x]), (12b)

[u.sup.c.sub.[SIGMA],n] = 8 cos [B.sub.x] cos [B.sub.y] cos [B.sub.z]. (12c)


[mathematical expression not reproducible]. (13)

The notations [S.sub.x], [S.sub.y], and [S.sub.z] in Equation (13) represent the directional cosines of the numerical wavenumber [??], hence:

[S.sup.2.sub.x] + [S.sup.2.sub.y] + [S.sup.2.sub.z] = 1. (14)

"Numerical" plane wave solutions are then evaluated on the grids and summed over three symmetry groups. The results are denoted by [u.sup.f.sub.[SIGMA],n], [u.sup.e.sub.[SIGMA],n] and [u.sup.c.sub.[SIGMA],n] respectively. To keep the spherical symmetry intact, it is better to choose [S.sub.x], [S.sub.y], and [S.sub.z] rather than [[theta].sub.k] and [[phi].sub.k] because the latter are associated with the spherical coordinate system which uses a z-axis as a referenced direction, breaking the directional symmetry. Substituting Equations (12a)-(12c) back into Equation (1) results in the following numerical dispersion equation for 3D discretized Helmholtz equation using a general 27-point FD-like scheme:

[mathematical expression not reproducible] (15)

The four W's {[W.sub.x], x = 0, f, e, c} in Equation (15) are functions of the normalized analytic wavenumber V, whereas the three [u.sub.[SIGMA],n]'s {[u.sup.x.sub.[SIGMA],n], x = f, e, c} are functions of B, [S.sub.x], [S.sub.y], and [S.sub.z]. Hence, the normalized numerical wavenumber B is an implicit function of V and [S.sub.x], [S.sub.y], [S.sub.z]. If the FD-like scheme is perfect, B would linearly depend on V and would be independent of [S.sub.x], [S.sub.y], and [S.sub.z]. However, except for the 1D case, it is impossible to find a perfect scheme in either 2D or 3D cases. Now we are ready to investigate the numerical dispersion equation (Equation (15)) for various FD-like formulae, including the temporal dispersion analysis and the spatial dispersion (numerical anisotropy) analysis.

3.2. Numerical Phase Velocity and Group Velocity

To quantitatively investigate dispersion characteristics, we consider the numerical phase velocity [] as defined below:

[] [??] [omega]/[kappa] = [ck x [DELTA]]/[[kappa] x [DELTA]] [equivalent to] [V/B] x c. (16)

Here, [omega] represents the angular frequency and c is the speed of light in vacuum. The numerical group velocity [[??]] defined below is, in general, anisotropic.

[mathematical expression not reproducible], (17)

where [omega] is considered as an implicit function of the normalized numerical wavenumber [??]. Since [omega] is linearly proportional to V which is an implicit function of [??], we have the following constructive equation for computing the magnitude of the numerical group velocity:

[mathematical expression not reproducible]. (18)

The relative error of [] and [] are defined as:

[mathematical expression not reproducible]. (19)

3.3. Surface of Non-degenerate Directions

For the 3D dispersion study, we need to first determine the minimum set of propagation directions in the k space (wavenumber space) which is illustrated in Fig. 4 and Fig. 5(a). Note that the labels and captions for the [k.sub.x], [k.sub.y] and [k.sub.z] axes of the k space are shorten as x, y and z to simplify the drawings of both figures. We refer to the surface of non-degenerate directions (NDD) as the surface of intersecting points between the NDD and the unit sphere. To find this surface of NDD we first apply the spatial symmetry to constrain the surface to the first octant. Furthermore by closer examining the three [u.sub.[SIGMA],n]'s in Equations (12a)-(12c), we see that we can exchange the placement orders of [S.sub.x], [S.sub.y], [S.sub.z] without changing the equations. Hence we have even symmetry with respect to both [k.sub.x] = [k.sub.y] and [k.sub.x] = [k.sub.z] planes. The equations that set the NDD domain are given by:

[mathematical expression not reproducible]. (20)

This surface of NDD is further illustrated by the red triangular-like patch in Fig. 5(a), where we show one of the NDD elements in the 3-D free space by the unit vector [??]. It is one sixth of the surface area of an octant and hence 1/48 of the entire unit sphere. This minimal set is constructed to simplify the analysis--every vector within this domain will occur 47 more times in the unit sphere. In addition, seven representative directions are chosen from the NDD set for our dispersion analysis. The definitions are given in Table 1 and illustrated in Fig. 5(b). Points Ai, A2 and A3 are intersections of vectors [[??].sup.1], [[??].sup.2] and [[??].sup.3] with the NDD surface. Each is parallel to the direction (1, 0, 0), (1, 1, 0) and (1, 1, 1) respectively. Three other directions marked by the intersection points [M.sub.1], [M.sub.2] and [M.sub.3] are also chosen to be, respectively, the midpoints of arc [mathematical expression not reproducible], and [??]. The addition of the centroidal point Q of this NDD surface completes the final selection.

Seven chosen propagation directions are listed in Table 1.


4.1. Numerical Study of Relative Phase/Group Velocity Dispersion Error

In general the B-V equation (Equation (15)) is highly nonlinear and is often very difficult to solve from either a given B to V or from a given V to B. In Section 4.2 we will apply the first-order analysis of the B-V relation for the LFE-27 stencil as we did previously for the 2D LFE-9 case [12]. We can only perform this type of analysis with highly accurate compact stencils or when B and V are both small. For greater values of B, V and for other low-order compact stencils we use numerical techniques, including trial and error methods, to solve the nonlinear B-V equation, Equation (15). This will be the main focus in this section.

4.1.1. FD2-7 Stencil

The relative phase and group veloctiy errors as functions of the sampling density for the classical FD2-7 stencil are plotted in the left and right part of Fig. 6. Each subplot shows seven curves computed numerically from the nonlinear [??]-V relations defined in Equations (16)-(17). They appear to be linear in log-log scale which is consistent with the [V.sup.2] convergent characteristics of the dispersion error. The wide spread between the lines is an indication of large numerical anisotropy of this classical formula. Note that the sampling density required to ensure less than one percent relative error is 13 and 22 for the phase and group velocity respectively.

4.1.2. FD2-27 Stencil

In Fig. 7 the relative phase and group velocity errors as a function of the sampling density for the FD2-27 stencil are ploted in the same manner as in Fig. 6. The errors curves also appear to be linear in log-log scale. The lack of spreading between the lines is an indication of trivial numerical anisotropy of this FD2-27 formula. Its performance is on par with FD2-7.

4.1.3. Sutmann's FD6-27 Stencil

The curves of relative phase and group veloctiy errors for the 6th-order accurate Sutmann's FD6-27 stencil are shown in Fig. 8. They also appear to be linear in log-log scale except at a much steeper slope due to its 6th-order accuracy. Some numerical anisotropy is observed for this FD6-27 formula. The same one percent relative error thresholds for the required sampling density are around 3 and 4 for the phase and group velocity respectively.

4.1.4. LFE-FC-7 (3-D RD-FD) Stencil

4.1.5. LFE-EC-13 Stencil

4.1.6. LFE-CR-9 Stencil

The relative phase and group velocity errors for the three second-order accurate LFE stencils are shown in Figs. 9-11. Due to the complex nature of LFE stencil coefficients, the relative error curves appear to be quasi-linear in log-log scale. We include them for the sake of being thorough since with the exception of LFE-FC-7 stencil, the other edge-centered or corner-based stencil can not be used individually as we reported in the accompanying paper [1]. One common fact shared by all three second-order accurate LFE stencils is that the most dispersive directions are along one of the Cartesian axis [mathematical expression not reproducible] and the least dispersive direction is along the direction that passes through the centroid of the NDD surface.

4.1.7. LFE-27 Stencil

The relative phase and group velocity errors for the sixth-order accurate LFE-27 stencil are plotted in Fig. 12. The summary plot comparing all four 3D compact stencils is shown in Fig. 13. From these curves we see that the 3D LFE-27 stencil possesses the best numerical quality in reducing phase and group velocity dispersion errors. Reading straight off the curves of Fig. 12, we see that the required sampling density to ensure less than one percent relative phase error is [N.sub.[lambda]] [greater than or equal to] 2.5 and [N.sub.[lambda]] [greater than or equal to] 2.7 for the group velocity. In order to be able to graph and compare the dispersion errors among all four stencils, we chose the worst possible propagation direction across all formulae in Fig. 13. The data shows that, with some calculation, second-order accuracy is observed for the classical FD2-7 and LFE-FC-7 stencils. We can also numerically verify that both FD6-27 and LFE-27 are sixth-order accurate. As for the discrepancy between the two, the most obvious difference is that LFE-27 stencil enjoys a more-than 12 dB relative error advantage over the FD6-27 formula from the log-log plot.

4.2. First-order Dispersion Analysis of LFE-27 Formula

Having reviewed the numerical studies of the dispersion characteristics of various 3D compact stencils, we will now look into the analytic characteristics of the LFE-27 formula. From numerical evidence, we know that the LFE-27 B-V relation is almost linear due its high-order accuracy. We wish to obtain an analytic formula for the small difference between the two normalized wavenumbers, i.e., we seek for the relative difference [[epsilon].sub.LFE-27] defined below:

[[epsilon].sub.LFE-27] = [B - V]/V. (21)

Hence, we may write the normalized numerical wavenumber B as V(1 + [[epsilon].sub.LFE-27]). Keeping only the first power of [[epsilon].sub.LFE-27] we shall refer the following analytic investigation as the first-order dispersion analysis.

Rewriting Equation (11) in terms of directional cosines, we have:

[mathematical expression not reproducible]. (22)

Using this new notation together with the trigonometry sum and difference identities, the three group sums resulting from a "numerical" plane wave solution are rewritten as:

[u.sup.f.sub.[SIGMA],n] = 2[cos (B[S.sub.x]) + cos (B[S.sub.y]) + cos (B[S.sub.z])], (23a)

[mathematical expression not reproducible], (23b)


[mathematical expression not reproducible]. (23c)

Without going into too much detail, in the remaining section we will follow the logic of deriving [[epsilon].sub.LFE-9], the first-order analytic dispersion formula [12, Equation (43)] for the 2D LFE-9 stencil. To further simplify Equations (23a)-(23c) we need to expand the cosine functions with a cosine function as the argument. For 3D analysis we consider the very similar expression given in Equation (56) of Stratton 1941, page 409 which can also be found in Equation (11)-(85) of Ishimaru [26].

[mathematical expression not reproducible]. (24)

Applying Equation (A6) to Equations (23a)-(23c), the three numerical group sums [u.sub.[SIGMA],n]'s can be explicitly written as:

[mathematical expression not reproducible], (25a)

[mathematical expression not reproducible], (25b)


[mathematical expression not reproducible]. (25c)

The significance of using [[GAMMA].sup.x.sub.2n] (referred to as the "directional" function) in above equations is that we want to separate the radial and the directional dependence of [u.sup.x.sub.[SIGMA],n] without making it more complex. [[GAMMA].sup.x.sub.2n] is a linear combination of Legendre polynomials whose arguments contain various linear combinations of directional cosines. The advantage of choosing directional cosines as the independent variables in "directional" functions is that it invokes only the regular Legendre polynomials. Otherwise, had we chosen the zenith and azimuthal angles of the spherical coordinate system we would have to deal with additional associated Legendre functions. In addition, we are surprised to find out that the three "directional" functions are almost identical except for the proportionality constants ([[alpha].sup.x.sub.2n], x = f, e, c) as shown below:

[[GAMMA].sup.x.sub.2n] = [[alpha].sup.x.sub.2n] x [[xi].sub.2n] ([S.sub.x], [S.sub.y]), x = f, e, c, n = 1, 2, .... (25d)


[mathematical expression not reproducible], (25e)


[mathematical expression not reproducible]. (25f)

This implies that the directional dispersion behavior of the three second-order accurate stencils LFE-FC-7, LFE-EC-13 and LFE-CR-9 are becoming increasingly indistinguishable as [N.sub.[lambda]] [right arrow] [infinity].

Omitting the details of complex algebraic manipulation, we show that by the principle of first-order analysis, the relative LFE-27 B-V error can be derived from the following intermediate form:

[mathematical expression not reproducible]. (26)

Here {[u.sup.x.sub.[SIGMA],a], x = f, e, c} stands for one of the three group sums of analytic plane wave defined as:

[mathematical expression not reproducible], (27a)

[mathematical expression not reproducible], (27b)


[mathematical expression not reproducible]. (27c)

Note that the specific form of the "analytical" [u.sup.x.sub.[SIGMA],a] a of Equations (27a)-(27c) is exactly the same as the "numerical" [u.sup.x.sub.[SIGMA],n] given in Equations (25a)-(25c) except that the normalized numerical wavenumber B in the argument of spherical Bessel functions is replaced by the normalized wavenumber V.

Substituting Equations (27a)-(27c) and Equation (10) into Equation (26), and after skipping half a dozen intermediate steps, the final first-order formula for the relative LFE-27 B-V error is given below:

[mathematical expression not reproducible]. (28)

We first note that [[epsilon].sub.LFE-27] is of the sixth power of V and it is multiplied by a small [10.sup.-6] factor. For small normalized wavenumbers the directional behavior is dictated by [[xi].sub.8], the eighth-order of the directional function given in Equation (A12). If we were to compare the first-order analytic dispersion error formula of 3D [[epsilon].sub.LFE-27] with the 2D [[epsilon].sub.LFE-9] (which is given below) we see that they are both of the 6th power in V. They have similar 8th-order directional dependence in the numerators, i.e., cos 8[theta] for the 2D and [[xi].sub.8] for the 3D LFE formulae. The most significant difference is that the error for 3D is about twice that of 2D.

[mathematical expression not reproducible]. (29)


Numerical errors are classified into two types: local and global errors. The dispersion error is considered to be the global error because it represents the deviation of the phase of a numerical plane wave propagating in an infinite grid space. The local error, however, includes deviation in both amplitude and phase from the exact analytic solution. The LFE local field error refers to the LFE-interpolation error of the 2D/3D field (sitting in the central location of a 2D/3D square/cube) from its 8/26 enclosing points using LFE-9/LFE-27 formula. The local truncation function [T.sub.LFE-27] for LFE-27 stencil is given by:

[T.sub.LFE-27] [??] [u.sub.0] - [u.sup.LFE-27.sub.0], (30)

where LFE-27 interpolated field [u.sup.LFE-27.sub.0] is given by:

[mathematical expression not reproducible]. (31)

The coefficients W's {[W.sub.x], x = f, e, c} are given in Equation (10). The exact local field can be written in terms of an infinite series of products of the spherical Bessel function and the spherical harmonics as shown below:

[mathematical expression not reproducible]. (32)

From the orthogonal equality of the spherical harmonics function, we may obtain the unknown coefficients with the following integrals:

[mathematical expression not reproducible], (33a)

[mathematical expression not reproducible], (33b)

[mathematical expression not reproducible]. (33c)

In Appendix A, we document the detailed derivation of [T.sub.LFE-27] based on the assumption that the local field is from either a propagating plane wave or some linear combination of propagating plane waves. The main result is given below:

[mathematical expression not reproducible]. (34)

We see that the leading term is proportional to [10.sup.-7] x [V.sup.8]. It is also multiplied by the weighted sum of three 8th-order spherical harmonics coefficients {[a.sup.k.sub.8], k = 0, 4, 8}. In order to further understand the behavior of high-order spherical harmonics coefficients, please see Appendix B, where we dive into further detail regarding the spherical Fourier-Bessel expansion of a plane wave. With these coefficients, we are able to numerically evaluate the magnitudes [a.sup.0.sub.8], 1680[a.sup.4.sub.8] and 2620800[a.sup.8.sub.8] revealing that these three terms are of the order of ten for a propagating plane wave in an arbitrary direction (see Fig. 14 and Equations (B4a)-(B4b)) Hence, the local truncation error of the LFE-27 interpolated field is approximately given by:

[T.sub.LFE-27] [approximately equal to] [10.sup.-6] x [V.sup.8]. (35)

In the worst-case scenario, the local truncation error for LFE-27 evaluated at the lowest sampling density [N.sub.[lambda]] = 2 (corresponding to a maximal V = [pi]) is [10.sup.-6] x [[pi].sup.8] [approximately equal to] 0.01. The error will, of course, increase for local fields containing evanescent waves when the EM field is inadequately sampled, i.e., when [N.sub.[lambda]] [much less than] 2.


In our accompanying paper [1] we analytically derived the LFE-27 stencil for the 3D Helmholtz operator in a homogeneous medium. We also computed, using LFE-27 stencil, 3-D Green's functions for accumulated amplitude and phase errors. In this paper, we first numerically investigate the global dispersion behavior of an arbitrary plane wave in the seven selected directions. We also derive, based on first-order analysis, [[epsilon].sub.LFE-27] which is an analytical expression for obtaining the normalized numerical wavenumber B from the normalized wavenumber V. Toward the end we derive the local truncation error function [T.sub.LFE-27] and show that LFE-27 interpolated fields from its 26 immediate surrounding fields is accurate to the sixth power of V.

Solving 2D Helmholtz problems with FD-FD-like methods are practical with modern computing resources. Solving 3-D Helmholtz problems by FD-FD methods require solving significantly larger matrix equations. If N is the number of unknowns in one dimension, a 3D FD-FD computation using a direct solver will respectively require storage and computational costs proportional to [N.sup.5] and [N.sup.7].

There are still many CLF related issues to be solved. We have considered the homogeneous case for both the 2D and 3D Helmholtz equation. These formulations extend naturally to continuous, smoothly-varying media. For modeling other general passive EM-optical waveguides/devices, we have to consider EM wave fields in media with discontinuous material properties. We need to develop special 2D and 3D LFE stencils for cells near a dielectric interface. We also have to consider a vector formulation for 3D structures with sudden changes of dielectric indices. The method of connected local fields is still in an early stage of development, and requires further research before it becomes a viable tool for modeling 2D and 3D complex dielectric devices.


In this paper, we investigate the dispersion characteristics and conduct local-error analysis of 3-D compact SFBS-based formulae of LFE-FC-7, LFE-EC-13, LFE-CR-9 and LFE-27. Comparing LFE-FC-7 with the classical FD2-7 equation, we found that the former produces about half the errors of the latter stencil. Although both are sixth-order accurate in V, when comparing LFE-27 with Sutmann's FD6-27 equation, we found that the LFE-27 stencil has a 12 dB advantage over FD6-27 in both relative phase and group velocity errors. Furthermore, the classical second-order FD formula requires more than twenty sampling points per wavelength to achieve less than 1% relative error in both phase and group velocities, whereas LFE-27 needs only three points per wavelength to match the same performance. A mere factor-of-seven reduction in the EM field sampling density leads to a whopping reduction by a factor of 800,000 in CPU computation times when running a 3D FD-FD computation with a direct linear equation solver.


We are grateful for the support of the National Science Council of the Republic of China under the contracts 101-2221-E-110-073.


The truncation functions for the three second-order accurate LFE stencils are given below in this appendix. The three O functions, {[O.sub.x]([V.sup.8]), x = f, e, c}, were given in Equations (33a)-(33c) of Ref. [1]. They contained only terms up to [V.sup.6] because the eighth-power terms were not needed to derive the 3D LFE-27 stencil. They are all given below with the [V.sup.8] terms in order to obtain the analytic local truncation function [T.sub.LFE-27] for the LFE-27 stencil.

A.1. LFE-FC-7 Truncation Function

[mathematical expression not reproducible] (A1)

[mathematical expression not reproducible], (A2)

[mathematical expression not reproducible]. (A3)

A.2. LFE-EC-13 Truncation Function

[mathematical expression not reproducible] (A4)

[mathematical expression not reproducible]. (A5)

[mathematical expression not reproducible]. (A6)

A.3. LFE-CR-9 Truncation Function

[mathematical expression not reproducible] (A7)

[mathematical expression not reproducible], (A8)

[mathematical expression not reproducible]. (A9)

A.4. LFE-27 Truncation Function

Equations (A3), (A6) and (A9) all contain [V.sup.8] terms for the three LFE truncation functions. Given these we are able to determine [O.sub.LFE-27], the LFE truncation function as defined below:

[mathematical expression not reproducible]. (A10)

The local truncation function [T.sub.LFE-27] for the LFE-27 stencil is given in Equation (30) and the numerical value [u.sup.LFE-27.sub.0] calculated by LFE-27 is given in Equation (31). The sixth-order accurate LFE-27 stencil was derived using four "special" weighting factors [W.sub.x] so that the LFE-27 truncation error function would not contain any [V.sup.4], [V.sup.6] terms. Hence, the leading term in [O.sub.LFE-27] will be proportional to [V.sup.8]. The exact equation for [O.sub.LFE-27] is given below:

[mathematical expression not reproducible], (A11)

where the W coefficients for the [O.sub.LFE-27] function are the same as the W coefficients for the LFE-27 formulae given by Equation (10). The W coefficients can be normalized by setting the coefficient of the central field [u.sub.0] to 1. The normalized W coefficients are given by:

[mathematical expression not reproducible] (A12)

The local truncation error function of the LFE-27 formula can be obtained by applying the three [O.sub.x] functions (given Equations (A3), (A6) and (A9)) to Equation (A11). The final result is given below:

[mathematical expression not reproducible]. (A13)


We begin with Equation (24) where a propagating plane wave along the z-axis is expanded in infinite SFBS of (Equation (32)). The extension to a plane wave in an arbitrary direction {[[theta].sub.k], [[phi].sub.k]} can be obtained by noting that (See Ref. [26, page 146, Equation (5)-(130)]):

[mathematical expression not reproducible]. (B1)

Hence, we have:

[mathematical expression not reproducible], (B2a)


[mathematical expression not reproducible]. (B2b)

From the two above equations, we can determine the explicit form of SFBS coefficients corresponding to a plane wave. Namely,

[mathematical expression not reproducible], (B3a)

[mathematical expression not reproducible], (B3b)

[mathematical expression not reproducible]. (B3c)

In Equations (B3b)-(B3c), the index m of the associated Legendre polynomial runs from 1 to l. Now we are ready to examine the magnitude of coefficients [a.sup.0.sub.8], [a.sup.4.sub.8] and [a.sup.8.sub.8] in the local truncation function of Equation (34).

[mathematical expression not reproducible], (B4a)

[mathematical expression not reproducible]. (B4b)


[1.] Chang, H.-W. and S.-Y. Mu, "Semi-analytical solutions of the 3D Homogeneous Helmholtz equation by the method of connected local fields," Progress In Electromagnetics Research, Vol. 142, 159-188, 2013.

[2.] Smith, G. D., Numerical Solution of Partial Differential Equations, 2nd Edition, Oxford University Press, 1978.

[3.] Hall, C. A. and T. A. Porsching, Numerical Analysis of Partial Differential Equations, Prentice-Hall, Englewood Cliffs, New Jersey, 1990

[4.] Jo, C.-H., C. Shin, and J. H. Suh, "An optimal 9-point, finite-difference, frequency-space, 2-D scalar wave extrapolator," Geophysics, Vol. 61, No. 2, 529-537, 1996.

[5.] Nehrbass, J. W., J. O. Jevtic, and R. Lee, "Reducing the phase error for finite-difference methods without increasing the order," IEEE Transactions on Antennas and Propagation, Vol. 46, No. 8, 1194-1201, 1998.

[6.] Singer, I. and E. Turkel, "High-order finite difference method for the Helmholtz equation," Computer Methods in Applied Mechanics and Engineering, Vol. 163, 343-358, 1998.

[7.] Singer, I. and E. Turkel, "Sixth order accurate finite difference schemes for the Helmholtz equation," Journal of Computational Acoustics, Vol. 14, 339-351, 2006.

[8.] Sutmann, G., "Compact finite difference schemes of sixth order for the Helmholtz equation," Journal of Computational and Applied Mathematics, Vol. 203, 15-31, 2007.

[9.] Hadley, G. R., "High-accuracy finite-difference equations for dielectric waveguide analysis I: Uniform regions and dielectric interfaces," Journal of Lightwave Technology, Vol. 20, No. 7, 1210-1218, 2002.

[10.] Hadley, G. R., "High-accuracy finite-difference equations for dielectric waveguide analysis II: Dielectric corners," Journal of Lightwave Technology, Vol. 20, No. 7, 1219-1231, 2002.

[11.] Chang, H.-W. and S.-Y. Mu, "Semi-analytical solutions of the 2D Homogeneous Helmholtz equation by the method of connected local fields," Progress In Electromagnetics Research, Vol. 109, 399-424, 2010.

[12.] Mu, S.-Y. and H.-W. Chang, "Theoretical foundation for the method of connected local fields," Progress In Electromagnetics Research, Vol. 114, 67-88, 2011.

[13.] Tsukerman, I., "Electromagnetic applications of a new finite-difference calculus" IEEE Transaction on Magnetics, Vol. 41, No. 7, 2206-2225, 2005

[14.] Fernandes, D. T. and A. F. D. Loula, "Quasi optimal finite difference method for Helmholtz problem on unstructured grids," Int. J. Numer. Meth. Engng., Vol. 82, 1244-1281, 2010.

[15.] Chang, H.-W. and Y.-H. Wu, "Analysis of perpendicular crossing dielectric waveguides with various typical index contrasts and intersection profiles," Progress In Electromagnetics Research, Vol. 108, 323-341, 2010.

[16.] Engquist, B. and A. Majda, "Absorbing boundary conditions for numerical simulation of waves," Applied Mathematical Science, Vol. 74, 1765-1766, 1977.

[17.] Chang, H.-W., W.-C. Cheng, and S.-M. Lu, "Layer-mode transparent boundary condition for the hybrid FD-FD method," Progress In Electromagnetics Research, Vol. 94, 175-195, 2009

[18.] Harari, I. and E. Turkel, "Accurate finite difference methods for time-harmonic wave propagation," Journal of Computational Physics, Vol. 119, No. 2, 252-270, 1995.

[19.] Trefethen, L. N., "Group velocity in finite difference schemes," SIAM Rev, Vol. 24, No. 2, 113-136, 1982.

[20.] Anne, L. and Q. H. Tran, "Dispersion and cost analysis of some finite difference schemes in one-parameter acoustic wave modeling," Computational Geosciences, Vol. 1, 1-33, 1997.

[21.] Peterson, A. F., S. L. Ray, and R. Mittra, Computational Method for Electromagnetics, IEEE Press, New York, 1998.

[22.] Rao, K. R., J. Nehrbass, and R. Lee, "Discretization errors in finite methods: Issues and possible solutions," Comput. Methods Appl. Mech. Engrg., Vol. 169, 219-236, 1999.

[23.] Taflove, A. and S. C. Hagness, Computational Electrodynamics: The Finite-difference Time-domain Method, 3rd Edition, Artech House, Norwood, MA, 2005.

[24.] Spotz, W. F. and G. F. Carey, "A high-order compact formulation for the 3D Poisson equation," Numerical Methods for Partial Differential Equations, Vol. 12, No. 2, 235-243, 1996.

[25.] Chang, H.-W. and S.-Y. Mu, "3-D LFE-27 formulae for the method of connected local fields," Optics & Photonics Taiwan, International Conference, Taipei, Taiwan, Dec. 6-8, 2012.

[26.] Ishimaru, A., Electromagnetic Propagation, Radiation, and Scattering, Prentice Hall, Englewood Cliffs, NJ, 1991.

Sin-Yuan Mu (1,2) and Hung-Wen Chang (1,2), *

(1) Institute of Electro-optical Engineering, National Sun Yat-sen University, Kaohsiung 80424, Taiwan, R.O.C.

(2) Department of Photonics, National Sun Yat-sen University, Kaohsiung 80424, Taiwan, R.O.C.

Received 1 September 2013, Accepted 30 October 2013, Scheduled 8 November 2013

* Corresponding author: Hung-Wen Chang (

Caption: Figure 1. Illustration of a compact 7-point stencil layout. The central point is shown in black while all face-centered points are shown in red.

Caption: Figure 2. Illustration of a compact 13-point stencil layout. The central point is shown in black while all edge-centered points are shown in green.

Caption: Figure 3. Illustration of a compact 9-point stencil layout. The blue points are all eight corners of a cubic cell.

Caption: Figure 4. Drawing of plane of symmetry used in the derivation of non-degenerate directions (in the k space). The space [S.sub.0] represents the spherical surface of the unit sphere. The colorless [E.sub.i], i = 1, 2, 3 plane restricts the domain to the first octant. The green [S.sub.1] plane and the orange [S.sub.2] intersect at the blue [k.sub.x] = [k.sub.y] = [k.sub.z] line. [E.sub.2], [S.sub.1], [S.sub.2] together with [S.sub.0] define the set of non-degenerate directions which is further illustrated in Fig. 5.

Caption: Figure 5. (a) Drawing of the set of non-degenerate directions (in the k space). (b) Intersection of the red surface in Fig. 5(a) with seven selected directions of interest listed in Table 1. These points are three vertices denoted by [A.sub.1,2] and [A.sub.3] and three midpoints [M.sub.1], [A.sub.2] and [A.sub.3] between {[A.sub.i]}. The centroid of this NDD surface is denoted by Q.

Caption: Figure 6. (a) Relative phase and (b) group velocity errors of the classical FD2-7 stencil. Each of the seven curves is ploted as a function of the sampling density for one of the selected directions as defined in Table 1.

Caption: Figure 7. (a) Relative phase and (b) group velocity errors of the classical FD2-27 stencil. Each curve is ploted as a function of the sampling density for a fixed direction (defined in Table 1).

Caption: Figure 8. (a) Relative phase and (b) group velocity errors of Sutmann's FD6-27 stencil.

Caption: Figure 9. (a) Relative phase and (b) group velocity errors of the 3D RD-FD (LFE-FC-7) stencil.

Caption: Figure 10. (a) Relative phase and (b) group velocity errors of the LFE-EC-13 stencil.

Caption: Figure 11. (a) Relative phase and (b) group velocity errors of the LFE-EC-13 stencil.

Caption: Figure 12. Relative phase and group velocity errors of the 3D LFE-27 stencil.

Caption: Figure 13. Phase velocity dispersion performance comparison of various 3D compact stencils.

Caption: Figure 14. Plot of three weighted 8th-order spherical harmonics coefficients in Equation (34).
Table 1. Seven selected directions for subsequent analysis.

Direction     [[theta].sub.k]     [[phi].sub.k]

   [??]         0[degrees]              X
   [??]         45[degrees]        0[degrees]
   [??]       54.74[degrees]       45[degrees]
   [??]        22.5[degrees]       0[degrees]
   [??]       47.63[degrees]     24.20[degrees]
   [??]       27.37[degrees]       45[degrees]
   [??]       31.65[degrees]     24.20[degrees]

Direction                    Property

   [??]                 (0, 0, 1) direction
   [??]                 (1, 0, 1) direction
   [??]                 (1, 1, 1) direction
   [??]        [M.sub.1] is the middle point of [??]
   [??]        [M.sub.2] is the middle point of [??]
   [??]        [M.sub.3] is the middle point of [??]
   [??]           Q is the intersection point of
                   [mathematical expression not
                      reproducible], and [??]
COPYRIGHT 2013 Electromagnetics Academy
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2013 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Author:Mu, Sin-Yuan; Chang, Hung-Wen
Publication:Progress In Electromagnetics Research
Date:Nov 1, 2013

Terms of use | Privacy policy | Copyright © 2021 Farlex, Inc. | Feedback | For webmasters |