1810 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 46, NO. 12, DECEMBER 1998 A Recursive Green s Function Method for Boundary Integral Analysis of Inhomogeneous Domains Michael A. Jensen, Member, IEEE, and Jim D. Freeze, Student Member, IEEE Abstract The recursive Green s function method (RGFM) for computation of fields scattered by two-dimensional (2-D) inhomogeneous dielectric bodies is presented. The algorithm efficiently constructs the Green s function for the inhomogeneous region by recursively combining known Green s functions from smaller subdomains. The fields on the scatterer surface are then computed using a boundary integral formulation. Proper implementation of the RGFM results in computational and storage complexities which scale as and , respectively, where is the total number of discrete cells in a domain. Comparisons of results obtained using the RGFM with those computed from moment method and exact solutions show the efficiency and accuracy of the technique. Index Terms Boundary integral equations, dielectric scatterer, integral equations, scattering. I. INTRODUCTION DETERMINING the behavior of electromagnetic fields in and near inhomogeneous dielectric domains finds significant practical application in such areas as inverse scattering and imaging [1], [2], target identification, waveguide analysis, biomedical research [3] [6], and antenna design [7], [8]. Consequently, considerable research concerning efficient computational methods for modeling inhomogeneous obstacles has appeared in the literature. The majority of these methods are based upon either the finite-difference technique, the finiteelement algorithm, or the method of moments (MoM), with each technique offering a different set of advantages and disadvantages. When applying these techniques to large inhomogeneous domains, the computational and storage costs of the algorithm become key considerations. This is particularly true for integral equation solutions via the MoM, which, for unknowns, typically have computational and storage complexities of and , respectively. Several methods have been developed that significantly reduce these costs. For example, the conjugate gradient fast Fourier transform (CG-FFT) algorithm [5] and the stacked spectral iterative technique (SIT) [6] exploit the convolutional properties of the volume integral equation to achieve an efficient solution for single-source configurations. For simulations requiring multiple source configurations, recursive techniques such as the nested equivalence principle algorithm (NEPAL) [9], [10] are highly effective. Manuscript received June 5, 1996; revised May 2, 1997. The authors are with the Department of Electrical and Computer Engineering, Brigham Young University, Provo, UT 84602 USA. Publisher Item Identifier S 0018-926X(98)09688-4. Fig. 1. Two-dimensional domain divided into unit sections. Recently, a new approach known as the recursive Green s function method (RGFM) has been introduced for analysis of one-dimensional (1-D) inhomogeneous domains [11] [13]. This method is similar to NEPAL for the volume integral equation and the nested dissection approach for finite-element analysis [14]. However, the RGFM is unique in that it constructs the Green s function for the inhomogeneous domain using an efficient recursive procedure. Boundary integral techniques can then be used to determine the field behavior around the obstacle for different illumination configurations. This 1-D RGFM has been applied to quasi two-dimensional (2-D) optical waveguide structures [15], but the rigorous generalization of the technique to 2-D geometries has not yet appeared. In this paper, we extend the RGFM to allow solution of the scalar wave equation for 2-D inhomogeneous domains. The approach is found to have computational and storage complexities of and , respectively, representing a considerable savings over traditional volume MoM solutions. Additionally, because the method produces the obstacle Green s function, it can be used to conveniently couple inhomogeneous domains to radiators and other scatterers. Computational examples are presented which demonstrate that for some geometries the RGFM offers improved accuracy over the volume MoM. Additional examples illustrate the flexibility and computational efficiency of the approach. II. THEORY Consider a 2-D space that contains an inhomogeneous region bounded by the contour and surrounded by a homogeneous region , as shown in Fig. 1. The space is characterized by the complex relative permittivity , where is constant for . Given this domain 0018 926X/98$10.00 1998 IEEE Authorized licensed use limited to: Brigham Young University. Downloaded on February 3, 2009 at 11:19 from IEEE Xplore. Restrictions apply. JENSEN AND FREEZE: RECURSIVE GREEN S FUNCTION METHOD FOR ANALYSIS OF INHOMOGENEOUS DOMAINS 1811 description, we seek to determine the solution , , to the scalar wave equation given as (1) where is the free-space wavenumber. Solution of (1) for a given incident field is often achieved using a MoM discretization of the volume integral equation [16]. The difficulty with this approach is that for very large domains, the required computational and storage resources can become overly intensive for modern computational platforms. An alternative approach to solving (1) involves first determining the Green s function for , which satisfies (2) where represents the Dirac delta function. Once is known, then coupled boundary integral equations for the fields on can be formulated and solved using a surface MoM. Such a technique offers the advantage of significantly reduced computational and storage costs for the field solution, but suffers from the fact that the Green s function construction is generally difficult and costly to perform. In the following, we illustrate the use of the RGFM to efficiently solve (2) and subsequently find in (1). A. Unit Section Green s Function The first step in the construction of involves subdividing the region into smaller rectangular regions , , referred to as unit sections, as implied in Fig. 1. It is assumed that within each subdomain, the permittivity remains constant at its local average value (3) where is the area of . This piecewise constant approximation of the permittivity is reasonably accurate provided that the unit cell dimensions remain small and that abrupt changes in permittivity coincide with the cell interfaces. Using this permittivity description, the unit cell Green s function on can be computed analytically. In general, any boundary condition for on can be chosen, although use of homogeneous Neumann boundary conditions simplifies the RGFM development as well as the later formulation of the boundary integral equations. If we assume that the domain has local coordinates where and , then can be expressed as (4) (5) for . The coefficient , while for . If , then the positions of and must simply be reversed in (4). Additionally, an equivalent form for Fig. 2. Two adjacent unit sections to be combined into a single composite domain. can be obtained by interchanging and , and , and and in these expressions. The Green s function possesses an interesting property which will be used to advantage in the RGFM formulation. A simple development using (2) can be performed to show that [17] (6) where . The expression in (6) is also valid if we interchange the roles of and as well as and . This condition is particularly interesting if we apply it to a scenario where the source point is on a boundary such as at . In this case, we must take the limit as and and apply the homogeneous Neumann boundary condition on for the term where the observation point reaches the boundary before the source point does. This procedure results in the following: (7) For the case where , this condition becomes (8) B. Composite Green s Function Consider now the region in Fig. 2 where the Green s functions and are known on two adjacent domains and where is the interface between the two sections. To form for the composite domain , we first define such that (9) where . Clearly, satisfies the same differential equation as , but obeys different conditions on . Therefore, the dependence of must differ from that of only by an additive homogeneous solution of (2) which satisfies the Neumann boundary condition on all sides except the interface . Similarly, , must be a homogeneous solution of (2). Proper construction of these solutions is facilitated by making two key observations. 1) Equation (4) is a homogeneous solution to the wave equation and it must be combined with the form for in order to become a particular solution. 2) In light of (7) and (8), (4) does not satisfy the Neumann condition on if the source point is also on . Authorized licensed use limited to: Brigham Young University. Downloaded on February 3, 2009 at 11:19 from IEEE Xplore. Restrictions apply. 1812 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 46, NO. 12, DECEMBER 1998 These observations allow us to write (10) (11) where the functions are unknown and . To determine the unknown functions , we must apply continuity of and its normal derivative at . Mathematically, these continuity conditions assume the form (12) (13) where and . Now, substituting (10) and (11) into (13) and using the homogeneous Neumann condition for the term outside the integrals results in the relation (14) where approaches from the side of . Now using (7) for and (8) for in (14) reveals that . Finally, placing this result into (10) (12) gives (15) The important feature of (15) is that all of the terms are known except the function appearing in the integrand. To solve for this unknown, we project the and dependence of and and the dependence of onto basis sets complete on either or and the dependence of onto a basis complete on . As discussed in the Appendix, it is computationally beneficial to use basis functions that are of compact support. In our case, piecewise constant or pulse functions will be used. This choice provides the added advantage that the basis representation on is a subset of the representation on . Such an expansion can be written as (16) (17) where the series coefficients can be written in matrix form as (18) (19) (20) Explicit forms for are provided in the Appendix. Given this representation, it is clearly advantageous to use orthogonal basis sets (such as pulse functions) in order to simplify computation of . We now let represent the matrix of elements for which and represent the matrix of elements for which and . Then solution of (15) results in (21) Inserting into discrete forms of (10) and (11) results in the matrix expressions (22) (23) where denotes a transpose. The relations (22) and (23) allow construction of the composite Green s function for any source and observation point from knowledge of the Green s functions of the two contributing domains. C. Two-Dimensional Recursion The above development indicates how to form a composite Green s function by combining two unit section Green s functions horizontally, as implied in Fig. 2. Naturally, the procedure can be duplicated to combine Green s functions vertically, with results similar to those in (22) and (23). With this in mind, a recursive procedure can be used to systematically construct the Green s function for an arbitrary domain. This procedure consists of dividing the unit sections in the domain into adjacent pairs which are combined horizontally using (22) and (23) as depicted in Fig. 3(a). The resulting regions are then paired and combined with adjacent domains in the vertical direction, as implied in Fig. 3(b). This horizontal/vertical combining is then repeatedly performed with the newly constructed Green s functions until the composite Green s function for the entire domain has been computed. III. FIELD EVALUATION Following construction of the Green s function for the region, the fields scattered by the obstacle for a given source can be determined by proper implementation of surface integral equations [16]. For the domain in Fig. 1 we can write PV (24) for the fields in , where denotes the outward normal derivative, is the domain perimeter, and the integral should be given a principal value interpretation. Additionally (25) is the unbounded-space Green s function. For the fields in , we use the fact that satisfies homogeneous Neumann boundary conditions to write (26) Authorized licensed use limited to: Brigham Young University. Downloaded on February 3, 2009 at 11:19 from IEEE Xplore. Restrictions apply. JENSEN AND FREEZE: RECURSIVE GREEN S FUNCTION METHOD FOR ANALYSIS OF INHOMOGENEOUS DOMAINS 1813 (a) (b) Fig. 3. (a) Horizontal and (b) vertical combining of domains to make a single composite domain using the RGFM equations. The coupled solution of (24) and (26) can be accomplished using a simple MoM discretization. To perform this, we use the basis function decompositions (27) (28) where should be the same form as used in (16). These representations along with the matrix obtained from the RGFM can be used to construct the coupled matrix equation (29) (30) where is defined in (20) and (31) (32) (33) (34) The function represents a weighting function which in this work is the point-matching operator. Solution of (29) and (30) results in solution for the fields on the boundary. Subsequent integration of these fields using (24) (with the half factor and principal value interpretation removed) allows determination of the scattered fields at any desired location within . IV. COMPUTATIONAL COMPLEXITY The computational requirements of the RGFM can be calculated as follows. Consider a square domain discretized into points, such that the number of points in each dimension is . The cost of computing all of the unit cells is . Now, consider the cost of combining four cells, with points in each dimension, as implied in Fig. 3. This will involve two matrix inversions for the matrices one of size and one of size with a total cost of operations assuming an inversion routine. Additionally, the computations performed to evaluate (22) and (23) for horizontal and vertical combining are dominated by the matrix products, at an estimated cost of operations. Since there are cells in the domain, a total of inversions and combinations are required. Also, the combining process must be repeated times. The cost of inverting the final matrix is . Summing over each recursive step in the algorithm as well as the final matrix solution results in the cost given by (35) (36) where we have used that at each step. The algorithm storage requirements are dominated by the matrix required in the surface MoM implementation, resulting in a storage cost of . It is noteworthy that the RGFM can also be implemented to allow computation of field values internal to the inhomogeneous domain. While the details of this procedure have not been demonstrated here, a similar analysis shows that the asymptotic computational cost involved with recursively updating Green s functions with observation points internal to the domain increases to . Additionally, the storage requirements for this procedure become . V. RESULTS In this section, we illustrate the performance of the RGFM in computing the fields scattered from various homogeneous and inhomogeneous dielectric cylinders. Comparisons are made with results from the surface MoM, volume MoM, and closedform expressions. Numerical investigations have shown that use of a five-point Gaussian quadrature integration for all Authorized licensed use limited to: Brigham Young University. Downloaded on February 3, 2009 at 11:19 from IEEE Xplore. Restrictions apply. 1814 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 46, NO. 12, DECEMBER 1998 Fig. 4. Bistatic scattering width for the square homogeneous cylinder with and computed using the surface MoM (SMOM), volume MoM (VMOM), and RGFM. Fig. 5. Bistatic scattering width for the square homogeneous cylinder with and computed using the surface MoM, volume MoM, and RGFM. required integrals provides relatively good accuracy. Pulsebasis functions with point matching are used in both the surface and volume MoM solutions. A. Homogeneous Square Cylinder To test the accuracy of the RGFM, we first apply it to compute the scattering from a homogeneous dielectric cylinder with a square cross section. The geometry is illustrated in the inset of Fig. 4. The cylinder has sides of length , where is the free-space wavelength, and relative permittivity . The spatial discretization uses an 8 8 grid of cells in order to provide ten cells per wavelength resolution. The curves in Fig. 4 compare the far-zone bistatic scattering width computed using the RGFM, the surface MoM, and the volume MoM for an incident plane wave as shown. Excellent agreement between the three results can be observed, with the best correlation existing between the RGFM and surface MoM results. This simple test shows that the RGFM faithfully constructs the Green s function for this domain. Fig. 5 illustrates a similar result for the square cylinder with side length and permittivity . In this case, a 16 16 grid is used to maintain the ten cells per wavelength resolution. Here, the surface MoM and RGFM results agree very well, but they differ from those obtained using the volume MoM. Such differences between the results of surface and Fig. 6. Bistatic scattering width for the two-layer circular cylinder shown computed using the eigenfunction expansion, volume MoM, and RGFM. The cylinder permittivity is for and for . volume integral equation solutions have been noted in the literature [18]. B. Two-Layer Circular Cylinder The discrepancies observed in Fig. 5 raise questions concerning the accuracy of the surface and volume formulations. To further investigate this issue, we examine scattering from a circular cylinder with two layers in the radial direction, as indicated in the inset of Fig. 6. The cylinder permittivity is for and for , where is the cylinder radial coordinate. The curves in Fig. 6 compare the bistatic scattering width from the eigenfunction series solution with that obtained using the RGFM and volume MoM. In the latter two cases, the cylinder is modeled using 16 16 cells with a stair-stepped approximation of the cylindrical interfaces. For this particular geometry, the RGFM results agree favorably with the exact solution, while the volume MoM solution shows some error. This improvement in accuracy offered by the RGFM is an important feature. The error in the RGFM in the forward and backscattered directions is most likely due to the block discretization of the cylindrical surface. C. Arm Model One currently important application of electromagnetic scattering from inhomogeneous objects involves microwave imaging of biological tissue. In light of this, consider applying the RGFM to the simple model of the human arm shown in Fig. 7. In this diagram, all dimensions are given in wavelengths at a frequency of 1.2 GHz. The model consists of a circular cylinder with representing the muscle and two smaller cylinders with representing the bone. The RGFM grid uses 32 32 cells in order to accurately model the fields in the high permittivity muscle. Fig. 8 illustrates the bistatic scattering width obtained for this computation. To test the accuracy of the RGFM in simulating this structure, the algorithm is also applied to the simple concentric cylinder shown in the inset of Fig. 6 with for and for . The comparison of this result with that obtained from the exact solution in Fig. 8 shows Authorized licensed use limited to: Brigham Young University. Downloaded on February 3, 2009 at 11:19 from IEEE Xplore. Restrictions apply. JENSEN AND FREEZE: RECURSIVE GREEN S FUNCTION METHOD FOR ANALYSIS OF INHOMOGENEOUS DOMAINS 1815 Fig. 7. Simple model of the arm used in the RGFM computation. All dimensions are in wavelengths at a frequency of 1.2 GHz. Fig. 8. Bistatic scattering width for a two-layer cylinder ( for and for ) computed using the exact and RGFM solutions. Also shown is the RGFM solution for the arm model of Fig. 7. Fig. 9. CPU time versus number of unit cells for the RGFM and volume MoM. that the RGFM provides a reasonably accurate solution. This fact provides some confidence concerning the accuracy of the RGFM result for the simple arm model. D. Computational Complexity Fig. 9 illustrates the computational time required versus the number of cells for the RGFM and volume MoM. It should be emphasized that represents the number of unit cells in either method, so that in the RGFM, the actual number of grid points is and the number of unknowns is (since there are two unknowns for each surface grid point). The computer platform used was an HP735 workstation with 128 MB of memory. In this comparison, a simple LU Decomposition with back substitution is used to perform the matrix solutions. As a result, the timing numbers could be improved using more elaborate linear system solution techniques, particularly for the volume MoM algorithm. However, the plot does illustrate that: 1) the RGFM requires significantly less computation time as compared to the volume MoM and 2) the RGFM computational complexity scales as . Also noteworthy is the fact that the volume MoM has a lower computational complexity for small numbers of unit cells since, in this case, it has fewer unknowns than the RGFM. Finally, the storage requirement of the RGFM allows solution of larger domains than the volume MoM which requires storage. As a result, the curve for the volume MoM in Fig. 9 stops at due to the memory constraints of the computational platform used. VI. CONCLUSIONS We have presented a novel numerical technique for solving the scalar wave equation in 2-D inhomogeneous dielectric domains. The methodology constructs the Green s function of the region by recursively combining known Green s functions from smaller subdomains. Boundary integral techniques are then used to determine the fields at the domain boundary for a given excitation. Numerical results show that the RGFM faithfully constructs the Green s function for inhomogeneous domains and provides highly accurate results for scattering from various cylindrical structures. Comparison with closedform solutions for canonical dielectric cylinders has revealed that the RGFM provides improved accuracy when compared with the volume MoM. Additionally, it allows simulation of domains for multiple source configurations with an asymptotic computational complexity of and a storage requirement of . This allows solution of larger problems with less computational time as compared to traditional MoM solutions.Work is currently under way to allow computation of fields internal to the dielectric domain and to simulate threedimensional geometries using the RGFM. Reports of these activities will appear in future communications. APPENDIX INTEGRATED GREEN S FUNCTIONS The choice of basis functions to use in (16) is important since it can significantly impact the convergence properties of (4) when used in (19). For example, if sinusoidal basis functions with global support are used, the resulting series obtained using term-by-term integration in (19) is divergent for points where and slowly convergent for other points. However, if simple piecewise constant or pulse functions are used having compact support, then the series resulting from (19) converges rapidly for all points except when and the boundary integrations occur in and (i.e., , or ). In this case, it is most convenient to use the form discussed immediately following (4) where the roles of and are interchanged. With this in mind, consider that the pulse functions have widths or depending on whether they span the coordinate or . In this case, if we let and Authorized licensed use limited to: Brigham Young University. Downloaded on February 3, 2009 at 11:19 from IEEE Xplore. Restrictions apply. 1816 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 46, NO. 12, DECEMBER 1998 , the integrated series assumes the form (37) where if or if or (38) REFERENCES [1] W. C. Chew and Y. M. Wang, Reconstruction of two-dimensional permittivity distribution using the distorted born iterative method, IEEE Trans. Med. Imaging, vol. 9, pp. 218 225, June 1990. [2] N. Joachimowicz, C. Pichot, and J. P. Hugonin, Inverse scattering: An iterative numerical method for electromagnetic imaging, IEEE Trans. Antennas Propagat., vol. 39, pp. 1742 1752, Dec. 1991. [3] D. M. Sullivan, O. P. Gandhi, and A. Taflove, Use of the finitedifference time-domain method for calculating EM absorption in man models, IEEE Trans. Biomed. Eng., vol. 35, pp. 179 186, Mar. 1988. [4] P. J. Dimbylow, FDTD calculations of the SAR for a dipole closely coupled to the head at 900 MHz and 1.9 GHz, Phys. Med. Biol., vol. 38, pp. 361 368, Feb. 1993. [5] D. T. Borup, D. M. Sullivan, and O. P. Gandhi, Comparison of the FFT conjugate gradient method and the finite-difference time-domain method for the 2-D absorption problem, IEEE Trans. Microwave Theory Tech., vol. MTT-35, pp. 383 395, Apr. 1987. [6] R. Kastner and R. Mittra, A new stacked two-dimensional spectral iterative technique (SIT) for analyzing microwave power deposition in biological media, IEEE Trans. Microwave Theory Tech., vol. MTT-31, pp. 898 904, Nov. 1983. [7] M. A. Jensen and Y. Rahmat-Samii, EM interaction of handset antennas and a human in personal communications, Proc. IEEE, vol. 83, pp. 7 17, Jan. 1995. [8] H. R. Chuang, Human operator coupling effects on radiation characteristics of a portable communication dipole antenna, IEEE Trans. Antennas Propagat., vol. 42, pp. 556 560, Apr. 1994. [9] W. C. Chew and C. C. Lu, The use of Huygens equivalence principle for solving the volume integral equation of scattering, IEEE Trans. Antennas Propagat., vol. 41, pp. 897 904, July 1993. [10] C. C. Lu and W. C. Chew, The use of Huygens equivalence principle for solving 3-D volume integral equation of scattering, IEEE Trans. Antennas Propagat., vol. 43, pp. 500 507, May 1995. [11] K. B. Kahen, Analysis of distributed-feedback lasers: A recursive Green s function approach, Appl. Phys. Lett., vol. 61, no. 17, pp. 2012 2014, Oct. 1992. [12] K. B. Kahen, Analysis of distributed-feedback lasers using a recursive Green s functional approach, IEEE J. Quantum Electron., vol. 29, pp. 368 373, Feb. 1993. [13] J. D. Freeze, M. A. Jensen, and R. H. Selfridge, A unified Green s function analysis of complicated DFB lasers, IEEE J. Quantum Electron., vol. 33, pp. 1253 1259, Aug. 1997. [14] A. George, Nested dissection of a regular finite element mesh, SIAM J. Numer. Anal., vol. 10, pp. 345 363, Apr. 1973. [15] K. B. Kahen, Recursive-Green s-function analysis of wave propagation in two-dimensional nonhomogeneous media, Physical Rev. E, vol. 47, pp. 2927 2933, Apr. 1993. [16] W. C. Chew, Waves and Fields in Inhomogeneous Media. New York: Van Nostrand Reinhold, 1990. [17] J. Mathews and R. L. Walker, Mathematical Methods of Physics. New York: Benjamin, 1970, pp. 272 273. [18] T. K. Sarkar, E. Arvas, and S. Ponnapalli, Electromagnetic scattering from dielectric bodies, IEEE Trans. Antennas Propagat., vol. 37, pp. 673 676, May 1989. Michael A. Jensen (S 93 M 95) received the B.S. (summa cum laude) and M.S. degrees in electrical engineering from Brigham Young University (BYU), Provo, UT, in 1990 and 1991, respectively, and the Ph.D. degree in electrical engineering at the University of California, Los Angeles (UCLA), in 1994. From 1989 to 1991, he was a Graduate Research Assistant in the Lasers and Optics Laboratory at BYU. From 1991 to 1994 he was a Graduate Student Researcher in the Antenna Laboratory at UCLA. Since 1994 he has been an Assistant Professor in the Electrical and Computer Engineering Department at Brigham Young University. His main research interests include radiation and propagation for personal communications, numerical electromagnetics, optical fiber communication, and implementation of finite-difference schemes on massively parallel computer architectures. Dr. Jensen received a National Science Foundation Graduate Fellowship in 1990. He is a member of Eta Kappa Nu and Tau Beta Pi. Jim D. Freeze (S 95) received the B.S. and M.S. degrees (electrical and computer engineering) from Brigham Young University, Provo, UT, in 1991. He is currently working toward the Ph.D. degree at the same university. He is currently with Lexmark International as the Director of Systems Research in the Consumer Products Division. His research interests include mathematical and numerical modeling of electromagnetic waves and development of efficient data compression algorithms. Authorized licensed use limited to: Brigham Young University. Downloaded on February 3, 2009 at 11:19 from IEEE Xplore. Restrictions apply.