- Email: [email protected]
An efficient numerical method for exterior and interior inverse problems of Helmholtz equation
387
WAVE MOTION 13 (1991) 387-399 ELSEVIER
AN EFFICIENT NUMERICAL METHOD INVERSE PROBLEMS OF HELMHOLTZ
FOR EXTERIOR EQUATION
AND INTERIOR
S.L. WANG and Y.M. CHEN Department of Applied Mathematics and Statistics, State Universiiy of New York at Stony Brook, Stony Brook, NY 117943600, USA
Received 17 April 1990, Revised 27 July 1990
The generalized pulse-spectrum technique (GPST), an efficient and versatile iterative inversion algorithm, is used with adaptive grids to solve both exterior (scattering) and interior (cavity) boundary-shape inverse problems of two-dimensional Helmholtx equation. Numerical simulations of nontrivial examples are carried out to test the feasibility and to study the general characteristics of GPST without the real measurement data. It is found that GPST does efficiently produce very good resuits.
1. Introduction The exterior inverse problems can be exemplified by the inverse scattering problems in which the shape, size and orientation of a scatterer are determined from the knowledge of incident and scattered waves in the exterior of the scatterer. The interior inverse problems can be exemplified by the inverse closed cavity problems in which the shape, size and orientation of the closed cavity are determined from the knowledge of incident and scattered waves in the interior of the cavity. In general, the knowledge or known data of the incident and scattered waves can come from a mixture of different locations and time or frequency intervals, e.g., at (x,, t,,,) or (x,, w,), I= 1,2,3,. . . , L, m = 1,2,3,. . . , A4. Both the exterior and interior inverse problems are difficult to solve due to the fact that they often are both ill-posed and nonlinear. Here only the shapes, sizes and orientations of the scatterers and cavities in the exterior and interior inverse problems of Helmholtz equation respectively are treated, and it is our intention to consider this paper as a detailed presentation of Refs. [l] and [2]. The inverse scattering problems of time-harmonic scalar waves with the known data came from a single frequency and many spatial locations lead to the inverse scattering problems of the Helmholtz equation. For the inverse scattering problems of the Helmholtz equation in the low to resonance frequency region, many methods are available, e.g., method of analytic continuation, explicit optimization methods, implicit optimization methods, etc. The method of analytic continuation [3] and [4], relied on the numerical analytic continuation, subjects to severe instability from the ill-posedness. In the explicit optimization methods [SJ, [6] (an excellent review) - [12], the inverse scattering problem is explicitly formulated as a nonlinear optimization problem first, then a nonlinear optimization technique with regularization to control instability is used to solve the problem; the regularization can either be applied to the full nonlinear problem or to each linear sub-problem arising from the linearization in the nonlinear optimization technique. In the implicit optimization methods, e.g., the generalized pulse-spectrum technique (GPST) [l] and [2], based on a combination of a Newton-like iterative method and the Tikhonov regularization method [13], is a class of versatile and efficient iterative algorithm, for solving multi-parameter inverse 01658641/91/$03.50
@ 1991 - Elsevier Science Publishers B.V.
388
S. L. Wang, Y. M. Chen / Exterior and interior inverse problems
problems of a system of nonlinear partial can be performed in either the space-time
differential equations. In GPST, the parameter identification (pulse) domain or the space-frequency (spectrum) domain or
both. There the inverse
is formulated
differential
equations
For the inverse
scattering instead
problem
of a nonlinear
closed cavity problems
optimization
as a perturbation
or iteration
problem
of partial
problem.
of scalar waves, there is no existing
research
can be found except
the well known
Kac’s inverse eigenvalue problem, “can one hear the shape of a drum” [ 14]? Their main difference is that the known data of the inverse cavity problem here are the incident and scattered waves in the interior of cavity and the known data of Kac’s inverse problem are all the normal modes of vibration. While there is only a partial analytic solution to Kac’s problem [ 151, GPST has been used to solve simple
inverse cavity problems cavity
problems
numerically
here are quite
with some success [2]. It should be pointed
different
from the inverse
open
out that the inverse
cavity problems
closed
such as the vocal tract
inverse problem [ 161. Since the incipience of GPST in 1974 [17], it has been generalized, improved and applied to many practical inverse problems. If only the shapes, sizes and orientations of the cavities in the exterior and interior inverse problems of Helmholtz equation respectively need to be determined, which is the case here, then one can make GPST dedicated to these special types of simple inverse problems by incorporating an adaptive body-fitted grid system. In this way, the computational efficiency of GPST is greatly improved. First, the dedicated GPST inversion algorithm for these types of inverse problems is presented. Then various computational aspects of GPST are described. Next, numerical simulations of simple but non-trivial examples are carried out on the adaptive grid system to test the feasibility and to study the general characteristics of this dedicated GPST without the real measurement data. Finally, a comprehensive discussion of the numerical results and the future improvement on GPST to attend higher efficiency is given.
2. Generalized pulse-spectrum For simplicity,
the following
in E* are cosidered (v’+
technique (GPST)
k2)U =f(x),
r”*(au/ar
exterior
and interior
boundary
value problems
of the Helmholtz
equation
respectively,
-iku)
x=(r,B)EE*-.f2,,
+ 0,
u(x)=O,
as r + 00, (radiation
xES,=R(e),
(1)
condition)
and (V’+k*)u=f(x),
XE&,
xeS,=R(O),
z&)=0,
(2)
where E* is the two-dimensional Euclidean space, 0, is a star-shaped finite sub-domain in E* with the boundary S, , 0, is a sub-domain in E * surrounding n, and has a circular outer boundary S, , R, is also a sub-domain in E* surrounding fin2 and has a circular outer boundary S, concentric with S2, x = 0 is located in 0, (Fig. l), and the source f(x) has a compact support in either E*-0, or 0,. Here the exterior and interior inverse problems are to determine R(B) or S, from the measurement data or prescribed design data u(x) = h(x),
x = q E a3
and 0,)
respectively,
1=1,2,3
,...,
L,
(3)
such that u(x) satisfies (1) and (2) respectively. The GPST inversion algorithm begins by setting u n+, = u, +&l,,
R ,,+,=R,+SR,,
n=0,1,2,3
,...,
(4)
389
S. L. Wang, Y.M. Chen / Exterior and interior inverse problems
Fig. 1. Scheme of adaptive computational
grid system.
where R,( 0) is the initial guess for R( 0) and S is a small number such that the b-terms are smaller than their corresponding non-&-terms in some norms for numerical convergence. For the exterior inverse problem, upon substituting (4) into (1) with u and R replaced by u,+, and R ?I+17 respectively, one obtains (V’+k*)(u,+Su,)=f(x),
XEE*-~,.,,
n=0,1,2,3
,...,
as r+oo.
r”*[(~u,/~r-iku,)+(&3u,/~r-ik&,)]+O,
(5)
Upon neglecting terms of O( S*) and higher, the system for u, same as (1) is obtained, (V*+k*)u,
=f(x),
XE E*-R,,,,
r”‘(au,/ar-iku,)+O,
u,=O,
XER,,
n = 0, 1,2,3, . . . ,
(6)
as r+m,
and a similar system for au,, (V*+k*)Su, =0,
XE E*-a,,,,
&I, = -(au,/ar)GR,,
XE R,,
n =O, 1,2,3,
(7)
as r+m.
r”*(a&/ar-iksu,)+O,
*. ,
Similarly, for the interior inverse problem one obtains the system for u, same as (2), XE n,,,,
(V’+ k*)u, =.0x),
u,=O,
XER,,
n = 0, 1,2,3, . . . ,
(8)
and a similar system for au,,, (V’+ k*)Su, =0,
XE f‘&,,
au,, = -(au,/ar)GR,,
XE R,,
n =0,1,2,3,.
...
(9)
By using the method of Green’s function, from (7) or (9) one obtains an integral relation between SR, and au,,
I R.
[aG,(x, x’)/av;][au,/ar’]t%R,(
0’) dl:, = SU,,
n = 0, 1,2,3,.
..,
(10)
S.L. Wang, Y.M. Chen / Exterior and interior inverse problems
390
where G,(x, x’) is the Green’s function of (7) or (9), v,, is the unit normal vector on R, pointing toward the domain of interest, and dl, is the incremental line segment of R,. For speeding up the convergence and eliminating the unknown u,+], in (10) one can set x to q, 1= 1,2,3,. . . , L, and replace u,+,(x,) by the known data h(x,), I = 1,2,3, . . . , L. Hence a system of Fredholm integral equations of the first kind for the unknown SR, (0) is obtained as
I
[aG,(x,, x’)/av~][au,/ar’]sR,(e’)
dl:, = h(x,) - u,(x,),
n =0, 1,2,3,.
. . , 1= 1,2,3,.
. . , L.
(11)
a,,
Now eqs. (4), (6) and (11) and eqs. (4), (8) and (11) form the basic theoretical structure of each iteration of GPST for solving the exterior and interior inverse problems of the Helmholtz equation respectively. Note that with some minor modifications, GPST can be used with equal efficacy to solve both exterior and interior inverse problems of the Helmholtz equation with Neumann or impedance boundary conditions and with measurement data containing derivatives of u(x).
3. Computational
consideration
Theoretically, the efficiency of a numerical algorithm depends on the intrinsic mathematical structure of the algorithm, but in practice it depends even more on how efficiently one can implement every minute computational step which is very much computer-architectureand compiler-dependent. Here an adaptive irregular polar computational grid system (Fig. 1) is used. In this grid system, the circular finite computational domain consists of three connected but non-overlapping sub-domains, 0, , L& and L$, where 0, represents the scatterer or cavity, the interface S, between 0, and L!* represents the shapes, sizes and orientations of the scatterer or cavity, and LZ2and 0, represent the region in the exterior of the scatterer; let the grid point be denoted by (ri, 6,) = (i, j). In 0, and L&, the locations of the grid points are re-adjusted along fixed constant j-lines after each iteration to accommodate the shape, size and orientation of the scatterer or cavity respectively; in this way, one can achieve high accuracy with a small number of grid points; the boundary of the scatterer or cavity S, should be a constant i-line. In L&, the grid system is a fixed circular cylindrical grid to avoid the unnecessary computational effort in generating grid points; however, due to the truncation of the infinite space domain, the radiation condition in (6) and (7) is replaced by the approximate radiation condition [ 181, dv/dr + (2r)-‘v
- ikv = 0,
XE&, n=0,1,2,3
v = u,, au,,
,....
(12)
To deal with the irregular incremental quadrilaterals, one considers the area integration of the Helmholtz equation on a typical quadrilateral Ai,j,, centered at (i, j) with (i - l/2, j - l/2), (i - l/2, j + l/2), (i + l/2, j-1/2) and (i+ l/2, j+ l/2) as the four vertices. From (6) or (8) and Green’s identity, the area integral becomes
I L.” I
au,/dv,
dl,, + k2 I,,,.Iu”du~=la,j”~fdu.. ..
(13)
n=O,1,2,3,..., .*
where Li,j,, is the boundary of Ai,j,, and V, is the unit vector normal to Li,j,,. Upon using the rectangle rule and the center difference scheme to approximate the integrals and the normal derivative respectively, one obtains a 9-point finite difference formula for u,(x), ifl
1
i+l
1
a-i-1 p=i-l
YmaU”,m,P=f;,j, i=l,2,3
,...,
Z,j=l,2,3
,...,
J,n=0,1,2,3
,....
(14)
S.L. Wang, Y.M. Chen / Exterior and interior inverse problems Hence
the above discretizations
391
of (6) with (12) and of (8) lead respectively to
B ‘,“.U,,=F,
n=0,1,2,3
,...,
(15)
Bz-U.=F,
n=0,1,2,3
,...,
(16)
and
where By and Bz both are known N x N symmetric sparse matrices with half bandwidth J, N = I X J, U, is the N x 1 matrix with all u,,~,~, i = 1,2,3, . . . , I, j = 1,2,3, . . . , J, as its components, and F is the N x 1 matrix with all _&, i = 1,2,3,. . . , Z, j = 1,2,3,. . . , J, as its components. Here, eqs. (15) and (16) can be solved by using an efficient LU-decomposition sparse matrix solver. Straightforward discretization and evaluation of the system of Fredholm integral equations of the first kind (11) and calculating G,(x, x’) require very large computational effort. A more efficient way is to consider the discretization of (7) and (9) by using the same 9-point discretization schemes for (6) and (12) and for (8) with the boundary values of Su, as equivalent sources. One obtains the similar linear algebraic systems respectively for Su,, B~-SUn=(Cn-SR,,O)T,
n=0,1,2,3
,...,
(17)
B~~SUn=(Cn*6R,,,0)T,
n=0,1,2,3
,...,
(18)
and
where MY,,is the N x 1 matrix with all a~,,~,~,i = 1,2,3, . . . , I, j = 1,2,3, . . . , J, as its components, 6R, is the J x 1 matrix with all cSR,~ = SR,(Bj), j = 1,2,3, . . . , J, as its components, and C,, is the known J x J diagonal matrix obtained from the boundary condition (a known relation between au, and SR, at the boundary R,). Equations (17) and (18) lead respectively to (B~)-‘~(Cn~SRn,0)T=8Un,
n=0,1,2,3
,...,
(19)
(B~)-‘*(Cn-SRn,O)T=SUn,
n=0,1,2,3
,....
(20)
and
Upon collecting all of the equations in (19) and (20) corresponding to the data points g, I = 1,2,3, . . . , L, only and replacing u,+,(x,)-u,(x,) by h(x,)-u,(x,), 1=1,2,3,. .., L, one obtains the compact linear algebraic systems respectively for the unknown SR,, 0:.
SR, = En,
n=0,1,2,3
,...,
(21)
n=0,1,2,3
,...,
(22)
and DF-SR,=E,,
where 0: and Df are Lx J full matrices and E, is a Lx 1 matrix with h(x,) - u.(x,), 1= 1,2,3,. . . , L, as its components. It is important to point out that in forming (21) or (22) only L rows of (B’,“)-’ or (B’,“)-’ corresponding to the data points q, 1= 1,2,3,. . . , L, are needed. Since the LU-decomposition of B’,” and B’,”have been already performed in solving (15) and (16), they are saved so that the L rows of (By)-’ or (Br)-’ can be obtained by performing back substitution L times. Hence a great saving in computational effort is achieved. It is also clear that (21) and (22) are the most compact forms of the discretization of the system of Fredholm integral equations of the first kind (11).
392
S. L. Wang, Y. M. Chen / Exterior and interior inverse problems
Here D’,” and Df are rectangular and very ill-conditioned matrices. Then the well-known Tikhonov regularization method is used to overcome this difficulty, i.e., instead of solving (21) and (22), one solves their regularized systems respectively,
[(D'n")T~D~+AZ]~6R,=(D~)T~En, n=0,1,2,3,...,
(23)
[(D'.")Td.l~+AZ]~6Rn=(D~)T~En, n=0,1,2,3,....
(24)
and
Finally, eqs. (4), (15) and (23) or (4), (16) and (24) form the basic computational structure of each iteration of the GPST inversion algorithm for solving the exterior or interior inverse problems of the Helmholtz equation respectively. It is extremely difficult to give a precise comparison of the computational efficiency of various numerical algorithms for solving inverse problems of partial differential equations, because the actual performance of a numerical algorithm to achieve a certain prescribed accuracy depends on many factors, e.g., the intrinsic computational complexity of the numerical algorithm, the compatibility between the numerical algorithm and the architecture of the computer, memory storage, peripheral hardwares and softwares, skill of programming by taking the special features of the computer architure, etc. However, a partial answer to this difficult question can be obtained by eliminating the influences of architecture of the computer, memory storage and programming skill so that only the intrinsic computational complexity of the numerical algorithm will be used to estimate the computational efficiency. Again, this is very much problem-dependent; nevertheless, it will provide some tractable partial answer to the above-mentioned problem. Since the easiest and simplest way to perform a computational complexity analysis for a given numerical algorithm is to count the number of arithmetic floating point operations (FLO) needed to complete a certain computational task, here we shall carry out an asymptotic FL0 count for the GPST inversion algorithm. Here, since the addition and multiplication will take about the same amount of CPU time on a modern computer, the FL0 count will include all types of FLO’s. In each GPST iteration, solving the discretized scattering or cavity problems (15) or (16) respectively by LU-decomposition needs approximately 2ZJ*(J + 2) FLOs. In forming (21) or (22), L back substitutions are needed so that the approximate FL0 count is 4JL*. Approximately 2JL(J+ 1) FLOs are needed to form (23) or (24). Lastly, solving (23) or (24) by LU-decomposition requires ( J3/3 + J*) FLOs. Hence the total asymptotic FL0 count per iteration for the GPST inversion algorithm applied to two-dimensional exterior and interior inverse problems is approximately, (2Z+1/3)J3+(4Z+2L+1)J2+2L(2L+1)J.
(25)
In the case of L= J, (25) has a simpler form, (2Z+19/3)J3+(4Z+3)J2.
(26)
For the corresponding three-dimensional inverse problems with grid points Xi,j,mz (ri, 0j, &,), i = M, and the boundary surface S, = R(8, c$), the above, Z, j = 1,2,3, . . . , J, m = 1,2,3, . . . , 1,2,3,. . . described GPST inversion algorithm will work in the same manner except that the computational effort will be greatly increased. This can be seen from the approximate total asymptotic FL0 count per iteration for the GPST inversion algorithm, (2Z+1/3)J3M3+(4Z+2L+1)J2M2+2L(2L+1)JM,
(27)
and similarly, if L= JM, eq. (27) becomes (2Z+19/3)J3M3+(4Z+3)J2M2.
(28)
S.L.
Wang,
Y.M.
Chen / Exterior
and interior
393
inverse problems
4. Numerical simulation In order to test the feasibility for solving
the exterior
and to study the general
and interior
boundary-shape
characteristics inverse
problems
of the GPST inversion of two-dimensional
equation without real measurement data, the following numerical simulation First, one chooses R*( 0) to represent the correct boundary of a scatterer
algorithm Helmholtz
procedure is carried out. or cavity. Then the exterior
or interior boundary value problem (1) or (2) with S, = R*( 0) is solved by using the same finite difference scheme described in the previous section to generate the supposedly measured or design data u(x,) = h(x,), 1= 1,2,3, . . . ) L. Next, &,(B) is chosen to start the GPST inversion algorithm. Hence upon solving eqs. (4), (15) and (23) or eqs. (4), (16) and (24) numerically, similar
manner.
truncation
One continues
and round-off
errors
this procedure
until
in both generating
RI(e)
a numerical the numerical
is obtained; limit
RN(e)
Rz(B) can be obtained is reached.
in a
Other than the R*(O), any norm
data and computing
of R*( 0) -RN(O) can be used as a criterion for evaluating the performance of the GPST inversion algorithm. As a matter of fact, due to the ever presence of these truncation and round-off errors, one can consider them to be equivalent to some kind of pseudo-noise in the numerically generated measurement data.
POINT
SOURCE AT
Fig. 2. Comparison
of the calculated
shape . . . . and the exact Saturn-like scatterer source at 0 = 0, and full-aperture measurement.
with the initial
R = 9
guess - - - -, one point
Due to economical reasons, the numerical simulation here is carried out only for the examples with k = 0.1 and ka in the range from 0.2 to 0.6 (with “a” being the typical .dimension of scatterer or cavity). Otherwise, more computational grid points are needed to provide sufficient accuracy in solving the boundary value problems (6) and (8) with larger ku. For all examples of the exterior inverse problem, the total computational grid points is 240 (I = 10, J = 24); for all examples of the interior inverse problem, the total computational grid points is 128 (I = 8, J = 16). Moreover, the number of iterations is five and the Tikhonov regularization parameter A is set to 10m3 for all examples. For the exterior inverse problem, the first scatterer is a non-convex Saturn-like object. To investigate the effect of the angular location of a single point source on reconstructing the shape of this scatterer, numerical simulations are carried out for the point source at (9,0), (9, n/4), and (9, n/2) with the scattered
394
S.L.. Wang, KM.
Chen / Exterior
and interior inverse problems
POINT
Y
SOURCE
AT R=
t
9
2
Fig. 3. Comparison of the calculated shape . . . . and the exact Saturn-like scatterer source at 0 = m/4, and full-aperture measurement.
\tr
POINT
with the initial guess - - - -, one point
SOURCE
AT R=
9
with tthe initial guess - - - -, one point Fig. 4. Comparison of the calculated shape . . . . and the exact Saturn-like scatterer source at 0 = n/2, and full-aperture measurement.
S.L.
Wang,
Y.M.
Chen / Exterior
and interior
395
inverse problems
Y
I’
, ,
)/H--\ ‘. .,
\ \ \
*.. ~-*.,,
\tr
0
tji POINT
. ..a
. ..--
‘, \
.
. .. . /
\
ATR’9
Fig. 5. Comparison
1 . . ... .
SOURCE
\
l. .
/
/
\t/ C._
3 ,3.4 /
x POINT AT
SOURCE R =
9
,/’
of the calculated
with the initial guess - - - -, two point shape . . . . and the exact Saturn-like scatterer sources at 0 = 0 and T, and full-aperture measurement.
POINT AT
Fig. 6. Comparison
41
of the calculated
SOURCE R = 9
with the initial guess - - - -, one point source at shape . . . . and the exact nose cone 0 = 0, and the full-aperture measurement.
wave measured at every j-line (full-aperture problem); their numerical results are plotted in Figs. 2, 3, and 4 respectively. To study the effect of the multiple point sources on the reconstruction, a numerical simulation is carried out for the full-aperture problem with two point sources at (9,0) and (9, m); its numerical result is plotted in Fig. 5. To examine the effect of the limited aperture on the reconstruction, numerical simulations are carried out for a single point source at (9,0) with the scattered wave measured in the ranges -7r/2 c 0 G 1r/2 (back-scattering measurement) and rr/2 s 8 c 31r/2 (forward-scattering measurement); their numerical results are very similar to those in Figs. 3 and 4 respectively. Finally, to investigate the effect of the limited aperture and multiple point sources on the reconstruction, a numerical
396
S.L. Wang, Y.M. Chen / Exterior and interior inverse problems
Fig. 7. Comparison of the calculated shape . . . . and the exact elliptic cavity with the initial guess - - - - inside the cavity and measurement performed at points ‘x’.
/rp.... /
/’ ,/ /
\
/
8 .- _
/--
/
; f i+
. .
\
\
\
*
+ \ :I *.....:
\
Fig. 8. Comparison of the calculated shape . . . . and the exact elliptic cavity with the initial guess - - - - outside the cavity and measurement performed at points ‘x’.
\
+
‘r :.
‘1,
\ I? :,
\..’
\
\
\
--
Fig. 9. Comparison
of the calculated
shape
__
0
/
/
_-I
. . . . and the exact triangular cavity performed
/
with the initial guess - - - - and measurement
at points ‘x’.
simulation is carried out with two point sources at (9,O) and (9, T) and the scattered wave measured in the range -1r/2 c 8 s 1r/2; its numerical result is very similar to that in Fig. 2. The second, third, and fourth scatterers are a nose cone, a flying saucer and a rounded slender rod respectively: The scattered wave is measured at every j-line for these three examples; single point source at (9,0) is used for the second and fourth examples and two point sources at (9,O) and (9,~) are used
S.L. Wang, Y.M. Chen / Exterior and interior inverse problems
Fig. 10. Comparison square cavity-with
of the calculated shape. . . . and the exact the initial guess - - - - and measurement performed at points ‘x’.
397
Fig. 11. Comparison of the calculated shape . . . * and the exact butterfly-like cavity-with initial guess - - - - and measurement performed at points ‘x’.
for the third example. The numerical result for the second example is plotted in Fig. 6 and the numerical results for the examples can be found in [l]. For all examples of the full-aperture exterior inverse problem, the CPU time on VAX 8600 is approximately 3.2 seconds which is approximately equivalent to 15 seconds on VAX 11/750 and 0.04 seconds on CRAY X-MP/4 [19]. For all examples of the limited-aperture exterior inverse problem, the CPU time on VAX 8600 is approximately 2.8 seconds. For the interior inverse problem, four different cavities, an ellipse, a triangle, a square and a buttertly, are chosen as examples. Here for all examples, the source is uniformly distributed on a circle of radius 0.5, i.e., f( r, 6) = 1, 0 d r C 0.5, 0 s 0 d 27r, and f( r, 6) = 0 otherwise, and the scattered wave is measured at every j-line. To examine the effect of the sizes of initial guess on the reconstruction, both the interior and exterior circles are used as the initial guesses for the elliptic cavity; their numerical results are plotted in Figs. 7 and 8 respectively. The numerical simulation results for the triangular, square, and butterfly-like cavities are plotted in Figs. 9-11 respectively. For all these examples, the CPU time on VAX 8600 is approximately 1.2 seconds.
5. Discussion
From the results of numerical simulation, it is clear that the GPST inversion algorithm with an adaptive grid system does give very good results in reconstructing the shape, size and orientation of a scatterer or cavity. However, for solving the exterior inverse problem, it does bristle with similar difficulties as those of other inversion methods. For the Saturn-like object of the exterior inverse problem, the results in Figs. 2-4 show that for a single source with the full aperture measurement, the reconstruction is relatively poor in the back or shadow side regardless the direction of incident wave; this is probably due to the fact that our calculation done on a relatively uniform grid system is less accurate in the shadow region. The result in Fig. 5 shows that the addition of another source in the opposite direction of the original source makes
398
S.L. Wang, KM. Chen / Exterior and interior inverse problems
the reconstruction much better, because there is no shadow region here. The numerical result for a single source with half aperture back-scattering measurement, the reconstruction is good in the lit side and poor in the shadow side; moreover, the reconstruction here is not as good as its corresponding full aperture case (Fig. 2). The obvious explanation for this is that the measurement data for the half aperture case are less sufficient than those for the full aperture case; as a matter of fact, for the case of half aperture eq. (21) is an under determined system and hence eq. (23) becomes more ill-conditioned. The numerical result for a single source with half aperture forward-scattering measurement, the reconstruction is uniformly poor. Even with the addition of another source, the numerical result shows that it will not help the case of half aperture back-scattering measurement very much. In view of the above discussion, the results of numerical simulation for the nose cone, flying saucer and rounded slender rod provide no surprise. For the examples of the interior inverse problem, the results of numerical simulation in Figs. 7-l 1 are also very good. It has been found that the limited aperture problem of the exterior inverse problems has very little effect on the performance of GPST for solving the interior inverse problems; moreover, whether the initial guess is inside or outside of the correct cavity does not seem to make much difference in reconstruction (Figs. 7 and 8). So far we have not made any comparison between the GPST inversion algorithm and other methods for solving the interior inverse problem, because we have not been able to find any related research in the literature yet. The relative efficiency of each inversion algorithm can be estimated by comparing its CPU time on a computer with the CPU times of others on the same computer for solving a similar class of examples; when different computers are used, the equivalent CPU time on any other computer can be estimated from the actual CPU time on a given computer with reasonable accuracy [ 191. GPST inversion algorithm with an adaptive grid system is found to be more efficient than other existing methods, even though GPST inversion algorithm itself is a very general numerical algorithm. Unfortunately, the present version of GPST inversion algorithm still is not efficient enough for many real-time applications. For example, the asymptotic FL0 count needed for solving the corresponding three-dimensional exterior inverse problem with total grid points 2880 (I = 10, J = 24, M = 12) can be estimated from (27) or (28) and the estimated CPU times on VAX 8600 and CRAY X-MP/4 are 6 x 10’ s. and 60 s. respectively [19]. Efficiency improvement of GPST inversion algorithm has been in the making at all levels. For examples, at the matrix level, the development of an efficient preconditioned conjugate gradient method for solving (15), (16), (23) and (24) has been underway; at other levels of the hierarchy of GPST, methods based upon the principle of the nearest neighborhood dominance, hierarchical multigrid strategy, and hierarchical domain-decomposition have been in the process of development. Their results will be reported in the future.
Acknowledgement
This research was supported in part by the National Science Foundation under grant no. DMS-8811331.
References [l]
Y.M. Chen and S.L. Wang, “An efficient numerical method for determination of shapes, sizes and orientations of flaws for Review of Progress in Quantitative Nondestructive Evaluation 4, eds. D.O. Thompson and D.E. nondestructive evaluation”, Chimenti; Plenum Press, New York, 543-549 (1984). [Z] S.L. Wang and Y.M. Chen, “GPST algorithm with adaptive grid for solving exterior and interior inverse problems of Helmholtz equation”, Proc. 12th IMACS World Congress on Scientifc Computation 2, Paris, 411-413 (1988).
S.L. Wang, V.M. Chen / Exterior and interior
inverse problems
399
bistatic scattering”, Canad. J. Phys. 47, [3] V.H. Weston and W.M. Boerner, “An inverse scattering technique for electromagnetic 1177-1184 (1969). IEEE Trans. Ant. & fiop. AP-18, inverse scattering problem”, [4] W.A. Imbriale and R. Mittra, “The two-dimensional 633-642 (1970). [5] A. Roger, “Newton-Kantorovich algorithm applied to an electromagnetic inverse problem”, IEEE Trans. Ant. & Prop. AP-29, 232-238 (1981). [6] D. Colton, “The inverse scattering problem for time-harmonic acoustic waves”, SIAM Reu. 26, 323-350 (1984). [7] D. Colton and P. Monk, “A novel method for solving the inverse scattering problem for time-harmonic acoustic waves in the resonance region”, SIAM J. Appl. Math. 45, 1039-1053 (1985). [S] G. Kristenson and C.R. Vogel, “Inverse problems for acoustic waves using the penalized likelihood method”, Inoerse Problems 2,461-479 (1986). [9] D. Colton and P. Monk, “The numerical solution of the three-dimensional inverse scattering problem for time-harmonic acoustic waves”, SIAM J. Sci. Stat. Comput. 8, 278-291 (1987). Boundary Elements IX Vol. 3, Fluid Flow [lo] A. Kirsch and R. Kress, “An optimization method in inverse acoustic scattering”, and Potential Applications, eds. C.A. Brebbia, W.L. Wendland and G. Kuhn, Springer, Berlin, 3-18 (1987). [l l] A. Kirsch, R. Kress, P. Monk and A. Zinn, “Two methods for solving the inverse acoustic scattering problem”, Inverse Problems 4, 749-770 (1988). [12] A. Zinn, “On an optimisation method for the full- and the limited-aperture problem in inverse acoustic scattering for a sound-soft obstacle”, Inverse Problems 5, 239-253 (1989). [13] A.N. Tikhonov and V.Y. Arsenin, Solution of III-Posed Problems, Wiley, New York (1977). [14] M. Kac, “Can one hear the shape of a drum ?” Amer. Math. Monthly 73, l-23 (1966). [15] M.H. Protter, “Can one hear the shape of a drum? revisited,” SIAM Rev. 29, 185-197 (1987). [16] M.M. Sondhi and J.R. Resnick, “The inverse problem for the vocal tract: numerical methods, acoustical experiments, and speech synthesis”, J. Acoust. Sot. Amer. 73, 985-1002 (1983). [17] D.S. Tsien and Y.M. Chen, “A numerical method for nonlinear inverse problems in fluid dynamics,” Computational Methods in Nonlinear Mechanics, Univ. of Texas, Austin, 935-943 (1974). [18] A. Bayliss, M. Gunzburger and E. Turkel, “Boundary conditions for the numerical solution of elliptic equations in the exterior region”, SIAM J. Appl. Math. 42, 707-725 (1982). [ 191 J.J. Dongarra, “Performance of various computers using standard linear equations software in a Fortran environment”, Tech. Mem. No. 23, Argonne National Laboratory (1988).
WAVE MOTION 13 (1991) 387-399 ELSEVIER
AN EFFICIENT NUMERICAL METHOD INVERSE PROBLEMS OF HELMHOLTZ
FOR EXTERIOR EQUATION
AND INTERIOR
S.L. WANG and Y.M. CHEN Department of Applied Mathematics and Statistics, State Universiiy of New York at Stony Brook, Stony Brook, NY 117943600, USA
Received 17 April 1990, Revised 27 July 1990
The generalized pulse-spectrum technique (GPST), an efficient and versatile iterative inversion algorithm, is used with adaptive grids to solve both exterior (scattering) and interior (cavity) boundary-shape inverse problems of two-dimensional Helmholtx equation. Numerical simulations of nontrivial examples are carried out to test the feasibility and to study the general characteristics of GPST without the real measurement data. It is found that GPST does efficiently produce very good resuits.
1. Introduction The exterior inverse problems can be exemplified by the inverse scattering problems in which the shape, size and orientation of a scatterer are determined from the knowledge of incident and scattered waves in the exterior of the scatterer. The interior inverse problems can be exemplified by the inverse closed cavity problems in which the shape, size and orientation of the closed cavity are determined from the knowledge of incident and scattered waves in the interior of the cavity. In general, the knowledge or known data of the incident and scattered waves can come from a mixture of different locations and time or frequency intervals, e.g., at (x,, t,,,) or (x,, w,), I= 1,2,3,. . . , L, m = 1,2,3,. . . , A4. Both the exterior and interior inverse problems are difficult to solve due to the fact that they often are both ill-posed and nonlinear. Here only the shapes, sizes and orientations of the scatterers and cavities in the exterior and interior inverse problems of Helmholtz equation respectively are treated, and it is our intention to consider this paper as a detailed presentation of Refs. [l] and [2]. The inverse scattering problems of time-harmonic scalar waves with the known data came from a single frequency and many spatial locations lead to the inverse scattering problems of the Helmholtz equation. For the inverse scattering problems of the Helmholtz equation in the low to resonance frequency region, many methods are available, e.g., method of analytic continuation, explicit optimization methods, implicit optimization methods, etc. The method of analytic continuation [3] and [4], relied on the numerical analytic continuation, subjects to severe instability from the ill-posedness. In the explicit optimization methods [SJ, [6] (an excellent review) - [12], the inverse scattering problem is explicitly formulated as a nonlinear optimization problem first, then a nonlinear optimization technique with regularization to control instability is used to solve the problem; the regularization can either be applied to the full nonlinear problem or to each linear sub-problem arising from the linearization in the nonlinear optimization technique. In the implicit optimization methods, e.g., the generalized pulse-spectrum technique (GPST) [l] and [2], based on a combination of a Newton-like iterative method and the Tikhonov regularization method [13], is a class of versatile and efficient iterative algorithm, for solving multi-parameter inverse 01658641/91/$03.50
@ 1991 - Elsevier Science Publishers B.V.
388
S. L. Wang, Y. M. Chen / Exterior and interior inverse problems
problems of a system of nonlinear partial can be performed in either the space-time
differential equations. In GPST, the parameter identification (pulse) domain or the space-frequency (spectrum) domain or
both. There the inverse
is formulated
differential
equations
For the inverse
scattering instead
problem
of a nonlinear
closed cavity problems
optimization
as a perturbation
or iteration
problem
of partial
problem.
of scalar waves, there is no existing
research
can be found except
the well known
Kac’s inverse eigenvalue problem, “can one hear the shape of a drum” [ 14]? Their main difference is that the known data of the inverse cavity problem here are the incident and scattered waves in the interior of cavity and the known data of Kac’s inverse problem are all the normal modes of vibration. While there is only a partial analytic solution to Kac’s problem [ 151, GPST has been used to solve simple
inverse cavity problems cavity
problems
numerically
here are quite
with some success [2]. It should be pointed
different
from the inverse
open
out that the inverse
cavity problems
closed
such as the vocal tract
inverse problem [ 161. Since the incipience of GPST in 1974 [17], it has been generalized, improved and applied to many practical inverse problems. If only the shapes, sizes and orientations of the cavities in the exterior and interior inverse problems of Helmholtz equation respectively need to be determined, which is the case here, then one can make GPST dedicated to these special types of simple inverse problems by incorporating an adaptive body-fitted grid system. In this way, the computational efficiency of GPST is greatly improved. First, the dedicated GPST inversion algorithm for these types of inverse problems is presented. Then various computational aspects of GPST are described. Next, numerical simulations of simple but non-trivial examples are carried out on the adaptive grid system to test the feasibility and to study the general characteristics of this dedicated GPST without the real measurement data. Finally, a comprehensive discussion of the numerical results and the future improvement on GPST to attend higher efficiency is given.
2. Generalized pulse-spectrum For simplicity,
the following
in E* are cosidered (v’+
technique (GPST)
k2)U =f(x),
r”*(au/ar
exterior
and interior
boundary
value problems
of the Helmholtz
equation
respectively,
-iku)
x=(r,B)EE*-.f2,,
+ 0,
u(x)=O,
as r + 00, (radiation
xES,=R(e),
(1)
condition)
and (V’+k*)u=f(x),
XE&,
xeS,=R(O),
z&)=0,
(2)
where E* is the two-dimensional Euclidean space, 0, is a star-shaped finite sub-domain in E* with the boundary S, , 0, is a sub-domain in E * surrounding n, and has a circular outer boundary S, , R, is also a sub-domain in E* surrounding fin2 and has a circular outer boundary S, concentric with S2, x = 0 is located in 0, (Fig. l), and the source f(x) has a compact support in either E*-0, or 0,. Here the exterior and interior inverse problems are to determine R(B) or S, from the measurement data or prescribed design data u(x) = h(x),
x = q E a3
and 0,)
respectively,
1=1,2,3
,...,
L,
(3)
such that u(x) satisfies (1) and (2) respectively. The GPST inversion algorithm begins by setting u n+, = u, +&l,,
R ,,+,=R,+SR,,
n=0,1,2,3
,...,
(4)
389
S. L. Wang, Y.M. Chen / Exterior and interior inverse problems
Fig. 1. Scheme of adaptive computational
grid system.
where R,( 0) is the initial guess for R( 0) and S is a small number such that the b-terms are smaller than their corresponding non-&-terms in some norms for numerical convergence. For the exterior inverse problem, upon substituting (4) into (1) with u and R replaced by u,+, and R ?I+17 respectively, one obtains (V’+k*)(u,+Su,)=f(x),
XEE*-~,.,,
n=0,1,2,3
,...,
as r+oo.
r”*[(~u,/~r-iku,)+(&3u,/~r-ik&,)]+O,
(5)
Upon neglecting terms of O( S*) and higher, the system for u, same as (1) is obtained, (V*+k*)u,
=f(x),
XE E*-R,,,,
r”‘(au,/ar-iku,)+O,
u,=O,
XER,,
n = 0, 1,2,3, . . . ,
(6)
as r+m,
and a similar system for au,, (V*+k*)Su, =0,
XE E*-a,,,,
&I, = -(au,/ar)GR,,
XE R,,
n =O, 1,2,3,
(7)
as r+m.
r”*(a&/ar-iksu,)+O,
*. ,
Similarly, for the interior inverse problem one obtains the system for u, same as (2), XE n,,,,
(V’+ k*)u, =.0x),
u,=O,
XER,,
n = 0, 1,2,3, . . . ,
(8)
and a similar system for au,,, (V’+ k*)Su, =0,
XE f‘&,,
au,, = -(au,/ar)GR,,
XE R,,
n =0,1,2,3,.
...
(9)
By using the method of Green’s function, from (7) or (9) one obtains an integral relation between SR, and au,,
I R.
[aG,(x, x’)/av;][au,/ar’]t%R,(
0’) dl:, = SU,,
n = 0, 1,2,3,.
..,
(10)
S.L. Wang, Y.M. Chen / Exterior and interior inverse problems
390
where G,(x, x’) is the Green’s function of (7) or (9), v,, is the unit normal vector on R, pointing toward the domain of interest, and dl, is the incremental line segment of R,. For speeding up the convergence and eliminating the unknown u,+], in (10) one can set x to q, 1= 1,2,3,. . . , L, and replace u,+,(x,) by the known data h(x,), I = 1,2,3, . . . , L. Hence a system of Fredholm integral equations of the first kind for the unknown SR, (0) is obtained as
I
[aG,(x,, x’)/av~][au,/ar’]sR,(e’)
dl:, = h(x,) - u,(x,),
n =0, 1,2,3,.
. . , 1= 1,2,3,.
. . , L.
(11)
a,,
Now eqs. (4), (6) and (11) and eqs. (4), (8) and (11) form the basic theoretical structure of each iteration of GPST for solving the exterior and interior inverse problems of the Helmholtz equation respectively. Note that with some minor modifications, GPST can be used with equal efficacy to solve both exterior and interior inverse problems of the Helmholtz equation with Neumann or impedance boundary conditions and with measurement data containing derivatives of u(x).
3. Computational
consideration
Theoretically, the efficiency of a numerical algorithm depends on the intrinsic mathematical structure of the algorithm, but in practice it depends even more on how efficiently one can implement every minute computational step which is very much computer-architectureand compiler-dependent. Here an adaptive irregular polar computational grid system (Fig. 1) is used. In this grid system, the circular finite computational domain consists of three connected but non-overlapping sub-domains, 0, , L& and L$, where 0, represents the scatterer or cavity, the interface S, between 0, and L!* represents the shapes, sizes and orientations of the scatterer or cavity, and LZ2and 0, represent the region in the exterior of the scatterer; let the grid point be denoted by (ri, 6,) = (i, j). In 0, and L&, the locations of the grid points are re-adjusted along fixed constant j-lines after each iteration to accommodate the shape, size and orientation of the scatterer or cavity respectively; in this way, one can achieve high accuracy with a small number of grid points; the boundary of the scatterer or cavity S, should be a constant i-line. In L&, the grid system is a fixed circular cylindrical grid to avoid the unnecessary computational effort in generating grid points; however, due to the truncation of the infinite space domain, the radiation condition in (6) and (7) is replaced by the approximate radiation condition [ 181, dv/dr + (2r)-‘v
- ikv = 0,
XE&, n=0,1,2,3
v = u,, au,,
,....
(12)
To deal with the irregular incremental quadrilaterals, one considers the area integration of the Helmholtz equation on a typical quadrilateral Ai,j,, centered at (i, j) with (i - l/2, j - l/2), (i - l/2, j + l/2), (i + l/2, j-1/2) and (i+ l/2, j+ l/2) as the four vertices. From (6) or (8) and Green’s identity, the area integral becomes
I L.” I
au,/dv,
dl,, + k2 I,,,.Iu”du~=la,j”~fdu.. ..
(13)
n=O,1,2,3,..., .*
where Li,j,, is the boundary of Ai,j,, and V, is the unit vector normal to Li,j,,. Upon using the rectangle rule and the center difference scheme to approximate the integrals and the normal derivative respectively, one obtains a 9-point finite difference formula for u,(x), ifl
1
i+l
1
a-i-1 p=i-l
YmaU”,m,P=f;,j, i=l,2,3
,...,
Z,j=l,2,3
,...,
J,n=0,1,2,3
,....
(14)
S.L. Wang, Y.M. Chen / Exterior and interior inverse problems Hence
the above discretizations
391
of (6) with (12) and of (8) lead respectively to
B ‘,“.U,,=F,
n=0,1,2,3
,...,
(15)
Bz-U.=F,
n=0,1,2,3
,...,
(16)
and
where By and Bz both are known N x N symmetric sparse matrices with half bandwidth J, N = I X J, U, is the N x 1 matrix with all u,,~,~, i = 1,2,3, . . . , I, j = 1,2,3, . . . , J, as its components, and F is the N x 1 matrix with all _&, i = 1,2,3,. . . , Z, j = 1,2,3,. . . , J, as its components. Here, eqs. (15) and (16) can be solved by using an efficient LU-decomposition sparse matrix solver. Straightforward discretization and evaluation of the system of Fredholm integral equations of the first kind (11) and calculating G,(x, x’) require very large computational effort. A more efficient way is to consider the discretization of (7) and (9) by using the same 9-point discretization schemes for (6) and (12) and for (8) with the boundary values of Su, as equivalent sources. One obtains the similar linear algebraic systems respectively for Su,, B~-SUn=(Cn-SR,,O)T,
n=0,1,2,3
,...,
(17)
B~~SUn=(Cn*6R,,,0)T,
n=0,1,2,3
,...,
(18)
and
where MY,,is the N x 1 matrix with all a~,,~,~,i = 1,2,3, . . . , I, j = 1,2,3, . . . , J, as its components, 6R, is the J x 1 matrix with all cSR,~ = SR,(Bj), j = 1,2,3, . . . , J, as its components, and C,, is the known J x J diagonal matrix obtained from the boundary condition (a known relation between au, and SR, at the boundary R,). Equations (17) and (18) lead respectively to (B~)-‘~(Cn~SRn,0)T=8Un,
n=0,1,2,3
,...,
(19)
(B~)-‘*(Cn-SRn,O)T=SUn,
n=0,1,2,3
,....
(20)
and
Upon collecting all of the equations in (19) and (20) corresponding to the data points g, I = 1,2,3, . . . , L, only and replacing u,+,(x,)-u,(x,) by h(x,)-u,(x,), 1=1,2,3,. .., L, one obtains the compact linear algebraic systems respectively for the unknown SR,, 0:.
SR, = En,
n=0,1,2,3
,...,
(21)
n=0,1,2,3
,...,
(22)
and DF-SR,=E,,
where 0: and Df are Lx J full matrices and E, is a Lx 1 matrix with h(x,) - u.(x,), 1= 1,2,3,. . . , L, as its components. It is important to point out that in forming (21) or (22) only L rows of (B’,“)-’ or (B’,“)-’ corresponding to the data points q, 1= 1,2,3,. . . , L, are needed. Since the LU-decomposition of B’,” and B’,”have been already performed in solving (15) and (16), they are saved so that the L rows of (By)-’ or (Br)-’ can be obtained by performing back substitution L times. Hence a great saving in computational effort is achieved. It is also clear that (21) and (22) are the most compact forms of the discretization of the system of Fredholm integral equations of the first kind (11).
392
S. L. Wang, Y. M. Chen / Exterior and interior inverse problems
Here D’,” and Df are rectangular and very ill-conditioned matrices. Then the well-known Tikhonov regularization method is used to overcome this difficulty, i.e., instead of solving (21) and (22), one solves their regularized systems respectively,
[(D'n")T~D~+AZ]~6R,=(D~)T~En, n=0,1,2,3,...,
(23)
[(D'.")Td.l~+AZ]~6Rn=(D~)T~En, n=0,1,2,3,....
(24)
and
Finally, eqs. (4), (15) and (23) or (4), (16) and (24) form the basic computational structure of each iteration of the GPST inversion algorithm for solving the exterior or interior inverse problems of the Helmholtz equation respectively. It is extremely difficult to give a precise comparison of the computational efficiency of various numerical algorithms for solving inverse problems of partial differential equations, because the actual performance of a numerical algorithm to achieve a certain prescribed accuracy depends on many factors, e.g., the intrinsic computational complexity of the numerical algorithm, the compatibility between the numerical algorithm and the architecture of the computer, memory storage, peripheral hardwares and softwares, skill of programming by taking the special features of the computer architure, etc. However, a partial answer to this difficult question can be obtained by eliminating the influences of architecture of the computer, memory storage and programming skill so that only the intrinsic computational complexity of the numerical algorithm will be used to estimate the computational efficiency. Again, this is very much problem-dependent; nevertheless, it will provide some tractable partial answer to the above-mentioned problem. Since the easiest and simplest way to perform a computational complexity analysis for a given numerical algorithm is to count the number of arithmetic floating point operations (FLO) needed to complete a certain computational task, here we shall carry out an asymptotic FL0 count for the GPST inversion algorithm. Here, since the addition and multiplication will take about the same amount of CPU time on a modern computer, the FL0 count will include all types of FLO’s. In each GPST iteration, solving the discretized scattering or cavity problems (15) or (16) respectively by LU-decomposition needs approximately 2ZJ*(J + 2) FLOs. In forming (21) or (22), L back substitutions are needed so that the approximate FL0 count is 4JL*. Approximately 2JL(J+ 1) FLOs are needed to form (23) or (24). Lastly, solving (23) or (24) by LU-decomposition requires ( J3/3 + J*) FLOs. Hence the total asymptotic FL0 count per iteration for the GPST inversion algorithm applied to two-dimensional exterior and interior inverse problems is approximately, (2Z+1/3)J3+(4Z+2L+1)J2+2L(2L+1)J.
(25)
In the case of L= J, (25) has a simpler form, (2Z+19/3)J3+(4Z+3)J2.
(26)
For the corresponding three-dimensional inverse problems with grid points Xi,j,mz (ri, 0j, &,), i = M, and the boundary surface S, = R(8, c$), the above, Z, j = 1,2,3, . . . , J, m = 1,2,3, . . . , 1,2,3,. . . described GPST inversion algorithm will work in the same manner except that the computational effort will be greatly increased. This can be seen from the approximate total asymptotic FL0 count per iteration for the GPST inversion algorithm, (2Z+1/3)J3M3+(4Z+2L+1)J2M2+2L(2L+1)JM,
(27)
and similarly, if L= JM, eq. (27) becomes (2Z+19/3)J3M3+(4Z+3)J2M2.
(28)
S.L.
Wang,
Y.M.
Chen / Exterior
and interior
393
inverse problems
4. Numerical simulation In order to test the feasibility for solving
the exterior
and to study the general
and interior
boundary-shape
characteristics inverse
problems
of the GPST inversion of two-dimensional
equation without real measurement data, the following numerical simulation First, one chooses R*( 0) to represent the correct boundary of a scatterer
algorithm Helmholtz
procedure is carried out. or cavity. Then the exterior
or interior boundary value problem (1) or (2) with S, = R*( 0) is solved by using the same finite difference scheme described in the previous section to generate the supposedly measured or design data u(x,) = h(x,), 1= 1,2,3, . . . ) L. Next, &,(B) is chosen to start the GPST inversion algorithm. Hence upon solving eqs. (4), (15) and (23) or eqs. (4), (16) and (24) numerically, similar
manner.
truncation
One continues
and round-off
errors
this procedure
until
in both generating
RI(e)
a numerical the numerical
is obtained; limit
RN(e)
Rz(B) can be obtained is reached.
in a
Other than the R*(O), any norm
data and computing
of R*( 0) -RN(O) can be used as a criterion for evaluating the performance of the GPST inversion algorithm. As a matter of fact, due to the ever presence of these truncation and round-off errors, one can consider them to be equivalent to some kind of pseudo-noise in the numerically generated measurement data.
POINT
SOURCE AT
Fig. 2. Comparison
of the calculated
shape . . . . and the exact Saturn-like scatterer source at 0 = 0, and full-aperture measurement.
with the initial
R = 9
guess - - - -, one point
Due to economical reasons, the numerical simulation here is carried out only for the examples with k = 0.1 and ka in the range from 0.2 to 0.6 (with “a” being the typical .dimension of scatterer or cavity). Otherwise, more computational grid points are needed to provide sufficient accuracy in solving the boundary value problems (6) and (8) with larger ku. For all examples of the exterior inverse problem, the total computational grid points is 240 (I = 10, J = 24); for all examples of the interior inverse problem, the total computational grid points is 128 (I = 8, J = 16). Moreover, the number of iterations is five and the Tikhonov regularization parameter A is set to 10m3 for all examples. For the exterior inverse problem, the first scatterer is a non-convex Saturn-like object. To investigate the effect of the angular location of a single point source on reconstructing the shape of this scatterer, numerical simulations are carried out for the point source at (9,0), (9, n/4), and (9, n/2) with the scattered
394
S.L.. Wang, KM.
Chen / Exterior
and interior inverse problems
POINT
Y
SOURCE
AT R=
t
9
2
Fig. 3. Comparison of the calculated shape . . . . and the exact Saturn-like scatterer source at 0 = m/4, and full-aperture measurement.
\tr
POINT
with the initial guess - - - -, one point
SOURCE
AT R=
9
with tthe initial guess - - - -, one point Fig. 4. Comparison of the calculated shape . . . . and the exact Saturn-like scatterer source at 0 = n/2, and full-aperture measurement.
S.L.
Wang,
Y.M.
Chen / Exterior
and interior
395
inverse problems
Y
I’
, ,
)/H--\ ‘. .,
\ \ \
*.. ~-*.,,
\tr
0
tji POINT
. ..a
. ..--
‘, \
.
. .. . /
\
ATR’9
Fig. 5. Comparison
1 . . ... .
SOURCE
\
l. .
/
/
\t/ C._
3 ,3.4 /
x POINT AT
SOURCE R =
9
,/’
of the calculated
with the initial guess - - - -, two point shape . . . . and the exact Saturn-like scatterer sources at 0 = 0 and T, and full-aperture measurement.
POINT AT
Fig. 6. Comparison
41
of the calculated
SOURCE R = 9
with the initial guess - - - -, one point source at shape . . . . and the exact nose cone 0 = 0, and the full-aperture measurement.
wave measured at every j-line (full-aperture problem); their numerical results are plotted in Figs. 2, 3, and 4 respectively. To study the effect of the multiple point sources on the reconstruction, a numerical simulation is carried out for the full-aperture problem with two point sources at (9,0) and (9, m); its numerical result is plotted in Fig. 5. To examine the effect of the limited aperture on the reconstruction, numerical simulations are carried out for a single point source at (9,0) with the scattered wave measured in the ranges -7r/2 c 0 G 1r/2 (back-scattering measurement) and rr/2 s 8 c 31r/2 (forward-scattering measurement); their numerical results are very similar to those in Figs. 3 and 4 respectively. Finally, to investigate the effect of the limited aperture and multiple point sources on the reconstruction, a numerical
396
S.L. Wang, Y.M. Chen / Exterior and interior inverse problems
Fig. 7. Comparison of the calculated shape . . . . and the exact elliptic cavity with the initial guess - - - - inside the cavity and measurement performed at points ‘x’.
/rp.... /
/’ ,/ /
\
/
8 .- _
/--
/
; f i+
. .
\
\
\
*
+ \ :I *.....:
\
Fig. 8. Comparison of the calculated shape . . . . and the exact elliptic cavity with the initial guess - - - - outside the cavity and measurement performed at points ‘x’.
\
+
‘r :.
‘1,
\ I? :,
\..’
\
\
\
--
Fig. 9. Comparison
of the calculated
shape
__
0
/
/
_-I
. . . . and the exact triangular cavity performed
/
with the initial guess - - - - and measurement
at points ‘x’.
simulation is carried out with two point sources at (9,O) and (9, T) and the scattered wave measured in the range -1r/2 c 8 s 1r/2; its numerical result is very similar to that in Fig. 2. The second, third, and fourth scatterers are a nose cone, a flying saucer and a rounded slender rod respectively: The scattered wave is measured at every j-line for these three examples; single point source at (9,0) is used for the second and fourth examples and two point sources at (9,O) and (9,~) are used
S.L. Wang, Y.M. Chen / Exterior and interior inverse problems
Fig. 10. Comparison square cavity-with
of the calculated shape. . . . and the exact the initial guess - - - - and measurement performed at points ‘x’.
397
Fig. 11. Comparison of the calculated shape . . . * and the exact butterfly-like cavity-with initial guess - - - - and measurement performed at points ‘x’.
for the third example. The numerical result for the second example is plotted in Fig. 6 and the numerical results for the examples can be found in [l]. For all examples of the full-aperture exterior inverse problem, the CPU time on VAX 8600 is approximately 3.2 seconds which is approximately equivalent to 15 seconds on VAX 11/750 and 0.04 seconds on CRAY X-MP/4 [19]. For all examples of the limited-aperture exterior inverse problem, the CPU time on VAX 8600 is approximately 2.8 seconds. For the interior inverse problem, four different cavities, an ellipse, a triangle, a square and a buttertly, are chosen as examples. Here for all examples, the source is uniformly distributed on a circle of radius 0.5, i.e., f( r, 6) = 1, 0 d r C 0.5, 0 s 0 d 27r, and f( r, 6) = 0 otherwise, and the scattered wave is measured at every j-line. To examine the effect of the sizes of initial guess on the reconstruction, both the interior and exterior circles are used as the initial guesses for the elliptic cavity; their numerical results are plotted in Figs. 7 and 8 respectively. The numerical simulation results for the triangular, square, and butterfly-like cavities are plotted in Figs. 9-11 respectively. For all these examples, the CPU time on VAX 8600 is approximately 1.2 seconds.
5. Discussion
From the results of numerical simulation, it is clear that the GPST inversion algorithm with an adaptive grid system does give very good results in reconstructing the shape, size and orientation of a scatterer or cavity. However, for solving the exterior inverse problem, it does bristle with similar difficulties as those of other inversion methods. For the Saturn-like object of the exterior inverse problem, the results in Figs. 2-4 show that for a single source with the full aperture measurement, the reconstruction is relatively poor in the back or shadow side regardless the direction of incident wave; this is probably due to the fact that our calculation done on a relatively uniform grid system is less accurate in the shadow region. The result in Fig. 5 shows that the addition of another source in the opposite direction of the original source makes
398
S.L. Wang, KM. Chen / Exterior and interior inverse problems
the reconstruction much better, because there is no shadow region here. The numerical result for a single source with half aperture back-scattering measurement, the reconstruction is good in the lit side and poor in the shadow side; moreover, the reconstruction here is not as good as its corresponding full aperture case (Fig. 2). The obvious explanation for this is that the measurement data for the half aperture case are less sufficient than those for the full aperture case; as a matter of fact, for the case of half aperture eq. (21) is an under determined system and hence eq. (23) becomes more ill-conditioned. The numerical result for a single source with half aperture forward-scattering measurement, the reconstruction is uniformly poor. Even with the addition of another source, the numerical result shows that it will not help the case of half aperture back-scattering measurement very much. In view of the above discussion, the results of numerical simulation for the nose cone, flying saucer and rounded slender rod provide no surprise. For the examples of the interior inverse problem, the results of numerical simulation in Figs. 7-l 1 are also very good. It has been found that the limited aperture problem of the exterior inverse problems has very little effect on the performance of GPST for solving the interior inverse problems; moreover, whether the initial guess is inside or outside of the correct cavity does not seem to make much difference in reconstruction (Figs. 7 and 8). So far we have not made any comparison between the GPST inversion algorithm and other methods for solving the interior inverse problem, because we have not been able to find any related research in the literature yet. The relative efficiency of each inversion algorithm can be estimated by comparing its CPU time on a computer with the CPU times of others on the same computer for solving a similar class of examples; when different computers are used, the equivalent CPU time on any other computer can be estimated from the actual CPU time on a given computer with reasonable accuracy [ 191. GPST inversion algorithm with an adaptive grid system is found to be more efficient than other existing methods, even though GPST inversion algorithm itself is a very general numerical algorithm. Unfortunately, the present version of GPST inversion algorithm still is not efficient enough for many real-time applications. For example, the asymptotic FL0 count needed for solving the corresponding three-dimensional exterior inverse problem with total grid points 2880 (I = 10, J = 24, M = 12) can be estimated from (27) or (28) and the estimated CPU times on VAX 8600 and CRAY X-MP/4 are 6 x 10’ s. and 60 s. respectively [19]. Efficiency improvement of GPST inversion algorithm has been in the making at all levels. For examples, at the matrix level, the development of an efficient preconditioned conjugate gradient method for solving (15), (16), (23) and (24) has been underway; at other levels of the hierarchy of GPST, methods based upon the principle of the nearest neighborhood dominance, hierarchical multigrid strategy, and hierarchical domain-decomposition have been in the process of development. Their results will be reported in the future.
Acknowledgement
This research was supported in part by the National Science Foundation under grant no. DMS-8811331.
References [l]
Y.M. Chen and S.L. Wang, “An efficient numerical method for determination of shapes, sizes and orientations of flaws for Review of Progress in Quantitative Nondestructive Evaluation 4, eds. D.O. Thompson and D.E. nondestructive evaluation”, Chimenti; Plenum Press, New York, 543-549 (1984). [Z] S.L. Wang and Y.M. Chen, “GPST algorithm with adaptive grid for solving exterior and interior inverse problems of Helmholtz equation”, Proc. 12th IMACS World Congress on Scientifc Computation 2, Paris, 411-413 (1988).
S.L. Wang, V.M. Chen / Exterior and interior
inverse problems
399
bistatic scattering”, Canad. J. Phys. 47, [3] V.H. Weston and W.M. Boerner, “An inverse scattering technique for electromagnetic 1177-1184 (1969). IEEE Trans. Ant. & fiop. AP-18, inverse scattering problem”, [4] W.A. Imbriale and R. Mittra, “The two-dimensional 633-642 (1970). [5] A. Roger, “Newton-Kantorovich algorithm applied to an electromagnetic inverse problem”, IEEE Trans. Ant. & Prop. AP-29, 232-238 (1981). [6] D. Colton, “The inverse scattering problem for time-harmonic acoustic waves”, SIAM Reu. 26, 323-350 (1984). [7] D. Colton and P. Monk, “A novel method for solving the inverse scattering problem for time-harmonic acoustic waves in the resonance region”, SIAM J. Appl. Math. 45, 1039-1053 (1985). [S] G. Kristenson and C.R. Vogel, “Inverse problems for acoustic waves using the penalized likelihood method”, Inoerse Problems 2,461-479 (1986). [9] D. Colton and P. Monk, “The numerical solution of the three-dimensional inverse scattering problem for time-harmonic acoustic waves”, SIAM J. Sci. Stat. Comput. 8, 278-291 (1987). Boundary Elements IX Vol. 3, Fluid Flow [lo] A. Kirsch and R. Kress, “An optimization method in inverse acoustic scattering”, and Potential Applications, eds. C.A. Brebbia, W.L. Wendland and G. Kuhn, Springer, Berlin, 3-18 (1987). [l l] A. Kirsch, R. Kress, P. Monk and A. Zinn, “Two methods for solving the inverse acoustic scattering problem”, Inverse Problems 4, 749-770 (1988). [12] A. Zinn, “On an optimisation method for the full- and the limited-aperture problem in inverse acoustic scattering for a sound-soft obstacle”, Inverse Problems 5, 239-253 (1989). [13] A.N. Tikhonov and V.Y. Arsenin, Solution of III-Posed Problems, Wiley, New York (1977). [14] M. Kac, “Can one hear the shape of a drum ?” Amer. Math. Monthly 73, l-23 (1966). [15] M.H. Protter, “Can one hear the shape of a drum? revisited,” SIAM Rev. 29, 185-197 (1987). [16] M.M. Sondhi and J.R. Resnick, “The inverse problem for the vocal tract: numerical methods, acoustical experiments, and speech synthesis”, J. Acoust. Sot. Amer. 73, 985-1002 (1983). [17] D.S. Tsien and Y.M. Chen, “A numerical method for nonlinear inverse problems in fluid dynamics,” Computational Methods in Nonlinear Mechanics, Univ. of Texas, Austin, 935-943 (1974). [18] A. Bayliss, M. Gunzburger and E. Turkel, “Boundary conditions for the numerical solution of elliptic equations in the exterior region”, SIAM J. Appl. Math. 42, 707-725 (1982). [ 191 J.J. Dongarra, “Performance of various computers using standard linear equations software in a Fortran environment”, Tech. Mem. No. 23, Argonne National Laboratory (1988).