NASA Contractor Repa 181654 

IWE REPORT NO. 88-23 

ICASE 

BOUNDARY SHAPE IDENTIFICATION PROBLEMS 

IN TWO-DIMENSIONAL DOMAINS RELATED TO 

THERPAL TESTING OF MATERIALS 

N A S A -CR - 1 8 1 6 5 4 ) B OIJ N3 AR Y S fT AP E NP8-231541 

I3ENT IFICATION PROEILc’YS IN 7 Ln-DIY F,WS ;Ubi .iL 

DOMAINS R9LATZI) TO TrIEPMAL TYSY’ZK; OF 

PlhTERIALS Final Report :NASA) 33 pC5i~ ILA I! Tic 1 a s 

\;3/59 0141225 

H. T. Ranks 

Fumio Kojima 

Contract No. NAS1-18107 

March 1988 

INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERING 

NASA Langley Research Center, Hampton, Virginia 23665 

Operated by the Universities Space Research Association 

National Aeronautics and 

Space Administration 

LMlgley Rerssrch Cemtar 

Hampton, Virginia 23665 
BOUNDARY SHAPE IDENTIFICATION PROBLEMS IN 

TWO-DIMENSIONAL DOMAINS RELATED TO THERMAL 

TESTING OF MATERIALS 

H. T. Banks 

Center for Control Sciences 

Division of Applied Mathematics 

B r own University 

Providence, RI 02912 

Fumio Kojima 

Institute for Computer Applications in Science and Engineering 

NASA Langley Research Center 

Hampton, VA 23665 

Abstract 

This paper is concerned with the identification of the geometrical structure of the 

system boundary for a two-dimensional diffusion system. The domain identification 

problem treated here is converted into an optimization problem based on a fit-to-data 

criterion and theoretical convergence results €or approximate identification techniques 

are discussed. Results of numerical experiments to demonstrate the efficacy of the 

theoretical ideas are reported. 

This research was supported by the National Aeronautics and Space Administration under NASA Contract No. NAS1-18107 while the authors were in residence at the Insitute for Computer Applications in 

Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. Also, research 

was supported for the first author in part by the National Science Foundation under NSF Grant No. MCS8504316, by the Air Force Office of Scientific Research under contract F49620-86-C-011, and the National 

Aeronautics and Space Administration under NASA Grant NO. NAG-1-517. 
I. INTRODUCTION 

Domain identification problems are important in the design of engineering systems and 

frequently such problems are treated as a branch of the calculus of variations which involves 

nonlinear optimization techniques, optimal control theory, partial differential equations (elliptic, parabolic, hyperbolic, etc.) and related numerical methods. Domain identification 

for elliptic systems has been studied theoretically and numerically by many authors (see 

e.g., [5],[7],[10],[13]). For parabolic systems, a couple of numerical methods for identifying the domain or boundary have been investigated in [14],[15]. Until recently, most 

investigations concentrated on the “optimal shape design problem” which is motivated by 

numerous applications to structural, engine, airplane, ship designs, etc. (see [lo] and the 

references therein). In this paper, our concern for domain identification is motivated by an 

application that is different from these shape design problems. However, as we shall see, 

the resulting theoretical aspects are closely related. Recently, associated with the use of 

fiber reinforced composite materials for aerospace structures, there is growing interest in 

the detection and characterization of large structural flaws which may not be detectable 

by visual inspection. One recent effort has focused on non-destructive evaluation methods 

(NDE) based on the measurement of thermal diffusivity in composite materials (see e.g., 

[SI). Motivated by these problems, we consider the domain identification for parabolic 

systems. 

To explain our approach, we restrict our attention to a 2-D domain identification problem. We consider the bounded domain G(q) in two-dimensional Euclidean space as follows: 

where 51 + t(z1,q) is some parameterized real function which is assumed to characterize 

the unknown part of the boundary and q is a constant parameterization vector to be 

identified among values in a given compact admissible parameter set Q. As depicted in 

1 
1.2 

1 .o 

0.8 

x" 0.6 

0.4 

0.2 

0.0 I1 I 

\ 

\ 

\ 

\ 

\ 

- ._-- 

I 

I; '.__I 

I 

I 

. 

_-_-' 

ac, 

0.0 0.4 0.8 1.2 

XI 

Fig. 1.1. The spatial domain G and its boundary aG1, dG2, aG,, dG4. 

Fig. 1.1, we assume the boundary of G(q) consists of the following components: 

dG1 = {X = (51,~2)10 < 21 < 1, ~2 = 0) 

{X = (x~,x~)Iz~ = 1, 0 < ~2 < e} 

{Z = (~1,~2)(0 < XI < 1, 52 = r(Zl,q)} 

{Z = (Z~,X~)IZ~ = 0, 0 < 52 < e} 

dG2 = 

aG, = 

aG4 = 

The measurement system is described by the following 2-D diffusion equation: 

au(t' - clAu(t, X) + cou(t, z) = 0 in T x G(q) at 

with the initial and boundary conditions, 

(1.la) 

U(O,X) = E,-,(X) on G(q) 

2 

(l.lb) 
aU 

a52 

-ci-+ h(u - f) = 0 on T x dG1 

au -=o 

dn 

on T xdG2 

on T xdG, 

on TxdG1 

(l.lc) 

(1.ld) 

(1.le) 

where c1, co and h are thermal diffusivity, radiation coefficient and heat transfer coefficient, 

respectively, which are given constants, and where T denotes the time interval (0, tf) during 

which the process is observed. In the above system, f is the known boundary input defined 

on T x aG1 and iiro is the given initial function defined on SZ where SZ is a known bounded 

domain in R2 such that f2 2 G(q) for any qcQ. The system output is assumed to be on a 

subset C of the boundary dG1, and mathematically, the observation is taken as 

From a physical point of view, the system state u = u(t,x) represents the temperature 

distribution at time t at location x = (x1,z2) and, the boundary input f and the output 

y correspond, respectively, to the thermal source (for example, by a laser beam) and 

the observation of the temperature distribution at the surface of the material (e.g., by an 

infrared imager) (see [SI for more details). Thus, the search for structural flaws in materials 

may be formulated as an inverse problem for a heat diffusion system. The problem treated 

here is that of identifying, from input and output data { f, EO, y} on (T x dG1) x G x (T x 

E), the constant parameter vector q in Q determining the geometrical structure of the 

boundary aG,. 

In Section 2, we formulate this problem in an abstract setting in a Hilbert space. In Section 3, for computational purposes, we approximate the Hilbert space by finite dimensional 

subspaces and we discuss the convergence analysis for the approximate identification problems. In Section 4, a practical optimization technique based on a finite element approach 

is outlined. Some numerical results for a simple example are given in Section 5. 
11. PROBLEM FORMULATION AND BASIC ASSUMPTIONS 

For the discussions here, we restrict the geometrical structure of the boundary dG, by 

imposing the following hypotheses: 

(H-0) The admissible parameter set Q is a compact subset of R"; 

(H-1) For each gcQ, we have r(g)EWk(O, 1); 

(H-3) There are constants PI and P2 satisfying 0 < PI < i! < P2 < 00 such that, for qcQ, 

~ we have 

P1 L r(€,q) 5 P2 a-e. in (w); 

I 

l and 

(H-4) There exists a constant M such that 

I 

Ir(E,q) - r(€,G)ll,m L MI9 - GI for Q, GEQ, €+A 11, 

I where I - Il,oo denotes the norm of WA(0,l). 

1 We make the following assumptions for the class of system inputs: 

(H-5) L2 (n) ; 

and 

It follows from results in ([9], Ch. 3) that, under the hypotheses (H-1)-(H-3), (H-5), and 

(H-6), for each fixed geQ, there exists for (1.1) a unique solution u in L2(T;H1(G(q))). 

I Following standard procedures in optimal shape design techniques ([lo], Ch. 8, p. 125), 

we introduce the affine mapping 

4 
given by 

Note that this is equivalent to z1 = 21, z2 = r(z1, q)z2/e. Under this coordinate change, the 

system domain G(q) is transformed into the fixed domain 5 

G(q) --+ G = (0,l) x (0, e) 

which is independent of the parameter q. Using this coordinate transformation, we obtain 

the system state ii given by 

this transformed state then satisfies the system equation 

(2.14 

with 

ii(0,z) = ~(z) 0 ~-'(q) on (2.lb) 

aii cy aii 

azl e a22 

c1- - -z2rf(l,q)- +he = 0 on T x dG2 

(2.14 

(2.ld) 
In this system the coefficients are given by 

a11 := c1 

where r' denotes dr/dzl, while aG, has been mapped into 

If we consider a variational formulation similar to that in [2],[9], the system dynamics can 

be described by the variational form: 

ii(0) = To 0 T-l(q) 

where the bracket < a,. > denotes the scalar product in L2(6) and where o(q)(-,.) and 

L(q)(.) denote, respectively, a sesquilinear form on @(e) x H1(e) and a linear functional 

on H1(e). Explicitly, a(q)(-,-) and L(q)(.) are given by for q5,t,fxH1(@ by 

6 
respectively. With some tedious calculations, one can readily establish the following useful 

conditions on the sesquilinear form Q. 

Theorem 1: Let I - IV and I - IH denote the norms in the Hilbert spaces V = H1(G) and 

H = L2(6). Then, under the hypotheses (H-0) to (H-4), the sesquilinear forms a(q)(.,-) 

satisfy the following inequalities: There exist positive constants ICl, A, IC2, and IC3 such that 

for $,$cV we have 

Q(d(4, 4) 2 kll4lC - Xl4l$ (2.10) 

Proof: We wish to show first that the sesquilinear form Q is coercive. From (2.2)-(2.4) 

and (2.8), the principal part of the differential operator becomes 

for (€1, (2) €R2. 

By simple calculations and from (H-3), we have 

(2.13) 

where 

This means the operator is strongly elliptic. For the coefficients bj(z,q)(j = 1,2), from 

(2.5) and (2.6), it follows that 

(2.15) 

(2.16) 

7 
respectively. We note that, from (H-0) to (H-4), 

SUP Ir'(z, q)l < R 

ZlC[O,11 

where R is some constant independent of q. Hence, we obtain 

(2.17) 

lbj(z,q)l 5 K2 < 00 for j = 1,2, (2.18) 

where (assuming R > 1) 

(2.19) 

For the last two boundary integrals in (2.8), by virtue of (H-3), the following inequality 

holds: 

(2.20) 

where ds denotes a line element on aE and ac = dG1 U aG2 U dG1. From (2.14), (2.18), 

and (2.20), we can derive the coercivity property of the sesquilinear form. Namely, the 

sesquilinear form satisfies 

(2.21) 

where 

K3 = min (* h) . 4P,2 ' 

Friedrich's second inequality ([I], p. 124) asserts that if i3c is a nontrivial subset of 82, 

then there exists a positive constant CY such that 

(2.22) 

8 
By applying this to the last parenthesis of (2.21), we can conclude that 

where 

respectively. 

To prove the boundedness of a(q), we note that 

(2.23) 

(2.24) 

The first t'mee integrals of RHS in (2.24) satis€y 

(2.26) 

(2.27) 

respectively. From (2.2) to (2.4), it follows that, under hypotheses (H-1) and (H-3) (see 

W7)), 

sup laij(z,q)( 5 K6 < 00 (2.28) 

i,j<2 

~€[O,~lX IO,4 

where 

Clt2 Ks = K(R2 + 1). 

From (2.18) and (2.25)-(2.28), we have 

9 
Furthermore, the boundary integral term satisfies 

which follows from the trace inequality ([l], p. 124) 

where 

Consequently, we can prove the boundedness property of a(q) (e, -). 

To establish the continuity property, we note that, for any q and GEQ, 

Under the hypotheses (H-1) and (H-3), we argue that 

10 

(2.30) 

(2.31) 

(2.32) 
and 

Applying these inequalities into (2.32), we have 

where 

2Cl.t R 2~1.t~ R 2RP2 C1 R 2cllR RP2 lh 

p," ma+---, P," K$j = --mmaz(-, 1)+ maz( - ,l)+-max(-, 1)+- Pl Pl Pt Pt P1 P1 

From the hypotheses (H-4), we can thus infer the continuity of the sesquilinear form 

a(q)(., a) with respect to the parameter q in Q. The proof has been completed. 

For the system (2.7), the output can be represented as the restriction of G(t) to a subset 

C C aG1 of positive measure, i.e, 

We assume (see (H-0)) throughout that the admissible parameter set Q is a given compact 

subset of R". The fundamental identification problem considered here is based on the 

fit-to-data functional (see [2]) given by 

(2.35) 

where F = L2(C), {yd(t)ItrT are given observed data, and y(t,q) is the solution of (2.7) 

corresponding to qEQ. Then our problem is stated as follows: 

(IDP) Find q'cQ which minimizes J(q) given in (2.35) subject to the system (2.7) and 

(2.34). 

In the next section, we consider a family of approximating identification problems associated with (IDP). 

11 
111. APPROXIMATE IDENTIFICATION PROBLEMS 

The approximation scheme we have employed is based on the use of a finite element 

Galerkin approach to construct a sequence of finite dimensional approximating identification problems. Let us choose U~=l(q5r}~1 as a set of basis functions in H1(G). That is, 

for all N, {@}~v=l are linearly independent and UN span{qhi NN }i=l is dense in the V norm 

in V = IT1(@. We choose the approximation subspaces as 

HN := span{#,#,..- ,&}. 

Then, we can define the approximate solution of Eq. (2.7) by 

where wy(t, q) are chosen such that for j = 1,2,. . - , N, 

and 

(3.1) 

(3.2~) 

(3.2b) 

Hence the system (2.7) and the output (2.34) can be approximated by solving the system 

N-N c w (t,Q) + AN(a)wN(t,q) = FN(t,q) (3.34 

wN(0) = 57: (3.3b) 

[J”(t,q)]i := L(t,q)(#) for i = 1,2,.-.,N 

12 
The approximating identification problems thus take the following form: 

(AIDP)N Find tNc& which minimizes 

subject to the approximating system (3.3) and (3.4). 

Our convergence results for the finite element schemes are summarized in the following 

two theorems. 

Theorem 2: Let {qM} C Q be a sequence such that qM 4 qcQ as M 4 00 and 

let GN(qM) and G(q) be the solutions of Eqs. (3.3) and (2.7) corresponding to qM and q, 

respectively. Then, under hypotheses (H-0) to (H-6), we have iiN(qM) 4 G(q) strongly in 

L~(T; ~(6)) as N, M -+ 00. 

Theorem 3: Let ijN be a solution of the problem (AIDP)N. Then the sequence {GN} 

admits a convergent subsequence {ijNk} with ijNk 4 as k -+ 00. Moreover, G is a solution 

of the problem (IDP). 

The proof of Theorem 2 follows from the general convergence framework for parameter 

identification problems given in [3] and [4]. To ensure the desired convergence, it suffices to 

show that the sesquilinear form a(q) (-, .) satisfies the continuity, coercivity and boundedness conditions as stated in [3], [4]. But this is a result of Theorem 1 under the hypotheses 

(H-0) to (H-4). 

The proof of Theorem 3 can be carried out by using Theorem 2 and the compactness 

of &. Since ijN is a solution of the problem (AIDP)N, it is clear that 

Thus, if we can argue that for any qM 3 q in Q, 

NM in L2(T;F) as N,M + 00, Y (q 1 -+ Yk.) 

13 
then, we can obtain the desired inequality 

J(G) I J(q) for bQ 

by taking limits in (3.6). But the needed arguments follow immediately from Theorem 2 

where K is independent of aNk and t. 

IV. OPTIMIZATION TECHNIQUES FOR THE APPROXIhlATE ESTIh 

TION PROBLEMS 

LLet aN be an optimal solution of the problem (AIDP)N. Then a necessary condition 

for aN to be optimal is characterized by 

V,JN(aN) - (q - aN) 2 0 for VqcQ (44 

where V, denotes the gradient of JN(q) with respect to q. From Eq. (3.5), we have for 

k = 1,2,.-. ,n 

NN [v,JN(q)]k = q))r(cb (t, q) - yf(t))dt 0 

Using the same procedure as in [9], we can evaluate the gradient vector by (for k = 

1,2,.. - ,n) 

where vN(t,q) is the solution of the adjoint equation, 

14 
W"(tf,Q) = 0. (4.3b) 

In Eq. (4.3a), the matrix A*"(q) is given by 

where a*(q)(., e) is the adjoint sesquilinear form of a(q)(-, e) defined by 

:= Jj& 

Consequently, the optimality condition (4.1) of the problem (AIDP)N can be characterized 

by 

In the sequel, we discuss computer implementation of numerical schemes for the problem (AIDP)N. Since we can evaluate the gradient of the cost function using (4.2), many 

optimization techniques for the constrained problems are readily applicable to our problem 

(see [ll] and the references therein). For ease in exposition, here the compact set Q C RN 

is assumed to be defined by 

where ll and q denote a given constraint matrix and vector, respectively. For the numerical results reported in this paper, we used the gradient projection method [12] which is 

15 
a particularly useful technique for optimization problems with the linear inequality constraints such as those given in (4.5). We use this method as presented in [12]; the iterative 

algorithm for finding ijN can thus be stated as follows: 

Step 0: Choose an initial value q(O) in Q and set i = 0. 

and proceed to Step S; otherwise, proceed to Step 2. 

Step 2: Compute the current direction by 

where 

P = I - l-I;(rIpn;)-lIIp 

and lIp includes the gradient of all currently active constraints associated with matrix l-I. 

If g(') # 0, proceed to Step S; otherwise, proceed to Step 4. 

Step 3: Compute Am,, ('1. satisfying 

where fi is the largest step that may be taken from q(') along g(') witllout violating any 

constraint. If Afin = i, then add the new contraints to the matrix l-Ip and proceed to Step 

4; otherwise, the new approximation to the solution is given by 

(;+I) = ,('I + A('). g('). min 

Replace i + 1 by i and return to Step I. 

Step 4: Compute the vector O(q) by 

e(q(')) = -(npnp) I -1 n,v,J(q(')). 

16 
If all components of 0 are nonnegative, then set 

and terminate the computation; otherwise, delete the column of lTp corresponding to the 

smallest component of O(q(i)), replace i + 1 by i and return to Step I. 

V. NUMERICAL PROCEDURES 

In a series of numerical experiments, we used a test example constructed as follows: 

We chose a function r(q), generated the corresponding solution numerically, added random 

noise, and then used this as "data" for our inverse algorithm. The parameter function 

r(E,q) to be identified is a piecewise cubic polynomial function (see [6] for more details). 

We denote the knot sequence for r by 

0 = .ro" < 71" < ..- ; 7, n : 7,+1 = 1 

and the unknown function r((, q) is given by 

r(€,.) = Pr(€) = %,i + U2,i(€ - T?) +as$(€ - T32/2 + Q,i(E - T33/6 

for r? 5 €5 i = O,l,---,n. 

The unknown parameter vector q = {qi}rZ1 is then given by 

Further, we assume 

Pop, q) = PJ1, q) = e (5.3) 

Pb(0, q) = PXl, q) = 0. (5.4) 

Substituting (5.2), (5.3), and (5.4) into (5.1), the coefficients {ak,i} can be determined 

uniquely and r(t, q) satisfies the hypotheses (H-1), (H-2), and (H-4). In order to guarantee 

the hypothesis (H-3), we impose the constraints 

17 
Hence, the matrix ll (2n x n) and the vector 

class Q (see (4.5)) is given by 

(2n x 1) defining the admissible parameter 

- 

0 

1 

n= 

1 

-1 

c 

1 

-1 

0 

- 

!I= 

- P2 - 

-A 

i 

P2 

~ -P1- 

1 -1 1 

To discretize the system model by the finite element method, the domain G is divided 

into a finite number of elements {ek}fZ1 (K 5 N) and a number of nodes defined by 

{zi = (z~,~~)}~~ are selected in G. For convenience of computations, we set A! = 1 in e. 

Each element is preassigned as an axiparallel rectangle with nodes at the vertices. The 

restriction of 4: to any element ek is given by the bilinear polynomial form, 

., 

for z = (21, z2)cek k = 1,2,-. - , K and i = 1,2,-. - , N. 

The coefficients {c!?} can be chosen such that each polynomial form (5.5) satisfies the 

properties of a piecewise bilinear basis function (see e.g., [l], Ch. 5). The integration of 

element matrices CN,AN(q),A*N(q), and Cf, and the element vectors FN and Y/ can 

be computed numerically by a Gauss-Legendre formula. Thus, the state model (3.3) and 

its adjoint system (4.3) can be solved numerically by an implicit scheme with respect to 

discrete time t = ih (i = 0,1, - - - , m) where h = tj/rn. The evaluation of cost functional JN 

and its gradient V,JN is the computationally expensive part of our algorithm since these 

involve the integration of the states wN(t, q) and the adjoint states wN(t, q) with respect 

to time t over T. This can be accomplished by using the two-point Gauss formula. 

The input data are preassigned as 

f(E) = 0, for EdG1. 

18 
The known parameters c1, CO, and h in Ea. (1.1) were set as 

~1 = 0.034, ~2 = 0.001, h = 0.1. 

The observed data {yd(t)} were generated by solving the finite element model (3.3). The 

number of finite elements and nodes in the numerical experiments were set as K = 256 (= 

16 x 16) and N = 289 (= 17 x 17), respectively. The final time and number of time 

divisions were taken as ti = 10 and m = 100. Random noise at various levels from 0% to 

50% was added to the numerical solution, thereby producing simulated noisy "data" for 

the algorithm. The set C relative to data acquisition was given by 

P 

C = u Nzj 

j=l 

where Nzj denotes a neighborhood of points xi at aG1, i.e., 

Nzj = (~j - E,X~ + E) for j'= ?,2,---,p. 

Using such data, the estimation algorithm given in Section 4 was tested. 

Example 1: In this example, the dimension of unknown vector was taken as n = 4 and the 

knot sequence {r~};"=+o' was given by 

r,? = i/5 for i = 1,2,.-.,5. 

The values of the true parameters were chosen as 

qi = r(r,?) = 0.8 for i = 1,2,3,4. 

The lower and upper bounds of the unknown parameter vector were taken as PI = 0.3 and 

P2 = 1.1, respectively. The initial guesses for the parameters were given by 

Q~O) = 1 for i = 1,2,3,4. 

The number of sensors was taken as p = 9. Table 1 shows the estimated parameter numerical results for the data with noise free, 5%, lo%, and 50% relative noise and Figure 5.1 

19 
shows the estimated parameter function r( 6, a”) and true function r( 6, q) which correspond 

to the estimated boundary shape and true boundary for the 10% noise case. 

G1 G2 G3 G4 

True Value 0.800 0.800 0.800 0.800 

Initial Guess 1.000 1.000 1.000 1.000 

Noise iteration 6 0.850 0.882 0.880 0.849 

Free 

5% 

Noise 

10% 

Noise 

50% 

Noise 

+Xi54 1% - &I2 , 

3.36 x lo-’ 

0.820 0.819 0.821 0.819 

0.810 0.785 0.810 0.806 

0.851 0.879 0.878 0.844 

0.828 0.821 0.818 0.820 

0.815 0.797 0.791 0.814 

iteration 13 

iteration 17 

9.92 x lo-’ 

5.39 x lo-’ 

3.25 x lo-’ 

1.09 x 

5.61 x 

itkration 6 

iteration 13 

iteration 17 

0.849 0.853 0.861 0.834 

0.787 0.847 0.835 0.750 

0.827 0.792 0.799 0.796 

0.922 0.820 0.748 0.808 

0.902 0.793 0.772 0.886 

0.813 0.783 0.727 0.794 

iteration 6 

iteration 14 

iteration 18 

iteration 7 

iteration 15 

iteration 19 

2.51 x 

1.95 x 

7.20 x lo-’ 

3.36 x lo-’ 

3.40 x 

1.90 x lo-’ 

Table 5.1. True Value and Estimated Values in Example 1. 

Example 2: We chose the same dimension of unknown parameter vector as in Example 1 

and we also used the same knot sequence. In this example, however, the values of the true 

parameters were preassigned as 

q1 = q4 = 0.9 

respectively. The lower and upper bounds, initial guess of unknown vector, and number 

of sensors were given by the same values as in Example 1. Table 5.2 shows the numerical 

results obtained here for the various sets of noisy data. Figures 5.2 and 5.3 represent the 

estimated parameter function for the case of 20% and 50% noisy observation case. 

20 
1.2 

1 .o 

0.8 

x" 0.6 

0.4 

0.2 

0.0 

T-. 

I1 U& Fufiction 

initial Guess _____-__---- 

Estimated Function - - - - - - 

1 I 1 I 

0.0 0.2 0.4 0.6 

X1 

0.8 1 .o 

Figure 5.1. True Function and Estimated Function in Example 1 (10% Noise). 

21 
True Value 

61 62 63 64 +cis4 1% - ail2 

0.900 0.600 0.600 0.900 

Initial Guess 

Noise iteration 16 

Free iteration 25 

5% iteration 16 

Noise iteration 26 

10% iteration 16 

Noise iteration 27 

20% iteration 16 

Noise iteration 27 

25% iteration 16 

Noise iteration 26 

50% iteration 16 

Noise iteration 28 ~ ~ ~ 

Table 5.2. True Value and Estimated Values in Example 2. 

1.000 1.000 1.000 1.000 

0.817 0.792 0.778 0.821 

0.894 0.607 0.602 0.893 

0.953 0.694 0.811 0.906 

0.896 0.604 0.605 0.898 

0.902 0.807 0.674 0.949 

0.908 0.594 0.581 0.896 

0.917 0.683 0.812 0.904 

0.902 0.603 0.610 0.887 

0.939 0.819 0.684 0.962 

0.868 0.563 0.563 0.874 

1.02 0.618 0.680 0.915 

0.951 0.574 0.599 0.953 

Example 3: In this example, we deal with a somewhat more difficult case as compared 

with Examples 1 and 2. We set the dimension of paraIzleter space as n = 8 and we chose 

the knot sequence as 

7.14 x 

2.85 x 10-3 

5.93 x 10-2 

1.93 x lo-’ 

5.62 x 

5.70 x 

4.19 x lo-’ 

6.15 x 

1.66 x 

3.64 x lod2 ‘ 

5.39 x 10-~ 

1.95 X 

{ri 89 }i=o r,? = i/9 for i = 0,1,2, * ,9. 

True parameter values were given by 

respectively. Figure 5.4 shows the corresponding boundary shape to be identified. The 

number of sensors was taken as p = 17. The bounds and initial guesses for the parameter 

22 
N 

X 

0.2 

0.4 1 

Initial Guess 

Estimated Function - - - - - - 

- - - - - - - - - - - -. 

- 

I I I I 

True Function 

Figure 5.2. True Function and Estimated Function in Example 2 (20% Noise). 

23 
1.2 1 

1 .o 

0.8 

0.6 

0.4 

0.2 

0.0 

True Function 

Initial Guess ____________ 

Estimated Function - - - - - - 

1 1 1 1 

0.4 0.6 0.8 1 .o 0.0 0.2 

X1 

Figure 5.3. True Function and Estimated Function in Example 2 (50% Noise). 

24 
1.2 

1 .o 

0.8 

x" 0.6 

0.4 

0.2 

0.0 

0.0 0.2 0.4 0.6 0.8 1 .0 

X1 

Figure 5.4. Unknown boundary shape in Example 3. 

vector were the same as in Examples 1 and 2. We ran numerical experiments for the case 

of noise free, 5%, lo%, 2O%, and 50% noisy observations. Table 5.3 shows the estimated 

parameter vector obtained here. Figure 5.5 represents the estimated boundary curve for 

the 10% noisy data. 

25 
Initial Guess 

Noise iteration 16 

Free iteration 23 

5% iteration 16 

Noise iteration 23 

10% iteration 16 

Noise iteration 24 

20% iteration 16 

Noise iteration 24 

50% iteration 16 

Noise iteration 24 

Throughout the numerical experiments, we checked the robustness of the algorithm with 

1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 

1.039 1.027 0.887 0.808 0.674 0.926 1.029 1.033 3.08 x 

1.014 1.007 0.943 0.624 0.615 0.943 1.003 1.017 7.27 x 

1.042 1.031 0.879 0.749 0.666 0.920 1.049 1.031 2.85 x 

1.023 0.987 0.948 0.626 0.628 0.937 1.002 1.020 7.90 x 

1.037 1.058 0.899 0.733 0.721 0.881 1.044 1.025 2.82 x 

1.010 1.009 0.955 0.586 0.578 0.948 0.994 1.000 6.20 x 

1.034 0.989 0.950 0.592 0.595 0.994 1.018 1.061 1.61 x 

1.022 0.985 0.915 0.616 0.610 0.941 0.993 1.061 1.06 x 

1.090 1.176 0.872 0.637 0.654 0.933 1.166 1.096 4.01 x 

1.030 1.033 0.939 0.566 0.603 0.933 1.046 1.052 1.47 x 

respect to noise in the observed data. Results in three examples indicated that the algorithm worked very well (i.e., as expected) for various noise levels. Furthermore, we checked 

the sensitivity of the algorithm with respect to the number of sensors. Specifically, we compared in Examples 2 and 3 the number of sensors (p) with the dimension of parameter 

space (n). In Example 2, for data with p = 5(> n = 4), the algorithm still yields an almost 

identical fit (to that for p = 9) even in 50% noise case while the fit could not be achieved 

under the reduced observation case p = 3(< n). Also, in Example 3, (where n = 8) the 

fit could not be obtained with p = 3 or p = 5, while the algorithm performed well with 

p = 9 (> n). Carrying out a large number of other numerical tests in addition to those 

reported for Examples 2 and 3, we suggest that the aIgorithm requires a number of sensors 

which is at least equal to the number of dimensions of parameter space, i.e., p 2 n. 

VI. CONCLUDING REMARKS 

In this paper, we have discussed techniques for estimating the system boundary shape 

in two dimensional parabolic systems. By using a simple coordinate transformation technique, the parabolic PDE defined on unknown spatially varying domain was converted 

26 
1.2 

1 .o 

0.8 

x" 0.6 

0.4 

0.2 

0.0 

True Function 

Initial Guess ------------ 

Estimated Function - - - - - - 

0.0 0.2 0.4 0.6 0.8 1 .o 

X1 

Figure 5.5. True Function and Estimated Function in Example 3 (10% Noise). 

27 
into the same type PDE with unknown coefficients defined on a fixed domain. Thus, our 

fundamental approach was placed within the theoretical framework for parameter identification problems given in [2],[3], and [4]. The practical utility of our algorithm is supported 

through a series of numerical experiments, a summary of which is given in Section 5. These 

simulations were carried out on the Sun Microsystems at ICASE, NASA Langley Research 

Center. For three different numerical examples, using data with no noise, the proposed 

algorithm yields an almost perfect fit, while, as expected, the fit degenerates significantly 

as noise in the observation becomes more pronounced. 

Although here we discuss only the case where the unknown boundary shape is represented by a simple function of one variable, our basic parameter estimation ideas and 

techniques can be readily extended to consider more general classes of geometrical structures for the system boundary. For example, we may also treat the case where the unknown 

boundary shape is characterized by 

r(q,q,z2) = o for (zl,zz)c~~. 

We are currently pursuing investigations for these cases. 

ACKNOWLEDGEMENT 

The authors would like to express their sincere appreciation to Dr. W. Winfree and 

Ms. Michele Heath (Instrument Research Division - Materials Characterization Instrumentation Section, NASA Langley Research Center) for numerous conversations and encouragement during the course of this research. Their suggestions and questions were an 

essential motivation of our efforts. 

28 
References 

[l] 0. Axelsson and V. A. Barker, Finite Element Solution of Boundary Value Problems, 

Academic Press, New York, 1984. 

[2] H. T. Banks, “On a variational approach to some parameter estimation problems,” 

Distributed Parameter Systems, Lecture Notes in Control and Information Sciences, 

Vol. 75 (1985), pp. 1-23, Springer-Verlag, New York. 

[3] H. T. Banks and K. Ito, “A theoretical framework for convergence and continuous 

dependence of estimates in inverse problems for distributed parameter systems,” 

LCDS/CCS Report No. 87-20, Brown University, March 1987; Appl. Math. Lett., 

Vol. 0 (1987), pp. 31-35. 

I41 H. T. Banks and K. Ito, “A unified framework for approximation in inverese problems 

for distributed parameter systems,” LCDS/CCS Report No. 87-42, Brown University, 

October 1987; Control Theory Adv. Tech., to appear. 

[5] D. Begis and R. Glowinski, “Application of the finite element method to the approximation of an optimum design problem,” Appl. Math. Optim., Vol. 2 (1975), pp. 

130-1 68. 

[6] C. de Boor, A Practical Guide to Splines, Applied Mathematical Science, Vol. 27, 

Springer-Verlag, New York, 1978. 

[7] D. Chenais, “On the existence of a solution in a domain identification problem,” 

J. Math. Anal. Appl., Vol. 52 (1975), pp. 189-219. 

[8] D. M. Heath, C. S. Welch, and W. P. Winfree, “Quantitative thermal diffusivity 

measurements of composites,” in Review of Progress in Quantitative Nondestructive 

Evaluation, (D. G. Thompson and D. E. Chimenti, eds.), Plenum Publ., Vol. 5B 

(1986), pp. 1125-1132. 

29 
[9] J. L. Lions, Optimal Control of Systems Governed by Partial Differential Equations, 

Springer-Verlag, New York, 1971. 

[ 101 0. Pironneau, Optimal Shape Design for Elliptic Systems, Springer-Verlag, New 

York, 1983. 

[ll] E. Polak, Computational Methods in Optimization, Academic Press, New York, 1971. 

[12] J. B. Rosen, “The gradient projection method for nonlinear programming, Part I: 

Linear constraints,” SIAM J. Appl. Math., Vol. 8 (1960), pp. 181-217. 

[13] J. Simon, “Differentiation with respect to the domain in boundary value problems,” 

Numer. Funct. Anal. Appl., Vol. 2 (1980), pp. 649-687. 

[14] Y. Sunahara, Sh. Aihara, and F. Kojima, “A method for spatial domain identification of distributed parameter systems under noisy observations, 

Proc. 9th IFAC World Congress, Budapest, Hungary, 1984, Pergamon Press, New 

York, 1984. 

[15] Y. Sunahara and F. Kojima, “Boundary identification for a two dimensional diffusion system under noisy observations,” Proc. 4th IFAC Symp. Control of Distributed 

Parameter Systems, UCLA, California, 1986, Pergamon Press, New York, 1986. 

30 
Iu/\sA MIWI ~craolcs am Report Documentation Page 

5oxe &wnrwstralm 

1. Report No. 

NASA CR- 18 1654 

ICASE Report No. 88-23 

2. Government Accession No. 

7. Authorls) 

H. T. Banks and Fumio Kojima 

17. Key Words (Suggested by Authork)) 

parameter estimation, distributed 

parameter systems, parabolic systems 

domain identification 

9. Performing Organization Name and Address 

Institute for Computer Applications in Science 

Mail Stop 132C, NASA Langley Research Center 

and Engineering 

VA 7- - 5775 

gency Name and Address 

National Aeronautics and Space Administration 

Lsngley Research Center 

Hampton, VA 23665-5225 

15. Supplementary Notes 

18. Distribution Statement 

59 - Mathematical and Computer Sciences 

64 - Numerical Analysis 

66 - Systems Analysis 

Unclassified - unlimited 

(General) 

3. Recipient's Catalog No. 

19. Security Classif. (of this report) 20. Security Classif. (of this page) 21. No. of pages 

Unclassified Unclassi € ied 32 

5. Report Date 

March 1988 

6. Performing Organization Code 

22. Price 

A0 3 

8. Performing Organization Report No. 

88-23 

10. Work Unit No. 

505-90-21-01 

11. Contract or Grant No. 

NAS1-18107 

13. Type of Report and Period Covered 

Contractor Report 

14. Sponsoring dgency Code 

Langley Technical Monitor: 

Richard W. Barnwell 

Submitted to Quart. Applied Math. 

Final Report 

16. Abstract 

This paper is concerned with the identification of the geometrical structure 

j of the system boundary for a two-dimensional diffusion system. The domain 

identification problem treated here is converted into an optimization problem 

based on a €it-to-data criterion and theoretical convergence results for 

approximate identification techniques are discussed. Results of numerical 

experiments to demonstrate the efficacy of the theoretical ideas are reported.