Miller, E.K. “Computational Electromagnetics” The Electrical Engineering Handbook Ed. Richard C. Dorf Boca Raton: CRC Press LLC, 2000 45 Computational Electromagnetics 45.1 Introduction 45.2 Background Discussion Modeling as a Transfer Function ? Some Issues Involved in Developing a Computer Model 45.3 Analytical Issues in Developing a Computer Model Selection of Solution Domain ? Selection of Field Propagator 45.4 Numerical Issues in Developing a Computer Model Sampling Functions ? The Method of Moments 45.5 Some Practical Considerations Integral Equation Modeling ? Differential Equation Modeling ? Discussion ? Sampling Requirements 45.6 Ways of Decreasing Computer Time 45.7 Validation, Error Checking, and Error Analysis Modeling Uncertainties ? Validation and Error Checking 45.8 Concluding Remarks 45.1 Introduction The continuing growth of computing resources is changing how we think about, formulate, solve, and interpret problems. In electromagnetics as elsewhere, computational techniques are complementing the more traditional approaches of measurement and analysis to vastly broaden the breadth and depth of problems that are now quantifiable. Computational electromagnetics (CEM) may be broadly defined as that branch of electromagnetics that intrinsically and routinely involves using a digital computer to obtain numerical results. With the evolu- tionary development of CEM during the past 20-plus years, the third tool of computational methods has been added to the two classical tools of experimental observation and mathematical analysis. This discussion reviews some of the basic issues involved in CEM and includes only the detail needed to illustrate the central ideas involved. The underlying principles that unify the various modeling approaches used in electromagnetics are emphasized while avoiding most of the specifics that make them different. Listed throughout are representative, but not exhaustive, numbers of references that deal with various specialty aspects of CEM. For readers interested in broader, more general expositions, the well-known book on the moment method by Harrington [1968]; the books edited by Mittra [1973, 1975], Uslenghi [1978], and Strait [1980]; the monographs by Stutzman and Thiele [1981], Popovic, et al. [1982], Moore and Pizer [1984], and Wang [1991]; and an IEEE Press reprint volume on the topic edited by Miller et al. [1991] are recommended, as is the article by Miller [1988] from which this material is excerpted. This chapter is excerpted from E. K. Miller, “A selective survey of computational electromagnetics,” IEEE Trans. Antennas Propagat., vol. AP-36, pp. 1281–1305, ?1988 IEEE. E.K. Miller Los Alamos National Laboratory ? 2000 by CRC Press LLC 45.2 Background Discussion Electromagnetics is the scientific discipline that deals with electric and magnetic sources and the fields these sources produce in specified environments. Maxwell’s equations provide the starting point for the study of electromagnetic problems, together with certain principles and theorems such as superposition, reciprocity, equivalence, induction, duality, linearity, and uniqueness, derived therefrom [Stratton, 1941; Harrington, 1961]. While a variety of specialized problems can be identified, a common ingredient of essentially all of them is that of establishing a quantitative relationship between a cause (forcing function or input) and its effect (the response or output), a relationship which we refer to as a field propagator, the computational characteristics of which are determined by the mathematical form used to describe it. Modeling as a Transfer Function The foregoing relationship may be viewed as a gener- alized transfer function (see Fig. 45.1) in which two basic problem types become apparent. For the anal- ysis or the direct problem, the input is known and the transfer function is derivable from the problem spec- ification, with the output or response to be deter- mined. For the case of the synthesis or inverse problem, two problem classes may be identified. The easier synthesis problem involves finding the input, given the output and transfer function, an example of which is that of determining the source voltages that produce an observed pattern for a known antenna array. The more difficult synthesis problem itself separates into two problems. One is that of finding the transfer function, given the input and output, an example of which is that of finding a source distribution that produces a given far field. The other and still more difficult is that of finding the object geometry that produces an observed scattered field from a known exciting field. The latter problem is the most difficult of the three synthesis problems to solve because it is intrinsically transcendental and nonlinear. Electromagnetic propagators are derived from a particular solution of Maxwell’s equations, as the cause mentioned above normally involves some specified or known excitation whose effect is to induce some to-be- determined response (e.g., a radar cross section, antenna radiation pattern). It therefore follows that the essence of electromagnetics is the study and determination of field propagators to thereby obtain an input–output transfer function for the problem of interest, and it follows that this is also the goal of CEM. Some Issues Involved in Developing a Computer Model We briefly consider here a classification of model types, the steps involved in developing a computer model, the desirable attributes of a computer model, and finally the role of approximation throughout the modeling process. Classification of Model Types It is convenient to classify solution techniques for electromagnetic modeling in terms of the field propagator that might be used, the anticipated application, and the problem type for which the model is intended to be used, as is outlined in Table 45.1. Selection of a field propagator in the form, for example, of the Maxwell curl equations, a Green’s function, modal or spectral expansions, or an optical description is a necessary first step in developing a solution to any electromagnetic problem. Development of a Computer Model Development of a computer model in electromagnetics or literally any other disciplinary activity can be decomposed into a small number of basic, generic steps. These steps might be described by different names but FIGURE 45.1The electromagnetic transfer function relates the input, output, and problem. ? 2000 by CRC Press LLC would include at a minimum those outlined in Table 45.2. Note that by its nature, validation is an open-ended process that cumulatively can absorb more effort than all the other steps together. The primary focus of the following discussion is on the issue of numerical implementation. Desirable Attributes of a Computer Model A computer model must have some minimum set of basic properties to be useful. From the long list of attributes that might be desired, we consider: (1) accuracy, (2) efficiency, and (3) utility the three most important as summarized in Table 45.3. Accuracy is put foremost because results of insufficient or unknown accuracy have uncertain value and may even be harmful. On the other hand, a code that produces accurate results but at unacceptable cost will have hardly any more value. Finally, a code’s applicability in terms of the depth and breadth of the problems for which it can be used determines its utility. The Role of Approximation As approximation is an intrinsic part of each step involved in developing a computer model, we summarize some of the more commonly used approximations in Table 45.4. We note that the distinction between an approximation at the conceptualization step and during the formulation is somewhat arbitrary, but choose to use the former category for those approximations that occur before the formulation itself. TABLE 45.1 Classification of Model Types in CEM Field Propagator Description Based on Integral operator Green’s function for infinite medium or special boundaries Differential operator Maxwell curl equations or their integral counterparts Modal expansions Solutions of Maxwell’s equations in a particular coordinate system and expansion Optical description Rays and diffraction coefficients Application Requires Radiation Determining the originating sources of a field and patterns they produce Propagation Obtaining the fields distant from a known source Scattering Determining the perturbing effects of medium inhomogeneities Problem type Characterized by Solution domain Time or frequency Solution space Configuration or wave number Dimensionality 1D, 2D, 3D Electrical properties of medium and/or boundary Dielectric, lossy, perfectly conducting, anisotropic, inhomogeneous, nonlinear, bianisotropic Boundary geometry Linear, curved, segmented, compound, arbitrary TABLE 45.2 Steps in Developing a Computer Model Step Activity Conceptualization Encapsulating observation and analysis in terms of elementary physical principles and their mathematical descriptions Formulation Fleshing out of the elementary description into a more complete, formally solved, mathematical representation Numerical implementation Transforming into a computer algorithm using various numerical techniques Computation Obtaining quantitative results Validation Determining the numerical and physical credibility of the computed results ? 2000 by CRC Press LLC 45.3 Analytical Issues in Developing a Computer Model Attention here is limited primarily to propagators that use either the Maxwell curl equations or source integrals which employ a Green’s function, although for completeness we briefly discuss modal and optical techniques as well. Selection of Solution Domain Either the integral equation (IE) or differential equation (DE) propagator can be formulated in the time domain, where time is treated as an independent variable, or in the frequency domain, where the harmonic TABLE 45.3 Desirable Attributes in a Computer Model Attribute Description Accuracy The quantitative degree to which the computed results conform to the mathematical and physical reality being modeled. Accuracy, preferably of known and, better yet, selectable value, is the single most important model attribute. It is determined by the physical modeling error (e P ) and numerical modeling error (e N ). Efficiency The relative cost of obtaining the needed results. It is determined by the human effort required to develop the computer input and interpret the output and by the associated computer cost of running the model. Utility The applicability of the computer model in terms of problem size and complexity. Utility also relates to ease of use, reliability of results obtained, etc. TABLE 45.4 Representative Approximations that Arise in Model Development Approximation Implementation/Implications Conceptualization Physical optics Surface sources given by tangential components of incident field, with fields subsequently propagated via a Green’s function. Best for backscatter and main-lobe region of reflector antennas, from resonance region (ka > 1) and up in frequency. Physical theory of diffraction Combines aspects of physical optics and geometrical theory of diffraction, primarily via use of edge-current corrections to utilize best features of each. Geometrical theory diffraction Fields propagated via a divergence factor with amplitude obtained from diffraction coefficient. Generally applicable for ka > 2–5. Can involve complicated ray tracing. Geometrical optics Ray tracing without diffraction. Improves with increasing frequency. Compensation theorem Solution obtained in terms of perturbation from a reference, known solution. Born–Rytov Approach used for low-contrast, penetrable objects where sources are estimated from incident field. Rayleigh Fields at surface of object represented in terms of only outward propagating components in a modal expansion. Formulation Surface impedance Reduces number of field quantities by assuming an impedance relation between tangential E and H at surface of penetrable object. May be used in connection with physical optics. Thin-wire Reduces surface integral on thin, wirelike object to a line integral by ignoring circumferential current and circumferential variation of longitudinal current, which is represented as a filament. Generally limited to ka < 1 where a is the wire radius. Numerical Implementation ?f /?x ? (f + – f – )/(x + – x – ) òf (x)dx ? ?f (x i )Dx i Differentiation and integration of continuous functions represented in terms of analytic operations on sampled approximations, for which polynomial or trigonometric functions are often used. Inherently a discretizing operation, for which typically Dx < l/2p for acceptable accuracy. Computation Deviation of numerical model from physical reality Affects solution accuracy and relatability to physical problem in ways that are difficult to predict and quantify. Nonconverged solution Discretized solutions usually converge globally in proportion to exp(–AN x ) with A determined by the problem. At least two solutions using different numbers of unknowns N x are needed to estimate A. ? 2000 by CRC Press LLC time variation exp( j w t ) is assumed. Whatever propagator and domain are chosen, the analytically formal solution can be numerically quantified via the method of moments (MoM) [Harrington, 1968], leading ultimately to a linear system of equations as a result of developing a discretized and sampled approximation to the continuous (generally) physical reality being modeled. Developing the approach that may be best suited to a particular problem involves making trade-offs among a variety of choices throughout the analytical formulation and numerical implementation, some aspects of which are now considered. Selection of Field Propagator We briefly discuss and compare the characteristics of the various propagator-based models in terms of their development and applicability. Integral Equation Model The basic starting point for developing an IE model in electromagnetics is selection of a Green’s function appropriate for the problem class of interest. While there are a variety of Green’s functions from which to choose, a typical starting point for most IE MoM models is that for an infinite medium. One of the more straightforward is based on the scalar Green’s function and Green’s theorem. This leads to the Kirchhoff integrals [Stratton, 1941, p. 464 et seq.], from which the fields in a given contiguous volume of space can be written in terms of integrals over the surfaces that bound it and volume integrals over those sources located within it. Analytical manipulation of a source integral that incorporates the selected Green’s function as part of its kernel function then follows, with the specific details depending on the particular formulation being used. Perhaps the simplest is that of boundary-condition matching wherein the behavior required of the electric and/or magnetic fields at specified surfaces that define the problem geometry is explicitly imposed. Alternative formulations, for example, the Rayleigh–Ritz variational method and Rumsey’s reaction concept, might be used instead, but as pointed out by Harrington [in Miller et al., 1991], from the viewpoint of a numerical imple- mentation any of these approaches lead to formally equivalent models. This analytical formulation leads to an integral operator, whose kernel can include differential operators as well, which acts on the unknown source or field. Although it would be more accurate to refer to this as an integrodifferential equation, it is usually called simply an integral equation. Two general kinds of integral equations are obtained. In the frequency domain, representative forms for a perfect electric conductor are (45.1a) (45.1b) where E and H are the electric and magnetic fields, respectively, r, r￠ are the spatial coordinate of the observation and source points, the superscript inc denotes incident-field quantities, and j(r,r￠) = exp[–jk*r – r￠*]/*r – r￠* is the free-space Green’s function. These equations are known respectively as Fredholm integral equations of the first and second kinds, differing by whether the unknown appears only under the integral or outside it as well [Poggio and Miller in Mittra, 1973]. Differential-Equation Model A DE MoM model, being based on the defining Maxwell’s equations, requires intrinsically less analytical manipulation than does derivation of an IE model. Numerical implementation of a DE model, however, can differ significantly from that used for an IE formulation in a number of ways for several reasons: nEr n nHrr,r n Er,r r,r r ′=′￠′ ￠￠ -￠× ￠￠? ￠￠? ò inc () { [ ()]( ) [ ( ) ( )} ; 1 4p wm j j j ds S S nHrnHr nnHr r,r r′=′+′￠′ ￠′ ￠? ￠￠? ò () () [ ()] ( )};2 1 2 inc ds p j S S ? 2000 by CRC Press LLC TABLE 45.5 Comparison of IE- and DE-Field Propagators and Their Numerical Treatment 1.The differential operator is a local rather than global one in contrast to the Green’s function upon which the integral operator is based. This means that the spatial variation of the fields must be developed from sampling in as many dimensions as possessed by the problem, rather than one less as the IE model permits if an appropriate Green’s function is available. 2.The integral operator includes an explicit radiation condition, whereas the DE does not. 3.The differential operator includes a capability to treat medium inhomogeneities, non-linearities, and time variations in a more straightforward manner than does the integral operator, for which an appro- priate Green’s function may not be available. These and other differences between development of IE and DE models are summarized in Table 45.5, with their modeling applicability compared in Table 45.6. Modal-Expansion Model Modal expansions are useful for propagating electromagnetic fields because the source-field relationship can be expressed in terms of well-known analytical functions as an alternate way of writing a Green’s function for special distributions of point sources. In two dimensions, for example, the propagator can be written in terms of circular harmonics and cylindrical Hankel functions. Corresponding expressions in three dimensions might involve spherical harmonics, spherical Hankel functions, and Legendre polynomials. Expansion in terms of analytical solutions to the wave equation in other coordinate systems can also be used but requires computation Differential Form Integral Form Field propagator Maxwell curl equations Green’s function Boundary treatment At infinity (radiation condition) Local or global “lookback” to approximate outward propagating wave Green’s function On object Appropriate field values specified on mesh boundaries to obtain stairstep, piecewise linear, or other approximation to the boundary Appropriate field values specified on object contour which can in principle be a general, curvilinear surface, although this possibility seems to be seldom used Sampling requirements No. of space samples N x μ (L/DL) D N x μ (L/DL) D–1 No. of time steps N t μ (L/DL) ? cT/dtN t μ (L/DL) ? cT/dt No. of excitations N rhs μ (L/DL) N rhs μ (L/DL) (right-hand sides) Linear system L is problem size D is no. of problem dimensions (1, 2, 3) T is observation time DL is spatial resolution dt is time resolution Sparse, but larger Dense, but smaller. In this comparison, note that we assume the IE permits a sampling of order one less than the problem dimension, i.e., inhomogeneous problems are excluded. Dependence of solution time on highest-order term in (L/DL) Frequency domain T w μ N x 2(D–1)/D+1 = (L/DL) 3D–2 T w μ N x 3 = (L/DL) 3(D–1) Time domain Explicit T t μ N x N t N rhs = (L/DL) D+1+r T t μ N x 2 N t N rhs = (L/DL) 2D–1+r ; 0 ￡ r ￡ 1 Implicit T t μ N x 2(D–1)/D+1 = (L/DL) 3D–2 , D = 2, 3; T t μ N x 3 = (L/DL) 3(D–1) μ N x N t N rhs = (L/DL) 2+r , D = 1; 0 ￡ r ￡ 1 Note that D is the number of spatial dimensions in the problem and is not necessarily the sampling dimensionality d. The distinction is important because when an appropriate Green’s function is available, the source integrals are usually one dimension less than the problem dimension, i.e., d = D – 1. An exception is an inhomogeneous, penetrable body where d = D when using an IE. We also assume for simplicity that matrix solution is achieved via factorization rather than iteration but that banded matrices are exploited for the DE approach where feasible. The solution-time dependencies given can thus be regarded as upper-bound estimates. See Table 45.10 for further discussion of linear-system solutions. ? 2000 by CRC Press LLC of special functions that are generally less easily evaluated, such as Mathieu functions for the two-dimensional solution in elliptical coordinates and spheroidal functions for the three-dimensional solution in oblate or prolate spheroidal coordinates. One implementation of modal propagators for numerical modeling is that due to Waterman [in Mittra, 1973], whose approach uses the extended boundary condition (EBC) whereby the required field behavior is satisfied away from the boundary surface on which the sources are located. This procedure, widely known as the T-matrix approach, has evidently been more widely used in optics and acoustics than in electromagnetics. In what amounts to a reciprocal application of EBC, the sources can be removed from the boundary surface on which the field-boundary conditions are applied. These modal techniques seem to offer some computational advantages for certain kinds of problems and might be regarded as using entire-domain basis and testing functions but nevertheless lead to linear systems of equations whose numerical solution is required. Fourier transform solution techniques might also be included in this category since they do involve modal expansions, but that is a specialized area that we do not pursue further here. Modal expansions are receiving increasing attention under the general name “fast multipole method,” which is motivated by the goal of systematically exploiting the reduced complexity of the source-field interactions as their separation increases. The objective is to reduce the significant interactions of a Green’s-function based matrix from being proportional to (N x ) 2 to of order N x log (N x ), thus offering the possibility of decreasing the operation count of iterative solutions. TABLE 45.6 Relative Applicability of IE- and DE-Based Computer Models Time Domain Frequency Domain DE IE Issue DE IE Medium ??Linear ?? ~ x Dispersive ? x Lossy ?? ? ~ Anisotropic ? x Inhomogeneous ? x ? x Nonlinear x x ? x Time-varying x x Object ~ ? Wire ~ ? ??Closed surface ?? ??Penetrable volume ~ ? Open surface ~ ? Boundary Conditions ??Interior problem ?? ~ ? Exterior problem ~ ? ??Linear ?? ??Nonlinear x x ??Time-varying x x ~ x Halfspace ~ ? Other Aspects ~ ~ Symmetry exploitation ?? ~ ? Far-field evaluation ~ ? x ~ Number of unknowns ~ ? ? ~ Length of code ~ x Suitability for Hybridizing with Other: ~ ? Numerical procedures ?? x ~ Analytical procedures ~ ? x ~ GTD x ? ? signifies highly suited or most advantageous. ~ signifies moderately suited or neutral. x signifies unsuited or least advantageous. ? 2000 by CRC Press LLC Geometrical-Optics Model Geometrical optics and the geometrical theory of diffraction (GTD) are high-frequency asymptotic techniques wherein the fields are propagated using such optical concepts as shadowing, ray tubes, and refraction and diffraction. Although conceptually straightforward, optical techniques are limited analytically by the unavail- ability of diffraction coefficients for various geometries and material bodies and numerically by the need to trace rays over complex surfaces. There is a vast literature on geometrical optics and GTD, as may be ascertained by examining the yearly and cumulative indexes of such publications as the Transactions of the IEEE Antennas and Propagation Society. 45.4 Numerical Issues in Developing a Computer Model Sampling Functions At the core of numerical analysis is the idea of polynomial approximation, an observation made by Arden and Astill [1970] in facetiously using the subtitle “Numerical Analysis or 1001 Applications of Taylor’s Series.” The basic idea is to approximate quantities of interest in terms of sampling functions, often polynomials, that are then substituted for these quantities in various analytical operations. Thus, integral operators are replaced by finite sums, and differential operators are similarly replaced by generalized finite differences. For example, use of a first-order difference to approximate a derivative of the function F(x) in terms of samples F(x + ) and F(x – ) leads to (45.2a) and implies a linear variation for F(x) between x + and x – as does use of the trapezoidal rule (45.2b) to approximate the integral of F(x), where h = x + – x – . The central-difference approximation for the second derivative, (45.2c) similarly implies a quadratic variation for F(x) around x 0 = x + – h/2 = x – + h/2, as does use of Simpson’s rule (45.2d) to approximate the integral. Other kinds of polynomials and function sampling can be used, as discussed in a large volume of literature, some examples of which are Abramowitz and Stegun [1964], Acton [1970], and Press et al. [1986]. It is interesting to see that numerical differentiation and integration can be accomplished using the same set of function samples and spacings, differing only in the signs and values of some of the associated weights. Note also that the added degrees of freedom that arise when the function samples can be unevenly spaced, as in Gaussian quadrature, produce a generally more accurate result (for well-behaved functions) for dF x dx Fx Fx h xxx () ( ) ( ) ? - ￡￡ +- -+ ; Fxdx h Fx Fx x x () [( ) ( )]?+ - + ò +- 2 dFx dx Fx Fx Fx h 2 2 0 2 2() [( ) ( ) ( ) ? -+ +- Fxdx h Fx Fx Fx x x () [( ) ( ) ( )]?++ - + ò +- 6 4 0 ? 2000 by CRC Press LLC a given number of samples. This suggests the benefits that might be derived from using unequal sample sizes in MoM modeling should a systematic way of determining the best nonuniform sampling scheme be developed. The Method of Moments Numerical implementation of the moment method is a relatively straightforward, and an intuitively logical, extension of these basic elements of numerical analysis, as described in the book by Harrington [1968] and discussed and used extensively in CEM [see, for example, Mittra, 1973, 1975; Strait, 1980; Poggio and Miller, 1988]. Whether it is an integral equation, a differential equation, or another approach that is being used for the numerical model, three essential sampling operations are involved in reducing the analytical formulation via the moment method to a computer algorithm as outlined in Table 45.7. We note that operator sampling can ultimately determine the sampling density needed to achieve a desired accuracy in the source–field rela- tionships involving integral operators, especially at and near the “self term,” where the observation and source points become coincident or nearly so and the integral becomes nearly singular. Whatever the method used for these sampling operations, they lead to a linear system of equations or matrix approximation of the original integral or differential operators. Because the operations and choices involved in developing this matrix descrip- tion are common to all moment-method models, we shall discuss them in somewhat more detail. When using IE techniques, the coefficient matrix in the linear system of equations that results is most often referred to as an impedance matrix because in the case of the E-field form, its multiplication of the vector of unknown currents equals a vector of electric fields or voltages. The inverse matrix similarly is often called an admittance matrix because its multiplication of the electric-field or voltage vector yields the unknown-current vector. In this discussion we instead use the terms direct matrix and solution matrix because they are more generic descriptions whatever the forms of the originating integral or differential equations. As illustrated in the following, development of the direct matrix and solution matrix dominates both the computer time and storage requirements of numerical modeling. In the particular case of an IE model, the coefficients of the direct or original matrix are the mutual impedances of the multiport representation which approximates the problem being modeled, and the coeffi- cients of its solution matrix (or equivalent thereof) are the mutual admittances. Depending on whether a subdomain or entire-domain basis has been used (see Basic Function Selection), these impedances and admittances represent either spatial or modal interactions among the N ports of the numerical model. In either case, these TABLE 45.7Sampling Operations Involved in MoM Modeling DE Model IE Model Equation L(s￠)f(s￠) = g(s￠) L(s,s￠)f(s￠) = g(s) Sampling of: Unknown via basis- functions b j (s￠) using f(s￠) ? ?a j b j (s￠) Subdomain bases usually of low order are used. Known as FD procedure when pulse basis is used, and as FE approach when bases are linear. Can use either subdomain or entire-domain bases. Use of latter is generally confined to bodies of rotation. Former is usually of low order, with piecewise linear or sinusoidal being the maximum variation used. Equation via weight functions w i (s) <w i (s),L(s,s￠)?a j b j (s￠)> = <w i (s),g(s) > to get Z ij a j = g i Pointwise matching is commonly employed, using a delta function. Pulse and linear matching are also used. Pointwise matching is commonly employed, using a delta function. For wires, pulse, linear, and sinusoidal testing is also used. Linear and sinusoidal testing is also used for surfaces. Operator Operator sampling for DE models is entwined with sampling the unknown in terms of the difference operators used. Sampling needed depends on the nature of the integral operator L(s,s￠). An important consideration whenever the field integrals cannot be evaluated in closed form. Solution of: Z ij a j = g i for the a j Interaction matrix is sparse. Time-domain approach may be explicit or implicit. In frequency domain, banded-matrix technique usually used. Interaction matrix is full. Solution via factorization or iteration. ? 2000 by CRC Press LLC coefficients possess a physical relatability to the problem being modeled and ultimately provide all the infor- mation available concerning any electromagnetic observables that are subsequently obtained. Similar observations might also be made regarding the coefficients of the DE models but whose multiport representations describe local rather than global interactions. Because the DE model almost always leads to a larger, albeit less dense, direct matrix, its inverse (or equivalent) is rarely computed. It is worth noting that there are two widely used approaches for DE modeling, finite-difference (FD) and finite-element (FE) methods. They differ primarily in how the differential operators are approximated and the differential equations are satisfied, i.e., in the order of the basis and weight functions, although the FE method commonly starts from a variational viewpoint, while the FD approach begins from the defining differential equations. The FE method is generally better suited for modeling problems with complicated boundaries to which it provides a piecewise linear or higher order approximation as opposed to the cruder stairstep approximation of FD. Factors Involved in Choosing Basis and Weight Functions Basis and weight function selection plays a critical role in determining the accuracy and efficiency of the resulting computer model. One goal of the basis and weight function selection is to minimize computer time while maximizing accuracy for the problem set to which the model is to be applied. Another, possibly conflicting, goal might be that of maximizing the collection of problem sets to which the model is applicable. A third might be to replicate the problem’s physical behavior with as few samples as possible. Some of the generic combinations of bases and weights that are used for MoM models are listed in Table 45.8 [Poggio and Miller from Mittra, 1973]. Basis Function Selection. We note that there are two classes of bases used in MoM modeling, subdomain and entire-domain functions. The former involves the use of bases that are applied in a repetitive fashion over subdomains or sections (segments for wires, patches for surfaces, cells for volumes) of the object being modeled. The simplest example of a subdomain basis is the single-term basis given by the pulse or stairstep function, which leads to a single, unknown constant for each subdomain. Multiterm bases involving two or more functions on each subdomain and an equivalent number of unknowns are more often used for subdomain expansions. The entire-domain basis, on the other hand, uses multiterm expansions extending over the entire object, for example, a circular harmonic expansion in azimuth for a body of revolution. As for subdomain expansions, an unknown is associated with each term in the expansion. Examples of hybrid bases can also be found, where subdomain and entire-domain bases are used on different parts of an object. Although subdomain bases are probably more flexible in terms of their applicability, they have a disadvantage generally not exhibited by the entire-domain form, which is the discontinuity that occurs at the domain boundaries. This discontinuity arises because an n s -term subdomain function can provide at most n s – 1 degrees of continuity to an adjacent basis of the unknown it represents, assuming one of the n s constants is reserved for the unknown itself. For example, the three-term or sinusoidal subdomain basis a i + b i sin(ks) + c i cos(ks) used for wire modeling can represent a current continuous at most up to its first derivative. This provides continuous charge density but produces a discontinuous first derivative in charge equivalent to a tripole charge at each junction. TABLE 45.8Examples of Generic Basis/Weight-Function Combinations Method jth Term of Basis ith Term of Weight Galerkin a j b j (r￠) w i (r) = b j (r) Least square a j b j (r￠) Q(r)?e(r)/?a i Point matching a j d(r – r j ) d(r – r i ) General collocation a j b j (r￠) d(r – r i ) Subsectional collocation U(r j )?a jk b k (r￠) d(r – r i ) Subsectional Galerkin U(r j )?a jk b k (r￠) U(r i )?b i (r) r￠ and r denote source and observation points respectively; a j , a jk are unknown constants associated with the jth basis function (entire domain) or the kth basis function of the jth subsection (subdomain); U(r k ) is the unit sampling function which equals 1 on the kth subdomain and is 0 elsewhere; b j (r￠) is the jth basis function; w i (r) is the ith testing function; d(r – r i ) is the Dirac delta function; Q(r) is a positive-definite function of position; and e(r) is the residual or equation error [from Poggio and Miller in Mitra (1973)]. ? 2000 by CRC Press LLC As additional terms are used to develop a subdomain basis, higher-order continuity can be achieved in the unknown that the basis represents, assuming still that one unknown is reserved for the amplitude of the multiterm basis function. In the general case of the n s -term subdomain basis, up to n s – 1 constants can be determined from continuity conditions, with the remainder reserved for the unknown. The kind of basis function employed ultimately determines the degree of fit that the numerical result can provide to the true behavior of the unknown for a given order of matrix. An important factor that should influence basis-function selection, then, is how closely a candidate function might resemble the physical behavior of the unknown it represents. Another consideration is whether a system of equations that is numerically easier to solve might result from a particular choice of basis and weight function, for example, by increasing its diagonal dominance so that an iterative technique will converge more rapidly and/or reduce the number of significant interactions. Various evolving approaches having names such as “impedance-matrix localization.” “fast multipole method,” “spatial decomposition,” and “multilevel matrix-decomposition” are being developed with these goals. Weight Function Selection. The simplest weight that might be used is a delta function which leads to a point- sampled system of equations, but point sampling of the field operators can be sensitive to any numerical anomalies that might arise as a result of basis function discontinuities. Distributed, multiterm weight functions can also be used on either a subdomain or an entire-domain basis to provide a further smoothing of the final equations to be solved. One example of this is the special case where the same functions are used for both the bases and weights, a procedure known as Galerkin’s method. The kind of testing function used ultimately determines the degree to which the equations can be matched for a given basis function and number of unknowns. Some specific examples of basis and weight function combinations used in electromagnetics are summarized in Table 45.9. Computing the Direct Matrix We observe that obtaining the coefficients of the direct matrix in IE modeling is generally a two-step process. The first step is that of evaluating the defining integral operator in which the unknown is replaced by the basis functions selected. The second step involves integration of this result multiplied by the weight function selected. When using delta-function weights, this second step is numerically trivial, but when using nondelta weights, such as the case in a Galerkin approach, this second step can be analytically and numerically challenging. Among the factors affecting the choice of the basis and weight functions, therefore, one of the most important is that of reducing the computational effort needed to obtain the coefficients of the direct matrix. This is one of the reasons, aside from their physical appeal, why sinusoidal bases are often used for wire problems. In this case, where piecewise linear, filamentary current sources are most often used in connection with the thin-wire approximation, field expressions are available in easily evaluated, analytical expressions. This is the case as well where Galerkin’s method is used. TABLE 45.9Examples of Specific Basis/Weight-Function Combinations Application jth Term of Basis ith Term of Weight 1D/wires Constant—a j U j (s) Delta function—d(s – s i ) 1D/wires Piecewise linear—a j1 (s – s j – d j /2) Piecewise linear—(s – s i – d i /2) + a j2 (s – s j + d j /2 ) + (s–s i + d i /2) 1D/wires 3-term sinusoidal—a j1 + a j2 sin[k(s – s j )] Delta function—d(s – s i ) + a j3 cos[k(s – s j )] 1D/wires Piecewise sinusoidal—a j1 sin[k(s – s j – d j /2)] Piecewise sinusoidal— + a j2 sin[k(s – s j + d j /2 )] sin[k(s – s i – d i /2)]+ sin[k(s – s i + d i /2)] 2D/surfaces Weighted delta function—a j d(s – s j )D j Delta function— d(s – s i ) 2D/rotational surfaces Piecewise linear axially, and exp(in?) Same (Galerkin’s method) azimuthally 2D/surfaces Piecewise linear Same (Galerkin’s method) 2D/surfaces Piecewise linear subdomain/Fourier series Same (Galerkin’s method) entire domain 3D/volumes Piecewise linear Same (Galerkin’s method) d k is the length of wire segment k; D k is the area of surface patch k. ? 2000 by CRC Press LLC Aside from such special cases, however, numerical evaluation of the direct-matrix coefficients will involve the equivalent of point sampling of whatever order is needed to achieve the desired accuracy as illustrated below. Using a wirelike one-dimensional problem to illustrate this point, we observe that at its most elementary level evaluation of the ijth matrix coefficient then involves evaluating integrals of the form (45.03) where K(s,s￠) is the IE kernel function, and s n and s￠ m are the nth and mth locations of the observation and source integration samples. Thus, the final, direct-matrix coefficients can be seen to be constructed from sums of the more elementary coefficients z(i,j,m,n) weighted by the quadrature coefficients p m and q n used in the numerical integration, which will be the case whenever analytical expressions are not available for the Z i,j . These elementary coefficients, given by w i (s n )b j (s￠ m )K(s n ,s￠ m ), can in turn be seen to be simply products of samples of the IE kernel or operator and sampled basis and testing functions. It should be apparent from this expanded expression for the direct-matrix coefficients that interchanging the basis and weight functions leaves the final problem description unchanged, although the added observation that two different IEs can yield identical matrices when using equivalent numerical treatments is less obvious. Computing the Solution Matrix Once the direct matrix has been computed, the solution can be obtained numerically using various approaches, ranging from inversion of the direct matrix to developing a solution via iteration as summarized in Table 45.10. A precautionary comment is in order with respect to the accuracy with which the solution matrix might be obtained. As computer speed and storage have increased, the number of unknowns used in modeling has also increased, from a few tens in earlier years to hundreds of thousands now when using IE models and millions of unknowns when using DE models. The increasing number of operations involved in solving these larger TABLE 45.10Summary of Operation Count for Solution of General Direct Matrix Having N x Unknowns To Obtain To Obtain Method Solution Matrix Solution Comments Cramer’s rule Expand in co-factors leading to ? ~N x ! Not an advisable procedure but useful to illustrate just how bad the problem could be Inversion N x 3 N x 2 Provides RHS-independent solution matrix Factorization N x 3 /3 N x 2 RHS-independent solution matrix Iteration General — N x 2 – N x 3 Each RHS requires separate solution With FFT — N x – N x 2 Same, plus applicability to arbitrary problems uncertain Symmetry Reflection (1 to 2 p )2(N x /2 p ) 3 N x 2 /2 p For p = 1 to 3 reflection planes Translation (Toeplitz) n x 3 [t(log 2 t) 2 ] N x 2 For n x unknowns per t sections of translation Rotation (circulant) log 2 (N x )N x 3 /n 2 N x For n rotation sectors and a complete solution m log 2 (N x )(N x /n) 3 m For m = 1 to n modes Banded General N x W 2 N x W For a bandwidth of W coefficients Toeplitz N x log 2 N x W 2 For a bandwidth of W coefficients ZwsbsKssdsd pqwsbsKss pqzijmnij N ij i C j C mninjm nm n Nij m Mij mn n N , () () (,)(,) () [()(,)] ()()(,) (,,,);, ,..., = ￠￠￠ ? ￠￠ == òò ?? == = rr 11 1 1 ( ,)(,)ij m Mij ?? =1 ? 2000 by CRC Press LLC matrices increases sensitivity of the results to roundoff errors. This is especially the case when the direct matrix is not well conditioned. It is therefore advisable to perform some sensitivity analyses to determine the direct-matrix condition number and to ascertain the possible need for performing some of the computations in double precision. Obtaining the Solution When a solution matrix has been developed using inversion or factorization, subsequently obtaining the solution (most often a current) is computationally straightforward, involving multiplication of the right-hand side (RHS) source vector by the solution matrix. When an iterative approach is used, a solution matrix is not computed, but the solution is instead developed from RHS-dependent manipulation of the direct matrix. Motivation for the latter comes from the possibility of reducing the N x 3 dependency of the direct procedure. As problem size increases, the computation cost will be increasingly dominated by the solution time. 45.5 Some Practical Considerations Although the overall solution effort has various cost components, perhaps the one most considered is the computer time and storage required to obtain the numerical results desired. With the increasing computer memories becoming available, where even microcomputers and workstations can directly address gigabytes, the memory costs of modeling are becoming generally less important than the time cost, with which we are primarily concerned here. For each model class considered, the computer-time dependence on the number of unknowns is presented in a generic formula followed by the highest-order (L/DL) term in that formula to demonstrate how computer time grows with increasing problem size. Integral Equation Modeling Frequency Domain If we consider an IE model specifically, we can show that, in general, the computer time associated with its application is dependent on the number of unknowns N x in the frequency domain as T IE,w ? A fill N x 2 + A solve N x 3 + A source N x 2 N rhs + A field N x N rhs N fields ~ (L/DL) 3(D–1) (45.4a) where the A’s are computer- and algorithm-dependent coefficients that account for computation of A fill , the direct (impedance) matrix; A solve , the solution (admittance) matrix (assuming inversion or factorization); A source , the source response (currents and charges) for one of N rhs different excitations or right-hand sides (the g term of Table 45.7); and A field , one of N field fields, where A field ￡ A fill , depending on whether a near-field (=) or far- field (<) value is obtained. D is the problem dimensionality (for a wire IE model, D = 2 except when used for wire-mesh approximations of surfaces in which case D = 3); L is a characteristic length of the object being modeled; and DL is the spatial resolution required, being proportional to the wavelength. Time Domain A similar relationship holds for a time-domain IE model which uses N t time steps, T IE,t ? A fill N x 2 + A solve N x 3 + A source N x 2 N t N rhs + A field N x N t N rhs N fields ~ (L/DL) 2(D–1)+1+r , explicit approach, 0 ￡ r ￡ 1 ~ (L/DL) 3(D–1) , implicit approach (45.4b) with the A’s accounting for computation of the time-domain terms equivalent to their frequency-domain counterparts (with different numerical values), and r accounting for the RHS dependency. Although a direct matrix may require solution initially before time-stepping the model, that is normally avoided by using dt ￡ ? 2000 by CRC Press LLC Dx/c, which yields an explicit solution, in which case A solve = 0. As can be appreciated from these expressions, the number of unknowns that is required for these computations to be acceptably accurate has a strong influence on the computer time eventually needed. Differential Equation Modeling Frequency Domain DE modeling is less commonly used in the frequency domain primarily because the order of the matrix that results depends on (L/DL) D rather than the usual (L/DL) D–1 dependency of an IE model. On the other hand, the matrix coefficients require less computation whether the DE model is based on a finite-difference or finite- element treatment. Furthermore, the matrix is very sparse because a differential operator is a local rather than a global one, as is the integral operator. Matrix fill time is therefore generally not of concern, and the overall computer time is given approximately by T DE,w ? A solve N x W 2 + A source N x WN rhs + A field N x (D–1)/D N rhs N fields ~ (L/DL) 3D–2 (45.4c) which exhibits a dominance by the matrix-solution term. Note that the banded nature of the DE direct matrix has been taken into account where the bandwidth W varies as N x 0 , N x 1/2 , and N x 2/3 respectively (N x [(D – 1)/D] ), in one, two, and three dimensions. Time Domain Time-domain DE modeling can use either implicit or explicit solution methods for developing the time variation of the solution. An explicit technique is one whereby the update at each time step is given in terms of solved- for past values of the unknowns and the present excitation, with no interaction permitted between unknowns within the same time step, and is the approach used in a technique known as finite-difference time domain (FDTD). An implicit technique, on the other hand, does allow for interaction of unknowns within the same time step but can therefore require the solution of a matrix equation. In spite of this disadvantage, implicit techniques are important because they are not subject to Courant instability when cdt > Dc as is an explicit approach. Books entirely devoted to FDTD and its applications are now becoming available, one of which is by Kunz and Luebbers (1993). The solution time for the explicit case is approximated by T DE,t ? A source N x N t N rhs + A field N x (D–1)/D N rhs N fields ~ (L/DL) D+1+r , explicit approach; 0 ￡ r ￡ 1 (45.4d) ? A solve N x W 2 + A source N x WN t N rhs + A field N x (D–1)/D N rhs N fields ~ (L/DL) 3D–2 , for D = 2, 3 and ~ (L/DL) 2+r , for D = 1, implicit approach; 0 ￡ r ￡ 1 assuming a banded matrix is used to solve the implicit direct matrix. Discussion It should be recognized that the above computer-time estimates assume solutions are obtained via matrix factorization, an N 3 process, and that iterative techniques when applicable should be expected to reduce the ? 2000 by CRC Press LLC maximum order of the (L/DL) dependency but at the cost, however, of requiring the computation to be repeated for each RHS. We also emphasize that these comparisons consider only problems involving homogeneous objects, thereby providing a more favorable situation for IE models because their sampling dimensionality d = D – 1 for a problem dimensionality of D but which increases to d = D when an inhomogeneous object is modeled. Because of these and other factors that can lead to many different combinations of formulation and numerical treatment, the foregoing results should be viewed as only generic guidelines, with the computational characteristics of each specific model requiring individual analysis to obtain numerical values for the various A x coefficients and their (L/DL) dependency. It is relevant to observe that the lowest-order size dependency for three-dimensional problems is exhibited by the DE explicit time-domain model which is on the order of (L/DL) 4 . An additional factor that should be considered when choosing among computer models is the information needed for a particular application relative to the information provided by the model. A time-domain model, for example, can intrinsically provide a frequency response over a band of frequencies from a single calculation, whereas a frequency-domain model requires repeated evaluation at each of the frequencies required to define the wideband response. Iterative solution of the direct matrix may be preferable for problems involving only one, or a few, excitations such as is the case for antenna modeling, to avoid computing all N x 2 admittances of the solution matrix when only a single column of that matrix is needed. A DE-based model necessarily provides the “near” fields throughout the region being modeled, while an IE-based model requires additional compu- tations essentially the same as those done in filling the impedance matrix once the sources have been obtained to evaluate the near fields. For applications that require modest computer time and storage, these considerations may be relatively less important than those that strain available computer resources. Clearly, the overall objective from an applications viewpoint is to obtain the needed information at the required level of accuracy for the minimum overall cost. Sampling Requirements We may estimate the number of samples needed to adequately model the spatial, temporal, and angular variation of the various quantities of interest in terms of an object characteristic length L and sampling dimension d. This may be done from knowledge of the typical spatial and temporal densities determined from computer experiments and/or from invocation of Nyquist-like sampling rates for field variations in angle as a function of aperture size. The resulting estimates are summarized in Table 45.11 and apply to both IE and DE models. These may be regarded as wavelength-driven sampling rates, in contrast with the geometry-driven sampling rates that can arise because of problem variations that are small in scale compared with l. Geometry-driven sampling would affect primarily N x , resulting in larger values than those indicated above. We note that the computer time is eventually dominated by computation of the solution matrix and can grow as (L/DL) 3 , (L/DL) 6 , and (L/DL) 9 (or f 3 , f 6 , and f 9 ), respectively, for wire, surface, and volume objects modeled using integral equations and matrix factorization or inversion. Thus, in spite of the fact that mainframe computer power has grown by a factor of about 10 6 from the UNIVAC-1 to the CRAY2, a growth that is anticipated to continue during the near future as shown in Fig. 45.2, the growth in problem size is much less TABLE 45.11Nominal Sampling Requirements for Various Field Quantities Quantity Value N x , total number of spatial samples (per scalar unknown) ~(L/DL) d = (2pL/l) d N t , number of time steps for time-domain model ~ (L/DL) = (2pL/l) N f , number of frequency steps to characterize spectral response from frequency-domain model ~ (L/2DL) = N t /2 N rhs , number of excitation sources for monostatic radar cross section in one plane a ~ (4L/DL) = 8pL/l N fields , number of far fields needed for bistatic pattern in one observation plane a ~N rhs = (4L/DL) l is the wavelength at the highest frequency of interest; DL is the spatial resolution being sought; L is object maximum object dimension or dimension in observation plane; d is the number of spatial dimensions being sampled and is not necessarily the problem dimensionality D. The distinction is important because when an appropriate Green’s function is available, the source integrals are usually one dimension less than the problem dimension, i.e., d = D – 1. An exception is an inhomogeneous, penetrable body where d = D when using an integral equation. a Assuming ~6 samples per lobe of the scattering pattern are needed. ? 2000 by CRC Press LLC impressive, as illustrated by Fig. 45.3. The curves on this graph demonstrate emphatically the need for finding faster ways of performing the model computations, a point further emphasized by the results shown in Fig. 45.4 where the computer time required to solve a reference problem using various standard models is plotted as a function of frequency. FIGURE 45.2Raw and smoothed FLOP (floating-point operation) rates of mainframe computers and smoothed rate-of- change in speed, at year of introduction from the UNIVAC-1 to the projected performance of an advanced parallel system at year 2000. Future growth is increasingly dependent on computer architecture, requiring increasing parallelism as improve- ments due to component performance reach physical speed limits. FIGURE 45.3Time development of IE-based modeling capability for one-dimensional (e.g., a wire), two-dimensional (e.g., a plate), and three-dimensional (e.g., a penetrable, inhomogeneous cube) sampling of a problem of characteristic dimension L in wavelengths and matrix order N solvable in 1 h of computer time using mainframe computers introduced in the years indicated. Linear-systems solution using LU decomposition (an N 3 dependency) is assumed with number of unknowns proportional to L, L 2 and L 3 , respectively, without any problem symmetry being exploited. These results should be viewed as upper bounds on solution time and might be substantially reduced by advances in linear-system solution procedures. ? 2000 by CRC Press LLC 45.6 Ways of Decreasing Computer Time The obvious drawback of direct moment-method models as N x increases with increasing problem size and/or complexity suggests the need for less computationally intensive alternatives. There are various alternatives for decreasing the computer cost associated with solving electromagnetic problems using the method of moments. The basic intent in any case is either to reduce the direct cost of performing a given modeling computation or to reduce the number of modeling computations needed to obtain a desired result. An example for achieving the latter is to employ lower-order models that can accurately enough represent the behavior of observables as a function of space, angle, frequency, or time so that the sampling density of the first-principles model can be reduced. A specific example is used of the rational functions to approximate frequency-domain transfer func- tions. These might include analytical, computational, and experimental approaches or combinations thereof, about which further discussion and references may be found in Miller [1988]. 45.7 Validation, Error Checking, and Error Analysis Modeling Uncertainties The process of proceeding from an original physical problem to computed results is one that is subject to numerous uncertainties caused by a variety of factors. Perhaps foremost among these factors is the degree of arbitrariness associated with many of the choices that are made by the code developer and/or modeler in the course of eventually obtaining numerical results. Whereas the numerical evaluation of classical boundary-value problems such as scattering from a sphere is numerically robust in the sense that different workers using different FIGURE 45.4 An illustration of the frequency dependence of the CRAY2 computer time required for some standard computer models applied to the reference problem of a perfectly conducting, space-shuttle-sized object having a surface area of 540 m 2 [Miller, 1988]. At a sampling density of 100/l 2 , a total of 6,000 surface samples is assumed for an IE model at 100 MHz for which LU decomposition of the direct (impedance) matrix requires about 1 h of CRAY2 time. The top, LU, curve has a slope of f 6 as discussed in the text. The next two curves have slopes of f 4 , the upper corresponding to use of a TD DE model (FDTD) as well as an iterative solution of the direct IE matrix, assuming acceptable convergence occurs in 100 iteration steps. The third curve is for a 10-step iterative solution of the IE matrix. The bottom three curves have f 2 slopes. The upper two of these are for 100- and 10-step iterative solutions used in connection with a near-neighbor approximation (NNA), wherein only the 100 and 10 largest interaction coefficients are retained in the matrix, respectively. The bottom curve is for the physical-optics approximation, in which the induced current is computed from the incident magnetic field. The effects of these different frequency slopes on the computer time can be seen to be extreme, emphasizing the need for developing more efficient solution procedures. ? 2000 by CRC Press LLC computers and different software can obtain results in agreement to essentially as many significant figures as they wish, the same observation cannot be made for moment-method modeling. Modeling uncertainties can be assigned to two basic error categories, a physical modeling error e P and a numerical modeling error e N as outlined in Table 45.12. The former is due to the fact that for most problems of practical interest varying degrees of approximation are needed in developing a simplified or idealized problem representation that will be compatible with the computer code to be used for the modeling computations. The latter is due to the fact that the numerical results obtained are almost invariably only approximate solutions to that idealized representation. We note that although an analytical expression may in principle represent a formally exact solution, the process of obtaining numerical results in that case is still one that inevitably involves finite-precision evaluation of the formal solution. By its very nature, the physical modeling error requires some kind of measurement for its determination, except for those few problems whose analytical solution in principle involves no physical idealization or subsequent numerical approximation. One example of such problems is that of determining the scattering or radiating properties of the perfectly conducting or dielectric sphere. The numerical modeling error is itself composed of two components in general, the determination of which would normally involve one or more kinds of computation. The first and generally more important of these components is the solution error that arises because the computer model used, even if solved exactly, would not provide an exact solution for the idealized problem representation. The solution error arises essentially because the computer model is solved using a finite number of unknowns. The other, generally less important contributor to the numerical modeling error is the equation error that arises because the numerical results obtained from the computer model used may not numerically satisfy the modeling equations. The equation error may be caused both by roundoff due to the computer word size as well as the solution algorithm used, as in the case of iteration, for example. The impact of equation error can be expected to increase with increasing condition number of the direct matrix. Validation and Error Checking One of the most time consuming and long lasting of the tasks associated with any model development is that of validation. Long after work on the model has been completed, questions will continue to arise about whether a given result is valid or whether the model can be applied to a given problem. There are essentially two kinds of validation procedures that can be considered to answer such questions: (1) internal validation, a check that can be made concerning solution validity within the model itself; and (2) external validation, a check that utilizes information from other sources which could be analytical, experimental, or numerical. Existing computer models often do not perform internal checks on the results they produce but instead leave that as an exercise for the user. It would be of extremely great potential value if a variety of such checks could be built into the code and exercised as desired by the modeler. The topic of error checking and validation is an active one in CEM and receives a great deal of ongoing attention, for which the technical literature provides a good point of departure for the reader interested in more detail. TABLE 45.12Error Types that Occur in Computational Electromagnetics Category Definition Physical modeling error, e P Arises because the numerical model used is normally an idealized mathematical representation of the actual physical reality Numerical modeling error, e N Arises because the numerical results obtained are only approximate solutions to that idealized representation and consists of two components: (1)Solution error—The difference that can exist between the computed results and an exact solution even were the linear system of equations to be solved exactly, due to using a finite number of unknowns (2)Equation error—The equation mismatch that can occur in the numerical solution because of roundoff due to finite-precision computations or when using an iterative technique because of limited solution convergence ? 2000 by CRC Press LLC 45.8 Concluding Remarks In the preceding discussion we have presented a selective survey of computational electromagnetics. Attention has been directed to radiation and scattering problems solved using the method of moments, a general procedure applicable to differential- and integral-equation formulations developed by either the frequency domain or the time domain. Beginning from the viewpoint of electromagnetics as a transfer-function process, we concluded that the basic problem is one of developing source-field relationships or field propagators. Of the various ways by which these propagators might be expressed, we briefly discussed the Maxwell curl equations and Green’s- function source integrals as providing the analytical basis for moment-method computer models. We then considered at more length some of the numerical issues involved in developing a computer model, including the idea of sampling functions used both to represent the unknowns to be solved for and to approximate the equations that they must satisfy. Some of the factors involved in choosing these sampling functions and their influence on the computational requirements were examined. Next, we discussed some ways of decreasing the needed computer time based on either analytical or numerical approaches. Some closing comments were directed to the important problems of validation, error checking, and error analysis. Throughout our discussion, emphasis has been given to implementation issues involved in developing and using computer models as opposed to exploring analytical details. Defining Terms Computer model: Based on a numerical solution of some appropriate analytical formulation that describes a physical phenomenon of interest. The model is realized in an algorithm or computer code that reduces the formulation to a series of operations suitable for computer solution. Field propagator: The analytical description of how electromagnetic fields are related to the sources that cause them. Common field propagators in electromagnetics are the defining Maxwell equations that lead to differential equation models, Green’s functions that produce integral equation models, optical prop- agators that lead to optics models, and multipole expansions that lead to modal models. Integral equation: An analytical relationship in which the quantity whose solution is sought (the unknown) appears under an integral sign. When this is the only place that the unknown appears, the integral equation is commonly called a first-kind equation, while if the unknown also appears outside the integral, it is a second-kind integral equation. Method of moments: A general technique for reducing integral, differential (including partial), and integrod- ifferential equations to a linear system of equations or matrix. The moment method involves discretizing, sampling, and approximating the defining equations using basis or expansion functions to replace the unknown and testing or weighting functions to satisfy the defining equations. The matrix that results may be full (all coefficients nonzero) or sparse (only a few per row are nonzero), depending on whether the model is an integral or differential equation. Modeling errors: In essentially all computer modeling, there are two basic kinds of errors. One, the physical modeling error, arises from replacing some real-world physical problems with some idealized mathemat- ical representation. The other, the numerical modeling error, comes from obtaining only an approximate solution to that idealized representation. Usually, the numerical modeling error can be reduced below the physical modeling error if enough unknowns, i.e., a large enough matrix, are used to model the problem of interest. Sampling: The process of replacing a continuous physical quantity by some sequence of sampled values. These values are associated with the analytical function used to approximate the behavior of the physical quantity whose solution is sought and are the unknowns of the moment-method matrix. Sampling is also involved in determining how well the defining equations are to be satisfied. A common approach for equation sampling is point sampling, where the equations are explicitly satisfied at a series of discrete points in some prescribed region of space. Unknown sampling can involve localized basis functions, an approach called subdomain sampling, while if the basis functions reside over the entire region occupied by the unknown, the approach is called entire-domain sampling. ? 2000 by CRC Press LLC Solution domain: Electromagnetic fields can be represented as a function of time, or a time-domain descrip- tion, or as a function of frequency using a (usually) Fourier transform, which produces a frequency- domain description. Related Topics 35.1 Maxwell Equations ? 44.3 Numerical Methods ? 44.4 Modern Design Environment References M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions, Applied Mathematics Series, vol. 55, Washington, D.C.: National Bureau of Standards, 1964. F.S. Acton, Numerical Methods that Work, New York: Harper and Row, 1970. B.W. Arden and K.N. Astill, Numerical Algorithms: Origins and Applications, Reading, Mass.: Addison-Wesley, 1970. R.F. Harrington, Time-Harmonic Electromagnetic Fields, New York: McGraw-Hill, 1961. R.F. Harrington, Field Computation by Moment Methods, New York: Macmillan, 1968. K.S. Kunz and R.J. Luebbers, The Finite Difference Time Domain Method for Electromagnetics, Boca Raton, Fla.: CRC Press, 1993. E. K. Miller, “A selective survey of computational electromagnetics,” IEEE Trans. Antennas Propagat., vol. AP- 36, pp. 1281–1305, 1988. E.K. Miller, “Solving bigger problems—by decreasing the operation count and increasing the computation band- width,” invited article in special issue of IEEE Proc. Electromagnets, vol. 79, no. 10, pp. 1493–1504, 1991. E.K. Miller, L. Medgyesi-Mitschang, and E.H. Newman, Computational Electromagnetics: Frequency-Domain Method of Moments, New York: IEEE Press, 1991. R. Mittra, ed., Computer Techniques for Electromagnetics, New York: Pergamon Press, 1973. R. Mittra, ed., Numerical and Asymptotic Techniques in Electromagnetics, New York: Springer-Verlag, 1975. J. Moore and R. Pizer, Moment Methods in Electromagnetics: Techniques and Applications, New York: Wiley, 1984. A.J. Poggio and E.K. Miller, “Low frequency analytical and numerical methods for antennas,” in Antenna Handbook, Y. T. Lo and S. W. Lee, eds., New York: Van Nostrand Reinhold, 1988. B.D. Popovic, M.B. Dragovic, and A.R. Djordjevic, Analysis and Synthesis of Wire Antennas, Letchworth, Hert- fordshire, England: Research Studies Press, 1982. W.H. Press, B.R. Flannery, S.A. Teukolsky, and W.T. Vettering, Numerical Recipes, London: Cambridge University Press, 1986. B.J. Strait, ed., Applications of the Method of Moments to Electromagnetic Fields, St. Cloud, Fla.: SCEEE Press, 1980. J.A. Stratton, Electromagnetic Theory, New York: McGraw-Hill, 1941. W. L. Stutzman and G. A. Thiele, Antenna Theory and Design, New York: John Wiley, 1981. P. L. E. Uslenghi, ed., Electromagnetic Scattering, New York: Academic Press, 1978. J.H. Wang, Generalized Moment Methods in Electromagnetics, New York: Wiley Interscience, 1991. Further Information The International Journal of Numerical Modeling, published by Wiley four times per year, includes numerous articles on modeling electronic networks, devices, and fields. Information concerning subscriptions should be addressed to Subscription Department, John Wiley & Sons Ltd., Baffins Lane, Chichester, Sussex PO19 1UD, England. The Journal of the Acoustical Society of America is published by the American Institute of Physics on a monthly basis. Most issues contain articles about the numerical solution of acoustics problems which have much in common with problems in electromagnetics. Information about the society and journal can be obtained from Acoustical Society of America, 500 Sunnyside Blvd., Woodbury, NY 11797. ? 2000 by CRC Press LLC The Journal of the Applied Computational Electromagnetics Society is published two or three times a year, accompanied by a newsletter published about four times per year. The focus of the society and journal is the application of computer models, their validation, information about available software, etc. Membership and subscription information can be obtained from Dr. R.W. Adler, Secretary, Applied Computational Electromag- netics Society, Naval Postgraduate School, Code ECAB, Monterey, CA 93943. The Journal of Electromagnetic Waves and Applications is published by VNU Science Press. It contains numerous articles dealing with the numerical solution of electromagnetic problems. Information about the journal can be obtained from its editor-in-chief, Professor J.A. Kong, Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139. The Proceedings of the IEEE, Transactions on Microwave Theory and Techniques of the IEEE, Transactions on Antennas and Propagation of the IEEE, and Transactions on Electromagnetic Compatibility of the IEEE all are periodicals published by the Institute of Electrical and Electronics Engineers, about which information can be obtained from IEEE Service Center, 445 Hoes Lane, PO Box 1331, Piscataway, NJ 08855-1331. ? 2000 by CRC Press LLC