|Computational electromagnetics A review
S. M. Rao and N.
Department of EE, Auburn University, Auburn, AL
*Department of Aerospace, Indian Institute of
Science, Bangalore 560 012, India
In this paper, we describe the important
computational methods available to successfully predict the radar cross-section of a given
object. Further, we evaluate these methods based on the applicability to general
structures, efficiency and accuracy. We confine our attention to the integral equation
(IE) based methods although some passing references were made to the differential equation
COMPUTATIONAL Electromagnetics (CEM) has evolved rapidly during the past decade
to a point where extremely accurate predictions can be made for very general scattering
and antenna structures. In general, all the available methods may be classified broadly
into two categories, viz. a) differential equation (DE) solution methods and b) integral
equation (IE) methods.
Although the Maxwell curl equations are usually
first encountered in the time domain (TD), i.e. with time as an explicit, independent
variable, until relatively recently, most electromagnetic instruction and research has
taken place in the frequency domain (FD) where time-harmonic behaviour is assumed. A
principal reason for favouring the FD over the TD in the pre-computer era had been that a
FD approach was generally more tractable analytically. Furthermore, the experimental
hardware available for making measurements in past years was largely confined to the FD.
The inferior position of TD electromagnetics (EM)
began to change with the arrival of the digital computer, which has not only profoundly
affected what can be done numerically (or computationally), but also experimentally. Since
the beginning of what has come to be called computational electromagnetics (CEM) in the
early 1960s, there has been a steady growth in both TD and FD modeling. This growth, which
began slowly at first, was primarily confined to integral-equation (IE) treatments, but it
has become almost explosive over the last 10 years as TD differential-equation (DE)
modeling has attracted wide attention. This presentation summarizes the status of
computational electromagnetic modeling and highlights some of the current research areas.
Modeling choices in CEM
In discussing CEM, it is appropriate to consider two
- What are the various alternative modeling approaches available for
- What are the advantages of one model relative to the other
To answer both questions, we observe that there are
four major, first-principles, models in CEM, given by,
- Time Domain Differential Equation (TDDE) models, the use of which has
increased tremendously over the past several years, primarily as a result of much larger
and faster computers.
- Time Domain Integral Equation (TDIE) models, although available for
well over 30 years, have gained increased attention in the last decade. The recent
advances in this area make these methods very attractive for a large variety of
- Frequency Domain Integral Equation (FDIE) models which remain the
most widely studied and used models, as they were the first to receive detailed
- Frequency Domain Differential Equation (FDDE) models whose use has
also increased considerably in recent years, although most work to date has emphasized low
These four choices can actually be narrowed down to
two choices, i.e. a) IE models and b) DE models, depending on the mathematical
formulation. The well-known method of moments (MoM) in general, involves IE modeling
whereas the well-known finite element method (FEM) uses DE formulation.
General aspects of CEM modeling
The formulation and numerical development of a CEM
model, in general, involves a number of basic steps whether a DE or an IE approach is
being followed. In the following, we discuss some of these considerations.
Model development. For any numerical solution,
it is necessary to develop the required equations and solve them on a computer. The
equations thus developed must include the physics of the problem as well as the
geometrical features. The following four steps are carried out in CEM problems:
- Develop integral equations using potential theory along with
appropriate boundary conditions or alternatively, begin with the time-dependent Maxwell
curl equations or their equivalent to develop methods such as FD, TD or FEM.
- Sample these equations in space, and also in time if it is a
time-dependent equation, utilizing an appropriate geometrical space grid and suitable
basis and testing functions. Note that, depending on the choice of formulation, the space
grid may cover the structure and/or the surrounding space.
- Develop a set of simultaneous equations relating known and unknown
quantities. Generally, the known and unknown quantities are the excitation field or its
derivatives and the radiated/scattered field or induced current and charge, respectively.
- Generate a computer solution of this system in space and time as an
Comparison of DE and IE models
The two primary choices for CEM modeling, as
indicated so far, are those based on DEs and IEs. Whatever the details of the specific
approach, any numerical method in its most general form, provides a way of solving IEs and
DEs, involving approximating integrals as finite sums, and derivatives as finite
differences in the generic forms
leading, after some additional manipulation, to a
linear system of equations, or system matrix. The process of discretizing and
quantifying DEs numerically is known by various names including finite-difference,
finite-area (or volume), and finite-element procedures. The term finite-element is
usually, but not necessarily, associated with a variational formulation, while use of a
designation other than finite difference usually refers to the use of more general basis
and testing functions. A numerical model based on an IE is also called a boundary
element method in structural dynamics and acoustics.
Most computer modeling involves replacing an
infinite domain, first principles, analytical description of a problem by a finite domain,
discretized, numerical one. The numerical model is finite in nature because only a limited
number of unknowns of limited precision can be used in the solution process. An analytical
model, however, entails, symbolically at least, an infinite dimensionality such as that
exhibited by a series expansion for a sphere. However, it is worth noting that, from a
practical viewpoint, the analytical model is finite in nature as well because the process
of quantitatively evaluating any analytical model is automatically subject to limited
precision and accuracy. This means that any observable of interest will exhibit no
quantitative change over some specified dynamic range after an appropriate number of terms
have been summed in its series solution. This will be the case whenever we deal with numerical
answers as opposed to analytical solutions.
Some basic differences between DE and IE models are
- In general, the differential equation methods generate a sparse
matrix, while the integral equation methods generate full matrices.
- Homogeneous/inhomogeneous/anisotropic materials can be handled in a
relatively simple manner, while the level of complexity for the integral equation methods
varies enormously for each of these cases.
- The code generation is straight forward for DE methods. This is
usually not the case for integral equation methods.
- For DE methods, the solution space includes the objects
surroundings, the radiation condition is not enforced in exact sense, thus leading
to certain error in the solution. For the IE solution, the solution space is confined to
the object and the radiation condition is automatically enforced.
- The IE solutions are generally more accurate and efficient.
- Spurious solutions exist in DE methods whereas such solutions are
absent in IE methods.
- For the DE solutions, developing numerical solution using parallel
computer architecture is easy. However, for IE solutions, this is possible only in the
time-domain. A lot of time and effort is needed to generate a parallel version for
the frequency-domain IE solution.
Integral equation solutions in CEM
Mathematically speaking, an equation involving the
integral of an unknown function of one or more variables is known as integral equation.
One of the most common integral equations encountered in electrical engineering is the convolution
integral given by
In eq. (2), we note that the response function Y(t)
and the system function H(t, t ) is known and we need to determine the input
X(t ). Of course, if X(t ) and H(t, t ) are known and we need
to determine Y(t), then eq. (2) simply represents a integral relationship
which can be performed in a straightforward manner. We further note that H(t,
t ) is also commonly known as impulse response if eq. (2) represents the system
response of a linear system. In general, in mathematics and in engineering literature, H(t,
t ) is known as Greens function or kernel function. We also acknowledge that,
for some other physical systems Y(t) and X(t) may represent
the driving force and response functions, respectively.
Next, we note that eq. (2) is known as integral
equation of first kind. We also have another type of integral equation given by
where C1 and C2
In eq. (3), we note that the unknown function X(t)
appears both inside and outside the integral sign. Such equation is known as the integral
equation of second kind. Further, we also see in electrical engineering yet another
type of integral equation given by
which is known as integro-differential equation.
It may be noted that for a limited number of kernel
and response functions, in eqs. (24), it is possible to obtain the solution using
analytical methods. However, for a majority of practical problems, these equations can be
solved using numerical methods only. Fortunately, in this day and age, we can obtain very
accurate numerical solutions owing to the availability of fast digital computers. In the
following section, we discuss a general numerical technique, popularly known as Method
of Moments, to solve the integral equations (2)(4).
Method of moments solution
The method of moments (MoM) solution procedure was
first applied to electromagnetic scattering problems by Harrington1. Consider a
linear operator equation given by
AX = Y, (5)
where A represents the integral operator, Y
is the known excitation function and X is the unknown response function to be
determined. Now, let X be represented by a set of known functions, termed as basis
functions or expansion functions (p1, p2, p3,...,)
in the domain of A as a linear combination:
where a is are scalar
constants to be determined. Substituting eq. (6) into eq. (5), and using the linearity of A,
where the equality is usually approximate. Let (q1,
q2, q3,...) define a set of testing functions in the
range of A. Now, multiplying eq. (7) with each qj and using the
linearity property of the inner product, we obtain
for j = 1, 2, ... , N. The
set of linear equations represented by eq. (8) may be solved using simple matrix methods
to obtain the unknown coefficients a i.
The simplicity of the method lies in choosing the
proper set of expansion and testing functions to solve the problem at hand. Further, the
method provides a most accurate result if properly applied. While applying the method of
moments to complex practical problems, the solution region, in general, is divided into
triangular subdomains, as shown in Figure 1. Then, one can define suitable basis and
testing functions and develop a general algorithm to solve the electromagnetic problem2.
In Figure 2, we present the current induced on the
aircraft shown in Figure 1 using the numerical procedure described in ref. 2, when
illuminated by a 300 MHz electromagnetic plane wave. The plane wave is polarized
along the length of the aircraft and travelling perpendicular to the body. Efficient
solutions have been obtained for very complex problems using these methods in
electromagnetics and acoustics3. It is possible that these methods found
applications in other areas of engineering. Lastly, using the method of moments, solutions
have been obtained for initial value problems also4.
Fast multipole method
One major problem with MoM is the generation of a
dense matrix and for certain problems, the dimension of this matrix can be prohibitively
large. Usually, for electromagnetic and acoustic scattering problems, it is necessary to
divide the solution region into small enough
subdomains in order to
obtain accurate results. By small enough, we mean about 200300
subdomains per square wavelength. In usual practice, we may typically solve for several
thousand unknowns for large, complex problems. Quickly, this requirement becomes
expensive, in terms of computational resources, and may even become impossible to handle.
Hence, we look for alternate schemes to reduce the computational resources.
The fast multiple method (FMM) dramatically reduces
the time and memory required to compute radar cross sections and antenna radiation
patterns compared to dense matrix techniques5. It is fairly simple to implement
the FMM in a method of moments program to compute the electromagnetic scattering from
large bodies of arbitrary shape. In the following, we describe the essential steps
involved in the FMM implementation.
In the method of moments solutions to boundary
integral equations, one is faced with solving large systems of equations of the form
ZI = V,
where Z is a dense matrix and I
and V are column vectors of length N, where N is the number of
current expansion functions. Eq. (9) can be solved by a number of iterative schemes. These
techniques involve the computation of the product of Z and a solution vector
one or more times for each iteration. This operation takes O(N2)
operations and usually dominates all other operations in the iterative loop.
In the FMM, one groups the N basis functions
into M groups so that the basis functions in each group have neighbouring support.
For the simplest single stage FMM, the optimum value for M is proportional to ÷ N.
Let the index m run over the groups and the index a refer to a basis function
within a particular group. The dense matrix Z is then replaced with the
Z = ZĘ + VTV+,
where ZĘ , V, and T
are sparse matrices. ZĘ are those components of the original Z
matrix for interactions between nearby regions of the target (typically within about one
wavelength). The approximation can be made arbitrarily precise by the appropriate choice
of FMM parameters in the computation of V and T.
The components of V are given by
where f are the basis functions. V
is evaluated at K = ÷ N angles of k needed for a
quadrature over the surface of a sphere.
The sparse matrix T is
where RmmĘ is the distance
between the centers of the group m and mĘ . The number of terms in the sum,
L, is chosen to give the desired accuracy in the FMM expansions.
Recently, the FMM algorithm was implemented on
massively parallel architectures such as Intel Paragon. These machines typically consist
of several hundred fast RISC microprocessors interconnected by a communication network.
Let T(1) be the time required to compute ZI with one processor and T(P)
be the time required with P processors. Ideally, one should have T(P)/T(1) = O(P).
This is often difficult or impossible to achieve due to inherently sequential portions of
the algorithm and communication costs.
For the FMM, the essential problem is to find an
optimal distribution of the data structures. This is done by assigning one group to each
A logical extension of FMM is the development of
multilevel FMM algorithm. Further, we have yet another technique to improve the
computational speed known as adaptive integral method (AIM). The details of these
algorithms may be obtained from refs 6, 7.
Sparse matrix methods
The generation of a sparse matrix in the method of
moment solution procedure may be achieved in two ways, viz. a) by defining a special set
of basis functions to represent the unknown quantity or b) by handling the influence of
the kernel function in a novel way. The usage of well-known wavelet-type basis functions
to provide the required sparsity belong to the former category8 and the
application of fast multipole method (FMM) belongs to the latter category5. So
far, the wavelet-type basis functions have been applied to integral equations with one
variable only and it remains to be seen how these functions can be utilized for two or
more variable case. In contrast, in the FMM scheme, the matrix-vector product is carried
out in a novel way and seem to work well for more complex problems.
There is also yet another scheme, known as impedance
matrix localization scheme (IML) which achieves modest sparsity for simple problems9.
Notice that the kernel function is, in general, a decaying function with respect to the
distance between the source and observation points. Thus, with increasing distances, the
influence of a given source becomes negligible at a sufficiently distant observation point
and may be actually set to zero. The IML scheme cleverly exploits this fact. However,
there is a certain degree of arbitrariness in this scheme and seems to work for simple
Recently a new method, known as generalized sparse
matrix reduction scheme (GSMR), is proposed which seem to improve on the IML method10.
The basic concept utilized in the GSMR technique may be qualitatively illustrated as
follows. Following similar procedures of the MoM, a moment matrix is also generated in the
GSMR method. However, in contrast to the conventional moment method where interaction is
computed from each and every cell on other cells, only the interaction from the self-cell
and few neighbouring cells is computed in the GSMR technique. In fact, for single variable
problems (wire scatterers and two-dimensional infinite cylinders) only the self-term and
two neighbouring terms on either side of the self-cell are generated in this technique.
However, this technique, although appears promising, needs to be validated for more
In this paper, we have described the use of integral
equation based methods in computational electromagnetics. In general, these methods are
applicable to scatterers whose characteristic dimensions are of the order of a wavelength.
The iterative methods and FMM are then called in to extend the IE methods to scatterers of
larger dimension. Since the integral equation methods are global in nature, these methods
work very well for scatterers with smooth geometries. The local nature of the DE
formulation may be exploited to cater for the sharp variations, thus creating hybrid
techniques. Both DE and IE, even with the present day super computers, are suited for
computing the RCS of a fighter aircraft at best up to 1 GHz. For higher frequencies,
the recourse is often taken to asymptotic methods like the geometrical theory of
diffraction (GTD), uniform theory of diffraction (UTD), and uniform asymptotic theory
(UAT). In essence, the CEM has grown today to be a mature tool for predicting precisely
the scattering and radiation characteristics of most complex structures encountered in
real life from human body to aircraft.
Harrington, R. F., Field Computation by Moment Methods,
Macmillan, New York, 1968.
Rao, S. M., Wilton, D. R. and Glisson, A. W., IEEE Trans. Antenna
Propagation, 1982, 30, 409418.
Miller, E. K., MedgyesiMitschang, L. and Newman, E. H., Computational
Electromagnetics Frequency-Domain Method of Moments, IEEE Press, New
Rao, S. M., Time Domain Electromagnetics, Academic Press, New
Coifman, R., Rokhlin, V. and Wandzura, S., IEEE AP-S Mag.,
1993, 35, 712.
Song, J. M. and Chew, W. C., Microwave Optical Technol. Lett.,
1995, 10, 1419.
Bleszynski, E., Bleszynski, M. and Jaroszewicz, T., IEEE AP-S
International Symposium Digest. Seattle, WA, June 1994, pp. 416419.
Steinberg, B. Z. and Leviatan, Y., IEEE Trans. Antenna Propagation,
1993, 41, 610619.
Canning, F. X., IEEE Trans. Antenna Propagation, 1993, 41,
Rao, S. M. and Gothard, G. K., Microwave Optical Technol. Lett.,
1998, 19, 271274.