HIGH FREQUENCY COMPUTATION IN WAVE EQUATIONS AND OPTIMAL
DESIGN FOR A CAVITY
By
Jun Lai

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

Applied Mathematics - Doctor of Philosophy
2013

ABSTRACT
HIGH FREQUENCY COMPUTATION IN WAVE EQUATIONS AND
OPTIMAL DESIGN FOR A CAVITY
By
Jun Lai
Two types of problems are studied in this thesis.
One part of the thesis is devoted to high frequency computation. Motivated by fast
multiscale Gaussian wavepacket transforms and multiscale Gaussian beam methods which
were originally designed for initial value problems of wave equations in the high frequency
regime, we develop fast multiscale Gaussian beam methods for wave equations in bounded
convex domains in the high frequency regime. To compute the wave propagation in bounded
convex domains, we have to take into account reflecting multiscale Gaussian beams, which
are accomplished by enforcing reflecting boundary conditions during beam propagation and
carrying out suitable reflecting beam summation. To propagate multiscale beams efficiently,
we prove that the ratio of the squared magnitude of beam amplitude and the beam width
is roughly conserved, and accordingly we propose an effective indicator to identify significant beams. We also prove that the resulting multiscale Gaussian beam methods converge
asymptotically. Numerical examples demonstrate the accuracy and efficiency of the method.
The second part of the thesis studies the reduction of backscatter radar cross section
(RCS) for a cavity embedded in the ground plane. One approach for RCS reduction is
through the coating material. Assume the bottom of the cavity is coated by a thin, multilayered radar absorbing material (RAM) with possibly different permittivities. The objective
is to minimize the backscatter RCS by the incidence of a plane wave over a single or a set
of incident angles and frequencies. By formulating the scattering problem as a Helmholtz

equation with artificial boundary condition, the gradient with respect to the material permittivities is determined efficiently by the adjoint state method, which is integrated into
a nonlinear optimization scheme. Numerical example shows the RCS may be significantly
reduced.
Another approach is through shape optimization. By introducing a transparent boundary
condition, the unbounded scattering problem from a cavity is reduced to a bounded domain
problem. RCS reduction for the cavity is formulated as a shape optimization problem involving the Helmholtz equation. The existence of the minimizer is proved under an appropriate
constraint. Descent directions of the objective function with respect to the boundary may
be found via the domain derivative. It is used in a gradient-based optimization scheme to
find the optimal shape of the cavity. Numerical examples show that the RCS is effectively
reduced at different incident frequencies.

ACKNOWLEDGMENTS

Foremost, I would like to express my deepest gratitude to my advisor Prof. Gang Bao for his
continuous support through out my Ph.D study, for his personality, enthusiasm, patience,
and immense knowledge. His guidance is of great help to me during the time of my research
and writing of this thesis .
My sincere thanks goes to Prof. Jianliang Qian for guiding my research for the past
several years and helping me develop my skills in both mathematical analysis and computer
programming.
Besides them, I would also like to thank the rest of my thesis committee: Prof. Zhengfang
Zhou, Prof. Andrew Christlieb, and Prof. Yingda Cheng, for their encouragement, insightful
comments, and valuable references.
I would like to give special thanks to my group members: Junshan Lin, Yuliang Wang,
Justin Droba, Hai Zhang, Eric Wolf, Yuqi Hong, Guanghui Hu and Xiang Xu for the inspiring
discussions and all the fun we have had in the past few years.
I would also like to thank my parents, my younger brother and younger sister. They were
always supporting me and encouraging me with their best wishes.
Finally, I would like to thank my wife, Jing Li. She was always there cheering me up and
standing by me through the good and bad times.

iv

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

vii

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii

Chapter 1 Introduction . . . . . .
1.1 Gaussian Beam Method . .
1.2 Optimal design for a cavity
1.2.1 Coating Material . . .
1.2.2 Shape Optimization .
1.3 Outline . . . . . . . . . . . . .

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

1
1
2
4
4
5

Chapter 2 Gaussian beam method applied to wave equations . .
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Gaussian beam methods for the wave equation . . . . . . . .
2.2.1 Gaussian beams for initial value problems . . . . . . .
2.2.2 Incident and reflected beams for the wave equation .
2.3 Fast Multiscale Gaussian beams . . . . . . . . . . . . . . . . .
2.3.1 Basic setup . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.2 Initialization and propagation for the wave equation
2.4 Numerical strategies for treating bounded domains . . . . .
2.4.1 Beam initial decomposition . . . . . . . . . . . . . . . .
2.4.2 Beam summation . . . . . . . . . . . . . . . . . . . . . .
2.4.3 Identifying significant beams . . . . . . . . . . . . . . .
2.4.4 Overall algorithm . . . . . . . . . . . . . . . . . . . . . .
2.5 Convergence results . . . . . . . . . . . . . . . . . . . . . . . . .
2.6 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6.1 One-dimensional domain . . . . . . . . . . . . . . . . . .
2.6.1.1 Case 1: Linear velocity . . . . . . . . . . . . . .
2.6.2 Two-dimensional rectangular domains . . . . . . . . .
2.6.2.1 Case 1: Linear velocity . . . . . . . . . . . . . .
2.6.2.2 Case 2: Sinusoidal velocity . . . . . . . . . . .
2.6.3 Two-dimensional circular domains . . . . . . . . . . . .
2.6.3.1 Case 1: Sinusoidal velocity . . . . . . . . . . .
2.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

6
6
9
9
13
17
18
21
24
24
25
27
31
32
38
39
39
40
40
41
41
42
42

Chapter 3 Radar cross section reduction of a cavity: TM case . . . . . . .
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

62
62

v

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

3.2 Problem formulation and objective function . . . . . . . . . .
3.2.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Objective function . . . . . . . . . . . . . . . . . . . . . .
3.3 Gradient by the adjoint state method . . . . . . . . . . . . . .
3.4 Optimization method: Sequential Quadratic Programming
3.5 Direct solver for the forward scattering problem . . . . . . .
3.6 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . . .
3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Chapter 4 Radar cross section reduction of a cavity: TE case
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Problem formulation and objective function . . . . . . . .
4.2.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . .
4.2.2 Objective function . . . . . . . . . . . . . . . . . . . .
4.3 Gradient by the adjoint state method . . . . . . . . . . . .
4.4 Optimization algorithm . . . . . . . . . . . . . . . . . . . . .
4.5 Direct solver for the forward scattering problem . . . . .
4.6 Numerical Experiments . . . . . . . . . . . . . . . . . . . . .
4.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

63
63
67
68
72
75
76
81

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

. 85
. 85
. 86
. 86
. 90
. 91
. 95
. 97
. 101
. 104

Chapter 5 Optimal shape design of a cavity for radar cross section reduction110
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
5.2 Problem description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
5.3 Existence of the minimizer . . . . . . . . . . . . . . . . . . . . . . . . . . 118
5.4 Domain Derivative . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
5.5 The optimization algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 129
5.6 Numerical experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
5.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
Chapter 6 Summary and future work . .
6.1 Summary . . . . . . . . . . . . . . . .
6.1.1 Gaussian beam method . . .
6.1.2 Optimal design for a cavity
6.1.2.1 Coating material . .
6.1.2.2 Shape optimization
6.2 Future work . . . . . . . . . . . . . .

BIBLIOGRAPHY

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

. . . . . . . . . . . .

vi

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

140
140
140
141
141
141
142

. . . . . . . . . . . . . . . . 144

LIST OF TABLES

Table 3.1

Table 3.2

Table 4.1

A comparison for the gradients obtained from the adjoint state method
(ASM) and the finite difference method (FDM). . . . . . . . . . . .

71

Relative permittivies after the optimization, number of iterations and
CPU time for all the three examples. . . . . . . . . . . . . . . . . .

80

Relative permittivities after the optimization, number of iterations
and CPU time for all the three examples. . . . . . . . . . . . . . . . 104

vii

LIST OF FIGURES

Figure 2.1

Odd periodic expansion for a circular domain.(For interpretation of
the references to color in this and all other figures, the reader is referred to the electronic version of this dissertation.) . . . . . . . . .

26

Figure 2.2

(a)Tail exceeds the boundary; (b)Image principle . . . . . . . . . . .

27

Figure 2.3

Example 1. Linear model in 1D. (a) Comparison between the beam
solution and exact solution at ω = 620; (b) Comparison between the
beam solution and exact solution at ω = 980; (c) Windowed comparison between the beam solution and exact solution at ω = 620; (d)
Windowed comparison between the beam solution and exact solution
at ω = 980. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

Figure 2.3

(cont’d) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

Figure 2.4

Relative L2 errors for linear velocity in Example 1 . . . . . . . . . .

46

Figure 2.5

Example 2. Linear model in 1D. (a) Comparison between the beam
solution and exact solution at ω = 620; (b) Comparison between the
beam solution and exact solution at ω = 980; (c) Windowed comparison between the beam solution and exact solution at ω = 620; (d)
Windowed comparison between the beam solution and exact solution
at ω = 980. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47

Figure 2.5

(cont’d)

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

48

Figure 2.6

Relative L2 errors for linear velocity in Example 2 . . . . . . . . . .

49

Figure 2.7

Example 3. Linear model in 2D. (a) The exact solution at ω = 180;
(b) The beam solution at ω = 180; (c) Comparison for the slice at
x = 0.25; (d) Comparison for the slice at y = 0.125; (e)Windowed
comparison for the slice at x = 0.25; (f )Windowed comparison for
the slice at y = 0.125. . . . . . . . . . . . . . . . . . . . . . . . . . .

50

Figure 2.7

(cont’d) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

Figure 2.7

(cont’d) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52

viii

Figure 2.8

Relative L2 errors for linear velocity in Example 3 . . . . . . . . . .

53

Figure 2.9

Example 4. Sinusoidal model in 2D. (a) The exact solution at ω =
220; (b) The beam solution at ω = 220; (c) Comparison for the slice
at x = 0.25; (d) Comparison for the slice at y = 0.125;(e) Windowed
comparison for the slice at x = 0.25; (f ) Windowed comparison for
the slice at y = 0.125. . . . . . . . . . . . . . . . . . . . . . . . . . .

54

Figure 2.9

(cont’d) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

55

Figure 2.9

(cont’d) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

56

Figure 2.10

Relative L2 errors for sinusoidal velocity in Example 4

. . . . . . .

57

Figure 2.11

Example 5. Sinusoidal model in 2D. (a) The exact solution at ω =
180; (b) The beam solution at ω = 180; (c) Comparison for the slice
at x = 0.45; (d) Comparison for the slice at y = 0.45; (e)Windowed
comparison for the slice at x = 0.45; (f ) Windowed comparison for
the slice at y = 0.45. . . . . . . . . . . . . . . . . . . . . . . . . . . .

58

Figure 2.11

(cont’d) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

59

Figure 2.11

(cont’d) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60

Figure 2.12

Relative L2 errors for sinusoidal velocity in Example 5

. . . . . . .

61

Figure 3.1

Geometry of the problem . . . . . . . . . . . . . . . . . . . . . . . .

64

Figure 3.2

RCS reduction for k0 = 16π. (a) Optimized at θ = π/2,(b) Optimized
at θ = 3π/5.(c) Optimized at θ = 5π/6. . . . . . . . . . . . . . . . .

78

Figure 3.2

(cont’d) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

79

Figure 3.3

RCS reduction for k0 = 32π. (a) Optimized at θ = π/2; (b) Optimized at θ = 3π/5; (c) Optimized at θ = 5π/6. . . . . . . . . . . . .

82

Figure 3.3

(cont’d) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

83

Figure 3.4

RCS reduction by the combination of θ = 3π/5 and θ = 5π/6. (a)
Optimized for k0 = 16π, (b) Optimized for k0 = 32π. . . . . . . . .

84

ix

Figure 4.1

Geometry of the problem . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4.2

RCS reduction for k0 = 16π. (a)Optimized at θ = 3π/5.(b)Optimized
at θ = 5π/6.(c)Optimized with the combination of θ = 3π/5 and
θ = 5π/6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

Figure 4.2

(cont’d)

Figure 4.3

RCS reduction for k0 = 32π. (a)Optimized at θ = 3π/5.(b)Optimized
at θ = 5π/6.(c)Optimized with the combination of θ = 3π/5 and
θ = 5π/6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

Figure 4.3

(cont’d)

Figure 4.4

RCS reduction for k between 16π and 32π. (a)Optimized at θ = π/2.
(b)Optimized at θ = 3π/5. . . . . . . . . . . . . . . . . . . . . . . . 109

Figure 5.1

Geometry of the problem . . . . . . . . . . . . . . . . . . . . . . . . 112

Figure 5.2

Constraint of the problem . . . . . . . . . . . . . . . . . . . . . . . . 117

Figure 5.3

Initial shape . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

Figure 5.4

(a)Optimal shape for k = 2π at normal incidence; (b) RCS comparison for k = 2π, ”- - -” for initial shape, ”—–” for optimal shape. . . 136

Figure 5.5

(a)Optimal shape for k = 8π at normal incidence; (b) RCS comparison for k = 8π, ”- - -” for initial shape, ”—–” for optimal shape. . . 137

Figure 5.6

(a)Optimal shape for k = 32π at normal incidence; (b) RCS comparison for k = 32π, ”- - -” for initial shape, ”—–” for optimal shape. . 138

Figure 5.7

87

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

(a)Optimal shape for k = 32π at the incidence angle θ = 120◦ ; (b)
RCS comparison for k = 32π, ”- - -” for initial shape, ”—–” for
optimal shape. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139

x

Chapter 1
Introduction
This thesis considers two types of problems: high frequency computation for wave equations
and optimal design for a cavity. These two problems are connected in the way that the
technique applied in the high frequency computation can be used to design a fast solver for
the forward problem in the cavity optimization, especially for large cavities. In this thesis,
we will discuss them separately and leave the combination as a future work.

1.1

Gaussian Beam Method

Wave equations in the high frequency regime lead to computational difficulties in the sense
that a large number of grid points are required when applying direct methods such as finitedifference or finite-element methods. The size of the linear system resulted from these
methods increases rapidly as the frequency increases, which often exceeds the capability
of the computational source. Consequently, alternative methods such as geometrical-optics
based asymptotic methods are sought to compute such high-frequency wave phenomena.
One of the powerful geometrical-optics methods is the Gaussian beam method [4, 42, 45, 51],
which is able to treat caustics automatically.
The idea underlying Gaussian beams is to construct a solution that asymptotically satisfies the partial differential equation by superposing the so called beam solutions[45]. Each
beam solution is concentrated on a ray and the amplitude exponentially decays away from

1

the ray, where a ray is defined as the projected characteristic curve resulted from a Hamiltonian system onto the physical domain. The pure mathematics community has known the
existence of such solutions since sometime in the 1960s [4]. In [45], these solutions have
been used to obtain results on propagation of singularities in hyperbolic PDEs. Numerical
applications for Gaussian beams have been successfully established in many equations, including Helmholtz [33] and Schr¨dinger equations [31, 43, 32]. Other applications includes
o
seismic wave modeling [15], migration [22] and mountain wave simulation [51]. It has been
shown that the numerical implementations based on ray-centered coordinates are sometimes
computationally inefficient [15, 22]. Therefore, a purely Eulerian computational approach
based on [45, 51] was proposed in [33] that overcomes some of these difficulties. Extension
includes both high frequency waves and semi-classical quantum mechanics [31, 32].
In this thesis, motivated by the work in [44] which designed fast multiscale Gaussian
wavepacket transforms and multiscale Gaussian beam methods for pure initial value problems
of wave equations, we propose to develop fast multiscale Gaussian beam methods for wave
equations in bounded convex domains. Several difficulties are encountered when dealing with
bounded domain problems. Details are discussed in Chapter 2 as well as the corresponding
solutions for those issues.

1.2

Optimal design for a cavity

Radar cross section(RCS) is an important measure for the delectability of a target by radar
systems. The RCS from a cavity is significant since the overall RCS of a target is often
dominated by some cavities, such as the jet inlet of an aircraft. Therefore, effective reduction
of the RCS from a cavity has been an important problem in wave propagation with many

2

practical applications. In this thesis, we focus on a 2-D rectangular cavity which is embedded
in the ground plane and illuminated by a time-harmonic plane wave.
A great deal of studies have been done on both direct scattering and inverse scattering
problems for such a cavity. Mathematical analysis for the scattering problem can be found
in [2, 1, 3]. In particular, questions on the existence and uniqueness for the solution of the
inverse scattering problem have been investigated in [5]. When the cavity is large and deep,
the solution is highly oscillatory, which presents both stability issue and computational challenge. The stability of the solution with the dependence on the wave number is studied in
[10] for a 2D rectangular cavity in the TM polarization. In terms of numerical simulation,
standard methods including finite difference [9], finite element with boundary integral [25]
and method of moment(MoM) [24] have been developed. Mode matching methods have also
been proposed in [6],[11]. Note that the large cavity problem is equivalent to a high wave
number problem which is known to be computationally challenging. Therefore, high frequency techniques, such as bounding and shooting rays [16] and Gaussian beam asymptotics
[14] were also applied to the large cavity scattering problem. For non-planar opening cavities
or so called overfilled cavities, mathematical analysis and numerical simulations can been
seen in [52, 23]
Several approaches can lead to the reduction of RCS, including coating material, shaping,
passive cancellation and active cancellation[28]. All the methods have pros and cons. In
practice, the most common way to achieve this goal is through coating material and shaping,
both of which are considered in this thesis.

3

1.2.1

Coating Material

In this approach, assume a thin, multilayered radar absorbing material(RAM) is coated
horizontally at the bottom of the cavity for the reduction of backscatter RCS. Our objective
is to determine the permittivity of the material such that RCS is minimized over a range of
angles or a set of frequencies.
Many methods have been proposed on RCS reduction with RAM [37],[40],[17]. A popular
one is based on genetic algorithm (GA), which is a gradient-free optimization method and
widely used in engineering community. The problem is, even for a small parameter space,
GA usually requires thousands of generations to obtain the best population, or the global
optimal. While designing a fast algorithm for the direct scattering problem already presents
a huge challenge especially when the cavity is large and deep [25],[35], it is computationally
unaffordable to employ GA to find the optimal synthesized RAM of a cavity.
We apply a gradient based optimization method called Sequential Quadratic Programming(SQP) to the design problem, where the gradient is obtained through the adjoint state
method. Numerical experiments show the method is very efficient for RCS reduction, both
in TM(transverse magnetic) polarization and TE(transverse electric) polarization.

1.2.2

Shape Optimization

Properly shaping the cavity is another approach to achieve RCS reduction. Shape optimization has been studied extensively in many problems[47, 41]. When considering the shape
optimization for a cavity, one needs to take into account some practical constraints of the
cavity, otherwise a bizarre shape may be obtained. In this thesis, we assume the cavity is
area invariant and the boundary of the cavity does not oscillate too much. Exact description

4

of these constraints can be seen in Chapter 5. The main idea is to find the domain derivative for RCS with respect to the boundary of the cavity. A simple optimization scheme is
proposed with the integration of domain derivative. Numerical examples are given at the
last section of Chapter 5 to show the effectiveness of the algorithm.

1.3

Outline

The thesis is outlined as follows. In Chapter 2, Gaussian beam method applied to wave
equations in the bounded domain has been discussed. Chapter 3 discusses the optimal
design of a cavity for RCS reduction through the coating material in TM polarization,
and Chapter 4 continues the study on TE polarization. Chapter 5 proposes an optimization
scheme to design the proper shape of a cavity. All the algorithms presented in these chapters
are illustrated with numerical experiments. Chapter 6 concludes the whole discussion and
presents some future work.

5

Chapter 2
Gaussian beam method applied to
wave equations

2.1

Introduction

We consider the following initial-boundary value problem (IBVP) of the wave equation,


 u − V 2 (x)∆u = 0, x ∈ D, t > 0
 tt








u|t=0 = f1 (x)












(2.1)

ut |t=0 = f1 (x)
u|x∈∂D = 0

where D is a convex bounded domain in Rd , and the velocity V (x) is smooth, positive and
bounded away from zero. Compared to the slowly changing velocity function, initial conditions f1 (x) ∈ H 1 (D) and f2 (x) ∈ L2 (D) are assumed to be highly oscillatory, compactly
supported functions.
Many papers [33, 51, 49, 31, 50, 32, 43, 44, 38, 46] are devoted to the Gaussian-beam
based numerical methods for pure initial value problems of Schr¨dinger equation and wave
o
equations, while efficient numerical Gaussian beam method developed for wave equations
and Schr¨dinger equations in bounded domains seems to be rare. Some essential difficulties
o

6

arise in developing numerical Gaussian beams in the case of wave equations in bounded
domains. Here is a list of them:
1. One needs to take care of reflected beams. The theory to construct reflecting beams
has been first addressed in [45] and further detailed in [12]. Numerical attempts to the
construction of reflected beams must consider the geometry of the domain.
2. A bounded domain of general geometry may give rise to diffraction phenomena or
gliding rays along the boundary in the sense of geometrical optics, although such
phenomena is usually weak for a regular domain.
3. The existing initialization algorithms are generally established under the assumption of unboundedness of the domain. For instance, Fourier-Bros-Iagolnitzer(FBI)
transform[32], direct decomposition[33], and the multiscale Gaussian wavepackets transform introduced in [44].
4. The effective area of a beam is unbounded despite its exponential decay. How to extract
the beam solution in the original bounded convex domain becomes an issue.
Those difficulties need to be overcome for the design of an efficient numerical Gaussian
beam method in bounded domains. We therefore provide the following solutions for those
issues:
1. Reflected beams are constructed by enforcing the homogeneous Dirichlet boundary
boundary condition and the tangential continuity for the phase.
2. Since geometrical-optics based approaches including Gaussian beams are not able to
capture those effects, other theories, such as geometrical theory of diffraction, FourierAiry integrals, or gliding beams, are needed. Therefore, to avoid those potential issues,
7

we will assume that the domain is strictly convex and non-grazing hypothesis [45] holds,
and the latter will be satisfied when the initial data is compactly supported away from
the boundary.
3. We apply the multiscale Gaussian wavepackets transform to the initialization due to its
efficiency. Since the original fast multiscale Gaussian wavepacket transform is designed
for periodic functions, we propose to first carry out periodic, odd extensions of the
initial data, then apply the multiscale Gaussian wavepacket transform to the resulting
extended data.
4. A method-of-images based approach has been proposed for the summation step.
In addition, we also address the issue on how to identify significant beams so that beam
propagation can be carried out more efficiently. We prove that the ratio of the squared
magnitude of beam amplitude and the beam width is roughly conserved for each individual
beam, and this ratio can be used as an indicator to identify significant beams. This way the
number of propagated beams is significantly reduced.
The rest of the chapter is organized as follows. Section 2.2 introduces the Gaussian Beam
method applied to the wave equation within a bounded domain in the setting of Lagrangian
formulation. In Section 2.3, we briefly review the fast multiscale Gaussian wavepacket transform and its application in wave equation, details of which can be found in [44]. Section 2.4
gives the strategy about how to select the significant beams and addresses some other numerical issues. Section 2.5 carries out the convergence proof. Numerical results are provided
in Section 2.6 as well as some special treatments around the boundary for a circular domain.
Section 2.7 concludes the whole discussion.

8

2.2
2.2.1

Gaussian beam methods for the wave equation
Gaussian beams for initial value problems

We start from the initial value problem of the scalar wave equation in Rd :

utt − V 2 (x)∆u = 0,

x ∈ Rd , t > 0,

(2.2)

where V (x) is smooth, positive and bounded away from zero. Initial conditions u(0, x) =
f1 (x) ∈ H 1 (Rd ) and ut (0, x) = f2 (x) ∈ L2 (Rd ) are highly oscillatory functions.
We are looking for asymptotic solutions of the wave equation in geometrical-optics form,

A(x, t)eıωτ (x,t) ,

where τ (x, t) is the phase function, A(x, t) the amplitude function, and ı =

(2.3)
√

−1. In the

ansatz (2.3), the frequency ω is a large parameter, and an asymptotic solution for the wave
equation is sought in the sense that the wave equation (2.2) and its associated initial conditions are satisfied approximately with a small error when ω is large. After substituting the
ansatz (2.3) into the wave equation (2.2) and considering the leading orders in inverse powers
of the large parameter ω, we end up with the following eikonal and transport equations:

2
τt − V 2 (x)|∇x τ (x, t)|2 = 0,

(2.4)

2At τt − 2V 2 ∇x A · ∇x τ + A(τtt − V 2 trace(τxx )) = 0.

(2.5)

9

Factorizing the eikonal equation (2.4) gives

±
τt + G± (x, ∇x τ (x, t)) = 0,

(2.6)

where G± (x, ∇x τ (x, t)) = ±V (x)|∇x τ (x, t)| correspond to two polarized wave modes in the
second-order wave equation. Accordingly, we define the Hamiltonians,

G± (x, p) = ±V (x)|p|,

where G± (x, p) is clearly homogeneous of degree one in the momentum variable p.
To construct asymptotic solutions for the wave equation, we are going to use Gaussian
beams [45, 36, 51]. Because the two polarized wave modes will be treated essentially in the
same way, we consider the following generic situation for the eikonal equation:

τt + G(x, ∇x τ (x, t)) = 0,

(2.7)

where G can be taken to be either G+ or G− and τ to be either τ + or τ − . According to the
Gaussian beam theory [45, 36, 51], a single Gaussian beam is an asymptotic solution to the
wave equation, and it is concentrated near a ray path which is the x-projection of a certain
bicharacteristic. To construct a bicharacteristic, we apply the method of characteristics to
the eikonal equation (2.7) to obtain the following Hamiltonian system:

dx
= Gp , x|t=0 = x0 ,
dt
dp
p =
˙
= −Gx , p|t=0 = p0 ,
dt

x =
˙

10

(2.8)
(2.9)

where t is time parameterizing bicharacteristics. Solving this system yields the bicharacteristic
{(x(t), p(t)) : t ≥ 0},
which emanates from the initial point (x0 , p0 ) in phase space at t = 0. The corresponding
ray path is γ = {(x(t), t) : t ≥ 0}, which is defined in the (x, t)-space. Notice that along the
ray path γ = {(x(t), t) : t ≥ 0}, we have by construction p(t) = τx (x(t), t) due to the method
of characteristics. Furthermore, the phase function τ (x(t), t) along the ray path satisfies

dτ (x(t), t)
= τt (x(t), t) + p(t) · Gp (x(t), p(t)) = τt (x(t), t) + G(x(t), τx (x(t), t)) = 0,
dt

which implies that the phase function τ (x(t), t) does not change along γ because the Hamiltonian G is homogeneous of degree one; we will take τ (x(t), t) = 0.
So far we have computed the phase function τ and its first-order derivative p(t) =
τx (x(t), t) along the ray path γ = {(x(t), t) : t ≥ 0}. To construct a second-order Taylor expansion for the phase function along the ray path, one needs to compute the Hessian of
the phase along the ray. Following [45, 51], we differentiate the eikonal equation (2.7) with
respect to t and x near the ray path γ:

τt,x (x, t) + Gx (x, τx (x, t)) + τxx (x, t)Gp (x, τx (x, t)) = 0,

(2.10)

τt,t (x, t) + Gp (x, τx (x, t)) · τx,t (x, t) = 0.

(2.11)

Differentiating equation (2.10) further with respect to x yields the following Riccati equation

11

for M (t) = τxx (x(t), t):
dM (t)
+ Gxx + M (t)Gxp + GT M (t) + M (t)Gpp M (t) = 0,
xp
dt

(2.12)

which is appended with an initial condition M |t=0 = M0 = ıϵI, where ϵ is a positive number
of order O(1).
Although the Riccati equation (2.12) does not admit a global smooth solution in general,
it turns out that complexifying the equation by specifying a complex initial value will guarantee that a global smooth solution exists because of the underlying symplectic structure
associated with the related Hamiltonian system; see [45, 36, 51] for theoretical justification.
Now with the Hessian of the phase function at our disposal, we may solve the transport
equation (2.5) for the amplitude A(t) = A(x(t), t) along the ray path γ. Since τt (x(t), t) =
−G(x(t), p(t)) along the ray path, the transport equation (2.5) is reduced to the following:
)
dA A(x(t), t) ( 2
+
V (x(t))trace(M (t)) − Gx · Gp − GT M (t)Gp = 0,
p
dt
2G

(2.13)

which is appended with a suitable initial condition A|t=0 = A0 .
At this stage, we are ready to construct a single Gaussian beam along the ray path γ by
defining the following two global, smooth approximate functions for the phase and amplitude:

1
τ (x, t) ≡ p(t) · (x − x(t)) + (x − x(t))T M (t)(x − x(t)),
2

(2.14)

A(x, t) ≡ A(x(t), t) = A(t),

(2.15)

which are accurate near the ray path γ = {(x(t), t) : t ≥ 0}. These two functions allow us

12

to construct a single-beam asymptotic solution

Φ(x, t) = A(x, t) exp(ıωτ (x, t)).

(2.16)

This beam solution is concentrated on a single smooth curve γ = {(x(t), t) : t ≥ 0}, which is
the x-projection of the bicharacteristic {(x(t), p(t)) : t ≥ 0} emanating from (x0 , p0 ) at t = 0.
1
Because the phase τ (x, t) has an imaginary part, Im(τ (x, t)) = 2 (x − x(t))T Im(M (t))(x −
x(t)), Φ(x, t) has a Gaussian profile of the form
( ω
)
T Im(M (t))(x − x(t)) ,
exp − (x − x(t))
2

which is concentrated on the smooth ray path γ.

2.2.2

Incident and reflected beams for the wave equation

So far we have discussed how to apply Gaussian beam methods to the initial value problem
of the wave equation, but the domain we are interested in is bounded. The boundary of the
domain requires us to construct a reflected beam when an incident beam hits the boundary.
The derivation here relies on the results of [45].
Before continuing our discussion, we have to assume the non-grazing hypothesis:

xs (t0 ) · ν(x(t0 )) > 0
˙

(2.17)

where xs (t0 ) denotes the direction of the incident ray, t0 is the time when the incident ray
˙
hits the boundary ∂D at location x(t0 ), and ν denotes the exterior normal vector to ∂D. In
other words, the ray will not propagate along the boundary.
13

We further denote the incident and reflected beams by:

us = As (x, t)eıωτs (x,t)

(2.18)

ur = Ar (x, t)eıωτr (x,t)

(2.19)

In order to satisfy the zero boundary condition, we require:

(us + ur )|(x(t ),t ) = 0
0 0

(2.20)

Substituting equations (2.18),(2.19) in (2.20) yields:

As (x, t0 )eıωτs (x,t0 ) = −Ar (x, t0 )eıωτr (x,t0 )

(2.21)

The independence of ω would force:

τs (x, t0 ) = τr (x, t0 )
As (x, t0 ) = −Ar (x, t0 )

(2.22)
(2.23)

From equation (2.22), we impose the condition that all of their tangential and time derivatives
should be continuous at (x(t0 ), t0 ). Therefore, differentiating with respect to t on both sides
of equation (2.22) and using the eikonal equation (2.7) give us:

V (x(t0 ))|ps (x(t0 ), t0 )| = V (x(t0 ))|pr (x(t0 ), t0 )|

14

(2.24)

In order to make the reflected beam incoming, we need pr ̸= ps . Hence we have

pr = (I − 2νν T )ps ,

(2.25)

where ν is the outward normal at the reflection point.
For example, in 1-D it implies:
pr = −ps ;

(2.26)

in 2-D, assuming that α is the angle between the tangential line of the boundary and the
positive x2 -axis if x = (x1 , x2 ), we have

 
 

2 α − cos2 α −2 sin α cos α
 p1,s 
p1,r  sin
,
=
·

2 α − sin2 α
−2 sin α cos α cos
p2,s
p2,r

(2.27)

where the outward normal at the reflection point is defined to be (cos β, sin β) with β = π−α.
The second-order derivatives of the phase function for the reflected beam can be determined by differentiating twice with respect to t in equation (2.22) and combining equations
(2.10) and (2.11), and it follows that

Gp (x(t0 ), ps ) · Gx (x(t0 ), ps ) + GT (x(t0 ), ps )Ms Gp (x(t0 ), ps )
p
= Gp (x(t0 ), pr ) · Gx (x(t0 ), pr ) + GT (x(t0 ), pr )Mr Gp (x(t0 ), pr )
p

(2.28)

Substituting all the related quantities into equation (2.28) plus the continuity of the tangential components of the second order derivatives of τ (x, t), we can determine the relation
between the Hessian matrix Ms of the incident beam and Mr of the reflected beam.
For example, we give those relations in three cases:
15

• 1-D interval :
Mr = Ms + 2ps

Vx
.
V

(2.29)

• 2-D rectangular domain: let x = (x1 , x2 ). Assume that two sides of the rectangle are
parallel to the x1 -axis and the other two sides are parallel to the x2 -axis:
– If the beam hits the boundary parallel to the x1 -axis, then


 
p2 +p2 Vx
1,s 2,s 1

M11,r , M12,r  M11,s + 2 p1,s
V , −M12,s 


=

M21,r , M22,r
−M21,s ,
M22,s

(2.30)

– If the beam hits the boundary parallel to the x2 -axis, then

 

M
,
−M12,s

M11,r , M12,r   11,s


=

2 +p2
p1,s 2,s Vx 
2
M21,r , M22,r
−M21,s , M22,s + 2 p
V
2,s

(2.31)

Here we simply ignore the situation when the beam hits the corner of the rectangle,
since it causes diffraction and formula (2.20) does not apply any more. Since the
Gaussian method is asymptotic, the numerical accuracy will not be degraded without
those beams as those diffractions have exponentially small effects.
• 2-D circular domain: Consider a unit disk with the boundary parameterized by angle
θ, then:





M11,s 
M11,r 







−1
−1 · K · 
M12,r  = Kr
s M12,s  + Kr · B








M22,r
M22,s

16

(2.32)

where



sin2 θ

Kr

Ks

B

cos2 θ

−2 sin θ cos θ






= p sin θ p sin θ − p cos θ −p cos θ
2,r
1,r
2,r
 1,r



2
2
p1,r
2p1,r p2,r
p2,r


2θ
2θ
−2 sin θ cos θ
cos

 sin




= p sin θ p sin θ − p cos θ −p cos θ
2,s
1,s
2,s

 1,s


2
2
p1,s
2p1,s p2,s
p2,s


(p1,r − p1,s ) cos θ + (p2,r − p2,s ) sin θ






=

0




Vx1
Vx2
2 + p2 ) · (
(p1,s
2,s
V (p1,s − p1,r ) + V (p2,s − p2,r ))

(2.33)

(2.34)

(2.35)

p1,s , p2,s , p1,r and p2,r are defined in (2.27). By symmetry, M21,r = M12,r , so all the
elements are determined.
From equations (2.29)-(2.32), one can prove the property that the imaginary part from
Ms to Mr is still symmetric positive definite [12]. We thus have all the initial components
for the reflected beams. The propagation of theses reflected beams follow the same equations
in Section 2.1.

2.3

Fast Multiscale Gaussian beams

In order to construct a solution for the wave equation, it is also necessary for the asymptotic solution to satisfy the initial condition. However, since the initial condition may not
have the form of a Gaussian wave packet, we have to decompose the general initial profile into a superposition of Gaussian wave packets. Here we apply fast multiscale Gaussian
17

wavepacket transforms to initialize the beam propagation for the wave equation, resulting
in fast multiscale Gaussian beams for wave equations in bounded domains.

2.3.1

Basic setup

We follow [44] closely. Let N be a sufficiently large positive integer in the sense that
[−N/2, N/2]d is enough to cover the spectra of the initial conditions f1 ∈ H 1 (D) and
f2 ∈ L2 (D) in the Fourier domain. For simplicity, we assume that the domain is D = [0, 1]d .
Without loss of generality, N is assumed to be the power of 2. We only consider the discrete
version of the transform here. Define the spatial grid and Fourier grid to be
)
(n n
1 , 2 , · · · , nd : 0 ≤ n , n , · · · , n < N, n , n , · · · , n ∈ Z}
1 2
1 2
d
d
N N
N
N
N
Ω = {ξ = (ξ1 , ξ2 , · · · , ξd ) : − ≤ ξ1 , ξ2 , · · · , ξd < , ξ1 , ξ2 , · · · , ξd ∈ Z}
2
2

X = {x =

Partition the Fourier domain Ω into Cartesian coronae Cl for l ≥ 1 as follows:

C1 = [−4, 4]d
Cl = {ξ = (ξ1 , ξ2 , · · · , ξd ) : max |ξs | ∈ [4l−1 , 4l ]}, l ≥ 2
1≤s≤d
Each corona Cl is further partitioned into boxes:

Bl,i =

d
∏

[2l · is , 2l · (is + 1)]
s=1

where i = (i1 , i2 , · · · , id ) is ranging over all possible choices which satisfy Bl,i ⊂ Cl . To each
box Bl,i , we associate a smooth and compactly supported function gl,i with size Ll = 2W l ,

18

where W l = 2l is the length of box Bl,i . The function is approximately defined as
(

gl,i (ξ) ≈ e

|ξ−ξl,i |
−
σl

)2
, ξ∈Ω

(2.36)

where σl = W l /2, ξl,i is the center of box Bl,i . Based on gl,i (ξ), we can define the conjugate
filter hl,i (ξ) for each Bl,i :
gl,i (ξ)
hl,i (ξ) = ∑
, ξ∈Ω
2
l,i gl,i (ξ)

(2.37)

It follows that the products of gl,i (ξ) and hl,i (ξ) form a partition of unity:
∑

gl,i (ξ)hl,i (ξ) = 1.

l,i
We now define two sets of functions φl,i,k (x) and ψl,i,k (x), which are both Gaussian
wavepackets. Their constructions are based on gl,i (ξ) and hl,i (ξ), respectively.
In the Fourier domain, they are defined by:
k·ξ
1 −2πı L
l gl,i (ξ), k ∈ {0, 1, · · · , Ll − 1},
e
φl,i,k (ξ) =
ˆ
d/2
L
l
k·ξ
1 −2πı L
ˆ
l hl,i (ξ), k ∈ {0, 1, · · · , Ll − 1}.
ψl,i,k (ξ) =
e
d/2
L
l

19

(2.38)

(2.39)

In the spatial domain, they can be numerically evaluated by:
k
∑ 2πı(x− L )ξ
l gl,i (ξ), k ∈ {0, 1, · · · , Ll − 1}
φl,i,k (x) =
e
(N Ll )d/2 ξ∈Ω
k
∑ 2πı(x− L )ξ
1
l hl,i (ξ), k ∈ {0, 1, · · · , Ll − 1}
e
ψl,i,k (x) =
(N Ll )d/2 ξ∈Ω
1

(2.40)

(2.41)

Here the subtraction in the spatial domain must be understood as modulus the periodic
domain [0, 1]n .
The forward multiscale Gaussian wavepacket transform for a given discretized function
f on X is defined by:
k·ξ
1 −2πı L
ˆ
ˆ
ˆ
l hl,i (ξ)f (ξ),
cl,i,k =< ψl,i,k , f >=< ψl,i,k , f >=
e
d/2
ξ∈Ω Ll
∑

(2.42)

ˆ
where f (ξ) is the discretized Fourier transform of f , and < ·, · > denotes the L2 inner
product. It is proved that f ∈ L2 (D) can be expressed as [44]:

f (x) =

∑

cl,i,k φl,i,k (x).

(2.43)

l,i,k
From the definition (2.40) of φl,i,k (x), it approximately equals:

φl,i,k ≈

(√

)d 2πı(x− k )ξ −σ 2 π 2 |x− k |2
π
Ll l,i e l
Ll , x ∈ X,
σ
e
N Ll l

(2.44)

k
which can be taken as a Gaussian wavepacket centered at L , with frequency ξl,i . Together
l
with (2.43), it gives us the way to decompose a general L2 (D) function into a supposition
of Gaussian wave packets. Throughout the decomposition, the main computational cost is
20

on Fourier transform which can be accomplished in O(N d log N ); see [44].

2.3.2

Initialization and propagation for the wave equation

We first assume that D is a Cartesian domain in Rd . With all the preparation above, we
are now ready to employ multiscale Gaussian wavepacket transform to initialize the function
f1 ∈ H 1 (D) and f2 ∈ L2 (D). Following [44], we decompose the initial data by the forward
transform (2.42):

f1 (x) =
f2 (x) =

∑

al,i,k φl,i,k (x)

l,i,k
∑

(2.45)

bl,i,k φl,i,k (x)

(2.46)

l,i,k
Assume that the global asymptotic solution for equation (2.1) has the form:

uGB (x, t) =

∑

(c+ Φ+ (x, t) + c− Φ− (x, t))
l,i,k l,i,k
l,i,k l,i,k

(2.47)

l,i,k
where ′ +′ and ′ −′ represent two different wave modes. If we suppress the superscript ′ ±′ ,
Φl,i,k (x, t) is a Gaussian beam propagating in the space-time domain with initial condition
φl,i,k (x). It can be obtained by solving the following equation:

x = Gp
˙

(2.48)

p = −Gx
˙

(2.49)

˙
M = −(Gxp )T M − M Gxp − M Gpp M − Gxx
)
A ( 2
˙
A=−
V (x)trace(M ) − Gx · Gp − GT M Gp
p
2G
21

(2.50)
(2.51)

with the initial conditions:

k
x|t=0 =
Ll

(2.52)

ξl,i
p|t=0 = 2π
|ξl,i |

(2.53)

2
σl
I
|ξl,i |
)d
(√
π
σ
A|t=0 =
N Ll l

M |t=0 = ı · 2π 2

(2.54)
(2.55)

where I means the identity matrix. According to Section 2.2, at the reflection point we have
to set:

xr |t=t = xs |t=t
0
0

(2.56)

pr |t=t = r1 (ps )
0

(2.57)

Mr |t=t = r2 (Ms )
0

(2.58)

Ar |t=t = −As |t=t .
0
0

(2.59)

Expressions of r1 (·) and r2 (·) depend on the geometry of the domain. In particular, the
explicit form of r1 (ps ) is given in (2.26),(2.27) and r2 (Ms ) is given in (2.29),(2.30), (2.32)
for three special cases.
The Gaussian Beam solution corresponding to φl,i,k (x) is given by
ı·|ξ |τ
(x,t)
Φl,i,k (x, t) = Al,i,k (x, t)e l,i l,i,k

22

(2.60)

with

1
τl,i,k (x, t) ≡ pl,i,k (t)(x − xl,i,k (t)) + (x − xl,i,k (t))T Ml,i,k (x − xl,i,k (t))
2
Al,i,k (x, t) ≡ Al,i,k (x(t), t) = Al,i,k (t)
Now we only have to determine the coefficients c± to complete the construction. From
l,i,k
(2.47) and (2.45), letting t = 0 would yield

c+ φl,i,k (x, t) + c− φl,i,k (x, t) = al,i,k φl,i,k (x)
l,i,k
l,i,k

(2.61)

Differentiate (2.47) with respect to t on both sides, use equation (2.13) and let t → 0:
(
)
D0
+ − c− )
+ (x, p(0)) φ
(c
+ ı · |ξl,i |G
l,i,k (x) = ut |t=0 = bl,i,k φl,i,k (x)
l,i,k
l,i,k 2G+ (x(0), p(0))
(2.62)
where
D0 = c2 M (0) − Gp (x(0), p(0))(Gx (x(0), p(0)) + M (0)Gp (x(0), p(0))).
Since |ξl,i | represents the frequency of f1 and f2 , compared to D0 , |ξl,i | is fairly large.
In this sense, D0 is negligible and G(x, p(0)) can also be approximated by G(x(0), p(0)) due
to the narrow support of φl,i,k (x). Therefore, it is reasonable to have
(
)
+ − c− ) ı · |ξ |G+ (x, p(0)) φ
(c
l,i
l,i,k (x) ≈ ut |t=0 = bl,i,k φl,i,k (x)
l,i,k
l,i,k

23

(2.63)

Solving equation (2.61) and (2.63), we get
(
)
bl,i,k
+ =1 a
−
c
l,i,k 2 l,i,k ı · G+ ( k , 2πξ )
l,i
Ll
(
)
bl,i,k
− =1 a
+
c
l,i,k 2 l,i,k ı · G+ ( k , 2πξ )
l,i
Ll

(2.64)

(2.65)

Once all the coefficients c± and the related quantity x± (t), p± (t), M ± (t), A± (t)
l,i,k
l,i,k
l,i,k
l,i,k
l,i,k
are available, the global asymptotic solution can thus be determined by formula (2.47).

2.4
2.4.1

Numerical strategies for treating bounded domains
Beam initial decomposition

Note that the multiscale Gaussian wavepacket transform is designed for periodic functions.
We can simply expand f1 and f2 periodically before the transform is applied, but doing this
would ’mix’ the function at the two sides of D. Regarding the zero boundary condition, a
more appropriate way is to carry out periodic, odd expansions for f1 and f2 .
Rectangular domains. For example, in 1-D, we assume f1 is only defined on domain
D = [0, 0.5]. We choose f = f1 (x) on [0, 0.5] and f = −f1 (1 − x) on [0.5, 1] and take a
multiscale Gaussian wavepacket transform on f . Then the beams near the boundary only
incorporate the information at one side, and we can recover the function f1 only using the
beams with centers contained in [0, 0.5] by the summation method discussed in Section 2.4.2.
It is also consistent with the analytic method for wave equation with constant velocity and
zero boundary condition. The initialization in two dimensional rectangular domain can be
accomplished accordingly.

24

Circular domains. There are several differences between rectangular domain and circular
domain. One difference is the boundary for circular domain is smooth, so we do not have
to drop the beams caused by the nonsmoothness as in the rectangular domain. However,
circular domain also introduces the additional difficulty for initialization, if we note that
the original Fast Gaussian Wavepacket transform is only designed for Cartesian domain. To
overcome this difficulty, we use the idea of odd periodic expansion again. More specifically,
assume the computational domain D is a disk centered at O = (0.5, 0.5) with radius r = 0.25
and the initial condition is compactly supported in D. We need to extend the initial condition
to the square [0, 1] × [0, 1] so that the Gaussian Wavepacket transform can be applied. As
shown in Figure 2.1, for each point P ∈ [0, 1] × [0, 1]\D, find point Q on OP such that
|OQ||OP | = r2 , which is also called the dual point of P with respect to ∂D, and define the
function value at P as the opposite value at Q. We then can apply the transformation and
recover the initial function with beams in D by the summation method discussed in Section
2.4.2.

2.4.2

Beam summation

Another problem raised by Gaussian Beam method is the beam summation. Note that
a beam is not only defined on the ray x(t), but also affects the region around the ray,
which is an advantage over the Geometric Optics from the caustic’s point of view. But
when we are considering a bounded domain, the disadvantage of a beam is how to make it
satisfies the boundary condition, especially for beam near the boundary. A ’tail’ of the beam
appears outside the domain and causes the artificial ’tail effect’. Figure 2.2(a) illustrates
this phenomenon in 1-D.
Rectangular domains. Motivated by the method of images, the remedy here is to imagine
25

1
P
Q
O

0

1

Figure 2.1: Odd periodic expansion for a circular domain.(For interpretation of the references
to color in this and all other figures, the reader is referred to the electronic version of this
dissertation.)
the boundary as a mirror and the tail of the beam is reflected back from the boundary.
Intuitively speaking, the beam can ’see’ it self through the boundary, but with opposite sign
and the field inside the domain is the superposition of these two kinds of beams, see 2.2(b).
We call it ’image principle’, and it’s similar to the concept of incident beams and reflected
beams in Section 2.2, which is also aimed for satisfying the zero boundary condition. Like
the method of images, the beam summation in this way is exact for wave equation with a
constant velocity. Due to the high frequency, the effective area of a beam is small compared to
the whole domain, and the velocity is almost constant in such a small area, so the algorithm
still works even for variable velocity.
Circular domains. When we sum up all the beams, if the effective area of any beam gets

26

0

0
(a)

1
(b)

Figure 2.2: (a)Tail exceeds the boundary; (b)Image principle
out of the computational domain, we need to reflect it back along the circular boundary based
on the same idea as in the case of rectangular domains, since the zero boundary condition
on the circle has to be satisfied. The way to do the reflection is similar to the procedure in
the initialization, but is done reversely. Namely, we find the area inside the computational
domain D, denoted by S ∗ , that is dual to the effective area outside D, denoted by S. For
any point in S ∗ , the beam value should be modified by adding an opposite value from its
dual point in S.

2.4.3

Identifying significant beams

One of the main concern in the algorithm is the number of beams to be launched, since it
dominates the computational time. The number of beams not only depends on the frequencies of the initial condition, but also the dimension of the domain. In particular, beams
increase rapidly in high dimension, so how to choose the significant beams becomes an im27

portant issue in numerical simulation. The most optimal situation should be using the least
needed beams while still achieving the required accuracy. One natural way to choose beam
is to drop the beams with amplitudes below a given threshold. However, as the following
example will show, the amplitudes of some beams could have exponential decay or growth,
which implies, the numerical solution will degrade if those exponentially growing beams are
dropped.
Assume that in 1-D, V (x) = 1+bx, the transport equation (2.13) and the bicharacteristic
equation for p can be simplified as:

dA b p
−
A=0
dt
2 |p|
dp
= b|p|
dt

with the initial condition:

A(0) = A0
p(0) = p0

which yields an analytic solution for A if p0 > 0:

A(t) = A0 e1/2bt

Hence, unless b = 0, the amplitude can be either exponentially decaying or growing. In
particular, if the amplitude is exponentially growing, dropping this ray would result in accu-

28

racy loss. On the other hand, if we keep the ray with exponentially decaying amplitude, the
computational resource is wasted. Although it is a simple example with a linear velocity in
1-D, it illustrates the potential problem in choosing beams based on amplitude only. Hence,
we need a new way to measure the contribution of each beam.
For a beam Φ(x, t) given in the form of (2.60) without subscript, consider the following
quantity:
CΦ |A(x, t)|2
E(Φ(x, t)) = √
det(ℑ(M (t)))
which combines the amplitude and the beam width; here CΦ is some scaling constant related
to the weight of the current beam. It can be shown that E(Φ(x, t)) does not change too
much in the propagation in the sense of the following theorem.
Theorem 2.4.1. Assume that V (x) is C 2 (D). There exist positive constants C1 and C2 ,
only depending on V (x), such that E(Φ(x, t)) satisfies

E(Φ(x, t))
C1 ≤
≤ C2
E(Φ(x, 0))
for 0 ≤ t ≤ T .
Proof. By using the definition of G(x, p), the transport equation can be rewritten as:

d ln A
= −(2G)−1 (V 2 (x(t))trace(M (t)) − Gx · Gp − GT M (t)Gp )
p
dt
1
= − trace(Gpp M (t) + Gxp ) + G−1 Gx · Gp
2
1
d ln V (x(t))
= − trace(Gpp M (t) + Gxp ) +
2
dt

29

(2.66)

The last step was obtained by:

Vx
p
d ln V (x(t))
=
· x = Vx ·
˙
= G−1 Gx · Gp
dt
V
|p|
Thus, we have:
|A(0)|2 |V (x(t))|2
|A(t)|2 =
exp
|V (x(0))|2

(
−

∫ t
0

)
trace(Gpp ℜ(M (s)) + Gxp )ds

(2.67)

The differential equation for det(ℑ(M (t))) is:
(
)
d ln det(ℑ(M (t)))
dℑ(M (t))
= trace
(ℑ(M (t)))−1
dt
dt
(
= −trace Gxp + ℜ(M (t))Gpp
(
)
+ ℑ(M (t)) GT + Gpp ℜ(M (t)) (ℑ(M (t)))−1
xp
(
)

(2.68)
)

= −2trace Gxp + ℜ(M (t))Gpp

(2.69)

The last equality was given by the trace invariance under similarity transformation. Hence:
(
det(ℑ(M (t))) = det(ℑ(M (0)))exp

−2

∫ t
0

)
trace(Gpp ℜ(M (s)) + Gxp )ds

(2.70)

Combining (2.67) and (2.70), the proof is done.
Based on Theorem 4.1, E(Φ(x, t)) is bounded during the propagation of a beam, even at
reflections, since both |A(t)| and ℑ(M (t)) are unchanged at a reflection point. In particular,
E(Φl,i,k (x, 0)) = C(d)cl,i,k |ξl,i |d/2 , where C(d) only depends on the dimension d and N .

30

We thus use E(Φ(x, t)) as the criteria to choose the significant beams, which is better than
using amplitude only.

2.4.4

Overall algorithm

After all those considerations, here is a sketch for the overall algorithm for IBVP (2.1) by
the multiscale Gaussian beam method:
• Decompose the initial condition into a summation of beams by Multiscale Gaussian
Wavepacket Transform (2.42) and determine the coefficients c+ and c− by equal,i,k
l,i,k
tions (2.64) and (2.65).
• Choose a small number ϵ, and for each beam Φl,i,k (x, 0) = φl,i,k (x), evaluate E(Φl,i,k (x, 0)).
If E(Φl,i,k (x, 0)) > ϵ, propagate the beam by solving the equation system (2.48) with
initial condition (2.52) and reflection condition (2.56), otherwise the beam is dropped.
• At the time T , sum up all the beams by image principle given in Section 4.2.
We now analyze the computational cost of this algorithm, which lies in three parts. First
is the initialization by multiscale Gaussian wavepacket transform, which has a complexity
O(N d log N ). The second part is the propagation of all the Gaussian beams. The computation for tracing a single beam over a finite time period can be accomplished in O(1)
steps, so the total cost at this part is proportional to the total number of beams. For most
of the applications, like point sources, plane waves, and curvilinear wavefronts, the number
of beams is supposed to be small at given accuracy ϵ. The final part is the summation
step. As the support of each Gaussian beam is of size O(N 1/2 ) in each dimension, each
beam at time T covers about O(N d/2 ) points. Overall, the computational complexity is

31

O(N d log N + C · N d/2 ), where C denotes the number of beams being traced. It is much
more efficient compared to the O(N d+1 ) cost of standard finite difference or finite element
methods.

2.5

Convergence results

The convergence result of Gaussian beam methods for pure initial value problems of wave
equations was discussed in many papers; see [38, 8]. In terms of wave equations in bounded
convex domains, the convergence results for single-scale Gaussian beam methods based on the
FBI transform was provided by [12]. Inspired by the convergence results in [8] of multiscale
Gaussian beam methods for pure initial value problems of wave equations, we prove the
convergence result following the idea in [8].
The main convergence result is given here:
Theorem 2.5.1. Assume that u(x, t) is the exact solution of the wave equation (2.1) and
uGB (x, t) is the solution based on Gaussian Beam method. T is a finite time. The initial
condition satisfies:

f1 (x) =
f2 (x) =

∑

al,i,k φl,i,k (x)

l,i,k
∑

(2.71)

bl,i,k φl,i,k (x)

(2.72)

l,i,k

32

Let ξmin = minl,i,k {|ξl,i |: |ξl,i | associated with φl,i,k } ≥ 4. We have

sup ∥u(·, t) − uGB (·, t)∥ 1
+ sup ∥∂t u(·, t) − ∂t uGB (·, t)∥ 2
H (D)
L (D)
t∈[0,T ]
t∈[0,T ]
≼

1
(∥f1 ∥ 1
+ ∥f2 ∥ 2
)
H (D)
L (D)
1/2
ξmin

(2.73)

(2.74)

Here ’≼’ denotes the upper bound up to a constant multiple C, which is independent
of the frequency ξmin . The assumption on the initial condition (2.71) (2.72) is based on
the fact that the error introduced by the Gaussian Wavepacket Transform is negligible [8].
Since we are focusing on the high frequency problem, the lower bound on the frequency is
also reasonable. To prove the theorem, we recall some lemmas. From the construction of
Gaussian Beam, let e(x, t) = u(x, t) − uGB (x, t), then e(x, t) satisfies:


 e − V 2 (x)∆e(x, t) = P u

GB (x, t)
 tt







e|t=0 = 0












(2.75)

et |t=0 = 0
e|x∈∂D = BuGB (x, t)

The two terms P uGB (x, t) and BuGB (x, t) represent the propagation error and the boundary error caused by the approximation of Gaussian Beam method, respectively. From classic
PDE theory for wave equations of IBVP, we have:
Lemma 2.5.2. Assume that u(x, t) is the exact solution of the wave equation (2.1) and
uGB (x, t) is the solution based on Gaussian Beam method. Define e(x, t) = u(x, t) −

33

uGB (x, t), then the following estimates holds

sup ∥e(·, t)∥ 1
+ sup ∥∂t e(·, t)∥ 2
H (D)
L (D)
t∈[0,T ]
t∈[0,T ]
≼

sup ∥P uGB ∥ 2
+ ∥BuGB ∥ 1
L (D)
H ([0,T ]×∂D)
t∈[0,T ]

With the assumption on the initial condition (2.71) (2.72), Bao et al [8] essentially proved
the convergence of the propagation error:
Theorem 2.5.3. supt∈[0,T ] ∥P uGB ∥ 2
+ ∥f2 ∥ 2
)
≼ 1 (∥f1 ∥ 1
1/2
H (D)
L (D)
L (D)
ξmin
Hence, we only need to prove:

1
(∥f1 ∥ 1
+ ∥f2 ∥ 2
)
∥BuGB ∥ 1
≼
H (D)
L (D)
H ([0,T ]×∂D)
1/2
ξmin

(2.76)

Before we proceed the proof, we need two important lemmas. The first one gives the relation
between the coefficients and the norm of the initial condition (2.45) (2.46).
Lemma 2.5.4. Let ξl,i , al,i,k , bl,i,k , c+ , and c− be defined in (2.44), (2.45), (2.46),
l,i,k
l,i,k
(2.64), and (2.65), respectively. Then
∑

2
ξl,i (|c+ |2 + |c− |2 ) ≼
l,i,k
l,i,k

l,i,k

∑
l,i,k

2
(ξl,i |al,i,k |2 + |bl,i,k |2 ) ≼ ∥f1 ∥ 1
+ ∥f2 ∥ 2
H (D)
L (D)
(2.77)

The proof is based on the Fourier transform and the boundness of the velocity V (x); see
[44]. In fact, these three quantities are equivalent under a suitable assumption. The second
lemma is:

34

Lemma 2.5.5. Assume that the phase function associated to any beam Φl,i,k (x, t) in (2.60)
satisfy the condition that the imaginary part of the Hessian matrix Ml,i,k is symmetric and
∑
positive definite. Let dl,i,k be complex numbers such that l,i,k |dl,i,k |2 < ∞, and dl,i,k = 0
for small l. Hl,i,k (x, t) is a differentiable function both in x and t. Assume that there exists
an integer n > 0 such that for each (l, i, k) we have Hl,i,k (x, t) = O(|x − xl,i,k (t)|n ). Then

∥

∑ n
∑
ξ 2 dl,i,k Φl,i,k (x, t) · Hl,i,k (x, t)∥2 2
≼
|dl,i,k |2
l,i
L (D)
l,i,k
l,i,k

(2.78)

The proof fully utilizes the exponential decay of the Gaussian Beam. Details can be
found in [8]. It shows that the interaction between different beams can be controlled. Now
we are ready to prove the main result.
Proof of the main theorem. Let

uGB (x, t) = u+ (x, t) + u− (x, t)
GB
GB
∑
∑
where u+ (x, t) = l,i,k c+ Φ+ (x, t) and u− (x, t) = l,i,k c− Φ− (x, t).
GB
GB
l,i,k l,i,k
l,i,k l,i,k
We discuss u+ (x, t) only since u− (x, t) can be treated in the same way, and for
GB
GB
simplicity, we drop the superscript ‘+’. For each beam Φl,i,k (x, t), during the propagation,
it will bounce back and forth in the domain up to finite times N0 . Without loss of generality,
we assume that each beam reflects N0 times, since we can always add zero beams if it is
p
(x, t) to denote the pth (0 ≤ p ≤ N0 ) reflection of the beam. It is easy to
not. Use Φ
l,i,k
N
p
(x, t) can be both incident beam
see that except Φ0 (x, t) and Φ 0 (x, t), for any p, Φ
l,i,k
l,i,k
l,i,k
and reflected beam and we add ’s’ and ’r’ to denote the difference. Now the boundary value

35

can be rewritten as:
N0 −1
∑ ∑
p,s
p+1,r
BuGB (x, t) =
(cl,i,k Φ
(x, t) + cl,i,k Φ
(x, t))|x∈∂D
l,i,k
l,i,k
l,i,k p=0
N0 −1
p,s
∑ ∑
ı·|ξl,i |τ
(x,t)
p,s
l,i,k
=
(cl,i,k A
(x, t)e
l,i,k
l,i,k p=0
p+1,r
ı·|ξl,i |τ
(x,t)
p+1,r
l,i,k
+ cl,i,k A
(x, t)e
)|x∈∂D
l,i,k

Let x

(2.79)

(2.80)

p,s
p,s
and t
denote the position and time that the pth reflected beam hits the
l,i,k
l,i,k

boundary. Following the construction of reflected beams in (2.22) (2.23), for each (l, i, k) we
have:

p,s p,s p,s
p,s
p,s
p,s
(x
,t
) + O(x − x
,t − t
)
(x, t)|x∈∂D = A
l,i,k l,i,k l,i,k
l,i,k
l,i,k
l,i,k

(2.81)

p+1,r
p+1,r p,s p,s
p,s
p,s
(x, t)|x∈∂D = A
(x
,t
) + O(x − x
,t − t
)
l,i,k
l,i,k
l,i,k l,i,k
l,i,k
l,i,k

(2.82)

p,s p,s p,s
p,s
p,s
= −A
(x
,t
) + O(x − x
,t − t
)
l,i,k l,i,k l,i,k
l,i,k
l,i,k

(2.83)

A
A

and

p,s p,s p,s
p,s
p,s
p,s
p,s
p,s
τ
(x
,∂ τ
(x, t)|x∈∂D = τ
,t
) + (∂x τ
) · (x − x
,t − t
) (2.84)
l,i,k
l,i,k l,i,k l,i,k
l,i,k t l,i,k
l,i,k
l,i,k

 

p,s
p,s
p,s
 ∂xx τl,i,k , ∂xt τl,i,k  x − xl,i,k 
p,s
p,s
·

)·
,t − t
+ (x − x
l,i,k 
l,i,k
p,s  
p,s 
p,s T
(∂xt τ
) , ∂tt τ
t−t
l,i,k
l,i,k
l,i,k
(2.85)


3
p,s
x − xl,i,k 

+ O

p,s 
t−t
l,i,k

(2.86)

36

p+1,r
p+1,r p,s p,s
p+1,r
p+1,r
p,s
p,s
τ
(x, t)|x∈∂D = τ
(x
,t
) + (∂x τ
, ∂t τ
) · (x − x
,t − t
)
l,i,k
l,i,k
l,i,k l,i,k
l,i,k
l,i,k
l,i,k
l,i,k
(2.87)

p,s
p+1,r
p+1,r
 ∂xx τl,i,k , ∂xt τl,i,k  x − xl,i,k 
p,s
p,s

·
+ (x − x
,t − t
)·
l,i,k
l,i,k 
p,s 
p+1,r T
p+1,r  
t−t
(∂xt τ
) , ∂tt τ
l,i,k
l,i,k
l,i,k


 

(2.88)


3
p,s
x − xl,i,k 


+ O
p,s 
t−t
l,i,k

(2.89)

Here x ∈ ∂D and ∂x means the partial derivative with respect to the tangential component
of the boundary. By the continuity of the tangential components of phase function τ up to
the second order, we get:


3
p,s
x − xl,i,k 
p,s
p+1,r

τ
(x, t)|x∈∂D − τ
(x, t)|x∈∂D = O

l,i,k
l,i,k
p,s 
t−t
l,i,k

(2.90)

Now combining (2.81) and (2.90) yields:
p,s
p+1,r
ı·|ξl,i |τ
(x,t)
ı·|ξl,i |τ
(x,t)
p,s
p+1,r
l,i,k
l,i,k
A
(x, t)e
+A
(x, t)e
|x∈∂D
l,i,k
l,i,k
p,s
(x,t)
ı·|ξl,i |τ
p,s
p+1,r
l,i,k
=(A
(x, t) + A
(x, t))e
|x∈∂D
l,i,k
l,i,k
p,s
p+1,r
p,s
(x,t))
(x,t)−τ
(x,t) ı·|ξl,i |(τ
ı·|ξl,i |τ
p+1,r
l,i,k
l,i,k
l,i,k
− 1)|x∈∂D
(e
(x, t)e
+A
l,i,k
p,s
(x,t)
ı·|ξl,i |τ
p,s
p,s
l,i,k
)|
,t − t
· O(x − x
=e
l,i,k x∈∂D
l,i,k
p,s
(x,t)
ı|ξl,i |τ
p,s 3
p,s
p+1,r
l,i,k
) |
,t − t
|ξl,i | · O(x − x
(x, t)e
+A
l,i,k x∈∂D
l,i,k
l,i,k

37

(2.91)
(2.92)
(2.93)
(2.94)
(2.95)

Recall the fact that with the non-grazing hypothesis (2.17) the Hessian of the beam on the
boundary [0, T ] × ∂D still has a symmetric positive definite imaginary part [12], so instead
of the space domain D, we can apply Lemma 5.5 to BuGB (x, t) on [0, T ] × ∂D. First
|ξl,i | 1/2
we multiply each term in (2.79) by (
)
, and then use the inequality (2.78), which
ξmin
implies:
1 ∑
∥BuGB (x, t)∥2 2
≼
|cl,i,k |2
L ([0,T ]×∂D) ξmin
l,i,k

(2.96)

Note that by differentiating BuGB (x, t), we gain an extra coefficient ξl,i for each beam and
the others still keep bounded. More importantly, the property that the beam has Gaussian
decay on the boundary [0, T ] × ∂D is unchanged, so result (2.96) can be adapted to the H 1
norm as follows:
1 ∑ 2
∥BuGB (x, t)∥2 1
≼
ξl,i |cl,i,k |2
H ([0,T ]×∂D) ξmin
l,i,k

(2.97)

Now apply Lemma 5.4, and proof is done.

2.6

Numerical results

In this section we will test our algorithm in 1D and 2D domain. The experiment will be
implemented both in constant and variable velocity. Note that the analytic solution may not
be available particularly with variable velocity. We obtain the ’exact’ solution by fourth order
finite difference method [18] with a dense grid (up to 8192 points on each dimension within
length of 0.5) over the computational domain. The CFL condition is given by CF L = 0.5,
which is small enough to prevent the dispersion of the solution. Also note we choose the cut
off threshold ϵ in Section 2.4 to be 10−3 . All the errors are measured in the discrete relative
38

L2 norm, instead of the norm defined in (2.73).

2.6.1

One-dimensional domain

In 1D, consider the domain to be [0,0.5] and N = 8192. The initial condition is given by:

f1 = 2 sin(ωπx) exp(−40(x − 0.25)2 );

(2.98)

f2 = 0

(2.99)

The velocity changes in different cases. In order to confirm the convergence rate obtained by
the analysis, let the main frequency ω of the initial condition vary. In particular, we consider
ω = 140, ω = 300, ω = 620 and ω = 980. Given the exact solution u(x) and Gaussian Beam
solution uGB on the grid, the relative L2 error is evaluated by the following way:
√
ek =

2.6.1.1

1 N
2
N Σi=1 (ui − uGBi )
√
1 N 2
N Σi=1 ui

(2.100)

Case 1: Linear velocity

Example 1: Our first example is the wave equation with linear velocity. In particular,
V (x) = 1 + 0.5x. Figure 2.3 gives the result for frequencies ω = 620 and ω = 980 at
t = 0.5. The exact solution and Gaussian Beam solution are overlapped in these figures.
We see a good match between the two solutions and the differences are almost negligible.
Figure 2.4 graphs the convergence rate by evaluating (2.100) along the frequency. The rate
is approximately 0.5, which is in agreement with our analysis.
Example 2: The second example has the same velocity and initial profile as Example 1,

39

except that the initial derivative is nonzero, which is defined as:

f2 = ω sin(ωπx);

Figure 2.5 gives the result for frequencies ω = 620 and ω = 980 at t = 0.5. The two solution
are almost identical to each other. Figure 2.6 shows the convergence rate, which agrees with
the analysis as well. This example demonstrates the effectiveness of Gaussian Beam method
for nonzero derivative.

2.6.2

Two-dimensional rectangular domains

Consider the wave equation in the domain [0, 0.5] × [0, 0.5] for different velocities. Reflection
condition for the Hessian is given by eqauation (2.30) when the beam hits the boundary.
The domain is uniformly discretized by N × N = 2048 × 2048. Initial condition for both of
the following examples is given by:

f1 = 2 sin(ωπ(x + y)) exp(−40((x − 0.25)2 + (y − 0.25)2 ))

(2.101)

f2 = 0;

(2.102)

ω will vary in order to verify the convergence rate. In particular, let ω to be 100,120,180,220,
respectively. The computation stops at t = 0.5.

2.6.2.1

Case 1: Linear velocity

Example 3: The velocity is given by V (x, y) = 1 + 0.5(x + y). Figure 2.7(a) shows the exact
solution for ω = 180, Gaussian Beam solution for ω = 180 is given in 2.7(b). Figure 2.7(c)

40

shows the comparison for the slice at x = 0.25. Figure 2.7(d) shows the comparison for the
slice at y = 0.125. Numerical solution shows the algorithm works well in two dimension.
Figure 2.8 demonstrate the relative L2 convergence rate, which agrees with the analysis.
There are 36736 beams propagating in the computational domain at ω = 180, which only
accounts for 0.4% of the total number of beams.

2.6.2.2

Case 2: Sinusoidal velocity

Example 4: Consider V (x, y) = 1 + 0.25 sin(2π(x + y)). Part of the result is shown in
Figure 2.9. Figure 2.9(a) shows the exact solution for ω = 220. Gaussian Beam solution
for ω = 220 is shown in figure 2.9(b).Figure 2.9(c), 2.9(d) provide the comparison of the
solutions at different slices. Figure 2.10 shows the convergence result based on the four
frequencies. It matches the analysis very well.

2.6.3

Two-dimensional circular domains

We proceed by considering a circular domain which is centered at (0.5, 0.5) with radius 0.25.
We discretize the domain uniformly by the grid N × N = 2048 × 2048. The computation
stops at t=0.5.The purpose of this experiment is test our scheme in more general domain.
Choose the initial condition to be:

f1 = 2 sin(ωπ(x − 0.5)) sin((ω + 20)π(y − 0.5)) exp(−40((x − 0.5)2 + (y − 0.5)2 )) (2.103)
f2 = 0

(2.104)

In order to check the convergence rate, we choose ω to be 60,100,120,180. Also note that in
order to get an accurate solution to compare, special treatment is needed at the boundary
41

for finite difference solution. Here we choose the Cartesian embedded boundary method to
accommodate the curved boundary. Details can be found in [29].

2.6.3.1

Case 1: Sinusoidal velocity

Example 5: Consider V (x, y) = 1+0.25 sin(2πx) cos(2πy). The result is shown in Figure 2.11.
Figure 2.11(a), 2.11(b) show the exact solution and Gaussian Beam solution for ω = 180
respectively. Figure 2.11(c) and 2.11(d) provide the comparison for the slice at x=0.45
and y=0.45. Figure 2.12 plots the convergence result based on the four frequencies. The
convergence rate matches the analysis result.

2.7

Conclusion

In this chapter, we discuss how to construct the global Gaussian Beam solution in a bounded
domain with multiple reflections. The initialization process is completed by Multiscale Gaussian wavepacket transform, which is quite fast and accurate in the representation of initial
profile. The theoretical analysis is provided to show the convergence rate of the method. We
also give the numerical examples to validate our algorithm, which yields a good match with
the compared solution. The convergence rate is justified in those numerical examples. We
also propose the criteria to choose significant beams after initialization, which turns out to
be fairly effective in practice.
Future work includes to extend the algorithm to discontinuous velocity, which requires to
consider the refraction at the discontinuity. The difficulty is to control the amount of beams
during the propagation since it increases every time at the discontinuity. Another future
work is to consider three dimensional bounded domain problem, especially high frequency

42

Maxwell’s equation.

43

2
Gaussian Beam solution
Exact solution

1.5
1
0.5
0
−0.5
−1
−1.5
−2
0

0.1

0.2

0.3

0.4

0.5

(a)
2
Gaussian Beam solution
Exact solution

1.5
1
0.5
0
−0.5
−1
−1.5
−2
0

0.1

0.2

0.3

0.4

0.5

(b)
Figure 2.3: Example 1. Linear model in 1D. (a) Comparison between the beam solution and
exact solution at ω = 620; (b) Comparison between the beam solution and exact solution at
ω = 980; (c) Windowed comparison between the beam solution and exact solution at ω = 620;
(d) Windowed comparison between the beam solution and exact solution at ω = 980.

44

Figure 2.3: (cont’d)

Gaussian Beam solution
Exact solution

1.5
1
0.5
0
−0.5
−1
−1.5
0.24

(c)

0.25

0.26

0.27

Gaussian Beam solution
Exact solution

1.5
1
0.5
0
−0.5
−1
−1.5

(d)

0.215

0.22

0.225

0.23

45

0.235

0.24

0.245

−1

10

Error vs Frequency
Line with slope −1/2

−2

10

−3

10

2.2

10

2.4

10

2.6

10

2.8

10

Figure 2.4: Relative L2 errors for linear velocity in Example 1

46

2
Gaussian Beam solution
Exact solution

1.5
1
0.5
0
−0.5
−1
−1.5
−2
0

0.1

0.2

0.3

0.4

0.5

(a)
2
Gaussian Beam solution
Exact solution

1.5
1
0.5
0
−0.5
−1
−1.5
−2
0

0.1

0.2

0.3

0.4

0.5

(b)
Figure 2.5: Example 2. Linear model in 1D. (a) Comparison between the beam solution and
exact solution at ω = 620; (b) Comparison between the beam solution and exact solution at
ω = 980; (c) Windowed comparison between the beam solution and exact solution at ω = 620;
(d) Windowed comparison between the beam solution and exact solution at ω = 980.

47

Figure 2.5: (cont’d)

Gaussian Beam solution
Exact solution

1.5
1
0.5
0
−0.5
−1
−1.5
0.205

(c)

0.21

0.215

0.22

0.225

0.23

0.235

Gaussian Beam solution
Exact solution

1.5
1
0.5
0
−0.5
−1
−1.5

(d)

0.24

0.25

48

0.26

0.27

−1

10

Error vs Frequency
Line with slope −1/2

−2

10

−3

10

2.2

10

2.4

10

2.6

10

2.8

10

Figure 2.6: Relative L2 errors for linear velocity in Example 2

49

0

1.5
1

0.1

0.5
0.2
0
0.3
−0.5
0.4

0.5
0

−1

0.1

0.2

0.3

0.4

0.5

−1.5

(a)
0
1

0.1

0.5
0.2
0
0.3
−0.5
0.4

0.5
0

−1

0.1

0.2

0.3

0.4

0.5

(b)
Figure 2.7: Example 3. Linear model in 2D. (a) The exact solution at ω = 180; (b) The beam
solution at ω = 180; (c) Comparison for the slice at x = 0.25; (d) Comparison for the slice
at y = 0.125; (e)Windowed comparison for the slice at x = 0.25; (f )Windowed comparison
for the slice at y = 0.125.

50

Figure 2.7: (cont’d)

1.5
Gaussian Beam solution
Exact solution

1
0.5
0
−0.5
−1

(c)

−1.5
0

0.1

0.2

0.3

0.4

0.5

1.5
Gaussian Beam solution
Exact solution

1
0.5
0
−0.5
−1

(d)

−1.5
0

0.1

0.2

0.3

51

0.4

0.5

Figure 2.7: (cont’d)

Gaussian Beam solution
Exact solution

1
0.5
0
−0.5
−1

0.22

(e)

0.24

0.26

0.28

0.3

Gaussian Beam solution
Exact solution

1
0.5
0
−0.5
−1

(f)

0.22

0.24

0.26

52

0.28

0.3

Error vs Frequency
Line with slope −1/2

−0.6

10

−0.7

10

−0.8

10

−0.9

10

1

10

2

10

3

10

Figure 2.8: Relative L2 errors for linear velocity in Example 3

53

0
3
0.1

2
1

0.2

0
0.3

−1
−2

0.4

−3
0.5
0

0.1

0.2

0.3

0.4

0.5

(a)
0

3
2

0.1

1
0.2
0
0.3

−1
−2

0.4

−3
0.5
0

0.1

0.2

0.3

0.4

0.5

(b)
Figure 2.9: Example 4. Sinusoidal model in 2D. (a) The exact solution at ω = 220; (b)
The beam solution at ω = 220; (c) Comparison for the slice at x = 0.25; (d) Comparison
for the slice at y = 0.125;(e) Windowed comparison for the slice at x = 0.25; (f ) Windowed
comparison for the slice at y = 0.125.

54

Figure 2.9: (cont’d)

1.5
Gaussian Beam solution
Exact solution

1
0.5
0
−0.5
−1

(c)

−1.5
0

0.1

0.2

0.3

0.4

0.5

3
Gaussian Beam solution
Exact solution

2
1
0
−1
−2

(d)

−3
0

0.1

0.2

0.3

55

0.4

0.5

Figure 2.9: (cont’d)

Gaussian Beam solution
Exact solution

1
0.5
0
−0.5
−1

0.22

(e)

0.24

0.26

0.28

0.3

Gaussian Beam solution
Exact solution

2
1
0
−1
−2

(f)

0.22

0.24

0.26

56

0.28

0.3

−0.3

10

Error vs Frequency
Line with slope −1/2

−0.4

10

−0.5

10

−0.6

10

−0.7

10

1

10

2

10

3

10

Figure 2.10: Relative L2 errors for sinusoidal velocity in Example 4

57

0.8

0.3

0.6
0.4

0.4

0.2
0.5

0
−0.2

0.6

−0.4
−0.6

0.7

−0.8
0.3

0.4

0.5

0.6

0.7

(a)

0.8

0.3

0.6
0.4

0.4

0.2
0.5

0
−0.2

0.6

−0.4
−0.6

0.7

−0.8
0.3

0.4

0.5

0.6

0.7

(b)
Figure 2.11: Example 5. Sinusoidal model in 2D. (a) The exact solution at ω = 180; (b)
The beam solution at ω = 180; (c) Comparison for the slice at x = 0.45; (d) Comparison
for the slice at y = 0.45; (e)Windowed comparison for the slice at x = 0.45; (f ) Windowed
comparison for the slice at y = 0.45.

58

Figure 2.11: (cont’d)

Gaussian Beam solution
Exact solution

0.15
0.1
0.05
0
−0.05
−0.1
−0.15
0.3

(c)

0.4

0.5

0.6

0.7

Gaussian Beam solution
Exact solution

0.4
0.3
0.2
0.1
0
−0.1
−0.2
−0.3
−0.4

(d)

0.3

0.4

0.5

59

0.6

0.7

Figure 2.11: (cont’d)

0.2

Gaussian Beam solution
Exact solution

0.15
0.1
0.05
0
−0.05
−0.1
−0.15
0.42

(e)

0.44

0.46

0.48

0.5

0.52

Gaussian Beam solution
Exact solution

0.4
0.3
0.2
0.1
0
−0.1
−0.2
−0.3
−0.4

(f)

0.48

0.5

0.52

60

0.54

0.56

0.58

Error vs Frequency
Line with slope −1/2
−0.3

10

−0.4

10

1

10

2

10

3

10

Figure 2.12: Relative L2 errors for sinusoidal velocity in Example 5

61

Chapter 3
Radar cross section reduction of a
cavity: TM case

3.1

Introduction

We consider a two dimensional rectangular cavity embedded in an infinite ground plane,
which is a simplified model for many important applications, like jet inlet and nozzle of an
aircraft. Such structure is generally considered as an identifier of an object on the radar
system due to its dominance of the overall RCS. Therefore, reducing the radar cross section(RCS) scattered from a cavity has practical importance. One of the effective approach
to reduce RCS is through coating material[28]. Assume a thin, multilayered RAM is coated
horizontally at the bottom of the cavity for the reduction of RCS, as shown in Figure 3.1.
The problem is how to determine the permittivity of the material so that RCS is minimized.
In this chapter, the problem is addressed by the gradient-based sequential quadratic
programming (SQP) method. SQP is a one of the most effective optimization methods
particularly for nonlinear constraint like partial differential equations (PDE)[39]. The basic
idea of this method is to generate steps by solving a sequence of quadratic subproblems.
The gradient involved in the quadratic subproblems is provided accurately by the continuous
adjoint state method. The accuracy plays an important role in the design especially when
the backscatter RCS is highly sensitive to the change of permittivities.
62

In the following section, the optimal design of a cavity with multilayered RAM is formulated as a minimization problem, with the constraint being formulated as a Helmholtz
equation with artificial boundary condition. The objective function and design variables are
also defined in this section. The continuous adjoint method is applied in Section 3.3 for the
gradient of the backscatter RCS. Section 3.4 introduces the SQP optimization technique to
the RCS reduction. Section 3.5 gives a brief introduction to the direct solver. Numerical
experiments are provided in Section 3.6 to show the application of the method. Section 3.7
concludes the whole chapter.

3.2
3.2.1

Problem formulation and objective function
Formulation

Consider a rectangular cavity Ω = [0, a] × [−b, 0] embedded in a ground plane illuminated
by a plane wave, as illustrated in Figure 3.1. The problem is in 2-D by assuming the cavity
and the materials are invariant in the z direction. Above the ground plane is empty space
with dielectric permittivity ε0 . The surface of the ground plane Γc and the boundary S of
the cavity are assumed to be perfect conductors. The cavity is filled with inhomogeneous
material in layered structure with permittivities εj , j = 1, 2, · · · , n, where n is the number
of layers and presumably small from the practical point of view. Assume the magnetic
permeability µ0 is the same everywhere.
In 2-D, there are two fundamental polarizations in the wave propagation, namely, TM
(transverse magnetic) polarization and TE (transverse electric) polarization. In this chapter,
we consider TM polarization only, which implies the magnetic field is transverse to the xy
plane and the electrical field is simply expressed as (0, 0, u(x, y)). The goal is to determine
63

the permittivities of the coating material within a small thickness h such that the backscatter
RCS is minimized.

Ei
y

Es
Er
Γ

x

Ω
εn

ε2
ε1

Γc

S

Figure 3.1: Geometry of the problem

Assume the incident field is given by a plane wave: ui = ei(αx−βy) , where α =
√
k0 cos θ, β = k0 sin θ, θ is the incident angle with respect to the x−axis. Let k = ω εµ0 be
√
the wavenumber in the space. In particular, k0 = ω ε0 µ0 is the wavenumber in the upper
half space R2+ . Denote u(x, y) and us (x, y) as the total electrical field and the scattered
field respectively. The total field u(x, y) can be expressed as :

u(x, y) = us (x, y) + ui − ei(αx+βy)

(3.1)

where −ei(αx+βy) is identified as the reflected field ur caused by the infinite ground plane.

64

In addition, u(x, y) satisfies the following equation:

∆u + k 2 u = 0 in Ω ∪ R2+
(3.2)
u = 0 on S ∪ Γc
together with the radiation condition satisfied by us :
(
lim
ρ→∞
√
where ρ =

∂us
− ik0 us
∂ρ

)
=0

(3.3)

x2 + y 2 .

In general, analyzing Eq. (3.2) directly is difficult both theoretically and numerically
due to the unboundedness of the domain. The general approach is to introduce some kind
of artificial boundaries so that the domain is reduced to be bounded. Here we introduce
a transparent boundary condition, which essentially is Drichlet-to-Neumann (DtN) map,
to truncate the unbounded domain. The idea of transparent boundary condition has been
applied to many wave scattering problems; for instance, see [19]. Following the formulation
in [2], we derive the boundary condition based on the Fourier Transform.
The exterior scattered field us (x, y) satisfies the equation:

2
∆us + k0 us = 0 in R2+
us = ψ(x, 0) on Γ

(3.4)

us = 0 on Γc

together with the radiation condition (3.3). Here we assume ψ(x, 0) is the known scattered
field at the aperture Γ.

65

Take the Fourier Transform of us with respect to variable x and denote it by us (ξ, y).
Eq. (3.4) becomes:
∂ 2 us
2
+ (k0 − ξ 2 )us = 0
∂y 2

(3.5)

us (ξ, 0) = ψ(ξ, 0) on Γ ∪ Γc

(3.6)

where ψ(x, 0) denotes the zero extension of ψ(x, 0) on y = 0. This equation admits a simple
solution:
i
us (ξ, y) = ψ(ξ, 0)e

√
2
k0 −ξ 2 y

(3.7)

We drop another linearly independent solution by the radiation condition (3.3). Taking an
inverse Fourier Transform of (3.7) yields the following solution for us :
√
∫
2
i k0 −ξ 2 y iξx
1
us (x, y) =
ψ(ξ, 0)e
e dξ
2π R

(3.8)

Differentiating both sides with respect to y at y = 0 gives us:
∫ √
∂us
i
2
k0 − ξ 2 ψ(ξ, 0)eiξx dξ
=
y=0+ 2π R
∂n

(3.9)

We use y = 0+ since the derivative is derived from y > 0. From the expression (3.1) for the
total field u and the fact ui + ur = 0 on y = 0, it yields:

∂u
= I(u) + g
∂n y=0+

66

(3.10)

where g = 2iβeiαx and I(·) is the boundary operator defined by:
∫ √
i
2
I(u) =
k0 − ξ 2 u(x, 0)eiξx dξ
2π R

(3.11)

By the continuity condition on the boundary Γ, boundary condition (3.10) is also satisfied
by u in Ω. We therefore reduce the unbounded domain problem for u to the following:
∆u + k 2 u = 0 in Ω
(3.12)

u = 0 on S
∂u
= I(u) + g(x) on Γ
∂n

Problem (3.12) is uniquely solvable for any incident plane wave. Proof can be found in [2].

3.2.2

Objective function

The RCS reduction problem may be formulated as minimizing the following cost function :
∫
2
4 k0
σ :=
sin θ
ue−ik0 x cos θ dx
k0 2
Γ

(3.13)

which is the backscatter RCS for TM polarization[26]. In the case of the RCS reduction over
a set of incident angles {θ1 , θ2 , · · · , θm }, the cost function may be changed to the following:

σt :=

m
∑

wl σ(θl )

(3.14)

l=1
where wl > 0, l = 1, 2, · · · , m are the weights for different angles. Furthermore, depending
on which angle is important, these weights may be chosen accordingly.

67

Since u depends on ε from the model (3.12), the optimization problem can be rewritten
as:
min σ(ε)
ε∈Λ

(3.15)

where the admissible set Λ for electrical permittivity ε satisfies the following two conditions:
• Layered medium: ε(y) in Λ is constant for each layer, either lossy or lossless. In other
′′
words, ε(y) = εj at the jth layer, where εj = ε′ + iεj , j = 1, 2, · · · , n, n is the number
j
of layers of the coating material. In particular, permittivity with nonzero imaginary
part corresponds to lossy material.
′′
′
′′
• Boundness: We assume there exist εmin = ε′
min + iε min , εmax = εmax + iε max
′′
′′
′′
such that ε′
< ε′ < ε′
max and ε min < εj < ε max for j = 1, 2, · · · , n. Namely,
min
j
the permittivity ε(y) is bounded between εmin and εmax .
With the well defined objective function and design variables, in order to employ a gradient
descent method, it is essential to calculate the gradient of the objective function efficiently,
which we will introduce in the next section.

3.3

Gradient by the adjoint state method

To use a gradient based optimization method for the backscatter RCS (3.13) with respect to
the permittivity, we need to evaluate the gradient along with the PDE constraint (3.12). Here
we adopt the adjoint state method to give the gradient for TM polarization. The adjoint state
method has a broad range of applications in nonlinear optimizations; see [20] for example.
Theoretically, the gradient obtained by the adjoint state method is more accurate compared
to the approximated gradient given by the finite difference. Numerically, it is also more
68

efficient since the evaluation only requires solving an extra adjoint problem regardless of the
number of variables. Here we give a formal derivation of the gradient.
Denote the far field coefficient:
∫
k0
Pθ (ε) =
sin θ
ue−ik0 x cos θ dx
2
Γ

(3.16)

Then
σ(ε) =

4
|P (ε)|2
k0 θ

Let δε be a ”small” perturbation to the permittivity ε. We denote the linearization of σ(ε)
with respect to δε by δσ, which is:
)
4
δP Pθ (ε)
δσ = 2Re
k0
(

(3.17)

where δP denotes the linearization of Pθ (ε) and Pθ is the complex conjugate of Pθ . From
(3.16), it is easy to calculate that
∫
k0
sin θ
δP =
δue−ik0 x cos θ dx
2
Γ

(3.18)

where δu solves the linearized problem
∆δu + k 2 δu = −(ω 2 µ0 δε)u in Ω
δu = 0 on S
∂δu
= I(δu) on Γ
∂n

69

(3.19)

Let u∗ solves the adjoint state equation:
∆u∗ + k 2 u∗ = 0 in Ω
u∗ = 0 on S

(3.20)

∂u∗
= I ∗ (u∗ ) + eik0 x cos θ Pθ (ε) on Γ
∂n
where I ∗ (u∗ ) is the adjoint operator of I(u), given by
∫ √
∗ (u∗ ) = −i
2
I
k0 − ξ 2 u∗ (x, 0)eiξx dξ
2π R
Given equations (3.19) and (3.20), an integration by parts yields that
∫
Ω

ωµ0 δε · uu∗ dx =

∫
Γ

δue−ik0 x cos θ ds · Pθ (ε)

(3.21)

Comparing with (3.17), we then have
(

∫

δσ = 4Re sin θ
Ω

)
2 µ δε · uu∗ dx
ω 0

(3.22)

Up to a constant mutiple, the gradient of the cost function σ(ε) is the function g(ε) such
that

∫ 0
δσ(ε) = Re

−b

δεg(ε)dy

(3.23)

Comparing (3.22) with (3.23), we arrive at the formula for the gradient:

g(ε) = 4 sin θω 2 µ0

70

∫ a
0

uu∗ dx

(3.24)

Re(ε1 )
Im(ε1 )
Re(ε2 )
Im(ε2 )
ASM -8.8402E-6 -3.2758E-5 -3.8514E-4 -1.5E-3
FDM -8.8332E-6 -3.2730E-5 -3.8514E-4 -1.5E-3
Ratio
1.0008
1.0008
1.0000
1.0000

CPU Time
2.24s
3.89s

Table 3.1: A comparison for the gradients obtained from the adjoint state method (ASM)
and the finite difference method (FDM).

where u∗ is the solution of (3.20).
For the cost function (3.14), since it is a linear combination of (3.13), the gradient can
be easily deduced from (3.24). It can be seen that the adjoint problem shares the same
structure as the original scattering problem, hence the same solver may be used for the
adjoint problem. This is a huge advantage given the level of difficulties for solving the
scattering problem from a cavity. Using the finite-dimensional counterpart of the gradient,
Table 3.1 gives a comparison between the gradient obtained by (3.24) and the numerical
gradient obtained by finite difference. Assume the dimension of the cavity is [0, 1] × [−0.3, 0]
with two layers of coating material placed at the bottom. The permittivities of the two
layers are (1 + i, 2 + 2i) with the thickness of each layer being 5.8480E − 4. Consider the
normal incidence of a plane wave with k0 = 16π on the cavity. The ratio between these
two gradients are close to one, but the computational time by (3.24) is much smaller than
the time spent by the finite difference. The time gap will enlarge as the number of design
variables increase. The efficiency and accuracy of the gradient obtained by the adjoint state
method provide a huge advantage for the optimization method.

71

3.4

Optimization method: Sequential Quadratic Programming

We adopt the Sequential Quadratic Programming(SQP) method to optimize the problem
(3.15) largely because of its outstanding performance on the nonlinear constraint problem.
For a sufficiently smooth problem, this method provides a superlinear convergence near
the region of optimal point under some suitable assumptions[39]. The principle idea is
to generate a sequence of quadratic programming(QP) subproblems that are based on the
approximation of Lagrangian function in a small region. For each QP subproblem, the
solution is used to form a search direction for a line search. In general the Hessian matrix
used in the QP subproblem is not obtained analytically, but instead is updated by some
quasi-Newton methods. Here we give three steps on how to apply SQP to our problem. One
resorts to [39] for a complete discussion on this method.
(
)
( )j=2n
Step 1: Let d = Re(ε1 ), Im(ε1 ), Re(ε2 ), · · · , Im(εn ) T = di j=1 . Define a sequence
j=4n
of functions {hj (d)}j=1 to be:
h4(j−1)+1 (d) = d2j−1 − ε′
max ,
h4(j−1)+2 (d) = ε′
min − d2j ,
h4(j−1)+3 (d) = d2j−1 − ε′′ max ,
h4(j−1)+4 (d) = ε′′ min − d2j ,
for j = 1, 2, · · · , n.

It is easily seen that hj (ε) ≤ 0 means ε is bounded below and above by the two constants

72

εmin and εmax . Hence, Problem (3.15) can be rewritten in the following form:

min σ(d)
d∈R2n

(3.25)

subject to: hj (d) ≤ 0,
for j = 1, 2, · · · , 4n.

Define the Lagrangian function for (3.25):

L(d, λ) = σ(d) +

4n
∑

λj hj (d),

(3.26)

j=1
where λj is the Lagrangian multiplier.
Step 2: It is difficult to solve (3.26) directly due to the nonlinearity of σ(d), so the idea
is to replace (3.26) by its quadratic approximation, which results in a QP subproblem, so
that many QP algorithms can be used. At the m−th iteration, the QP subproblem of SQP
when d = d(m) is given by:

1
min g(d(m) )T x + xT Hm x
2
x∈R2n

(3.27)

subject to: hj (d(m) + x) ≤ 0,
for j = 1, 2, · · · , 4n,

where g(d(m) ) is the gradient of σ(d) at d(m) and Hm is the Hessian of the Lagrangian
function (3.26), which is in fact the Hessian of σ(d) by the linearity of the constraint. It is
usually approximated by some quasi-Newton methods to keep the Hessian positive definite.
Here we choose the wildly used Broyden-Fletcher-Goldfarb-Shanno(BFGS) formula, which
73

is given by:
Hm+1 = Hm −

T
T
H m xT xm H m y m y m
m
+
,
T
xT H m x m
y m xm
m

(3.28)

where

xm = d(m+1) − d(m) ,

(3.29)

ym = g(d(m+1) ) − g(d(m) ).

(3.30)

The Hessian keeps being positive definite whenever the initial Hessian matrix H0 is positive
T
T
definite and ym xm > 0. Issues can happen when ym xm < 0. In that case we modify ym
T
component-wise such that ym xm > 0 is satisfied.
Step 3: The solution of the subproblem (3.27) yields a new step:

d(m+1) = d(m) + αx

(3.31)

The step length α ∈ (0, 1) can be determined by an appropriate defined merit function. Here
we use the following:
φ(ε) = σ(ε) +

4n
∑

µj max(0, hj (ε))

(3.32)

j
for sufficiently large µj , which gives penalty to solutions outside the admissible set.
We are now ready to state the minimization algorithm:
1. Given the convergence tolorence ϵ > 0, the initial value d(0) and H0 = I, where I is the
identity matrix, let m = 0.
2. Repeat, until ∥x∥ < ϵ
Evaluate g(d(m) ), Hm ,

74

Solve (3.27) to obtain x,
Set d(m+1) = d(m) + αx,
where α is determined by the merit function (3.32).
end
Note that the solution obtained by this way is by no means the global minimum. What
we expect is only local minimum. Depending on different applications, we either expect
there is a good initial value available or try different initial values in the admissible set and
choose the most ”optimal” one after running the optimization.

3.5

Direct solver for the forward scattering problem

Numerical solution of the electromagnetic scattering from a cavity has been long known as
a difficult problem [35], especially when the size of the cavity is large or the wavenumber
is large. Here we choose the algorithm developed in [9] as the direct solver, which is a
finite difference method. The idea of this algorithm is to use the discrete sine transform in
the horizontal direction and a Gaussian elimination in the vertical direction. The resulting
system becomes an M×M (M is number of mesh points at the opening of the cavity) linear
system, which is much smaller than the original one. It then can be solved by any linear
solver with an appropriate preconditioner. Since the algorithm highly utilizes the special
structure of the cavity problem, it is extremely fast and efficient, which fits exactly into the
requirement for a computation engine of the optimization method. The evaluation of the
boundary integral operator I in (3.11) also needs special care. We refer the readers to [48]
for details.

75

3.6

Numerical Experiments

In this section, we apply our algorithm to the RAM design through three numerical examples.
The goal is to show that RCS can be reduced with appropriate coating material inside the
cavity. Throughout, we assume the dimension of the cavity is 1m in width and 0.3m in
depth. The cavity is uniformly partitioned by a 512 × 512 mesh. The thickness of the
coating material for each layer is 0.6mm, up to 8 layers. For each layer, the initial value of
the relative permittivity is set to be 5 + 5i, and assume it is continuously changing between
1 + 0i and 100 + 100i. The artificial material is used only for illustration purpose. All the
computation is carried out on a laptop with Intel dual core 2.1Ghz and memory 4GB. For
physical interest, the RCS is expressed in term of decibel(dB), which is

RCS: = 10 log10 σ dB.

Example 1: Consider the scattering from the cavity by the incidence of a plane wave
with wavenumber k0 = 16π. Figure 3.2(a) gives the result for the optimization at normal
incidence. RCS are reduced by more than 100dB at θ = π/2. However, the reduction
immediately disappears after θ = π/2 and starts to oscillate. The overall effect of the coating
material is very limited. Figure 3.2(b) shows the result for RCS optimization at θ = 3π/5.
The deep well in the graph shows a large amount of reduction at θ = 3π/5. Unlike the
optimization at θ = π/2, the RCS reduction is not only confined around θ = 3π/5. It also
extends to the other angular sector with an average reduction around 20dB, which can be
considered as a large improvement. Similar effect has been shown in Figure 3.2(c), which
gives the RCS when the optimization is conducted at θ = 5π/6. As expected, large amount
of reduction appears at θ = 5π/6. Meanwhile, the RCS at the other angular sector also gets
76

reduced. The resulted relative permittivities are listed in Table 3.2 with the corresponding
number of iterations and CPU time. All the numerical experiments are finished in 121s and
the convergence is achieved in less than 35 steps, which shows the algorithm is very effective
and efficient.
The different performances between θ = π/2 and the other angles are caused by the
geometry of the cavity and the way the coating material works. Mathematical analysis
for such differences is very difficult, so we only give a possible explanation based on the
ray tracing theory[34]. Because of the thin thickness, the abosrbing effect of the coating
material is very limited. The main reason that RCS gets reduced is because the propagating
directions of the reflected rays are changed by the coating material. While the change may
not succeed for all incident angles, it usually will result in some unchanged angles or even
enhanced angles. Therefore, the optimized coating for normal angle may not work for some
oblique angles and vice versa.
Example 2: Consider the incidence of a plane wave with wavenumber k0 = 32π on the
cavity. Similar to Example 1, the optimization are performed at θ = π/2, 3π/5, 5π/6.
Figure 3.3(a) gives the result for the optimization at normal incidence. RCS are reduced
by more than 25dB at θ = π/2. However, the result is also similar to Example 1, since
the reduction appears only in a small interval around θ = π/2 and tend to disppear beyond
θ = 5π/9. Figure 3.3(b) shows the result for RCS optimization at θ = 3π/5. The deep well
in the graph shows a large amount of reduction at θ = 3π/5. The RCS reduction is not
only limited around θ = 3π/5. It also extends to the other angular sector with an average
reduction around 10dB. Figure 3.3(c) gives the optimized RCS when the optimization is
conducted at θ = 5π/6. As expected, a large amount of reduction appears at θ = 5π/6.
Meanwhile, the RCS at the other angular sector also gets reduced. Table 3.2 contains all the
77

60
Without coating
With coating

40

RCS (dB)

20
0
−20
−40
−60
−80

100

120

140
θ (Degree)

160

180

(a)

50

Without coating
With coating

RCS (dB)

0

−50

−100

−150

100

120

140
θ (Degree)

160

180

(b)
Figure 3.2: RCS reduction for k0 = 16π. (a) Optimized at θ = π/2,(b) Optimized at
θ = 3π/5.(c) Optimized at θ = 5π/6.
78

Figure 3.2: (cont’d)

60
Without coating
With coating

40
20

RCS (dB)

0
−20
−40
−60
−80
−100
−120
(c)

100

120

140
θ (Degree)

79

160

180

resulted relative permittivities with the corresponding number of iterations and CPU time.
It again shows the algorithm is very efficient and effective.
Example 3: In this example, we try to define the objective function as a combination of
RCS at θ = 3π/5 and 5π/6. The formula is given in (3.14) with equal weights for the two
angles. Results are shown in Figure 3.4 for both k = 16π and k = 32π and the corresponding
CPU time are given in Table 3.2. It takes much longer time for the method to converge
compared to the previous two examples. The combination of these two angles provides a
smoother reduction instead of a sharp decrease at some particular angle. Consequently, it
may have more practical applications.

k0
θ
1st(Re)
1st(Im)
2nd(Re)
2nd(Im)
3rd(Re)
3rd(Im)
4th(Re)
4th(Im)
5th(Re)
5th(Im)
6th(Re)
6th(Im)
7th(Re)
7th(Im)
8th(Re)
8th(Im)
Iter.
CPU

Example 1
16π
16π
π/2
3π/5
1.786 11.206
0
5.387
4.142 11.788
0
6.016
8.067 10.085
0
0
13.563 18.199
0
0
20.638 24.967
0
0
29.312 46.344
0
0
39.624 65.949
0
0
51.624 46.428
0
25.258
8
31
26.2s 112.5s

16π
5π/6
11.558
6.326
13.151
9.729
12.831
4.895
21.972
7.066
27.592
5.880
42.611
5.786
47.110
3.951
58.850
4.013
24
120.3s

Example 2
32π
32π
π/2
3π/5
1
10.588
0
6.431
1
9.329
0
10.143
1
4.457
0
5.722
1
7.485
0
8.257
1
5.863
0
7.032
1
12.887
0
5.986
1
9.514
0
1.904
37.788 13.725
0
1.800
107
15
299.4s 59.3s

32π
5π/6
10.542
6.283
9.162
9.563
4.169
4.461
7.133
6.088
5.653
3.844
13.136
1.784
10.602
0
15.910
0
13
59.7s

Example 3
16π
32π
3π/5, 5π/6 3π/5, 5π/6
100
13.985
30.978
40.326
100
59.028
67.296
100
100
1
0
0
100
1
0
0
90.459
1
0
0
1
1
0
0
1
1
0
0
1
37.949
0
0
246
177
2306s
1816s

Table 3.2: Relative permittivies after the optimization, number of iterations and CPU time
for all the three examples.

80

3.7

Conclusion

In this chapter, the design of a multilayered RAM for reducing RCS of a cavity is formulated
as a minimization problem. The descent direction for the cost function is evaluated through
the adjoint state method. Subsequently, SQP is integrated with the gradient to obtain the
optimal absorbing materials for the cavity. The algorithm is implemented with a fast and
accurate direct solver for the scattering problem. Numerical results show that RCS is reduced
significantly with RAM coated at the bottom of the cavity. It is also observed that different
objective function generally leads to different RAM design and RCS reduction.

81

30
Without coating
With coating

20

RCS (dB)

10
0
−10
−20
−30
−40

100

120

140
θ (Degree)

160

180

(a)

50
Without coating
With coating

RCS (dB)

0

−50

−100

−150

100

120

140
θ (Degree)

160

180

(b)
Figure 3.3: RCS reduction for k0 = 32π. (a) Optimized at θ = π/2; (b) Optimized at
θ = 3π/5; (c) Optimized at θ = 5π/6.
82

Figure 3.3: (cont’d)

40
Without coating
With coating

20
0
RCS (dB)

−20
−40
−60
−80
−100
−120

(c)

100

120

140
θ (Degree)

83

160

180

60
Without coating
With coating

40

RCS (dB)

20
0
−20
−40
−60
100

120

140
θ (Degree)

160

180

(a)

30
Without coating
With coating

20

RCS (dB)

10
0
−10
−20
−30
−40
100

120

140
θ (Degree)

160

180

(b)
Figure 3.4: RCS reduction by the combination of θ = 3π/5 and θ = 5π/6. (a) Optimized for
k0 = 16π, (b) Optimized for k0 = 32π.
84

Chapter 4
Radar cross section reduction of a
cavity: TE case

4.1

Introduction

In this chapter we continue our study on the RCS reduction for a cavity embedded in the
ground plane. Assume the layered coating material with different permittivities is placed
at the bottom of the cavity. We are trying to determine the permittivity of the coating
material such that the RCS is minimized. The optimization problem is again addressed by
the gradient-based sequential quadratic programming(SQP) method. However, instead of
TM(transverse magnetic) polarization, we focus on the RCS reduction for TE(transverse
electric) polarization, which implies the electric field is transverse to the xy plane and the
magnetic field is pointing along the z direction, as shown in Figure 4.1.
Two difficulties arise when one is dealing with TE polarization:
1. The equation governing the TE polarization is a generalized Helmholtz equation, where
the design variable, or the permittivity of the coating material is included in a differential operator. Special care must be taken when using the adjoint state method to
find the gradient of the objective function.
2. The fast algorithm for solving the generalized Helmholtz equation in [9] does not apply
85

to TE polarization with nondifferentiable permittivity.
Both of these issues will be addressed in this chapter. In particular, we extend the
algorithm in [9] to efficiently solve the forward problem with permittivity being piecewise
constant. The outline of this chapter is as follows: in Section 4.2, the optimal design of a
cavity with multilayered radar absorbing material is formulated as a minimization problem,
with the constraint being formulated as a generalized Helmholtz equation with artificial
boundary condition. The objective function and design variables are also defined in this
section. The continuous adjoint method is applied in Section 4.3 for the gradient of the
backscatter RCS. In Section 4.4, SQP optimization technique is applied to the RCS reduction.
Section 4.5 introduces the extended algorithm for solving the generalized Helmholtz equation.
Numerical experiments are provided in Section 4.6 to show the efficiency of the method.
Section 4.7 concludes the whole discussion.

4.2
4.2.1

Problem formulation and objective function
Formulation

Consider a rectangular cavity Ω = [0, a]×[−b, 0] embedded in the ground plane illuminated by
a time harmonic plane wave, as illustrated in Figure 4.1. The problem is in 2-D by assuming
the cavity and the material are invariant in the z direction. Above the ground plane, the
medium is assumed to be homogeneous with dielectric permittivity ε0 . The surface of the
ground plane Γc and the boundary S of the cavity are assumed to be perfect conductors.
The bottom of the cavity is coated by a thin, inhomogeneous material in layered structure
with permittivities given by εj , j = 1, 2, · · · , L, where L is the number of layers. Above the

86

coating material, the cavity is filled with the same material as the upper half space. For
√
convenience, denote ε as the permittivity function in terms of variable y and k = ω εµ0
as the wavenumber in the whole space. In this paper, we consider TE(transverse electric)
polarization and assume the magnetic permeability µ0 is the same everywhere. The goal
is to determine the permittivities of the coating material within a thickness h, which is
presumably small, such that the backscatter RCS is minimized.

ui
y

z

(0,0)

x

S

us

ur
Γ

Ω
εL

Γc

ε2
ε1

Figure 4.1: Geometry of the problem
Assume the incident field is given by a plane wave: ui = ei(αx−βy) , where α =
√
k0 cos θ, β = k0 sin θ, θ is the incident angle with respect to the x-axis and k0 = ω ε0 µ0
is the wavenumber in the upper half space R2+ . Denote u(x, y) and us (x, y) as the total
magnetic field and the scattered field respectively. Following the formulation in [26], the
total field u(x, y) satisfies:

u(x, y) = us (x, y) + ui + ei(αx+βy)

(4.1)

where ei(αx+βy) is identified as the reflected field ur . We have the following generalized
87

Helmholtz equation for u(x, y), which can be easily derived from the time harmonic Maxwell’s
equation:

∇·(

1
∇u) + u = 0 in Ω ∪ R2+
k2

∂u
= 0 on Γc ∪ S
∂n

(4.2)
(4.3)

We now employ Fourier Transform method to derive a transparent boundary condition, by
which Eq (4.2) is reduced to a bounded domain problem. Ideas of transparent boundary
condition have been used in many wave propagation problems. One resorts to [19] and
reference therein for the applications.
The exterior scattered field us (x, y) satisfies the equation:

2
∆us + k0 us = 0 in R2+
∂us
= ϕ(x) on Γ
∂n
∂us
= 0 on Γc
∂n

(4.4)

together with the radiation condition:
(
lim
ρ→∞

∂us
− ik0 us
∂ρ

)
=0

(4.5)

√
Here ρ = x2 + y 2 and we assume ϕ(x) is known. Let ϕ(x) denote the zero extension of
ϕ(x) on y = 0. Taking a Fourier Transform for us with respect to x and denote it by us (ξ, y),

88

the equation (4.4) becomes:
∂ 2 us
2
+ (k0 − ξ 2 )us = 0
∂y 2
∂ us
= ϕ(ξ) on y = 0
∂n

(4.6)
(4.7)

Using the outgoing radiation condition (4.5), we easily obtain the solution for (4.6):
√
2
i k0 −ξ 2 y
ϕ(ξ)
us (ξ, y) = √
e
,
2 − ξ2
i k0

if ξ ̸= k0

(4.8)

Taking an inverse Fourier Transform with respect to ξ and letting y = 0 yields:
∫
1
s (x, 0) = 1
√
u
ϕ(ξ)eiξx dξ
2πi R k 2 − ξ 2
0

(4.9)

s
Observing that ϕ(x) = ∂u = ∂u on Γ, we therefore obtain the transparent boundary
∂n
∂n
condition for u:

(

∂u
u(x, 0) = T
∂n

)
+ g,

on Γ

(4.10)

where
(

∂u
T
∂n

)

∫
1
1
∂u iξx
√
=
e dξ,
2πi R k 2 − ξ 2 ∂n
0

g = 2eiαx ,

and ∂u denotes the zero extension of ∂u on y = 0.
∂n
∂n
Equipped with the transparent boundary condition, Eq. (4.2) is reduced to a bounded

89

domain problem:

∇·(

1
∇u) + u = 0,
k2

in Ω

∂u
= 0,
∂n

on S
( )
∂u
u(x, 0) = T
+ g,
∂n

(4.11)
on Γ

Eq. (4.11) is uniquely solvable and the proof can be found in [2]. The normal derivative of
u at the aperture Γ can be found through the numerical solution of (4.11) and then used to
evaluate the backscatter RCS. The formula of RCS is given in the next section.

4.2.2

Objective function

The RCS reduction problem may be formulated as minimizing the following cost function :
∫
∂u −ik x cos θ 2
4 1
0
e
dx
σ :=
k0 2 Γ ∂n

(4.12)

which is the backscatter RCS for TE polarization[26]. In the case of the RCS reduction over
a set of incident angles {θ1 , θ2 , · · · , θr }, the cost function may be changed to the following:

σtotal :=

r
∑

wl σ(θl )

(4.13)

l=1
where wl > 0, l = 1, 2, · · · , r are the weights for different angles. Furthermore, depending
on which angle is important, these weights may be chosen accordingly.
Since u depends on ε from the model (4.11), the optimization problem can be rewritten

90

as:
min σ(ε)
ε∈Λ

(4.14)

where the admissible set Λ for electrical permittivity ε satisfies the two conditions given
in Section 3.2.2. We intend to use a gradient based optimization algorithm for the RCS
reduction. It is therefore very important to evaluate the gradient of the objective function
accurately and efficiently.

4.3

Gradient by the adjoint state method

To determine the descent direction of the backscatter RCS (4.12) with respect to the permittivity, we need to evaluate the gradient along with the PDE constraint (4.11). It can be
approximated by the numerical derivative. However, given the sensitive dependence of RCS
on the permittivity and the amount of computational cost to evaluate (4.11), it is generally
undesirable to simply use the numerical approximation. Our approach here is to derive the
gradient for TE polarization based on the adjoint state method. This method is well-known
for its effectiveness in evaluating the gradient when the state variable comes from a solution
of a forward problem. Lots of nonlinear optimization problems have been studied based on
this method; for example, see [20]. Here we provide a former derivation for the gradient
based on the adjoint state method.
For convenience, let us first define the far field coefficient:
∫
1
∂u −ik x cos θ
0
Pθ (ε) =
e
dx,
2 Γ ∂n

91

(4.15)

Comparing with (4.12), we have

σ(ε) =

4
|P (ε)|2
k0 θ

Let δε be a ”small” perturbation to the permittivity ε. We denote the linearization of σ(ε)
with respect to δε by δσ, which is:
(

)
4
δσ = 2Re
δP Pθ (ε)
k0

(4.16)

where δP denotes the linearization of Pθ (ε) and Pθ is the complex conjugate of Pθ . From
(4.15), it is easy to calculate that
∫
∂δu −ik x cos θ
1
0
δP =
e
dx
2 Γ ∂n

(4.17)

where δu solves the linearized problem
(
∇·

)
( 2
)
ω µ0 δε
1
∇δu + δu = ∇ ·
∇u
in Ω
k2
k4

∂δu
= 0 on S
∂n
(
)
∂δu
δu = T
on Γ
∂n

(4.18)

Let u∗ solve the adjoint state equation:
(
∇·

)
1
∗ + u∗ = 0 in Ω
∇u
k2

∂u∗
= 0 on S
∂n
( ∗)
1
∗ = T ∗ ∂u
u
+ k 2 eik0 x cos θ Pθ (ε) on Γ
∂n
2
92

(4.19)

where k 2 is the complex conjugate of k 2 and T ∗ (·) is the adjoint operator of T (·), given by
( ∗)
∫
∂u iξx
1
1
∗ ∂u
T
=−
e dξ
√
∂n
2πi R
2 − ξ 2 ∂n
k0
Given equations (4.18) and (4.19), Green’s Theorem implies that
∫
Ω

u∗ ∇ ·

( 2
)
∫
∫
ω µ0 δε
1 ∗ ∂δu
1 ∂u∗
∇u dx =
u
dx −
δudx
∂n
k4
∂Ω k 2
∂Ω k 2 ∂n

(4.20)

Applying the integration by parts yields:
( 2
)
( 2
)
( 2
)
∫
∫
∗ ∇· ω µ0 δε ∇u dx = −
∗ · ω µ0 δε ∇u dx+
∗ ω µ0 δε ∂u dx (4.21)
u
∇u
u
k4
k4
k 4 ∂n
Ω
Ω
∂Ω

∫

Recalling the homogeneous Neumann boundary condition, the second term on the right
hand side of (4.21) is zero since the permittivity is unchange near the aperture Γ. Hence,
combining (4.20) and (4.21) together with the property of adjoint operator, we obtain:
∫
−

Ω

∇u∗ ·

( 2
)
∫
ω µ0 δε
1
∂δu −ik x cos θ
0
∇u dx = Pθ (ε)
e
dx
4
2
∂n
k
Γ

(4.22)

Comparing with (4.16), we then have
) )
( 2
∫
4
∗ · ω µ0 δε ∇u dx
δσ = −2Re
∇u
k0 Ω
k4
(

(4.23)

Up to a constant multiple, the gradient of the cost function σ(ε) is the function g(ε) such
that

∫ 0
δσ(ε) = Re

−b
93

δεg(ε)dy

(4.24)

Comparing (4.23) with (4.24), we arrive at the formula for the gradient:
(
) ( )
∫ a
8 2
1
∗ · ∇u dx
g(ε) = − ω µ0
∇u
k0
0 k4

(4.25)

In other words, to find the gradient at a particular point ε ∈ Λ, we need to solve the forward
scattering problem (4.11) as well as the adjoint state problem (4.19).
Remark 4.3.1. The gradient given in (4.25) is in fact L2 gradient, which is different from
the usual vector form. In the numerical simulation, we need the finite dimensional counterpart of the gradient, which points to the maximal increase of RCS in the admissible set.
Since the change δε is a constant for each layer, the finite discretization for (4.24) implies
the gradient can be found through the following vector:
∫ −b+h
∫ −b+2h
∫ −b+Lh
(
g(ε)dy,
g(ε)dy, · · · ,
g(ε)dy, )
−b
−b+h
−b+(L−1)h

(4.26)

where h is the thickness of each layer and L is the number of layers. In practice, we can
also split each element in (4.26) into real and imaginary part, which gives the corresponding
change for the real and imaginary part of δε.
For the cost function (4.13), since it is a linear combination of (4.12), the gradient can
be easily deduced from (4.25). It can be seen that the adjoint problem shares the same
structure as the original scattering problem, hence the same solver may be used for the
adjoint problem. This is a huge advantage given the level of difficulties for solving the
scattering problem from a cavity.

94

4.4

Optimization algorithm

Similar to Chapter 3, we propose an optimization algorithm based on SQP. For completeness,
we restate the algorithm here. Define the Lagrangian function:

L(ε, λ) = σ(ε) +

4L
∑

λj hj (ε),

(4.27)

j=1
where λj is the Lagrangian multiplier, hj (ε) gives the constraint on the permittivity:
h4(j−1)+1 = −Re(εj ) + ε′
min
h4(j−1)+2 = Re(εj ) − ε′
max
h4(j−1)+3 = −Im(εj ) + ε′′ min

(4.28)

h4(j−1)+4 = Im(εj ) − ε′′ max
j = 1, 2, · · · , L
Here we consider the real part and imaginary part of the permittivity separately, so that the
boundness of ε can be simply expressed as

hj (ε) ≤ 0

(4.29)

The m-th subproblem of SQP during the iteration is given by:

1
min σ(εm ) + g(εm )T (d) + dT Hm d
2
d∈CL

(4.30)

with the linear constraint (4.28) satisfied by εm + d, where g(εm ) is the gradient of σ(ε) at
εm , Hm is the Hessian of the Lagrangian function (4.27), which is approximated by BFGS
95

method.
Solution of the subproblem (4.30) provides an update for ε along the direction d

εm+1 = εm + αd

(4.31)

The step length α ∈ (0, 1) can be determined by an appropriate defined merit function. We
use the following merit function

φ(ε) = σ(ε) +

4L
∑

µj max(0, hj (ε))

(4.32)

j
for sufficiently large µj > 0.
Now we are ready to state the minimization algorithm:
1. Given initial permittivity ε0 , tolerance ϵ, Let m = 0, α = 1.
2. Repeat, until α > ϵ
Evaluate g(εm ), Hm , h(εm ).
Solve (4.30) to obtain d.
Set εm+1 ← εm + αd.
While φ(εm+1 ) > φ(εm ), do
α ← α/2.
end
Remark 4.4.1. Since the problem (4.14) we consider is generally not convex with respect to
ε, the solution from the algorithm above may not converge to the global one. In practice, we
either expect a good initial guess is provided, or try different initial values in the parameter
space and choose the most optimal one. More complicated strategy is possible, but we will

96

not pursue it here.
Since our algorithm is programmed in MATLAB, lots of existing functions can be called.
Basically, we only have to provide the objective function, the constraint, the gradient and
a direct solver for the forward scattering problem (4.11). We already discussed how to find
the gradient in Section 4.3. In the next section, we present an accurate and efficient direct
solver for the optimization algorithm.

4.5

Direct solver for the forward scattering problem

Similar to Chapter 3, we extend the finite difference method developed in [9] to be the direct
solver. The idea to solve the resulting linear system is to use the discrete cosine transform
in the horizontal direction and a Gaussian elimination process in the vertical direction. The
resulting system becomes an M×M (M is number of mesh points at the opening of the cavity)
linear system, which is much smaller than the original one. Since the original algorithm in
[9] requires the permittivity ε of the medium inside the cavity to be differentiable in the
vertical direction for TE polarization, which is not satisfied for layered medium. Therefore,
we modify the algorithm to make it work for our case.
Assume Ω is composed by L layers with each layer denoted by Ωl , l = 1, · · · , L and
the associated permittivity is given by εl . Similar to [9], the interior and the aperture
M +1,N +1
of the cavity is uniformly partitioned by the mesh points {xm , yn }m=1,n=1 with hx =
xm+1 − xm and hy = yn+1 − yn . We further assume there exist mesh points on each
interface of two layers with n = n1 , n2 , · · · , nL−1 . The algorithm will not lose the generality
by this assumption since the extension to nonuniform mesh in y direction is straightforward.
Define um,n as the corresponding numerical solution at (xm , yn ). Since the material is
97

homogeneous in each layer, the second order finite difference scheme in the l-th layer is
simply given by the following:

um+1,n − 2um,n + um−1,n um,n+1 − 2um,n + um,n−1
+
+ ω 2 µ0 εl um,n = 0,
h2
h2
x
y
where (xm , yn ) ∈ Ωl (4.33)

Equations at different layers are coupled together via the continuity condition:

u+ = u − ,

(4.34)

1 ∂u+
1 ∂u−
= −
,
ε+ ∂y
ε ∂y

(4.35)

where ’+’ and ’−’ denote the upper and lower side of an interface, respectively.
Condition (4.34) implies it is unnecessary to distinguish the value on the two sides of an
interface. Condition (4.35) can be approximated by the finite difference with points on the
the same vertical line. Here we use the second order approximation in order to be consistent
with (4.33). For n = n1 , n2 , · · · , nL−1 , the approximation is given by:
1 1/2um,n+2 − 2um,n+1 + 3/2um,n
1 −1/2um,n−2 + 2um,n−1 − 3/2um,n
=
εl+1
hy
εl
hy
(4.36)
For simplicity, the first order approximation is applied for the homogenous Neumann
boundary condition at ∂Ω. It is also used to discretize the normal derivative and the improper

98

integral at the aperture Γ, which results in the following discretization[26]:

um,N +1 =

N
∑
n=1

(u
Gm,n

m,N +1 − um,N
hy

)
+ gm ,

m = 1, 2, · · · , M

(4.37)

where

 [
i
 1 + 2i ln(0.1638k h )]h

0 x x
π
2
Gm,n =

 i (2)
 H (k |x − x |)h

m x
0 n
2 0

if n = m
(4.38)
if n ̸= m

gm = 2eiαxm .

(4.39)

For compactness, the scheme can be rewritten in the matrix form:

∗
(Ax ⊗ IN + IM ⊗ A∗ + IM ⊗ D∗ )u + (IM ⊗ aN +1 )uN +1 = 0
y
GuN + (hy IM − G)uN +1 = hy g

(4.40)
(4.41)

where ⊗ denotes the tensor product(Kronecker product), IM is the M × M identity matrix,
∗
IN is similar to the N ×N identity matrix but the entry at (nl , nl ) is zero for l = 1, 2, · · · , L−
1, and G = (Gm,n )M
m=1.n=1 .

99

For the other matrices, their forms are:













−1 1
0
 g1 


 




 



 1 −2 1
0
 g 


 
 2 



1 . 
1 


 . 
 .. 
.. .. ..
Ax =
g=
 aN +1 = 2  . 

.
.
.
. 



h2 
hy  
x

 




 



0
gM −1 
1 −2 1 


 




 


1 −1
1
gM
M ×M
(N +1)×1
(M +1)×1
(4.42)

A∗ =
y

1
h2
y


−1 1


 1 −2
1



..
..
..
.
.
.




1
−2




1 εl+1 −2 εl+1

εl
2 εl


























1



(
εl+1 )
 · · · n -th row
3 1+
−2 1

εl
l
2
2



−1
−2 1



.. .. ..

.
.
.



1 −2 1 


1 −2
N ×N
(4.43)

100

and

D∗ = diag(ε1 , · · · , ε1 , 0, ε2 , · · · , εL−1 , 0, εL , · · · , εL ),

(4.44)

u = (u1,1 , · · · , u1,N , u2,1 , · · · , u2,N , · · · , uM,1 , · · · , uM,N )T ,
un = (u1,n , · · · , uM,n )T
where D∗ is a diagonal matrix with zero appearing at the nl -th row, l = 1, 2, · · · , L.
We then can apply the discrete cosine transform to Ax and the rest of the algorithm
is completely analogy to the method introduced in [9], which we will omit and refer the
readers to [9] for details. Since the algorithm highly utilizes the special structure of the
cavity problem, it is extremely fast and efficient, which exactly fits the requirement for a
computation engine of the optimization method.

4.6

Numerical Experiments

In this section, we apply our algorithm to the RAM design through three numerical examples.
Throughout, we assume the dimension of the cavity is 1m in width and 0.3m in depth. In
other words, Ω = [0, 1] × [−0.3, 0]. We partition Ω uniformly by a 512 × 512 mesh. The
thickness of the coating material for each layer is 2.3mm, with four layers in total. Assume
the relative permittivity is continuously changing between 1 + 0i and 100 + 100i. Note that
these settings are only used for illustration purpose. Our algorithm applies to any finite
number of layers and any finite bounds of the permittivity. All the numerical experiments
are carried out on a laptop with 2.1Ghz Intel dual core and 4GB memory. For physical

101

interest, the RCS is expressed in terms of decibel(dB), which is

RCS: = 10 log10 σ dB.

Initially, we assume the layered RAM has a uniform permittivity, namely, 5 + 5i. As we
mentioned before, the resulted permittivity of the material after the optimization may highly
depend on the initial value. We certainly have tried other initial values, but no significant
improvement was observed in terms of RCS reduction. Therefore, all the examples shown
here are based on the same initial value.
Example 1: Consider the scattering from a cavity by the incidence of a plane wave with
wavenumber k0 = 16π. Figure 4.2(a) gives the result for the optimization at θ = 3π/5. RCS
is reduced by more than 100dB at θ = 3π/5. The reduction is not only confined around
θ = 3π/5. It also extends to the other angular sector with an average reduction around
10dB, which can be considered as a large improvement. Similar effect has been shown in
Figure 4.2(b), which gives the RCS when the optimization is conducted at θ = 5π/6. The
deep well in the graph shows a large amount of reduction at θ = 5π/6. Meanwhile, the
RCS at the other angular sector also gets reduced. Figure 4.2(c) shows the RCS reduction
when the optimization is conducted by the combination of θ = 3π/5 and θ = 5π/6. The
weights for these two angles are set to be equal. Compared to the optimization at a single
angle, optimization with multiple angles avoids the sharp reduction at a particular angle
and provides a smoother reduction for all the angles. The resulted relative permittivities are
listed in Table 4.1 with the corresponding number of iterations and CPU time. Generally,
the time and steps for optimization with two angles are much longer than the time spent for
a single angle.

102

Example 2: Consider the incidence of a plane wave with wavenumber k0 = 32π on the
cavity. Similar to Example 1, the optimization are performed at 3π/5, 5π/6 and both.
Figure 4.3(a) gives the result for the optimization at θ = 3π/5. RCS are reduced by more
than 100dB at θ = 3π/5. At the other angular sectors, RCS also gets reduced, although
the performance is not as good as θ = 3π/5. Figure 4.3(b) shows the result for RCS
optimization at θ = 5π/6. The deep well in the graph shows a large amount of reduction at
θ = 5π/6. Meanwhile, the reduction is also occurred at the other angles. Figure 4.3(c) gives
the optimized RCS when the optimization is conducted by the combination of θ = 3π/5 and
3π/5. The sharp reduction at a particular angle disappears but the overall reduction gets
improved. Table 4.1 contains all the resulted relative permittivities with the corresponding
number of iterations and CPU time. It again shows more time is needed for combination of
multiple angles.
Example 3: In the last example, we test our algorithm by showing a RCS reduction
over a range of frequencies. In particular, we consider the reduction between 16π and 32π,
which implies the frequency ranges from 2.4Ghz to 4.8Ghz. Take the objective function as
a combination of RCS at k = 16π and k = 32π. The formula is similar to (4.13) with equal
weights for the two wavenumber. Results are shown in Figure 4.4 for both θ = π/2 and
θ = 3π/5 and the corresponding CPU time are given in Table 4.1. From Figure 4.4(a), we
see some reduction occurs at the two end. However, the overall effect is merely some kind
of horizontal shift of RCS. Generally it implies a reduction at one frequency accompanies an
enhancement at another frequency for normal incidence. The situation is changed for oblique
incidence, as shown in 4.4(b), which gives an overall reduction for all the frequencies.

103

k
θ
1st(Re)
1st(Im)
2nd(Re)
2nd(Im)
3rd(Re)
3rd(Im)
4th(Re)
4th(Im)
Iter.
CPU

Example 1
16π
16π
16π
3π/5
5π/6 3π/5, 5π/6
4.669 5.374
1.000
6.737 5.066
3.293
6.786 7.005
1.000
6.341 4.837
2.662
11.107 10.176
29.271
6.767 4.931
0.000
16.573 13.985
3.937
9.981 6.613
10.149
15
16
100
48.9s
52.7s
1113s

Example 2
Example 3
32π
32π
32π
16π-32π
16π-32π
3π/5 5π/6 3π/5, 5π/6
π/2
3π/5
4.809 5.004
1.000
1.000
93.098
4.925 4.936
0.000
0.009
62.367
3.659 4.793
54.316
1.000
20.880
4.370 4.775
25.262
0.000
18.116
2.892 4.684
1.000
1.000
6.769
2.732 4.384
5.669
0.000
0.042
4.788 5.141
16.503
8.587
8.338
1.017 4.035
0.506
0.177
0.000
10
11
129
24
81
50.5s 49.1s
2404s
33.8s
1420s

Table 4.1: Relative permittivities after the optimization, number of iterations and CPU time
for all the three examples.

4.7

Conclusion

In this paper, the design of a multilayered RAM for reducing RCS of a cavity is formulated
as a minimization problem. The descent direction for the cost function is evaluated through
the adjoint state method. Subsequently, SQP is integrated with the gradient to obtain the
optimal absorbing materials for the cavity. The algorithm is implemented with a fast and
accurate direct solver for the scattering problem. Numerical results show that RCS is reduced
significantly with RAM coated at the bottom of the cavity. It is also observed that different
weights in the cost function generally lead to different RAM design and RCS reduction.
Future work includes considering RCS reduction for three dimensional cavity.

104

50
Without coating
With coating

RCS (dB)

0

−50

−100

100

120

140
θ (Degree)

160

180

(a)

40
Without coating
With coating

20
0
RCS (dB)

−20
−40
−60
−80
−100
−120

100

120

140
θ (Degree)

160

180

(b)
Figure 4.2: RCS reduction for k0 = 16π. (a)Optimized at θ = 3π/5.(b)Optimized at θ =
5π/6.(c)Optimized with the combination of θ = 3π/5 and θ = 5π/6.
105

Figure 4.2: (cont’d)

30
Without coating
With coating

20

RCS (dB)

10
0
−10
−20
−30
(c)

100

120

140
θ (Degree)

106

160

180

50
Without coating
With coating

RCS (dB)

0

−50

−100

−150

100

120

140
θ (Degree)

160

180

(a)

50
Without coating
With coating

RCS (dB)

0

−50

−100

−150

100

120

140
θ (Degree)

160

180

(b)
Figure 4.3: RCS reduction for k0 = 32π. (a)Optimized at θ = 3π/5.(b)Optimized at θ =
5π/6.(c)Optimized with the combination of θ = 3π/5 and θ = 5π/6.
107

Figure 4.3: (cont’d)

30
Without coating
With coating

20

RCS (dB)

10
0
−10
−20
−30
−40
(c)

100

120

140
θ (Degree)

108

160

180

40
Without coating
With coating

30

RCS (dB)

20
10
0
−10
−20
−30
−40

2.5

3

3.5
4
Frequency (GHz)

4.5

20
Without coating
With coating

10

RCS (dB)

0
−10
−20
−30
−40

2.5

3

3.5
4
Frequency (GHz)

Figure 4.4: RCS reduction for k between 16π and 32π.
(b)Optimized at θ = 3π/5.

109

4.5

(a)Optimized at θ = π/2.

Chapter 5
Optimal shape design of a cavity for
radar cross section reduction

5.1

Introduction

Our goal in this chapter is to propose an optimization method to reduce the RCS by a
proper shaping of the cavity. It is observed in [16] that using an appropriate longitudinal
bending, the RCS of a cylindrical cavity was reduced significantly over a set of incident
angles. In this paper, the RCS reduction is formulated as a minimization problem involving
the Helmholtz equation. To reduce the open cavity problem into a bounded domain problem,
an artificial boundary is introduced, which is essentially a Dirichlet to Neumann map. The
minimization problem is solved by an iterative scheme. During each iteration, a forward
scattering problem and an adjoint problem are solved. The solutions are then combined
together to give a descent direction of the backscatter RCS. Since only the boundary data
is needed, both of them are solved by the boundary integral method. The boundary of the
cavity is parameterized by a summation of cosine functions, and the shape of the cavity
is controlled by the amplitude of those functions. The relation between the perturbation
of the amplitudes and the descent direction is obtained by solving a least square problem.
Numerical examples show the algorithm is effective in finding the optimal shape with reduced
RCS at different frequencies.
110

The rest of the chapter is organized as follows. Section 5.2 gives the formulation of the
problem and define the objective function. In Section 5.3, the existence of the minimizer
for the optimization problem is established. Descent direction of the objective function is
given in Section 5.4 via the domain derivative. Section 5.5 gives the parametrization of
the boundary and the optimization algorithm that leads to the optimal shape of the cavity.
Numerical examples are provided in Section 5.6 to illustrate the effectiveness of the algorithm.
Section 5.7 concludes the whole discussion.

5.2

Problem description

Consider a cavity Ω embedded in an infinite ground plane Γc , as shown in Fig. 5.1. The
opening of the cavity is denoted as Γ. The cavity is assumed to be filled with the same
material as the upper half space, which is vacuum. Thus the dielectric coefficient ε and
magnetic permeability µ are the same everywhere. We focus on the two dimensional problem
by assuming that the cavity has no variance in the x3 direction. There are two fundamental
polarizations for the EM wave propagation: TM (Transverse magnetic) polarization and
TE (Transverse electric) polarization. For simplicity, we consider the TM polarization only,
which is the case when the magnetic field is transverse to the invariant direction and the
electrical field E = (0, 0, u(x1 , x2 )). The boundary of the cavity Ω is composed by a piecewise
smooth curve. Denote S1 , S2 and S3 as the left, right and bottom boundary, respectively
and let S = S1 ∪ S2 ∪ S3 . The ground plane Γc and the boundary S are assumed to be
perfect conductor, so the total field u satisfies the following boundary condition:

u = 0 on S ∪ Γc .

111

(5.1)

Assume the incident electrical field is a time harmonic plane wave

ui
Γc

x3
x3

us

θ

ur

Γ

S1

ε, µ

Ω

x1

S2
S3

Figure 5.1: Geometry of the problem

ui = ei(αx1 −βx2 ) ,

where α = k cos θ, β = k sin θ, k 2 = ω 2 εµ is the wavenumber of the empty space with
frequency ω and θ is the incident angle with respect the positive x1 direction. Following the
formulation in [2], the total field u after the illumination takes the form:

u = us + ui + ur ,

(5.2)

where us is the scattered field and ur = −ei(αx1 +βx2 ) is the reflected field caused by the
2
ground plane. Denote R+ as the upper half space. The equation for the total field is:

112

2
∆u + k 2 u = 0 in Ω ∪ R+ ,
(5.3)
u = 0 on S ∪ Γc ,

together with the Sommerfeld radiation condition:
(

∂us
lim
r
− ikus
r→∞
∂r
√

√
where r =

)
= 0,

(5.4)

x2 + x2 .
1
2

Since to analyze the problem (5.3) directly is difficult due to the unboundedness of
the domain, we introduce a transparent boundary condition on Γ to reduce the unbounded
problem into a bounded one. The idea of transparent boundary condition, which is essentially
the Dirichlet to Neumann map, has been studied in many papers [5],[2],[3]. The boundary
condition can be derived either by Green’s function method or by the Fourier transform [9].
Here we give the derivation based on the Fourier transform.
Since ui + ur = 0 on Γ ∪ Γc , the scattered field satisfies the following equation in the
upper half space:

2
∆us + k 2 us = 0 in R+ ,
us = 0 on Γc ,
us = φ on Γ,

along with the radiation condition (5.4), where φ is assumed to be known.
113

(5.5)

Take the Fourier transform of the Helmholtz equation in (5.5) with respect to x1 and
denote the Fourier Transform of us as us

2
∂x us + (k 2 − ξ 2 )us = 0 for x2 > 0.
2

(5.6)

It admits two linearly independent solutions. We choose the one that will satisfy the radiation
condition and represent the ”outgoing” wave
√
i k 2 −ξ 2 x2 s
us = e
u (ξ, 0).

(5.7)

After the inverse Fourier Transform, we obtain the identity

us =

∫

1
2π R

√
i k 2 −ξ 2 x2 s
u (ξ, 0)eiξx dξ.
e

(5.8)

Taking a normal derivative of both sides at x2 = 0 yields the transparent boundary condition
for the scattered field
∂us
|
= I(us ),
∂n x2 =0

(5.9)

where the normal direction n towards the upper half space and the boundary operator is
defined by

∫ √
i
k 2 − ξ 2 f (ξ, 0)eiξx1 dξ.
I(f ) =
2π R

(5.10)

Using the identity (5.2), we can derive the corresponding boundary condition for the total
field u(x1 , x2 ). Since on the ground plane it holds ui + ur = 0, the same boundary operator
I can be applied to the total field u at x2 = 0. The resulting boundary condition for u

114

satisfies:
∂u
= I(u) + g
∂n

on Γ,

(5.11)

where g = −2iβeiαx1 is the summation of the normal derivative of the incident and reflected
fields. Now the equation for the total field u(x1 , x2 ) becomes:

∆u + k 2 u = 0 in Ω,
u = 0 on S,
∂u
= I(u) + g
∂n

(5.12)
on Γ.

For further discussions on the properties of the boundary operator I, we define the space

1/2
H00 (Γ) = {u ∈ H 1/2 (Γ) : ∃˜ ∈ H 1/2 (R) such that u = 0 on Γc and u = u|Γ }
u
˜
˜
In this setting, we call u the zero extension of u to H 1/2 (R).
˜
˜1
The scattering problem (5.12) has an equivalent weak formulation: Find u ∈ H0 (Ω) =
1/2
{ϕ ∈ H 1 (Ω), ϕ = 0 on S, ϕ|Γ ∈ H00 (Γ)} such that
∫
a(u, υ) =
Γ

˜1
gυds ∀υ ∈ H0 (Ω).

(5.13)

Here the bilinear form is defined by
∫
a(u, υ) =
Ω

∫
∇u · ∇υ −

Ω

k 2 uυ −

∫
I(˜)υds.
u

(5.14)

Γ

The proof of uniqueness and existence for the solution of (5.13) can be found in [2]. In

115

this chapter, we focus on the RCS reduction for the cavity. In the two dimension, the RCS
is defined by:
σ(φ) = lim 2πr
r→∞
√
where r =

|us (r, φ)|2
,
|ui |2

(5.15)

x2 + x2 and φ is the observation angle.
1
2

If the observation angle coincides with the incident angle, it is called backscatter RCS
or monostatic RCS. Since most radar systems are monostatic [28], we only focus on the
reduction of backscatter RCS in this chapter. In order to evaluate the RCS, we need to
find the far field pattern of the scattered field. Since we assume the whole upper space is
vacuum, the far field pattern can be evaluated through the field at the opening of the cavity.
In particular, for TM polarization, the backscatter RCS can be computed by:
∫
2
4 k
σ :=
sin θ
ueikx1 cos θ dx .
k 2
Γ

(5.16)

We next introduce the design problem. The idea is to change the shape of the cavity so
that the backscatter RCS is minimized. It is easy to see that the design problem is trivial if
no constraint is imposed, since the RCS is zero for any angle and any frequency if no cavity
exists. Therefore, we impose the following constraints on the cavity:
• The opening Γ of Ω is fixed;
• The area of Ω is invariant;
• Ω is bounded in a set D.
In particular, as shown in Fig. 5.2, there exist two fixed regions K and D as the lower bound
and upper bound of the cavity. In order to satisfy the constraints, we choose to design the
cavity in the following way:
116

K

Ω
D

Figure 5.2: Constraint of the problem
The left boundary S1 of Ω can be parameterized by a C 2 function with respect to x2 .
The right boundary S2 is kept to be parallel to S1 . To satisfy the area invariance, the depth
of the cavity is assumed to be fixed, and the bottom of Ω is simply a line segment with the
end points connecting S1 and S2 .
More specifically, S1 is parameterized by f (x2 ), where f (x2 ) is a C 2 function chosen
from a compact subset with respect to C 1,β , 0 ≤ β < 1. Denote the admissible set as Λ:

Λ = {f ∈ C 2 ([−h, 0])| f (0) = 0, (f (x2 ), x2 ) ∈ D,
(f (x2 ), x2 ) ∈ K, ∥f ∥ 2 ≤ C}.
/
C

(5.17)

where h is the depth of the cavity, C is a fixed positive constant and the norm ∥ · ∥ 2 is the
C

117

usual maximal norm up to the second derivative, namely:

∥f ∥ 2 =
C

2
∑

max |f i (x)|.
i=0 x∈[−h,0]

(5.18)

The compactness of the admissible set is used to ensure the solution would converge. We
also denote the width of the opening Γ as a.
Given the incidence angle, our goal is to solve the following minimization problem:

min σ(S1 ).
S1 ∈Λ

(5.19)

The existence of the minimizer for the objective functions is followed from the theorem:
Theorem 5.2.1. There exists at least one minimizer for the problem (5.19) in Λ.
The theorem will be proved in the following section. Note that we only give the existence of the minimizer. Whether the minimizer is unique or not is still open, nevertheless,
numerical experiments show that there are generally more than one minimizer.

5.3

Existence of the minimizer

The existence of the minimizer is based on the continuous dependence of the solution u to
(5.13) on the boundary S1 . We first define the mapping

F : S1 → u|Γ .

(5.20)

As the proof for the well-posedness of the solution u, the property of the boundary operator
I(·) plays an important role in the continuity argument. We begin by recalling a lemma in
118

[2].
1/2
Lemma 5.3.1. The boundary operator I: H00 (Γ) → H −1/2 (Γ) is continuous, where
1/2
H −1/2 (Γ) is denoted as the dual space of H00 (Γ).
The following theorem shows the mapping F is continuous with respect to S1 . A similar
result for obstacle scattering problem can be found in [30].
1/2
Theorem 5.3.2. The mapping F : S1 → u|Γ is continuous from C 1 into H00 (Γ).
1/2
Proof. For u|Γ ∈ H00 (Γ), a direct calculation shows:
∫
−Re

∫ √
I(˜)u = Im
u
Γ

R

ˆ
k 2 − ξ 2 |u(ξ)|2 dξ ≥ 0.
˜

(5.21)

Therefore, there exist positive constants C1 , C2 such that
∫
∫
∫
u
Re{a(u, u)} = Ω ∇u · ∇u − Ω k 2 uu − Re Γ I(˜)uds
(5.22)
≥ C1 ∥∇u∥2 2
− C2 ∥u∥2 2
.
L (Ω)
L (Ω)

Hence, by the Lax-Milgram Lemma and the Fredholm alternative, the uniqueness of the the
scattering problem [2] implies the existence and the continuous dependence of the solution
˜1
on the right-hand side. Define the operator B : H0 (Ω) → H −1 (Ω), where H −1 (Ω) is the
˜1
dual space of H0 (Ω), such that:

< Bu, υ >= a(u, υ).

˜1
where < ·, · > is the dual system between H −1 (Ω) and H0 (Ω).
119

(5.23)

From the discussion above, the operator B is invertible and holds

∥u|Γ ∥ 1/2
≤ ∥u∥ ˜ 1
≤ ∥B −1 ∥∥g∥ −1/2 .
H0 (Ω)
H
(Γ)
H00 (Γ)

(5.24)

Define the transformation:

(y1 (x), y2 (x)) = (x1 − f (x2 ), x2 ) for (x1 , x2 ) ∈ Ω,

(5.25)

where f (x2 ) is the parametrization of S1 . It maps the cavity Ω onto the rectangle R =
[0, a] × [−h, 0]. Then the bilinear form can be rewritten as:

a(u, υ) =

∫ { ∑
2
R

j,k=1

}
∫
∂u ∂υ
2 uυ J −
−k
bjk
I(˜)υds,
u
∂yj ∂yk
Γ

(5.26)

where J denotes the Jacobian of the transformation and the coefficients bjk are given by
2
∑ ∂yj ∂y
k.
bjk =
∂xi ∂xi
i=1

(5.27)

Since f (x2 ) ∈ C 1 , the coefficients in (5.27) are well defined. Denote the C 1 norm of f (x2 )
as ∥f ∥. Then bjk depends continuously on ∥f (x2 )∥ and the Jacobian J is identically one, so
the bilinear form a also depends continuously on the C 1 norm of f (x2 ). Let Bf denote the
operator B that depends on f . Using the perturbation argument on the Neumann series, a
small perturbation δf to f results in:

∥B −1 − B −1 ∥ ≤ C1 (f )∥δf ∥,
f
f +δf

120

(5.28)

where C1 is a constant depending on f . Therefore, from (5.24), there exists another constant
C2 , such that
∥uf |Γ − uf +δf |Γ ∥ 1/2
≤ C2 (f )∥δf ∥.
H
(Γ)

(5.29)

The proof is now complete.
Now we are ready to give the proof of Theorem 5.2.1 for the existence of the minimizer.
Proof. Since the admissible set Λ is a compact set with respect to C 1,β , 0 ≤ β < 1, the
existence result immediately follows from Theorem 5.3.2 and the fact that σ depends continuously on u|Γ .
We remark that the proof is given for the specific admissible set. For a more general
admissible set, like nonsmooth curve for S1 , the existence of the minimizer will be more
involved. We next proceed to give the derivation of the domain derivative, which leads to
the descent directions for the objective function.

5.4

Domain Derivative

Domain derivatives have been studied extensively in various problems, including shape optimization problems[47], [41], inverse scattering problems[27] and the image of local surface
displacement[7]. The general idea is based on the transformation of coordinates. The perturbed domain is mapped back to the original domain, so that the difference at the boundary
is transformed into the difference at the coefficients in the variational equation. Rigorous
calculation on the coefficients leads to the equation for the domain derivative.
2
Let V (x) be a vector field defined on S1 such that V (x) = (a(x2 ), 0) ∈ C0 (S1 ; R2 ),
2
where C0 (S1 ; R2 ) denotes the twice differentiable vector fields on S1 with a(0) = 0. For
121

such a vector field V (x), define

δ
2
S1 = {x + δV (x)|x ∈ S1 , V (x) ∈ C0 (S; R2 )}

(5.30)

as the perturbation of S1 with respect to V (x). Since S2 is parallel to S1 , V (x) is also
defined on S2 by a horizontal shift and the same perturbation is also applied to S2 . The
domain derivative of the forward mapping F with respect to the direction V (x) is defined
as:
δ
′ (S ) := lim F(S1 ) − F(S1 )
F 1
δ
δ→0

(5.31)

2
Theorem 5.4.1. Let u be the solution of (5.3). If S1 is C 2 and V (x) ∈ C0 (S1 ; R2 ), then
the domain derivative F ′ (S1 ) exists. Moreover, F ′ (S1 ) = u′ |Γ , where u′ solves:
∆u′ + k 2 u′ = 0

in

Ω,

u′ = −(V · n)

∂u
∂n

on S1 ∪ S2 ,
(5.32)

u′ = 0

on S3 ,

∂u′
= I(u′ ) on
∂n

Γ.

Proof. Assume V (x1 , x2 ) = (a(x2 ), 0). For a small number δ, let Ωδ denote the perturbed
cavity with respect to Ω along the direction V . The total field uδ satisfies the following

122

equation:

∆uδ + k 2 uδ = 0 in Ωδ ,
uδ = 0 on S,
∂uδ
∂n

(5.33)

= I(uδ ) + g

on Γ,

with the weak form:
∫

∫
Ωδ

∇u · ∇υ −

k 2 uυ −

Ωδ

∫

∫
I(˜)υds =
u

Γ

gvds,
Γ

˜1
∀v ∈ H0 (Ωδ ).

(5.34)

Define the mapping q:
(y1 , y2 ) = (x1 − δa(x2 ), x2 ),

(5.35)

which maps Ωδ back to Ω. Let uδ = uδ ◦ q −1 and υ = υ ◦ q −1 . The weak form in (5.34)
˜
˜
can be rewritten as:
∫ ( ∑
2
Ω

j,k=1

)
∫
∫
˜
∂ uδ ∂ υ
˜
2 uδ υ J −
δ )˜ =
−k ˜ ˜
υ
bjk
I(u υ
˜
g˜ds,
∂yj ∂yk
Γ
Γ

(5.36)

where bjk has the same definition as (5.27). Here, the J is the Jacobian for the transformation
q, which is one. Define the new bilinear form

aδ (˜, υ) =
u

∫
Ω

2
∑
j,k=1

∫
∂ u ∂υ
˜
2 uυ −
˜
−k ˜
I(u)υ,
˜
bjk
∂yj ∂yk
Γ

123

(5.37)

˜1
where u, υ ∈ H0 (Ω). In particular, uδ in (5.36) is the unique solution of
˜
˜
aδ (˜δ , υ) =
u

∫
gυds,
Γ

˜1
∀υ ∈ H0 (Ω).

(5.38)

Assume u is the solution of (5.13), it holds

a(˜δ − u, v) = a(˜δ , v) − aδ (˜δ , v)
u
u
u
∫
=
Ω

∇˜δ · ∇v −
u

∫
Ω

2
∑
j,k=1

∂ uδ ∂v
˜
bjk
.
∂yj ∂yk

(5.39)

Extend the definition of V (x) to Ω by the following:

V (x) = (a(x2 ), 0) for (x1 , x2 ) ∈ Ω.

(5.40)

It is easy to see that the matrix

(bjk )2×2 = I − δ(∇V + (∇V )T ) + O(δ 2 ),

(5.41)

where I is the 2 × 2 identity matrix. By the continuous dependence result in Theorem 3.2,
if δ → 0, we have

∫
uδ − u
˜
, v) →
∇uT (∇V + (∇V )T )∇v.
a(
δ
Ω

(5.42)

˜δ
Let u0 = limδ→0 u −u . Then u0 solves the following problem
δ
∫
a(u0 , v) =

Ω

∇uT (∇V + (∇V )T )∇v.

124

(5.43)

2
˜1
The standard regularity for elliptic equations shows u ∈ H0 (Ω) ∩ Hloc (Ω) [21], so let v ∈
2
˜1
H0 (Ω) ∩ Hloc (Ω). According to the vector identity and the fact that ∇ · V = 0, it holds
∇uT (∇V + (∇V )T )∇v = ∇(V · ∇v) · ∇u + ∇(V · ∇u) · ∇v − ∇ · [(∇u · ∇v)V ].

Using Green’s formula and the fact that V = 0 on Γ, we have
∫
a(u0 , v) =

∫

Ω

−(V · ∇v)∆u + ∇(V · ∇u) · ∇v +

S

(V · ∇v)

∂u
− (∇u · ∇v)V · nds.
∂n

Since u = v = 0 on S, it follows

(V · ∇v)

∂u
− (∇u · ∇v)V · n = 0.
∂n

Based on the equation of u, a(u0 , v) can be rewritten as
∫
a(u0 , v) =

Ω

k 2 u(V · ∇v) + ∇(V · ∇u) · ∇v.

The divergence theorem and ∇ · V = 0 yield that
∫
a(u0 , v) =

Ω

k 2 u(V · ∇v) + ∇(V · ∇u) · ∇v + k 2 (∇ · V )uv
(5.44)

∫
=
Ω

∇(V · ∇u) · ∇v − k 2 (V · ∇u)v.

125

The above identity implies u0 satisfies

∆u0 + k 2 u0 = (∆ + k 2 )(V · ∇u) in Ω,
u0 = 0 on S,

(5.45)

∂u0
= I(u0 ) on Γ.
∂n
Let u′ = u0 − V · ∇u. Then u′ solves (5.32). Note that V = 0 on Γ, so
u′ |Γ = u0 |Γ = lim
δ→0

uδ − u
˜
|Γ = F ′ (S).
δ

(5.46)

The proof is complete.
Remark: The domain derivative can be obtained for more general vector fields (i.e.
nonzero component in the vertical direction) following the same idea.
As the following theorem shows, the domain derivative offers a way to find the descent
direction of the objective function (5.19), which is crucial for the optimization algorithm.
Theorem 5.4.2. For S1 ∈ C 2 , let u be the solution of the forward scattering problem (5.12)
and v be the solution of the following boundary value problem:

∆v + k 2 v = 0

in

Ω,

v = 0

on

S,

(5.47)
∫

′
∂v
ueikx cos θ ds′
= I ∗ (v) + 2k sin2 θe−ikx cos θ
∂n
Γ

126

on

Γ,

2
where I ∗ is the adjoint operator of I. If V (x) ∈ C0 (S1 ; R2 ) satisfies that
∫

then

∫
∂u ∂v
∂u ∂v
(V (x) · n)Re(
(V (x) · n)Re(
)ds +
)ds < 0,
∂n ∂n
∂n ∂n
S1
S2

(5.48)

δ
σ(S1 )−σ(S1 )
dσ(S1 )
|δ=0 = limδ→0
< 0. Here V (x) is also defined on S2 by a
dδ
δ

horizontal shift.
Proof. Assume u′ solves the equation (5.32). Apply Green’s theorem
∫

v∆u′ + k 2 vu′ =

∫

∫
∫
∂u′
∂v ′
v
ds −
u ds +
v∆u′ + k 2 vu′ .
∂Ω ∂n
∂Ω ∂n
Ω

Ω

(5.49)

Since both u′ and v satisfy the Helmholtz equation, the above equation implies
∫

∫
∂u′
∂v ′
v
ds =
u ds.
∂Ω ∂n
∂Ω ∂n

(5.50)

Recalling the boundary condition for u′ and v, we have
∫

∫
∂u′
˜
v
ds =
vI(u′ )ds,
∂n
∂Ω
Γ

and
∫

∂v ′
u ds =
∂Ω ∂n

∫

∫
∫
2 θeikx cos θ u′ ueikx′ cos θ ds′ ds
∗ (˜)u′ ds +
I v
2k sin
Γ
Γ
Γ
∫
∂u ∂v
−
(V (x) · n)
ds.
∂n ∂n
S1 +S2

127

(5.51)

The property of the adjoint operator implies
∫

˜
vI(u′ )ds =

∫

Γ

I ∗ (˜)u′ ds.
v

(5.52)

Γ

Combining (5.50) and (5.52), we obtain
∫

2k sin2 θeikx cos θ u′

Γ

∫

′
ueikx cos θ ds′ ds =

Γ

∫
S1 +S2

(V (x) · n)

∂u ∂v
ds.
∂n ∂n

(5.53)

The conclusion follows immediately from the fact that the real part of the left hand side is
exactly

dσ(S1 )
|
.
dδ δ=0

For abbreviation, let us denote

φ1 (x) = Re( ∂u ∂v )|S ,
∂n ∂n 1
φ2 (x) = Re( ∂u ∂v )|S ,
∂n ∂n 2

where x = (f (x2 ), x2 ). Since S1 and S2 are parallel to each other, according to the theorem,
the descent direction can be chosen to satisfy:

V (x) · n = φ2 (x) − φ1 (x),

(5.54)

f ′ (x2 )
−1
,√
) is the unit out normal for S1 . Then the RCS may
′2 (x )+1
′2 (x )+1
f
f
2
2
be reduced after a small perturbation to the cavity along the direction V (x). However, such

where n = ( √

2
a V (x) may not be in C0 (S1 ; R2 ). In order to keep the perturbed boundary in the admissible
set, our idea is to project V (x) into a given vector space by solving a least square problem,
and then move the boundary along the projected direction and expect the RCS to decrease
128

along that direction. Details will be given in the following section.

5.5

The optimization algorithm

Without loss of generality, assume the depth of the cavity h = π. We choose to represent
the boundary by the finite Fourier series
M
∑

f (x2 ) =

βj (cos jx2 − 1),

x2 ∈ [−π, 0],

(5.55)

j=1
where β belongs to the admissible set UM

UM = {β ∈ RM : ρ1j ≤ βj ≤ ρ2j ,

j = 1, 2, · · · , M }.

(5.56)

Here we drop the sine modes since cosine modes are already sufficient to be a base for
a function f ∈ C 2 [−π, 0] when M → ∞ if we add a technical assumption that f ′ (0) =
f ′ (−π) = 0. In that case, f can be extended to a 2π periodic C 2 function by an even
extension over [−π, π], and all the sine modes will vanish in the Fourier series. In (5.55),
the additional term ”-1” for each cosine function is used to guarantee f (0) = 0. This kind
of parametrization can be seen in [16] for the design of cylindrical cavity using longitudinal
bending with M = 1. It is also a common parametrization for the boundary of obstacle in
the setting of inverse scattering problems[30]. Compared to directly controlling the value
on f with uniform discretization, controlling the coefficients can avoid high oscillations for
the output boundary. It is easy to see that by imposing narrow changes for large j in the
admissible set UM . we can obtain a smooth curve with low oscillation, which is more useful
from the practical point of view.
129

Assume we perturb the coefficients β by adding δβ = {δβ1 , δβ2 , · · · , δβM }. By linearity,
the corresponding change to f (x2 ) is:

δf =

M
∑

δβj (cos jx2 − 1).

(5.57)

j=1
From the discussion above, δf can be chosen based on the equation (5.54) in order to get
the RCS reduced. We therefore solve the least square problem
∫ 0 (

−1

min
δf (x2 ) √
− (φ2 (f (x2 ), x2 ) − φ1 (f (x2 ), x2 ))
δβ∈RM −h
f ′2 (x2 ) + 1

)2
dx2 .

(5.58)

The equation can be solved by any linear solver after finite discretization. If β new is the
perturbed coefficients based on δβ, it is necessary for βnew to satisfy the admissible condition. This can be accomplished by projecting β new back into the admissible set UM . Let
β = β old + αδβ for α ∈ (0, 1]. Then β new = P (β), which is defined as the following


 ρ
 1j



new =
βj
 ρ
 2j



 β
j

if

βj < ρ1j ,

if

βj > ρ2j ,

for j = 1, 2, · · · , M.

otherwise,

Here we do not try to find the optimal step length α due to the computational cost.
Instead, we use the bisection search on the interval [0, 1] to find a step length that can
reduce the RCS, and then we continue to the next iteration. This algorithm follows the
general idea in [20] for the optimal design of periodic antireflective structures.
We now state the minimization algorithm:
• Choose an initial value β 0 ∈ UM , set m = 0.
130

• for m=1,2,. . . , until stop
– Solve the forward problem (5.13), adjoint problem (5.47) and evaluate φ1 (x) and
φ2 (x). Solve the least square problem (5.58) to obtain δβ.
– if ∥δβ∥ < ϵ, stop;
– else let α = 1, while α > ϵ do
β m = P (β m−1 + αδβ);
if σ(S1 (β m )) < σ(S1 (β m−1 )), stop while;
α ← α/2;
– end.
• end.
To solve the scattering problem for a cavity, since only the boundary data is needed during the
optimization, we thus choose the boundary integral method. In particular, for the forward
scattering problem, we solve the equation:
1 u(x) = − ∫ G∂ u| + ∫ ∂ Gu| − ∫ GI(u| ) − ∫ Gg,
n s
Γ
Γ
S
Γ n
Γ
Γ
2
∫

∫
∫
∫
G∂n u|s = Γ ∂n Gu|Γ − Γ GI(u|Γ ) − Γ Gg
S

∀x ∈ Γ
(5.59)

∀x ∈ S

i 1
where G(ρ, ρ′ ) = − 4 H0 (k0 |ρ − ρ′ |). For the adjoint problem, we simply replace g by
2k sin2 θeikx cos θ

∫
Γ

.
131

′
ueikx cos θ ds′

Special care must be taken in dealing with the singularity of G and the corner point on
the boundary S. Effective numerical methods for (5.59) can be seen in [13].

5.6

Numerical experiment

In this section, we present three numerical examples to illustrate the effectiveness of the
algorithm and show that the backscatter RCS can be effectively reduced by appropriate
shaping. Throughout all the examples, we assume the admissible set UM for the coefficient
is {β ∈ R20 : −1/j 1.5 ≤ βj ≤ 1/j 1.5 ,

j = 1, 2, · · · , 20}. As we mentioned before,

imposing narrower change for larger j results in a less oscillatory shape. However, for specific
engineering problems, those design variables should be chosen accordingly. The initial cavity
has a rectangular shape as shown in Fig. 5.3, with the height h = 0.3m and width of the
opening a = 1m. The boundary is obtained by simply assuming βj = 0, for j = 1, 2, · · · , 20.
The forward scattering problem and the adjoint problem are solved by the boundary
integral method given in (5.59). The opening of the cavity is discretized uniformly by 512
subintervals as well as the bottom. For the curved boundary S1 and S2 , since both of
them are parameterized by a function on [−h, 0], we simply discretize [−h, 0] uniformly by
512 subintervals and assume S1 and S2 are piecewise constant on each subinterval. Note
that higher order elements will result in a more accurate design, and our approximation is
only given for illustration. The nonlocal boundary integral in (5.10) is transformed into a
hypersingular integral by the Fourier transform:
k0 ∫
I(u) = − 2j =
Γ

1 H (1) (k |x − x′ |)u(x′ , 0)dx′ ,
0
|x−x′ | 1

132

∫
where = ·dx denotes the Hadamard principal value (or finite part) integral. We refer readers
to [9] for details. We numerically evaluate the hypersingular integral by a Newton-Cotes
formula proposed by Wei in [48]. Although these two integrals are equivalent mathematically,
the hypersingular integral form does not need to truncate into finite terms compared to the
Fourier type integral and therefore is more accurate numerically.
In practice, the unit of RCS is in dB and the value is expressed as the logarithm of σ:

RCS = 10 log10 σ dB.

For the first example, we focus on reducing the backscatter RCS at the three wave numbers

Figure 5.3: Initial shape

k = 2π, k = 8π and k = 32π. We consider the minimization of the backscatter RCS at

133

the normal incidence. It corresponds to a low frequency wave(ω = 300MHz) for k = 2π.
After 500 iterations, we obtain the optimized shape shown in Fig. 5.4(a). Fortunately,
for such a low frequency case, the backscatter RCS is constantly decreasing during the
iteration. The reduction finally becomes extremely small and we take it to be convergent.
The comparison between the initial shape and the optimal shape for the backscatter RCS is
shown in Fig. 5.4(b). From the graph we can see the an average of 15dB reduction is achieved
on all the incident angles, although the cavity is only designed to reduce the backscatter RCS
at normal incidence. This speciality can only be found in the low frequency case, where there
is almost no oscillation in the backscatter RCS across all the incident angles.
For k = 8π, the optimal shape is given in Fig. 5.5. The RCS has been decreased by
more than 15dB at the normal incidence. As shown in Fig. 5.4, the range of incident angles
with RCS reduced is from 90◦ to 110◦ , which is smaller than the interval for k = 2π. For
k = 32π, the incident frequency is considerably high(4.8GHZ). the optimal shape is given in
Fig. 5.6. Similar to k = 8π, although the RCS has been decreased at the normal incidence
as expected, the reduction does not occur for all the incident angles. This is understandable
since roughly, high frequency scattering can be viewed as ”billiard-ball-bouncing scattering”,
where the scattered field can be approximated by the ray-tracing technique[16]. In such a
case, reduction of RCS at one incident angle always accompanies an increase at some other
angles. Nevertheless, it is still very useful since the angles near the normal is considered as
the ”threat sector” during the detection of an object[28].
For the second example, we try to reduce the backscatter RCS with the incidence angle
θ = 120◦ for k = 32π. Although the angle is less important than the normal incidence from
the practical point of view, it is given for illustration purpose. The result is shown in Fig. 5.7.
We see the RCS is decreased successfully around θ = 120◦ , but increased around θ = 100◦ .
134

This confirms the conclusion that RCS reduction through shaping can not be expected for
all angles for high frequency incidence[28].

5.7

Conclusion

In this chapter, we consider the optimal design of an open cavity for RCS reduction. The
problem is formulated as a minimization problem with the Helmholtz equation as the constraint. In addition, we make several assumptions on the shape of the cavity in order to make
the problem nontrivial. The Helmholtz equation satisfied by the total field is reduced to a
bounded domain problem by introducing a transparent boundary condition. The existence
of the minimizer has been proved based on the continuous dependence of the boundary. The
descent direction with respect to the boundary is also derived based on the domain derivative. A simple gradient-based optimization method is introduced with the integration of the
descent direction. We also test our algorithm by several numerical examples. The result
shows our approach is fairly effective in the RCS reduction.

135

(a)

Backscatter RCS(dB)

20

0

−20

−40

−60
100

120
140
Incident angle

160

180

(b)
Figure 5.4: (a)Optimal shape for k = 2π at normal incidence; (b) RCS comparison for
k = 2π, ”- - -” for initial shape, ”—–” for optimal shape.
136

(a)

Backscatter RCS(dB)

20

0

−20

−40

−60
100

120
140
Incident angle

160

180

(b)
Figure 5.5: (a)Optimal shape for k = 8π at normal incidence; (b) RCS comparison for
k = 8π, ”- - -” for initial shape, ”—–” for optimal shape.
137

(a)

Backscatter RCS(dB)

20

0

−20

−40

−60

100

120
140
Incident angle

160

180

(b)
Figure 5.6: (a)Optimal shape for k = 32π at normal incidence; (b) RCS comparison for
k = 32π, ”- - -” for initial shape, ”—–” for optimal shape.
138

(a)

20
0
−20
−40
−60
−80
−100

100

120

140

160

180

(b)
Figure 5.7: (a)Optimal shape for k = 32π at the incidence angle θ = 120◦ ; (b) RCS comparison for k = 32π, ”- - -” for initial shape, ”—–” for optimal shape.
139

Chapter 6
Summary and future work

6.1

Summary

Here is a brief summary of this thesis.

6.1.1

Gaussian beam method

A global Gaussian beam solution for wave equations in the bounded domain is constructed
through several steps. Beams are first initialized by the Multiscale Gaussian wavepacket
transform with odd, periodic extension of the initial condition. Propagation of these beams
follows the general Gaussian beam method presented in [45]. Boundary condition is satisfied
by the superposition of incident beams and reflected beams. At the final time, beam solutions
are summed up with the method-of-images to form a global solution of the wave equation.
A criteria to identify significant beams is proposed in order to reduce the computational
cost. Theoretical convergence rate is analyzed for the method and justified by the numerical examples. In particular, we test our algorithm in different bounded domains, namely,
one dimensional interval, two dimensional rectangular domain and two dimensional circular
domain. Different velocities, including constant case, linear case and sinusoidal case are all
considered in the numerical tests. The successful implementation of these tests validate the
effectiveness of the Gaussian beam method in bounded domain problems.

140

6.1.2

Optimal design for a cavity

6.1.2.1

Coating material

We propose an effective way to design a multilayered RAM for reducing RCS of a cavity
both in TM and TE polarizations. The problem is formulated as a minimization problem.
The descent direction for the cost function is evaluated through the adjoint state method.
Subsequently, SQP is integrated with the gradient to obtain the optimal absorbing materials
for the cavity. The algorithm is implemented with a fast and accurate direct solver for the
scattering problem. Numerical results show that RCS is reduced significantly with RAM
coated at the bottom of the cavity. It is also observed that different objective function
generally leads to different RAM design and RCS reduction.

6.1.2.2

Shape optimization

The optimal shape design of an open cavity for RCS reduction is considered in this thesis.
We formulate the design problem as a minimization problem with a Helmholtz equation
as the constraint. To avoid trivialness, several geometrical constraints are imposed on the
shape of the cavity. Similar to the optimization problem in coating material, a transparent
boundary condition is introduced to the Helmholtz equation in order to reduce the domain
of interest into a bounded one. The existence of the minimizer has been proved based on the
continuous dependence of the boundary. The descent direction with respect to the boundary
is also derived based on the domain derivative. A simple gradient-based optimization method
is introduced with the integration of the descent direction. The algorithm is tested through
several numerical examples. It shows the approach is fairly effective in the RCS reduction.

141

6.2

Future work

Lots of work needs to be done in the future. For high frequency problem, we consider
generalizing the algorithm to the case when the velocity term is discontinuous. Compared to
the smooth velocity, discontinuous velocity will generate a reflected beam and a transmitted
beam for each incident beam at the discontinuity. Computational complexity is increased
along the increase of the number of beams. New idea must be introduced to keep the number
of beams under control. We also consider implementing the Gaussian beam method by
parallel computation on GPU, since the beams are propagating independently in the domain
and the special structure of GPU allows each beam being evaluated by a single processor.
Significant acceleration are expected to be achieved when hundreds of such processors are
running simultaneously. Another future work is to consider three dimensional bounded
domain problem, especially high frequency Maxwell’s equation.
For the optimal design problem, future work includes combining coating material and
the shape optimization together or extending the method to three dimensional cavities.
Instead of Helmholtz equation, the three dimensional problem is formulated as Maxwell’s
equation with much more freedom at the incident angle. Although theoretical result can be
derived accordingly, numerical simulation are more difficult than the two dimensional case.
In fact, computation for large cavities has been considered as one of most challenging work
in electromagnetic community. Many methods have been proposed over the last few decades
[9, 25]. In terms of optimization, the method requires an efficient and accurate solver, which
can be extremely challenging given all the level of difficulties.
As we mentioned at the beginning, developing a fast solver for high frequency problem
and optimal design of a cavity should not be treated as a separate work. The long term

142

goal is to design a fast algorithm to solve the scattering problem for large cavities, which is
essentially a high frequency problem by rescaling and then used in the optimization.

143

BIBLIOGRAPHY

144

BIBLIOGRAPHY

[1] H. Ammari, G. Bao, and A. W. Wood. An integral equation method for the electromagnetic scattering from cavities. Math. Meth. Appl. Sci., 23:1057–1072, 2000.
[2] H. Ammari, G. Bao, and A. W. Wood. Analysis of the electromagnetic scattering from
a cavity. Japan J. Indust. Appl. Math., 19(2):301–310, 2002.
[3] H. Ammari, G. Bao, and A. W. Wood. A cavity problem for maxwells equation. Math.
Meth. Appl. Sci., 9:249–260, 2002.
[4] V. M. Babich and V. S. Buldyrev. Asymptotic methods in short wave diffraction problems
(in Russian). Nauka, Moscow, 1972.
[5] G. Bao, J. Gao, and P. Li. Analysis of direct and inverse cavity scattering problems.
Numer. Math. Theor. Meth. Appl., 4:419–442, 2011.
[6] G. Bao, J. Gao, J. Lin, and W. Zhang. Mode matching for the electromagnetic scattering
from three-dimensional large cavities. IEEE Antennas Wireless Propagat., 60:2004–
2010, 2012.
[7] G. Bao and J. Lin. Imaging of local surface displacement on inifinte ground plane: the
multiple frequency case. SIAM J. Appl. Math., 71:1733–1752, 2011.
[8] G. Bao, J. Qian, L. Ying, and H. Zhang. A convergent multiscale gaussian-beam
parametrix for the wave equation. Commun. PDEs, 38:1–43, 2013.
[9] G. Bao and W. Sun. A fast algorithm for the electromagnetic scattering from a large
cavity. SIAM J. Sci. Comput., 27:553–574, 2005.
[10] G. Bao, K. Yun, and Z. Zhou. Stability of the scattering from a large electromagnetic
cavity in two dimensions. SIAM J. Math. Anal., 44(1):383–404, 2012.
[11] G. Bao and W. Zhang. An improved mode-matching method for large cavities. IEEE
Antennas Wireless Propagat. Lett., 27:393–396, 2005.
[12] S. Bougacha, J. Akian, and R. Alexandre. Gaussian beams summation for the wave
equation in a convex domain. Commun. Math. Sci., 7:973–1008, 2009.
145

[13] C. Brebbia and J. Dominguez. Boundary Elements: An Introductory Course, Second
edition. WIT Press/Computational Mechanics Publications, Southampton, 1998.
[14] R. Burkholder and P. Pathak. Analysis of em penetration into and scattering by electrically large open waveguide cavities using gaussian beam shooting. Proc. IEEE, 79:1401–
1412, 1991.
[15] V. Cerveny, M. Popov, and I. Psencik. Computation of wave fields in inhomogeneous
media-gaussian beam approach. Geophys. J. R. Astr. Soc., 70:109–128, 1982.
[16] R. Chou. Reduction of the Radar Cross Section of Arbitrarily Shaped Cavity Structures.
PhD thesis, University of Illinois, Champaign, 1987.
[17] R. Chou and S. Lee. Modal attenuation in multilayered coating waveguide. IEEE Trans.
Microwave Theory Tech., 36:1167–1176, 1988.
[18] G. Cohen. Higher-order Numerical Methods for Transient Wave Equations. Springer,
2001.
[19] D. Colton and R. Kress. Inverse Acoustic and Electromagnetic Scattering Theory, Applied Mathematical Sciences 93. Springer-Verlag, Berlin, 1998.
[20] D. C. Dobson. Optimal design of periodic antireflective structures for the helmholtz
equation. Euro. J. Appl. Math., 4:321–340, 1993.
[21] L. Evans. Partial Differential Equations, Second edition. American Mathematical Society, Providence, 2010.
[22] N. Hill. Gaussian beam migration. Geophysics, 55:1416–1428, 1990.
[23] J. Huang and A. Wood. Numerical simulation of electromagnetic scattering induced by
an overfilled cavity in the ground plane. IEEE Trans. Antennas Propagat., 4:224–228,
2005.
[24] P. Huddleston. Scattering from conducting finite cylinders with thin coatings. IEEE
Trans. Antennas Propagat., 35:1128–1136, 1987.
[25] J. Jin. A finite element-boundary integral formulation for scattering by threedimensional cavity-backed apertures. IEEE Trans. Antennas Propagat., 39:97–104,
1991.

146

[26] J. Jin. The Finite Element Method in Electromagnetics, 2nd edition. Wiley, New York,
2002.
[27] A. Kirsch. The domain derivative and two applications in inverse scattering theory.
Inverse Problem, 9:81–96, 1993.
[28] E. Knott, J. Shaeffer, and M. Tuley. Radar Cross Section, Second edition. Scitech
Publishing Inc, Releigh, NC, 2004.
[29] H. Kreiss and N. Petersson. A second-order accurate embedded boundary method for
the wave equation with dirichlet data. SIAM J. Sci. Comput., 27:1141–1167, 2006.
[30] R. Kress and Z. Axel. On the numerical solution of the three dimensional inverse
obstacle scattering problem. J. Comput. Appl. Math., 42:49–62, 1992.
[31] S. Leung and J. Qian. Eulerian gaussian beam methods for schr¨dinger equations in
o
the semi-classical regime. J. Comput. Phys., 228:2951–2977, 2009.
[32] S. Leung and J. Qian. The backward phase flow and fbi-transform-based eulerian gaussian beams for the schrodinger equation. J. Comput. Phys., 229:8888–8917, 2010.
[33] S. Leung, J. Qian, and B. R. Eulerian gaussian beams for high frequency wave propagation. Geophysics, 72:SM61–SM76, 2007.
[34] H. Ling, R. Chou, and S. Lee. Shooting and bouncing rays: Calculating the rcs of an
arbitrarily shaped cavity. IEEE Trans. Antennas Propagat., 37:194–205, 1989.
[35] J. Liu and J. Jin. A special higher order finite-element method for scattering by deep
cavities. IEEE Trans. Antennas Propagat., 48:694–703, 2000.
[36] V. P. Maslov. The Complex WKB Method for Nonlinear Equations I: Linear theory.
Birkhauser Verlag, Basel, 1994.
[37] H. Mosallaei and Y. Rahmat-Samii. Rcs reduction of canonical targets using genetic
algorithm synthesized ram. IEEE Trans. Antennas Propagat., 48(10):1594–1606, 10
2000.
[38] M. Motamed and O. Runborg. Taylor expansion and discretization errors in gaussian
beam superposition. Wave Motion, 47:421–439, 2010.

147

[39] J. Nocedal and S. J. Wright. Numerical Optimization, Second Edition, Springer Series
in Operations Research. Springer Verlag, Berlin, 2006.
[40] S. Ohnuki and T. Hinata. Rcs of material partially loaded parallel-plate waveguide
cavities. IEEE Trans. Antennas Propagat., 51:337–344, 2003.
[41] O. Pironneau. Optimal Shape Design for Elliptic Systems. Springer-Verlag, New York,
1984.
[42] M. M. Popov. A new method of computation of wave fields using gaussian beams. Wave
Motion, 4:85–97, 1982.
[43] J. Qian and L. Ying. Fast gaussian wavepacket transforms and gaussian beams for the
schr¨dinger equation. J. Comput. Phys., 229:7848–7873, 2010.
o
[44] J. Qian and L. Ying. Fast multiscale gaussian wavepacket transforms and multiscale
gaussian beams for the wave equation. SIAM J. Multiscale Modeling and Simulation,
8:1803–1837, 2010.
[45] J. Ralston. Gaussian beams and the propagation of singularities. studies in partial
differential equations. maa studies in mathematics, 23. edited by w. littman. pages
206–248, 1983.
[46] P. Smereka and G. Russo. The gaussian wave packet transform: Efficient computation
of the semi-classical limit of the schr¨dinger equation. part 1 formulation and the one
o
dimensional case. J. Comput. Phys., 233:192–209, 2013.
[47] J. Sokolowski and J. Zolesio. Introduction to Shape Optimization: Shape Sensitivity
Analysis. Springer-Verlag, Berlin, 1992.
[48] W. Sun and J. M. Wu. Newton-cotes formulae for the numerical evaluation of certain
hypersingular integral. Computing, 75:297–309, 2005.
[49] N. Tanushev. Superpositions and higher order gaussian beams. Comm. Math. Sci.,
6:449–475, 2008.
[50] N. Tanushev, B. Engquist, and R. Tsai. Gaussian beam decomposition of high frequency
wave fields. J. Comput. Phys., 228:8856–8871, 2009.
[51] N. Tanushev, J. Qian, and J. Ralston. Mountain waves and gaussian beams. SIAM J.
Multiscale Modeling and Simulation, 6:688–709, 2007.
148

[52] T. Van and A. Wood. Analysis of transient electromagnetic scattering from overfilled
cavities. SIAM J. Appl. Math., 64(2):688–708, 2004.

149