IDENTIFICATION OF A CORE FROM BOUNDARY DATA
IDENTIFICATION OF A CORE FROM BOUNDARY DATA WOLFGANG RINGy
Abstract. The problem of determining the interface separating regions of constant density from boundary data of a solution of the corresponding potential equation is considered. An equivalent formulation as a nonlinear integral equation is obtained. Fourier methods are used to analyze and implement the problem. Numerical experiments based on a regularized least-squares method are presented. Key words. inverse source problem, nonlinear integral equation, Fourier transform, Tikhonov
regularization
AMS subject classi cations. 31B20, 35J05, 35R25, 35R30
1. Introduction. Suppose (r; ') denote polar coordinates in the plane and r = a(') is a continuous function, with graph of a contained in some domain R2 .
The subject of this work is the following problem:
Determine the interface r = ( a(') in u = fa (r; ') = 1 for 0 r < a(') 2 for a(') r 2 on R from known boundary data z1 = uj@ and z2 = @u @ n j@ .
(1)
We consider only the situation where = U , the open unit disk in R2 . The following introductory remarks, however, do not only apply to U , but to any simply connected open subset of the plane with smooth boundary. Problem (1) occurs in the context of recovering an unknown mass distribution inside a planar body from measurements of the potential at the boundary. Suppose a mass distribution consisting of two layers is given on : a region of constant density 1 inside (a core) and an annular region (a mantle) of constant density 2 enclosing the core, and suppose that the interface between the two regions is given by the curve r = a('). The gravitational potential u generated by this mass distribution is then a solution of the Poisson equation u = fa . Therefore (1) can be interpreted as the problem of nding the curve which separates the regions of constant density from measurements of the potential at the boundary of . Closely related to (1) is the inverse domain recognition problem in potential theory: Find the domain ~ in known g : @ ! R.
Z
~
E (x; ) d = g(x); x 2 @ , from
(2)
Here ~ denotes some open subset of and E (x; ) = 41 log jx , j2 : This work was supported by Fonds zur F orderung der Wissenschaftlichen Forschung, Austria, under P 7869-Phy and Bundesministerium fur Wissenschaft und Forschung, under project `Mathematische Optimierung und Inverse Prozesse' y Institut f ur Mathematik, Technische Universitat Graz, Kopernikusgasse 24, A-8010 Graz, Austria 1
2
W. RING
In fact, (1) and (2) are equivalent in the following sense. Suppose (1) holds for
a, z1 and z2 . Using Green's formula and well known identities from the theory of
boundary integral equations (c.f. Chen and Zhou [4, p 223, (6.43) and p 224, (6.45)]), we obtain Z Z Z fa() E (x; ) d = z1 (2x) , z1 () @@n E (x; ) d + z2 ()E (x; ) d
@
@
(3) and hence, with a denoting the core enclosed by a, we nd Z Z n E (x; ) d = ,1 z1 (2x) , z1 () @@n E (x; ) d 1 2 @
a (4)
Z
+
@
z2()E (x; ) d , 2
Z
E (x ) d
o
for x 2 @ . This shows that a (and hence a) can be obtained from z1 and z2 as a solution of (2) with g given by the right hand side of (4).R On the other hand, the Newtonian potential w(x) = a E (x; ) d satis es w = fa for 1 = 1, 2 = 0. We set z1 := g = wj@ and z2 := @w @ n j@ . Starting with (2), only z1 but not z2 is known. However, if we apply Green's formula and well known integral identities [4, p 223, (6.43) and p 224, (6.45)] to w() and E (x; ), we obtain the integral equation Z Z z2 () E (x; ) d = z1 (2x) + z1 () @@n E (x; ) d for x 2 @ : (5) @
@
The right-hand side of (5) de nes a linear operator L1 : H , 12 (@ ) ! H 21 (@ ) acting on z2 , which is Fredholm of index 0. In case = U , the kernel on L1 consists of the constant functions on @ , hence z2 is determined by (5) up to an additive constant c. (c.f. [4] for the facts mentioned above.) Using Gauss' divergence theorem, we obtain the additional relation Z
Z
Z @w z2 () d = () d = fa () d; @
@ @n
from which the constant c can be uniquely determined. Hence, we arrive at (1) with known boundary data z1 and z2 and unknown a. The domain recognition problem (2) was studied extensively in the past. We mention Novikov [14], Prilepko [15] [16], Cabayan and Belford [3], Brodsky [2], Strakhov and Brodsky [19], the book by Isakov [9] and the literature cited there. The starting point of yet another equivalent formulation is the boundary value problem (
u = fa on : u j@ = z 1
(6)
With u = u(a; z1) denoting the solution of (6), we can de ne the operator and we could solve
G(a; z1 ) = @@n u(a; z1) @
G(a; z1 ) = z2
(7)
IDENTIFICATION OF A CORE
3
for a. A disadvantage of (7) arises if we have to deal with noisy data. Suppose z1 and z2 are replaced by inaccurate measurements z1i and z2i . Then (7) becomes G(a; z1i ) = z2i : This means that we have to deal not only with perturbed data on the right-hand side but also with a perturbed operator Gi = G( ; z1i ), which makes the analysis of the problem rather complicated. A completely di erent, optimization theoretic approach to a problem similar to (1) is presented in Kunisch and Pan [11]. There the augmented Lagrangian method was used to obtain a formulation of the identi cation problem as a constrained minimization problem. In this paper we present an approach to (1) which is very similar to the considerations which led to (4), with the di erence that we use the Green's function for the Dirichlet problem on U instead of the fundamental solution E (x; ). As in (4) we obtain a nonlinear operator equation F ( a) = z (8) where F is an integral operator over the variable domain a and z = z (z1; z2 ). This method is restricted to those domains for which a Green's function is known. The loss of generality is however rewarded by a very compact form of the resulting integral equation for a. Moreover the choices = U and a = a(') make it possible to use Fourier analysis, which will turn out to be a very powerful tool for analyzing and implementing the problem. We show that (8) can be reduced to an in nite system of nonlinear functionals of a, attaining prescribed values which depend on the data z1 and z2 . We derive results on smoothness, di erentiability, injectivity and ill-posedness for (8). To deal with the ill-posedness of the problem, we approximate (8) by a regularized least-squares problem, and we check conditions which ensure stability and convergence for the approximating problem. Finally we present numerical results for di erent interfaces and di erent choices of the regularization term. 2. Preliminaries. We now x the terminology for the subsequent sections and we provide some basic technicalities which are needed henceforth. The notation (r; ') and (; ) is used for polar coordinates of points x 2 R2 and 2 R2 respectively. U = B (0; 1) denotes the open unit disk in the plane R2 . H 1 (U ) is the Sobolev space @u 2 L2 (U ) for i = 1; 2, endowed with the inner of all functions u 2 L2 (U ) for which @x i product
@u @v
@u @v hu; viH 1 (U ) = hu; viL2 (U ) + @x ; @x L2 (U ) + @x ; 2 : 1 1 2 @x2 L (U )
Moreover we de ne the Hilbert space H (; U ) = fv 2 H 1 (U ) : v 2 L2(U )g; (9) with inner product hu; viH (;U ) = hu; viH 1 (U ) + hu; viL2 (U ) : (10) Points on the boundary of U are identi ed with the corresponding angle in polar coordinates, i.e. we set @U = f' + 2Z : ' 2 Rg = R=2Z =: T. For l 0, the Sobolev space H l (T) is de ned by X H l (T) = ff 2 L2 (T) : (1 + n2 )l jf^(n)j2 < 1g; (11) n2Z
4 where
W. RING Z 2 1 ^ f ( ) e,in d f (n) = 2
0
(12)
denotes the Fourier transform of f . H l (T) is a Hilbert space with respect to the inner product X hf; giH l (T) = (1 + n2 )l f^(n)^g(n): (13) n2 Z
By H ,l (T) (l 0) we denote the dual of H l (T). With the Fourier transform de ned on H ,l (T) by f^(n) := hf; e,in' iH ,l (T);H l(T) we nd X hf; giH ,l (T);H l(T) = f^(n) g^(,n) (14) n2Z
for the duality pairing h; iH ,l (T);H l(T) and for f 2 H ,l (T) and g 2 H l (T). Moreover, H ,l (T) is characterized by (11) and the inner product on H ,l (T) is given by (14), in both cases with l replaced by ,l. We have X f = f^(n) ein' (15) n2Z
for f 2 H l (T), l 2 R, where the series (15) converges in H l (T). The embeddings H l (T) ,! H m(T) and C (T) ,! H s (T) are compact and the images are dense for l > m and s > 21 . We de ne the di erential operator D : H l (T) ! H l,1 (T) by X Df = in f^(n) ein' : (16) n2Z
It is a bounded linear operator with ker D given by the constant functions on T. The space C 1 (T) of all in nitely di erentiable functions on T is characterized by C 1 (T) = ff : T ! C : jnjk jf^(n)j ! 0 as n ! 1 for all k 2 N g: (17) The preceding facts follow as a simple special case from theorems on Sobolev spaces on smooth compact manifolds as presented e.g. in Wloka [20]. 12 1 The trace map 0 : H (U ) ! H (T) is the unique continuous extension of the operator 0 : C 1 (U ) ! C 1 (T) de ned by 0 u =1 uj, from C 1 (U ) to H 1 (U ). 0 is a bounded linear surjection from H 1 (U ) on H 2 (T). Analogously we de ne 1 : 12 , 1 1 H (; U ) ! H (T) to be the unique extension of 1 : C (U ) ! C (T), where @u
1 u = @ n T, and n denotes the unit outward normal vector to the boundary T. If there is no danger of misinterpretation, we write 0 u = u and 1 u = @u @ n . (c.f. [20, p 130, Satz 8.7] and Dautray and Lions [5, p 381, Lemma 1] for these facts). The Poisson kernel P : U ! C is de ned by
P (r; ') =
X
rjnj ein' :
n2Z
(18)
The solution u 2 H 1 (U ) of the Dirichlet problem (
u = 0 on U
0 u = f 2 H 12 (T)
(19)
5
IDENTIFICATION OF A CORE
is given by the Poisson integral
u(r; ') = 21
2
Z
0
P (r; ' , ) f ( ) d = P
X
rjnj f^(n) ein' :
(20)
n2Z rjnj f^(n) ein' ! u(r; ') in H 1 (U )
For a proof of this, observe that uN (r; ') := jnjN P and 0 uN = jnjN f^(n) ein' ! f in H 21 (T) as N ! 1. Hence 0 u = f and obviously u = 0 on U . This proves (20). Since uN = u = 0 on U , it is also clearPthat uN ! u in H (; U ) as N1 ! 1. By the trace theorem, we nd
1 uN = jnjN jnj f^(n) ein' ! 1 u in H , 2 (T). We can therefore state the following Proposition 2.1. Suppose given u 2 H 1 (U ), f 8 > :
1 u = g Then
g=
X
jnj f^(n) ein' ;
n2Z
and hence g^(n) = jnjf^(n) for all n 2 Z. Remark 2.2. For f 2 H l (T), l 2 R we de ne
f~ =
X
,i sign(n) f^(n) ein' :
n2Z
f~ is called the function conjugate to f , or the Hilbert transform of f . The Hilbert transform is a bounded linear operator on H l (T) with kernel given by the constants. (c.f. Edwards [7] for further details on conjugate functions.) A short computation shows that X jnj f^(n) ein' = (Df ) ; (21) n2Z
where D is the (tangential) derivative operator de ned in (16). With this notation we can write the conclusion of Proposition 2.1 in the form
1 u = (D 0 u) :
(22)
3. An Alternative Approach. Let us return to problem (1) with = U and given data z1 2 H 21 (T) and z2 2 H , 12 (T). We assume that i > 0 for i = 1; 2, 1 6= 2 and a 2 H 1 (T) with 0 a(') 1 for all ' 2 T. Our goal is to formulate the identi cation problem for a as an inverse problem F (a) = z; where z is a function on T depending on the data z1 and z2 . However, this cannot be done in a naive way since we have to take care of the fact that there is an in nite
6
W. RING
variety of possible boundary values (z1 ; z2) for a,solution u(a) of u = fa , for one and the same a. Hence a mapping a 7! u(a) 7! 0 u(a); 1 u(a) would be, without additional assumptions, a multivalued function. We therefore seek to construct an invariant z depending on z1 and z2 which is the same function for all pairs (z1 ; z2 ) within the set of all possible boundary values of u(a), i.e. z should not depend on the special choice of z1 and z2 but only on the interface a, and we must be able to construct z from any given pair (z1 ; z2 ). This can be done by means of the following Proposition 3.1. Suppose a 2 H 1 (T) with 0 a(') 1 for all ' 2 T and let u; v 2 H 1 (U ) ful ll u = v = fa on U . Then
1 u , (D 0 u) , 22 = 1 v , (D 0 v) , 22
(23)
in H , 21 (T). Proof. Since w = u , v 2 H 1 (U ) and w = 0 on U , we can apply Proposition 2.1 to obtain D 0 (u , v) = 1 (u , v);
,
and consequently
1 u , (D 0 u) , 22 = 1 v , (D 0 v) , 22
in H , 21 (T). Proposition 3.1 implies that we can take any solution u(a) of u = fa and de ne the mapping
, (24) a 7! F (a) := 1 u(a) , D 0 u(a) , 22 since this function is independent of the special choice of u(a). We can therefore
express (1) in the form:
Solve F (a) = z2 , (Dz1 ) , 2 for a: (25) 2 The 22 -term on the right-hand side of (25) could of course be brought to the left and be included in the de nition of F . We prefer, however, to split it from F in order to obtain a more compact representation of F , as will become apparent later. In order to derive an integral representation for the operator F , we use the (unique) solution ua 2 H 1 (U ) of the Dirichlet problem (
u = fa on U :
0 u = 0
(26)
From the de nition (24) we obtain
F (a) = 1 ua , 22 :
(27)
Note that the boundary values in (26) were chosen such that the (D 0 u) -term vanishes. Our aim is now to calculate the Neumann trace of ua . We de ne PN (; ) =
7
IDENTIFICATION OF A CORE P
jnj in jnjN e for N
2 N and we use the notation PN (; ') = PN (; ' , ). Operators acting on functions with respect to the variable are indicated by the corresponding subscript. Green's formula Z
U
ua () PN (; ') , ua() PN (; ') d = Z
@ u () P (; ') , u () @ P (; ') d a N a @n N @U @ n together with ua = fa , PN (; ') = 0 and 0 ua = 0 implies (1 , 2 ) and hence
(28)
Z
0
2Z a( ) 0
PN (; , ) d d + 2 =
2
Z
0
1 ua ( )
X
jnjN
ein(', ) d ;
1 , 2 Z 2Z a( ) P (; , ) d d N 2 0 0 Z 2 X 1 =
1 ua ( ) e,in d ein' , 22 : 2 jnjN
0
The sum on the right-hand side of (28) is a truncated Fourier series which converges to 1 ua in1 H , 21 (T). (Recall that ua 2 H (; U ) and by the trace theorem 1 ua 2 H , 2 (T).) PN (; ') converges with respect to pointwise on U towards the Poisson kernel P (; '), and it is bounded by the integrable function Q'(; ) = 4=j1 , ei(', ) j. It follows that limit and integral can be transposed in (28) for N ! 1. We can therefore combine (28), (27) and (25) to obtain and
F (a) = 1 2, 2
Z
0
2Z a( ) 0
P (; , ) d d
(29)
1 , 2 Z 2Z a( ) P (; ' , ) d d = z , (Dz ) , 2 , (30) 2 1 2 0 0 2 as an equation for a which is equivalent to (1). Let us return to (28). Carrying out the integration with respect to and taking the limit N ! 1, we nd X 1 , 2 Z 2 a( )jnj+2 F (a) = e,in d ein' : (31) 2 j n j + 2 0 n2Z
From this representation we immediately obtain Z 2 , a( )jnj+2 e,in d ; , 1 2 b F (a) (n) = 2 jnj + 2 0 and therefore, following (30),
1 , 2 Z 2 a( )jnj+2 e,in d = z^ (n) , jnj z^ (n) , 2 2 1 2 0 jnj + 2 2 0;n
(32)
(33)
8
W. RING
for all n 2 Z, where 0;n = 1 if n = 0 and 0;n = 0 else. This is a system of in nitely many (nonlinear) functionals which we have to solve simultaneously for a. Remark 3.2. In view of generalization to domains of arbitrary geometry we shall now give a brief explanation of the above formula (30) from a di erent perspective. Let G(x; ) denote the Green's function for the Dirichlet problem (6) on U , i.e. assume that the solution of (6) can be written in the form
u(x) =
Z
U
G(x; ) fa () d +
Z
@U
u() @@n G(x; ) d
(34)
for all x 2 U . We (formally) take now rst order traces on both sides of (34), where we di erentiate under the integrals on the right-hand side and evaluate at x 2 @U . A rigorous justi cation of these operations will be given at some other place. We obtain Z
@ G(x; ) f () d = @ u(x) , Z u() @ 2 G(x; ) d a @n @ nx @ n U @ nx @U for all x 2 @U . It is well known (c.f. Myint-U [13, p 289]) that G is given by
(35)
G(r; '; ; ) = 41 log[r2 + 2 , 2r cos(' , )] , 41 log[1 + r2 2 , 2r cos(' , )]:
A short calculation then gives
@ G(x; ) @ G(r; '; 1; ) = 1 1 , r2 1 P (r; ' , ): = = 2 @ n 2 1 + r , 2r cos(' , ) 2 2@U @ (36)
For the last equality we refer to Rudin [17, p 233, (2)]. From the symmetry of G it follows that 1 P (; ' , ): @ G(x; ) = @ nx x2@U 2 With this (35) has the form 1 Z 2Z 1 P (; ' , ) f (; ) d d a 2 0 0 Z 2 @ P (r; ' , ) d = 1 u(') , 21
0 u( ) @r r=1 0 Z 2 X = 1 u(') + 21 jnj
0 u( ) e,in d ein' 0 n2N (37) = 1 u(') , (D 0 u) sptilde('): Since the left-hand side of (37) can be written as
1 , 2 Z 2Z a( ) P (; ' , ) d d + 2 ; 2 0 0 2
we have derived (30) from the relation (35). Obviously (35) also makes sense if U is replaced by some arbitrary R2 , hence we can consider (35) as a generalization of (30) to domains di erent from the unit disk, provided that the Green's function for the Dirichlet problem on is known.
9
IDENTIFICATION OF A CORE
4. Smoothness, Di erentiability and Uniqueness. In this section we consider certain properties of the operator F as de ned in (29). We investigate smoothness properties of F (a) (') and calculate the Frechet derivative F 0 (a). The domain for F is de ned as the set dom F = fa 2 H 1 (T) : 0 a(') 1 for all ' 2 Tg:
(38)
Frequently we shall consider the interior (dom F ) of dom F in H 1 (T), given by (dom F ) = fa 2 H 1 (T) : 0 < a(') < 1 for all ' 2 Tg:
(39)
Proposition 4.1. Let F be as de ned in (29) and suppose that dom F and (dom F ) are given as in (38) and (39) respectively. Then
F (a) 2 H l (T) for every a 2 dom F and l < 12 . If a 2 (dom F ) , then F (a) 2 C 1 (T): Moreover, for every a 2 (dom F ) , there exists a neighbourhood U (a) of a in H 1 (T) such that
F : U (a) H 1 (T) ! H l (T) is Lipschitz continuous for every l 2 R, with Lipschitz constant given by k = j1 ,2 j c~, where c~ is some constant depending on U (a) and l. Proof. From (14) and (29) follows
(40) (41)
2
l
Z 2 2 X a( )jnj+2 e,in d , 1 2 2 l (1 + n ) H (T) = 2 jnj + 2 0 n2Z 2l X j1 , 2 j2 ((1jn+j +n2))2 : n2Z
kF (a)k2 l
2
n) l 2l,2 for all l 0 and jnj 1, the series (40) converges if Since ((1+ jnj+2)2 2 jnj 2l , 2 < ,1, i.e. if l < 12 . Let us now suppose that a 2 (dom F ) . Since a 2 H 1 (T) C (T), it follows that sup'2T ja(')j = kak1 < 1. We have
Z
2
jnj+2
a( ) e,in d jnjk F (a) b(n) = jnjk 1 2, 2 jnj + 2 0 k (42) jnjjn+j 2 j1 , 2 j kakj1nj+2 ! 0 as jnj ! 1 since kak1 < 1, and hence F (a) 2 C 1 (T), due to (17). The boundedness of the embedding H 1 (T) ,! C (T) implies that we can nd a neighbourhood U (a) of a 2 (dom F ) in H 1 (T) such that m := supf kbk1 : b 2 U (a) g < 1. We consider F : U (a) H 1 (T) ! C 1 (T) H l (T) ,
10
W. RING
with l 2 R. Then kF (a1 ) , F (a2 )kH l (T) =
X
n2Z
2
Z
2
0
a1 ( )jnj+2 , a2 ( )jnj+2 e,in d jnj + 2
2 1 2
j+1 (1 + n2 )l Z 2 ja ( ) , a ( )j jnX k j n j +1,k d 2 a ( ) a ( ) 1 2 1 2 2 0 n2Z (jnj + 2) k=0 X 1 j1 , 2 j (1 + n2 )l m2(jnj+1) 2 ka1 , a2 k1 k ka1 , a2 kH 1 (T);
j1 2, 2 j (43)
(1 + n2 )l 1 , 2
X
12
n2Z
where
k := j1 , 2 j
X
(1 + n2 )l m2(jnj+1)
n2Z
21
c1 < 1
with c1 such that kf k1 c1 kf kH 1 (T) for all f 2 H 1 (T). The convergence of the series in (43) is easily seen by the quotient criterion. Hence we have proved that F : U (a) H 1 (T) ! H l (T) is locally Lipschitz continuous for every l 2 R. Remark 4.2. If (25) has a solution a 2 (dom F ) , then Proposition 4.1 and (25) imply that although each of the boundary data z2 and (Dz2 ) is only in H , 12 (T) (and, in general, will not be in any space H l (T) with l > , 21 ) their di erence is in C 1 (T). We also point out that in Proposition 4.1 we have also proved the Lipschitz continuity of F with respect to the C (T)-norm on (dom F ) . Proposition 4.3. F : (dom F ) H 1 (T) ! H l (T), with l 2 R, is Frechet di erentiable at every point a 2 (dom F ) with derivative given by [F 0 (a) h](') = 1 2, 2
2
Z
0
a( ) P (a( ); ' , ) h( ) d
(44)
for all h 2 H 1 (T). The Fourier transform of F 0 (a) h is given by
[F 0 (a) h]b(n) = 1 , 2
2
Z
2
0
a( )jnj+1 h( ) e,in d
(45)
for h 2 H 1 (T) and n 2 N . Proof. First of all we show that the Fourier transform of (44) is in fact given by (45). Since for a 2 (dom F ) , 2 P a( ); ' , (46) 1 , kak1 < 1; we can exchange summation and integral in (44) and obtain
F 0 (a)h = 1 2, 2
XZ
n2Z 0
2
a( )jnj+1 h( ) e,in d ein'
11
IDENTIFICATION OF A CORE
which proves that (45) is correct. To prove Frechet di erentiability it is sucient (c.f. Berger [1, p 68, Th 2.1.13]) to ensure that 1. F 0 (a) : H 1 (T) ! H l (T) is a bounded linear Operator. 2. F 0 depends locally Lipschitz continuously on a, i.e. there exists an L > 0 (depending on a) such that
0
F (a) , F 0 (b) (47) L(H 1 (T);H l(T)) Lka , bkH 1 (T) for all b in some neighbourhood of a. (This is even more than required in [1].) 3. F 0 (a) is the Gateaux derivative of F at a, i.e. lim kF (a + #h) , F (a) , # F 0 (a) hkH l (T) = 0; #!0
where # 2 [0; 1) for every h 2 H l (T). We start with 1. and show that F 0 as de ned in (44) is a bounded linear operator from H 1 (T) into H l (T) for some arbitrary but xed l 2 R. Using (45) we nd
kF 0 (a) hkH l (T) =
X
(1 + n2 )l 1 , 2
2
n2Z
X
(1 , 2 ) c
,X
2
Z
0
2 12 j n j +1 , in a( ) h( ) e d
jnj+1) jh^ (n)j2 (1 + n2 )l kak2( 1
12
n2Z
1
jh^ (n)j2 2 = c khkL2(T)
n2Z
where c is chosen large enough such that (1 , 2 ) (1 + n2 ) 2l kakj1nj+1 c for all n 2 Z. (Recall that a 2 (dom F ) and therefore kak1 < 1.) Hence F 0 (a) : H 1 (T) L2 (T) ! H l (T) is a bounded linear operator. We note that this is a stronger result than required in 1. and that F 0 (a) can be uniquely extended to a bounded linear operator from L2 (T) into H l (T). For the proof of 2. we choose a neighbourhood V (a) of a in H 1 (T) such that kbk1 m < 1 for all b 2 V (a). Then we have kF 0 (a) h , F 0 (b) hkH l (T) =
X
(1 + n2 )l 1 , 2
n2Z
(1 , 2 )
2
X
2 ,
Z
0
a( ) , b( )
jnj
X
k=0
a( )k b( )jnj,k h( ) e,in
d
2 12
1
(1 + n2 )l (jnj + 1)2 m2jnj jh^ (n)j2 2 ka , bk1:
n2Z
If we choose constants c1 , c2 such that ka , bk1 c1 ka , bkH 1(T) for all b 2 V (a) and for all n 2 Z, we obtain
(1 + n2 )l,1 (jnj + 1)2 m2jnj c22 X
kF 0(a) h , F 0 (b) hkH l (T) (1 , 2 ) c1 c2 ka , bkH 1 (T) = L ka , bkH 1 (T) khkH 1 (T);
(1 + n2 ) jh^ (n)j2
n2Z
12
12
W. RING
with
L = (1 , 2 ) c1 c2 and consequently
(48)
F 0 (a) , F 0 (b) L(H 1 (T);H l(T)) L ka , bkH 1 (T)
for all b 2 V (a). It remains to prove that F 0 is the Gateaux derivative of F . We nd
kF (a + #h) , F (a) , # F 0 (a) hkH l (T) X 2 = (1 + n2 )l 1 , 2 n2Z Z
2
a( ) + #h( ) jnj+2 , a( )jnj+2 jnj+1 h( ) e,in d , # a ( ) jnj + 2
2 ,
0
= 1 2, 2 X
n2Z
# (1 , 2 )
2
Z
(1 + n2 )l X
0
j+2 , i 1 h jnX jnj+2#k h( )k a( )jnj+2,k e,in d jnj + 2 k=2 k
(1 + n2 )l
n2Z
2 12
1
a + p#jhj
jnj+2 2 1 jnj + 2
12
2 12
;
k where we p used #k,1 # 2 for all k 2 and # 2 (0; 1). The last series above converges if ka + # jhj k1 < 1 by the quotient criterion, hence the whole expression tends to 0 if # ! 0. We now address the problem of uniqueness of a solution a 2 dom F of F (a) = z for z 2 rg F . In Section 1. we proved that every solution a 2 C (T) with 0 a(') 1 is also a solution of the equivalent inverse potential problem (3). Hence, to prove injectivity of the operator F , it is sucient to prove unique solvability of (3). The question of uniqueness for (3) is a well investigated topic in the literature. It was rst shown by Novikov [14] that ~ in (2) is uniquely determined by the potential uj@ in the class of star-shaped open sets with smooth boundary. Note that a = f(r; ') : r < a(')g is star shaped with respect to the origin. A proof of a unique-solvability result for (2) in Rd (d 2) for star-shaped regions can be found in [9, p 55, Th 3.1.1]. The following theorem is an immediate consequence of this result. Theorem 4.4 (Uniqueness). Suppose that a; b 2 dom F , satisfying a(') > 0 or b(') > 0 for all ' 2 T and suppose that
F (a) = F (b) holds in H l (T) with l < 21 . Then
a=b in H 1 (T). Next we shall prove the injectivity of F 0 .
13
IDENTIFICATION OF A CORE
Theorem 4.5. The operator F 0 (a) : H 1 (T) ! H l (T) with l 2 R is injective for all a 2 (dom F ) . Proof. Suppose that f 2 ker F 0 (a) H 1 (T), i.e. (c.f. (45) ) Z
0
2
a(')jnj+1 f (') ein' d' = 0 for all n 2 Z:
(49)
We shall rst prove that (49) implies 2
Z
0
,
a(') f (') h a('); ' d' = 0
(50)
for all functions h which are harmonic on some neighbourhood V of a = (r; ') : r a(') . For this purpose x h harmonic on V and let hc : V ! R be a harmonic conjugate of h, i.e. suppose hc is chosen such that the function H (x + iy) = h(x; y) + ihc(x; y) is analytic on V . Here we have changed for the moment from the usual polar coordinates (r; ') to Cartesian coordinates (x; y). For a proof of the existence of a harmonic conjugate to a given harmonic function see Kellogg [10, p 345, Th III]. H can be approximated by polynomials
Qn (z ) =
mX (n) k=0
ck;n z k ;
with z = x + iy, ck;n 2 C and m(n) 2 N , uniformly on a by Runge's theorem (c.f. [17, p 271, Th 13.7]). Then
We set and
(n) m X , Re Qn (z ) , H (z ) = 12 ck;n z k + ck;n z k , h(z ): k=0
(51)
d0;n = 21 (c0;n + c0;n ) dk;n
(
1 2 ck;n 1 ck;n 2
for k 1 : for k ,1
With z = r ein' we can write (51) as Re (Qn , H ) =
mX (n) k=,m(n)
dk;n rjkj eik' , h(r; '):
(52)
Since jRe(z )j jz j for every complex number z , (52) implies that h can be approximated by a sequence of `harmonic polynomials'
qn (r; ') =
mX (n) k=,m(n)
dk;n rjkj eik'
14
W. RING
uniformly on a . Hence we obtain Z
2
0
,
a(') f (') h a('); ' d'
Z
2
0
(53)
,
,
a(') f (') h a('); ' , qn a('); ' d'
+
mX (n)
2
Z
k=,m(n)
dk;n
0
a(')jnj+1 f (') ein' d' :
The second term on the right-hand side of (53) is 0 due to (49) and the rst one tends to 0 since pn ! h uniformly on @ a . Hence the left-hand side of (53) must be 0 and (50) is proved. We consider now the Dirichlet problem (
w = 0 on
a , w @ a = w a('); ' = f (') :
(54)
since a 2 (dom F ) , there exists > 0 satisfying a(') for all ' 2 T. This together with f 2 H 1 (T) C (T) implies that f : @ a ! C in (54) is a continuous function on @ a . From this and from geometric properties of the boundary @ a (c.f. [6, p 336, Th 1 and p 340, Prop 4]) the existence of a classical solution w for (54) follows, that is to say w = 0 on a and lim w(x) = f (z ) x!z ,
for x 2 a and z 2 @ a. We de ne wk (r; ') = w ,(1 , k1 )r; ' for k 2 N . From the properties of a classical solution it follows ,that wk a('); ' ! f (') pointwise on T and by the maximum modulus principle w a('); ' is bounded by jf (')j. Hence, by ,k the dominated convergence theorem, wk a('); ' ! f (') in L1 (T). Moreover, wk is harmonic on Vk = (1 + k1 ) a . Employing (50) and the last arguments we obtain 0=
2
Z
0
,
a(') f (') wk a('); ' d' !
2
Z
0
a(') jf (')j2 d' 2 kf k2L2(T);
where 0 < a(') for all ' 2 T. Hence we nd f = 0 in H 1 (T). R Remark 4.6. The fact that 02 a(') f (') a(')jnj ein' d' = 0 for all n 2 Z implies that f = 0 in H 1 (T) is equivalent to the assertion that the functions fa(')jnj ein' gn2Z are dense in H 1 (T) (and hence in L2 (T)). It would be interesting to prove this consequence of Theorem 4.5 directly and not by solving boundary value problems which seems somewhat unnatural in this context. 5. Ill-Posedness, Least-Squares Formulation, and Regularization. In this section we consider the problem F (a) = z; (55) within the concept of least-squares methods. Solutions of (55) do not depend continuously on the data z , not even if we replace the parameter space H 1 (T) by L2(T) as the following example shows. We consider the sequence of functions fan gn2N dom F de ned by an (') = 12 + 2 cos n' = 12 + (ein' + e,in')
15
IDENTIFICATION OF A CORE
with 0 < < 14 and the constant function
r
a(') = 14 + 2 2: It is easy to see that r
2 kan , ak2L2(T) = 21 , 14 + 2 2 + 2 2 = const > 0; and hence an 6! a for n ! 1 in L2 (T). On the other hand we nd Z 2 h , , 1 1 + (ein + e,in )2 , 1 , 2 2i d = 0 1 2 b F (an ) , F (a) (0) = 2 2 2 4
0
for all n 1, and (recall that F (a) is constant) ,
Z
F (an ) , F (a) b(k) = 1 2, 2
2
0
1 1 + (ein + e,in )jkj+2 e,ik d jk j + 2 2
for n 1. The term 12 + (ein + e,in ) jkj+2 can be written in the form P all k 6= 0 and imn with coecients cm;n 2 C , cm;n 6= 0 only for nitely many m 2 Z. m2Zcm;n e , ik 2 Hence the orthogonality of the functions e in L (T) implies that F (an ) b(k) = 0 for all k 2 Z with 1 jkj < n. For jkj n we have
,
,
F (a
Z 2 , 1 1 + 2 jkj+2 d = 1 , 2 qjkj+2 1 2 k 2 jk j + 2 0 jk j + 2 2
n ) b(
)
where q = 21 + 2 < 1. We therefore obtain
kF (an ) , F (a)k2H l (T) =
X
,
,
k2Z
(1 , 2 )
X
jkjn
jkjn
jkj+2
(1 + k2 )l jqkj + 2 :
Using the quotient criterion it is easy to see that the series converges. Thus we have X
(1 + k2 )l F (an ) b(k) , F (a) b(k)
2 l qjkj+2 k2Z(1 + k ) jkj+2
P
jkj+2
(1 + k2 )l jqkj + 2 ! 0
for n ! 1 and consequently F (an ) ! F (a) in H l (T). Since the kth derivative Dk is a continuous operator from H k+1 (T) into H 1 (T) it follows that Dk F (an ) ! Dk F (a) in H 1 (T) and due to the compact embedding H 1 (T) ,! C (T) we have
Dk F (an ) ! Dk F (a) uniformly on T for all k 2 N . From this we can conclude that F (an ) ! F (a) even in the topology on C 1(T).
16
W. RING
The lack of continuous dependence of the interface a on the observation z in (55) is the reason why we consider the inverse problem of determining a from z as `illposed' (in the sense of Hadamard). An ecient and well investigated tool to deal with ill-posed problems is Tikhonov's regularization method. (c.f. Engl, Kunisch and Neubauer [8], Seidman and Vogel [18], the work [12].) The method consists in approximating solutions of (55) by minimizers of the Tikhonov functional n
o
min kF (a) , z k2H , 12 (T) + kD(a , a )k2L2 (T) ; a2dom F
(56)
where D : H 1 (T) ! L2 (T) and is some (small) regularization parameter. The function a in the regularization term plays the r^ole of an `origin' with respect to which the semidistance kD( , a )kL2 (T) is measured. It should be chosen from some a-priori information close to the actual solution of (55). In the following we assume that a 2 dom F H 1 (T) and z 2 H , 12 (T). Although we could use any H l (T)-space with l 2 (,1; 12 ) as output-space (c.f. Proposition 4.1), we consider only the case where l = , 21 . This is motivated by the fact that the data on the right-hand side of (55) have the form z = z2 , (Dz1 ) , 22 , where the functions z2 and (Dz1 ) are in H , 12 (T), only their di erence is smoother, e.g. in C 1 (T) if a 2 (dom F ) . If (Dz1 ) and z2 result from real data, we cannot expect the di erence z2 , (Dz1 ) to be reasonable accurate in the C 1 (T)-topology or in any H l (T)-norm with large l but only in the topology which is `natural' to both terms 12 , of the di erence, which is H (T1). Moreover, we assume that z is given by 1some inaccurate measurement z 2 H , 2 (T) of `unperturbed' data z0 2 rg F H , 2 (T), with kz0 , z kH , 21 (T) ; (57) where we suppose that a0 2 (dom F ) exists, satisfying F (a0 ) = z0 . By Theorem 4.4, a0 is then uniquely determined by z0 . In [12] conditions were obtained which ensure solvability, stability with respect to noisy data, and convergence of solutions of (56) | with a certain convergence rate | towards a solution of (55) if the noise level and the parameter go to 0. Applying the results [12] leads to the following 1 Theorem 5.1. Assume that F : dom F H 1 (T) ! H , 2 (T) is given as in (29) and (38) respectively, and a 2 H 1 (T), z 2 H , 12 (T) and > 0. Then there exists a solution a = a(1 ; z ) 2 dom F of (56). Solutions of (56) are stable in the sense that for zn ! z in H , 2 (T), the sequence an = a( ; zn ) has a H 1 (T)-convergent subsequence, and the limit of every weakly convergent subsequence of fang is solution of n
o
min kF (a) , z k2H , 12 (T) + kD(a , a )k2L2 (T) : a2dom F If z0 = F (a0 ) and kz0 , z kH , 12 (T) and if we choose = () then a( (); z ) ! a0 in H 1 (T). Moreover, if there exists ! 2 H 12 (T) such that
,
D2 (a0 , a ) (
or equivalently
X
Z 2 , , 1 2 ) = 2 a0 ( ) P a0 ( ); ' , 0
n2 (a0 , a )b(n) ein = (1 , 2 )
n2Z
X
!(') d'
!^ (n) a0 ( )jnj+1 ein :
n2Z
(58) (59)
17
IDENTIFICATION OF A CORE
and if k!kH 12 (T) is bounded by a certain constant M depending on a0 , then
p ka( (); z ) , a0 kH 1 (T) = O( ) and kF a( (); z ) , z0 kH , 12 (T) = O():
Remark 5.2. A detailed description of the dependence of the constant M on the solution a0 is ,given inthe Appendix. Condition (58) is in fact the explicit form of D D(a0 , a ) = F 0 (a0 ) !. We demonstrate that (58) includes certain smoothness , requirements on the di erence a , a0 . From (59) we see that the range of F 0 (a0 ) consists of functions whichPcan be expressed as series in terms of a0 ( )jnj+1 ein with coecients cn satisfying n2Zjnj jcn j2 < 1. Using the constant function 1 ,1 2 we , nd that a0 is always in the range of F 0 (a0 ) . It is easily, seen that the series 0 on the right-hand side of (59) converges uniformly, therefore F (a0 ) g 2 C (T) for all g 2 H 12 (T). Thus, if (58) is ful lled, we have D2 (a0 , a ) 2 C (T) and hence a0 , a 2 C 2 (T). This, however, is only the minimal smoothness requirement as the following example shows. Suppose that a0 = q < 1 is constant. Then X F 0 (a0 ) g = (1 , 2 ) g^(n) qjnj+1 ein 2 C 1 (T)
,
n2Z
for all g 2 H 21 (T). Hence in this case we have rg F 0 (a0 ) C 1 (T). 6. Numerical Results. As a starting point for the numerical implementation we use the system of equations (33). We introduce the set SN of real-valued trigonometric polynomials of degree less or equal than d = N2,1 given by SN = f N =
d X
j =,d
,
cNj eij' : cN0 2 R; cNj 2 C ; cNj = cN,j for j 6= 0
with odd N 2 N . We use the abbreviations 'Nk := 2N k and fkN := f N ('Nk ) for f N 2 SN and k = 0; ; N , 1. Given f N 2 SN there is a one-to-one correspondence between the Fourier-coecients fcN0 ; ; cNd g and the interpolation points ff0N ; ; fNN,1 g via the discrete Fourier transform (c.f. [5, pp 60, 61]), hence we can use ff0N ; ; fNN,1g as coordinate-vector for f N 2 SN instead of fcN,d; ; cNd g. In the following it is assumed that F acts on SN , i.e. we consider F (aN ) and F (aN ) b(n) for aN 2 SN . Since the exact evaluation of the integral in (33) for aN 2 SN and large values of jnj is very complicated, we used a numerical quadrature by means of the rectangular rule for the computation of the integral F (aN ) b(n) ' 1 , 2
N
NX ,1 (aN )jnj+2 k e,in'Nk : j n j + 2 k=0
(60)
Note that the fact that aN (0) = aN (2) implies that the rectangular rule and the trapezoid rule give the same numerical quadrature formula (60). We also mention that the dimension N of SN also de nes the mesh-width and therefore also the accuracy of the numerical integration (60). It follows from the de nition that the operator F depends linearly on the di erence = 1 , 2 , hence (56) can be transformed into h
i
min 2 kFe (a) , ze k2H , 12 (T) + e kD(a , a )k2L2 (T) ; a2dom F
18
W. RING
where Fe = 1 F , ze = 1 z and e = 2 are independent of . Therefore, if 6= 0, the optimization problem (56) depends on only as a scaling factor. This justi es to x
1 = 2 and 2 = 1 for all numerical examples. In order to obtain numerical data for the right-hand side of (33) from some known interface function a, we also used the numerical integration (60). We set NX d,1 (aNd )jnj+2 ,in'Nk d =: z^Nd (n); k e d k=0 jnj + 2 for all jnj Nd2,1 , and z^Nd (n) = 0 for jnj > Nd2,1 . For this (forward) calculation large values for Nd (usually Nd = 3003) were used to obtain reasonable accurate `real' data, while for the calculation of faNk gNk=0 from the data z^Nd , mostly N = 91 was chosen. The accuracy of the data was then arti cially destroyed by adding a random noise zNd , where
z^(n) ' 1N, 2
zNd ('Nk d ) = Rand(0; 1) , 12 h
i
with Rand(0; 1) a random number in [0; 1] and the noise level is chosen as one of the values in f0:05; 0:01; 0:005; 0:001; 0:0005g. All interface functions under consideration were centered around r = 0:8. This gave data z with kz k1 0:32. Hence the noise level varied between 0.075% and 7.5% with respect to the kk1 -norm. Moreover, the random function zNd was the same for all interface functions. A di erent possibility to construct data z from a known interface a would be to solve the forward problem u = fa by means of the fundamental solution F (x; ) = , 21 log jx , jRand to evaluate both Dirichlet and Neumann boundary data of the solution u(x) = F (x; ) fa() d at the points xk = (cos 'Nk d ; sin 'Nk d ). We could then use the discrete Fourier transform to obtain the Fourier transform of the data z = z2 , (Dz1 ) , 22 . With these approximations for the operator F and the data z , we can enter (56) and calculate the interface a as the minimizer of the approximated Tikhonov functional
J N (aN ) = (61)
jnj M2,1
+
, 2 (1 + n2 ), 12 1 N
X
X
jnj M2,1
NX ,1 (aN )jnj+2 ,in'Nk k e k=0 jnj + 2
2
, z^(n)
jnj2 a^N (n) , a^ (n) 2 :
The rst thing to come to mind is perhaps to choose M = N in (61). With this choice it was observed in the numerical computations that large values of M increased the computation time considerably. We were therefore bound to x M somewhere below 100. However, it turned out that (for xed M ) the value of N (the number of interpolation points) could be increased beyond M without increasing the computational error and without loosing too much computation time. This may be due to
19
IDENTIFICATION OF A CORE
the fact that [F (a)]b(n) decreases exponentially as jnj ! 1 (c.f. (42) with k = 0), and hence the in uence of the terms for large values of jnj to the minimum of (61) (the stabilized problem!) is small. For this reason we used `more variables than equations', in fact we xed N = 91 and M = 45 for all computations. Moreover for 0 < n M2,1 we have F N (aN ) b(,n) = F N (aN ) b(n), z^N (,n) = z^N (n) and (aN , a )b(,n) = (aN , a )b(n); hence each of the terms
,1 (aN )jnj+2 1 , 2 NX ,in'Nk , z^(n) 2 k e N k=0 jnj + 2
and N ^ (
a n) , a^ (n) 2
occurs twice in (61). If we identify aN 2 S N with its coordinate vector (aN0 ; aNN ,1 ) 2 RN and if we ignore the constraints 0 aN 1 we obtain for (61)
min
aN 2RN
+2 (62)
M2,1 X
n=1
,1 jaN j2 1 , 2 NX k , z^(0) 2 N k=0 2
, 2 (1 + n2 ), 12 1 N
M2,1
+2
X
n=1
NX ,1 jaN jn+2 ,in'Nk k e k=0 n + 2
jnj2 a^N (n) , a^ (n) 2
2
, z^(n)
:
This optimization problem was implemented using the numerical package Matlab and the optimization routine leastsq provided there, which is based on the Levenberg-Marquard algorithm for nonlinear least-squares problems. The Fourier transform (aN , a )b(n) occurring in the regularization term was computed by means of Matlab's fast Fourier transform routine fft. Formulation (62) also allows to calculate an analytic gradient which was used in the numerical computations. As a rst example we consider the interface given by the function
ac (') = 0:8 + 0:1 cos ';
(63)
which we named `Cosine'. The a-priori guess a waschosen 0. In Figure 1 the true ,1 of (62) as small interface ac as a solid line, and the solution-vector aNk (z; ) Nk=0 crosses, are plotted vs the angle '. The noise level was chosen to be 0 and the regularization parameter = 10,11 was (experimentally) determined such that the error
kaN (z; ) , ackL2 (T) '
X (
jnjN
21
aN (z; ))b(n) , a^c (n) 2
is minimal. This fact is indicated in the plot by the word `optimal'. Also the absolute errors with respect to the L2 (T)- and the H 1 (T)-norm are given. Figure 2 shows the result for the same interface function ac with perturbed data z N + zN , where = 0:05 (i.e. 7.5% noise). As in the previous example the regularization parameter
20
W. RING true interface (-) and solution of the regularized problem (+) 0.95
0.9 +++++++ 0.85
r
0.8
0.75
0.7
0.65
+++++ ++ ++ ++ ++ + ++ ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ++ + ++ ++ + guess=0 ++ + +++ ++ ++++++++
number of points=91 number of equations=45 regularization parameter alpha=1e-11, optimal H(-0.5)-noise=0, L2-noise=0 L2-error=2.907e-05, H1-error=9.476e-05
0.6 0
1
2
3
4
5
6
phi
Fig. 1. Cosine without Noise true interface (-) and solution of the regularized problem (+) 0.95
0.9 ++++++ 0.85
r
0.8
0.75
0.7
0.65
++ +++ ++ ++ ++ ++ + ++ ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ++ + ++ ++ ++ + guess=0 + ++ ++++++++++++
number of points=91 number of equations=45 regularization parameter alpha=0.002, optimal H(-0.5)-noise=0.001021, L2-noise=0.0143 L2-error=0.001338, H1-error=0.003781
0.6 0
1
2
3
4
5
6
phi
Fig. 2. Cosine with Noise ( = 0:05)
is optimal. Note that the optimal regularization parameter increases for growing noiselevel. It is also worth mentioning that the di erence between the L2 (T)-error and the H 1 (T)-error is very small which can be attributed to the fact that the representation of ac as a Fourier series has only two modes (ac 2 S3 ). An error in the higher modes is therefore not likely to occur since in (62) they are penalized by the regularization term with comparably large . The four plots in Figure 3 show numerical experiments with the interface function (
0:8 for ' 2 [0; , 1 ] [ [ + 1 ; 2] ah(') = 0:8 , 0:27 exp , 2 1 2 for ' 2 ( , 1 ; + 1 ) : 1,4(',)2 2 2
(64)
21
IDENTIFICATION OF A CORE
We called this function `Hole' and we investigated the behaviour of solutions of (62) for di erent noise-levels = 0; = 0:0005; = 0:005 and = 0:05. According to Theorem 5.1 the regularization parameter was chosen of order ; i.e. we choose = 10,8; = 10,7 and = 10,6 for the values 6= 0, where the value = 10,6 for = 0:05 is optimal. For this interface the optimal parameter is considerably smaller than for Cosine (with the same noise level), which is due to the fact that higher frequencies are needed to represent ah , hence it is not good to suppress them by a strong regularization term. true interface (-) and solution of the regularized problem (+)
true interface (-) and solution of the regularized problem (+)
0.95
0.95
0.9
0.9
0.85
0.85
0.8 ++++++++++++++++++++++++++++++++++++++++
+++++++++++++++++++++++++++++++++++++++
r
+
0.75
+
0.65
0.6 0
2
+
0.7
0.65
3
4
5
6
0.6 0
1
2
3
true interface (-) and solution of the regularized problem (+)
0.9
0.9
0.85
0.85
0.8 +++++++++++++++++++++++++++++++++++++++
++++++++++++++++++++++++++++++++++++++ +
r
+
0.6 0
2
++ ++ +++++++++++++ ++++++++++++++++++ ++ + + +
+ +
0.75
3
0.7
0.65
4
5
6
0.6 0
+
+ +
+
+ guess=0 + + + ++ number of points=91 number of equations=45 regularization parameter alpha=9e-07 H(-0.5)-noise=0.0001021, L2-noise=0.00143 L2-error=0.002288, H1-error=0.03717
1
++++ + +
0.8 ++++++++++++++++++++++++++++++++++
+
+
+
6
true interface (-) and solution of the regularized problem (+) 0.95
+
5
Noise ( = 0:0005)
0.95
0.65
4
phi
no Noise
0.7
+
guess=0 + + + + ++ number of points=91 number of equations=45 regularization parameter alpha=9e-08, optimal H(-0.5)-noise=1.021e-05, L2-noise=0.000143 L2-error=0.001425, H1-error=0.02805
phi
0.75
+
+
0.75
+
guess=0 + + + + ++ number of points=91 number of equations=45 regularization parameter alpha=1e-11 H(-0.5)-noise=0, L2-noise=0 L2-error=0.0004521, H1-error=0.01477 1
++++++++++++++++++++++++++++++++++++++ + +
+
+
+
0.7
0.8 ++++++++++++++++++++++++++++++++++++++++
+
+
+ + guess=0 + + ++ number of points=91 number of equations=45 regularization parameter alpha=9e-06 H(-0.5)-noise=0.001021, L2-noise=0.0143 L2-error=0.005702, H1-error=0.05715
1
2
phi
3
4
5
6
phi
Noise ( = 0:005)
Noise ( = 0:05)
Fig. 3. Hole with Increasing Noise and
In Table 1 and Table 2 we listed the H , 21 (T)-noise , the regularization parameter 2, the H 21 (T)- and the L2 (T)-error and ~ of the computed interface, and the rates and ~ for the preceding examples Cosine and Hole. The noise-levels in the H , 12 (T)-norm correspond to values for 2 f0p :05; 0:01; 0:005; 0:001; 0:0005g. Note that there seems to be a convergence rate for Cosine, while for Hole we can not observe this convergence rate neither in H 1 (T) nor in, L2 (T).Most likely this can be attributed to a lack of containment of D Dah in rg F 0 (ah ) (c.f. Theorem 5.1
22
W. RING
and Appendix, condition (R4) ).
H 21 -noise Reg. parameter H 1 -error 2 = 0.001021 0.0002043 0.0001021 0.00002043 0.00001021
0.002 0.0004 0.0002 0.00004 0.00002
0.003781 0.001556 0.0009800 0.0003477 0.0002492
0.0140 0.0119 0.0094 0.0059 0.0061
L2 -error ~
0.001338 0.0003865 0.0002247 0.00006623 0.00004551
~2 =
0.0018 0.00073 0.00049 0.00021 0.00020
Table 1
Error Rates for Cosine
H , 12 -noise Reg. parameter H 1 -error 2 = 0.001021 0.0002043 0.0001021 0.00002043 0.00001021
0.000009 0.0000018 0.0000009 0.0000018 0.00000009
0.05715 0.04042 0.03717 0.03089 0.02804
L2-error ~ ~2 =
3.199 7.997 13.532 46.705 77.007
0.005702 0.002681 0.002288 0.001685 0.001425
0.0318 0.0352 0.0513 0.1390 0.1989
Table 2
Error Rates for Hole
true interface (-) and solution of the regularized problem (+) 0.95
0.9
+ +++ + ++ + +++ + + + +++ ++ + +++ + 0.85 +++ +++ + ++ + + ++ +++ + +++ + 0.8 + +
r
+++ +++ ++ +++ + + ++ ++ +++ ++ + +++ + +++ ++ + +++ +++ + ++ ++
0.75
0.7
0.65
0.6 0
guess=0 number of points=91 number of equations=45 regularization parameter alpha=1e-11, optimal H(-0.5)-noise=0, L2-noise=0 L2-error=0.0007609, H1-error=0.0248 1
2
3
4
5
6
phi
Fig. 4. Zigzag without Noise
For the plots in Figures 4{7 we considered the piecewise linear interface function
23
IDENTIFICATION OF A CORE true interface (-) and solution of the regularized problem (+) 0.95
0.9
r
++++++ ++ ++ + ++++ + ++ + +++ + + 0.85 ++ ++ + ++ + ++ ++++++ + + ++ 0.8 ++ +
++ ++ ++++ +++ ++ + + +++ ++++ + ++++ +++ + + + + + + + + +++++
+
0.75
0.7
0.65
guess=0 number of points=91 number of equations=45 regularization parameter alpha=7e-05, optimal H(-0.5)-noise=0.001021, L2-noise=0.0143 L2-error=0.004186, H1-error=0.04132
0.6 0
1
2
3
4
5
6
phi
Fig. 5. Zigzag with Noise ( = 0:05)
`Zigzag', de ned by
8 4 >
Abstract. The problem of determining the interface separating regions of constant density from boundary data of a solution of the corresponding potential equation is considered. An equivalent formulation as a nonlinear integral equation is obtained. Fourier methods are used to analyze and implement the problem. Numerical experiments based on a regularized least-squares method are presented. Key words. inverse source problem, nonlinear integral equation, Fourier transform, Tikhonov
regularization
AMS subject classi cations. 31B20, 35J05, 35R25, 35R30
1. Introduction. Suppose (r; ') denote polar coordinates in the plane and r = a(') is a continuous function, with graph of a contained in some domain R2 .
The subject of this work is the following problem:
Determine the interface r = ( a(') in u = fa (r; ') = 1 for 0 r < a(') 2 for a(') r 2 on R from known boundary data z1 = uj@ and z2 = @u @ n j@ .
(1)
We consider only the situation where = U , the open unit disk in R2 . The following introductory remarks, however, do not only apply to U , but to any simply connected open subset of the plane with smooth boundary. Problem (1) occurs in the context of recovering an unknown mass distribution inside a planar body from measurements of the potential at the boundary. Suppose a mass distribution consisting of two layers is given on : a region of constant density 1 inside (a core) and an annular region (a mantle) of constant density 2 enclosing the core, and suppose that the interface between the two regions is given by the curve r = a('). The gravitational potential u generated by this mass distribution is then a solution of the Poisson equation u = fa . Therefore (1) can be interpreted as the problem of nding the curve which separates the regions of constant density from measurements of the potential at the boundary of . Closely related to (1) is the inverse domain recognition problem in potential theory: Find the domain ~ in known g : @ ! R.
Z
~
E (x; ) d = g(x); x 2 @ , from
(2)
Here ~ denotes some open subset of and E (x; ) = 41 log jx , j2 : This work was supported by Fonds zur F orderung der Wissenschaftlichen Forschung, Austria, under P 7869-Phy and Bundesministerium fur Wissenschaft und Forschung, under project `Mathematische Optimierung und Inverse Prozesse' y Institut f ur Mathematik, Technische Universitat Graz, Kopernikusgasse 24, A-8010 Graz, Austria 1
2
W. RING
In fact, (1) and (2) are equivalent in the following sense. Suppose (1) holds for
a, z1 and z2 . Using Green's formula and well known identities from the theory of
boundary integral equations (c.f. Chen and Zhou [4, p 223, (6.43) and p 224, (6.45)]), we obtain Z Z Z fa() E (x; ) d = z1 (2x) , z1 () @@n E (x; ) d + z2 ()E (x; ) d
@
@
(3) and hence, with a denoting the core enclosed by a, we nd Z Z n E (x; ) d = ,1 z1 (2x) , z1 () @@n E (x; ) d 1 2 @
a (4)
Z
+
@
z2()E (x; ) d , 2
Z
E (x ) d
o
for x 2 @ . This shows that a (and hence a) can be obtained from z1 and z2 as a solution of (2) with g given by the right hand side of (4).R On the other hand, the Newtonian potential w(x) = a E (x; ) d satis es w = fa for 1 = 1, 2 = 0. We set z1 := g = wj@ and z2 := @w @ n j@ . Starting with (2), only z1 but not z2 is known. However, if we apply Green's formula and well known integral identities [4, p 223, (6.43) and p 224, (6.45)] to w() and E (x; ), we obtain the integral equation Z Z z2 () E (x; ) d = z1 (2x) + z1 () @@n E (x; ) d for x 2 @ : (5) @
@
The right-hand side of (5) de nes a linear operator L1 : H , 12 (@ ) ! H 21 (@ ) acting on z2 , which is Fredholm of index 0. In case = U , the kernel on L1 consists of the constant functions on @ , hence z2 is determined by (5) up to an additive constant c. (c.f. [4] for the facts mentioned above.) Using Gauss' divergence theorem, we obtain the additional relation Z
Z
Z @w z2 () d = () d = fa () d; @
@ @n
from which the constant c can be uniquely determined. Hence, we arrive at (1) with known boundary data z1 and z2 and unknown a. The domain recognition problem (2) was studied extensively in the past. We mention Novikov [14], Prilepko [15] [16], Cabayan and Belford [3], Brodsky [2], Strakhov and Brodsky [19], the book by Isakov [9] and the literature cited there. The starting point of yet another equivalent formulation is the boundary value problem (
u = fa on : u j@ = z 1
(6)
With u = u(a; z1) denoting the solution of (6), we can de ne the operator and we could solve
G(a; z1 ) = @@n u(a; z1) @
G(a; z1 ) = z2
(7)
IDENTIFICATION OF A CORE
3
for a. A disadvantage of (7) arises if we have to deal with noisy data. Suppose z1 and z2 are replaced by inaccurate measurements z1i and z2i . Then (7) becomes G(a; z1i ) = z2i : This means that we have to deal not only with perturbed data on the right-hand side but also with a perturbed operator Gi = G( ; z1i ), which makes the analysis of the problem rather complicated. A completely di erent, optimization theoretic approach to a problem similar to (1) is presented in Kunisch and Pan [11]. There the augmented Lagrangian method was used to obtain a formulation of the identi cation problem as a constrained minimization problem. In this paper we present an approach to (1) which is very similar to the considerations which led to (4), with the di erence that we use the Green's function for the Dirichlet problem on U instead of the fundamental solution E (x; ). As in (4) we obtain a nonlinear operator equation F ( a) = z (8) where F is an integral operator over the variable domain a and z = z (z1; z2 ). This method is restricted to those domains for which a Green's function is known. The loss of generality is however rewarded by a very compact form of the resulting integral equation for a. Moreover the choices = U and a = a(') make it possible to use Fourier analysis, which will turn out to be a very powerful tool for analyzing and implementing the problem. We show that (8) can be reduced to an in nite system of nonlinear functionals of a, attaining prescribed values which depend on the data z1 and z2 . We derive results on smoothness, di erentiability, injectivity and ill-posedness for (8). To deal with the ill-posedness of the problem, we approximate (8) by a regularized least-squares problem, and we check conditions which ensure stability and convergence for the approximating problem. Finally we present numerical results for di erent interfaces and di erent choices of the regularization term. 2. Preliminaries. We now x the terminology for the subsequent sections and we provide some basic technicalities which are needed henceforth. The notation (r; ') and (; ) is used for polar coordinates of points x 2 R2 and 2 R2 respectively. U = B (0; 1) denotes the open unit disk in the plane R2 . H 1 (U ) is the Sobolev space @u 2 L2 (U ) for i = 1; 2, endowed with the inner of all functions u 2 L2 (U ) for which @x i product
@u @v
@u @v hu; viH 1 (U ) = hu; viL2 (U ) + @x ; @x L2 (U ) + @x ; 2 : 1 1 2 @x2 L (U )
Moreover we de ne the Hilbert space H (; U ) = fv 2 H 1 (U ) : v 2 L2(U )g; (9) with inner product hu; viH (;U ) = hu; viH 1 (U ) + hu; viL2 (U ) : (10) Points on the boundary of U are identi ed with the corresponding angle in polar coordinates, i.e. we set @U = f' + 2Z : ' 2 Rg = R=2Z =: T. For l 0, the Sobolev space H l (T) is de ned by X H l (T) = ff 2 L2 (T) : (1 + n2 )l jf^(n)j2 < 1g; (11) n2Z
4 where
W. RING Z 2 1 ^ f ( ) e,in d f (n) = 2
0
(12)
denotes the Fourier transform of f . H l (T) is a Hilbert space with respect to the inner product X hf; giH l (T) = (1 + n2 )l f^(n)^g(n): (13) n2 Z
By H ,l (T) (l 0) we denote the dual of H l (T). With the Fourier transform de ned on H ,l (T) by f^(n) := hf; e,in' iH ,l (T);H l(T) we nd X hf; giH ,l (T);H l(T) = f^(n) g^(,n) (14) n2Z
for the duality pairing h; iH ,l (T);H l(T) and for f 2 H ,l (T) and g 2 H l (T). Moreover, H ,l (T) is characterized by (11) and the inner product on H ,l (T) is given by (14), in both cases with l replaced by ,l. We have X f = f^(n) ein' (15) n2Z
for f 2 H l (T), l 2 R, where the series (15) converges in H l (T). The embeddings H l (T) ,! H m(T) and C (T) ,! H s (T) are compact and the images are dense for l > m and s > 21 . We de ne the di erential operator D : H l (T) ! H l,1 (T) by X Df = in f^(n) ein' : (16) n2Z
It is a bounded linear operator with ker D given by the constant functions on T. The space C 1 (T) of all in nitely di erentiable functions on T is characterized by C 1 (T) = ff : T ! C : jnjk jf^(n)j ! 0 as n ! 1 for all k 2 N g: (17) The preceding facts follow as a simple special case from theorems on Sobolev spaces on smooth compact manifolds as presented e.g. in Wloka [20]. 12 1 The trace map 0 : H (U ) ! H (T) is the unique continuous extension of the operator 0 : C 1 (U ) ! C 1 (T) de ned by 0 u =1 uj, from C 1 (U ) to H 1 (U ). 0 is a bounded linear surjection from H 1 (U ) on H 2 (T). Analogously we de ne 1 : 12 , 1 1 H (; U ) ! H (T) to be the unique extension of 1 : C (U ) ! C (T), where @u
1 u = @ n T, and n denotes the unit outward normal vector to the boundary T. If there is no danger of misinterpretation, we write 0 u = u and 1 u = @u @ n . (c.f. [20, p 130, Satz 8.7] and Dautray and Lions [5, p 381, Lemma 1] for these facts). The Poisson kernel P : U ! C is de ned by
P (r; ') =
X
rjnj ein' :
n2Z
(18)
The solution u 2 H 1 (U ) of the Dirichlet problem (
u = 0 on U
0 u = f 2 H 12 (T)
(19)
5
IDENTIFICATION OF A CORE
is given by the Poisson integral
u(r; ') = 21
2
Z
0
P (r; ' , ) f ( ) d = P
X
rjnj f^(n) ein' :
(20)
n2Z rjnj f^(n) ein' ! u(r; ') in H 1 (U )
For a proof of this, observe that uN (r; ') := jnjN P and 0 uN = jnjN f^(n) ein' ! f in H 21 (T) as N ! 1. Hence 0 u = f and obviously u = 0 on U . This proves (20). Since uN = u = 0 on U , it is also clearPthat uN ! u in H (; U ) as N1 ! 1. By the trace theorem, we nd
1 uN = jnjN jnj f^(n) ein' ! 1 u in H , 2 (T). We can therefore state the following Proposition 2.1. Suppose given u 2 H 1 (U ), f 8 > :
1 u = g Then
g=
X
jnj f^(n) ein' ;
n2Z
and hence g^(n) = jnjf^(n) for all n 2 Z. Remark 2.2. For f 2 H l (T), l 2 R we de ne
f~ =
X
,i sign(n) f^(n) ein' :
n2Z
f~ is called the function conjugate to f , or the Hilbert transform of f . The Hilbert transform is a bounded linear operator on H l (T) with kernel given by the constants. (c.f. Edwards [7] for further details on conjugate functions.) A short computation shows that X jnj f^(n) ein' = (Df ) ; (21) n2Z
where D is the (tangential) derivative operator de ned in (16). With this notation we can write the conclusion of Proposition 2.1 in the form
1 u = (D 0 u) :
(22)
3. An Alternative Approach. Let us return to problem (1) with = U and given data z1 2 H 21 (T) and z2 2 H , 12 (T). We assume that i > 0 for i = 1; 2, 1 6= 2 and a 2 H 1 (T) with 0 a(') 1 for all ' 2 T. Our goal is to formulate the identi cation problem for a as an inverse problem F (a) = z; where z is a function on T depending on the data z1 and z2 . However, this cannot be done in a naive way since we have to take care of the fact that there is an in nite
6
W. RING
variety of possible boundary values (z1 ; z2) for a,solution u(a) of u = fa , for one and the same a. Hence a mapping a 7! u(a) 7! 0 u(a); 1 u(a) would be, without additional assumptions, a multivalued function. We therefore seek to construct an invariant z depending on z1 and z2 which is the same function for all pairs (z1 ; z2 ) within the set of all possible boundary values of u(a), i.e. z should not depend on the special choice of z1 and z2 but only on the interface a, and we must be able to construct z from any given pair (z1 ; z2 ). This can be done by means of the following Proposition 3.1. Suppose a 2 H 1 (T) with 0 a(') 1 for all ' 2 T and let u; v 2 H 1 (U ) ful ll u = v = fa on U . Then
1 u , (D 0 u) , 22 = 1 v , (D 0 v) , 22
(23)
in H , 21 (T). Proof. Since w = u , v 2 H 1 (U ) and w = 0 on U , we can apply Proposition 2.1 to obtain D 0 (u , v) = 1 (u , v);
,
and consequently
1 u , (D 0 u) , 22 = 1 v , (D 0 v) , 22
in H , 21 (T). Proposition 3.1 implies that we can take any solution u(a) of u = fa and de ne the mapping
, (24) a 7! F (a) := 1 u(a) , D 0 u(a) , 22 since this function is independent of the special choice of u(a). We can therefore
express (1) in the form:
Solve F (a) = z2 , (Dz1 ) , 2 for a: (25) 2 The 22 -term on the right-hand side of (25) could of course be brought to the left and be included in the de nition of F . We prefer, however, to split it from F in order to obtain a more compact representation of F , as will become apparent later. In order to derive an integral representation for the operator F , we use the (unique) solution ua 2 H 1 (U ) of the Dirichlet problem (
u = fa on U :
0 u = 0
(26)
From the de nition (24) we obtain
F (a) = 1 ua , 22 :
(27)
Note that the boundary values in (26) were chosen such that the (D 0 u) -term vanishes. Our aim is now to calculate the Neumann trace of ua . We de ne PN (; ) =
7
IDENTIFICATION OF A CORE P
jnj in jnjN e for N
2 N and we use the notation PN (; ') = PN (; ' , ). Operators acting on functions with respect to the variable are indicated by the corresponding subscript. Green's formula Z
U
ua () PN (; ') , ua() PN (; ') d = Z
@ u () P (; ') , u () @ P (; ') d a N a @n N @U @ n together with ua = fa , PN (; ') = 0 and 0 ua = 0 implies (1 , 2 ) and hence
(28)
Z
0
2Z a( ) 0
PN (; , ) d d + 2 =
2
Z
0
1 ua ( )
X
jnjN
ein(', ) d ;
1 , 2 Z 2Z a( ) P (; , ) d d N 2 0 0 Z 2 X 1 =
1 ua ( ) e,in d ein' , 22 : 2 jnjN
0
The sum on the right-hand side of (28) is a truncated Fourier series which converges to 1 ua in1 H , 21 (T). (Recall that ua 2 H (; U ) and by the trace theorem 1 ua 2 H , 2 (T).) PN (; ') converges with respect to pointwise on U towards the Poisson kernel P (; '), and it is bounded by the integrable function Q'(; ) = 4=j1 , ei(', ) j. It follows that limit and integral can be transposed in (28) for N ! 1. We can therefore combine (28), (27) and (25) to obtain and
F (a) = 1 2, 2
Z
0
2Z a( ) 0
P (; , ) d d
(29)
1 , 2 Z 2Z a( ) P (; ' , ) d d = z , (Dz ) , 2 , (30) 2 1 2 0 0 2 as an equation for a which is equivalent to (1). Let us return to (28). Carrying out the integration with respect to and taking the limit N ! 1, we nd X 1 , 2 Z 2 a( )jnj+2 F (a) = e,in d ein' : (31) 2 j n j + 2 0 n2Z
From this representation we immediately obtain Z 2 , a( )jnj+2 e,in d ; , 1 2 b F (a) (n) = 2 jnj + 2 0 and therefore, following (30),
1 , 2 Z 2 a( )jnj+2 e,in d = z^ (n) , jnj z^ (n) , 2 2 1 2 0 jnj + 2 2 0;n
(32)
(33)
8
W. RING
for all n 2 Z, where 0;n = 1 if n = 0 and 0;n = 0 else. This is a system of in nitely many (nonlinear) functionals which we have to solve simultaneously for a. Remark 3.2. In view of generalization to domains of arbitrary geometry we shall now give a brief explanation of the above formula (30) from a di erent perspective. Let G(x; ) denote the Green's function for the Dirichlet problem (6) on U , i.e. assume that the solution of (6) can be written in the form
u(x) =
Z
U
G(x; ) fa () d +
Z
@U
u() @@n G(x; ) d
(34)
for all x 2 U . We (formally) take now rst order traces on both sides of (34), where we di erentiate under the integrals on the right-hand side and evaluate at x 2 @U . A rigorous justi cation of these operations will be given at some other place. We obtain Z
@ G(x; ) f () d = @ u(x) , Z u() @ 2 G(x; ) d a @n @ nx @ n U @ nx @U for all x 2 @U . It is well known (c.f. Myint-U [13, p 289]) that G is given by
(35)
G(r; '; ; ) = 41 log[r2 + 2 , 2r cos(' , )] , 41 log[1 + r2 2 , 2r cos(' , )]:
A short calculation then gives
@ G(x; ) @ G(r; '; 1; ) = 1 1 , r2 1 P (r; ' , ): = = 2 @ n 2 1 + r , 2r cos(' , ) 2 2@U @ (36)
For the last equality we refer to Rudin [17, p 233, (2)]. From the symmetry of G it follows that 1 P (; ' , ): @ G(x; ) = @ nx x2@U 2 With this (35) has the form 1 Z 2Z 1 P (; ' , ) f (; ) d d a 2 0 0 Z 2 @ P (r; ' , ) d = 1 u(') , 21
0 u( ) @r r=1 0 Z 2 X = 1 u(') + 21 jnj
0 u( ) e,in d ein' 0 n2N (37) = 1 u(') , (D 0 u) sptilde('): Since the left-hand side of (37) can be written as
1 , 2 Z 2Z a( ) P (; ' , ) d d + 2 ; 2 0 0 2
we have derived (30) from the relation (35). Obviously (35) also makes sense if U is replaced by some arbitrary R2 , hence we can consider (35) as a generalization of (30) to domains di erent from the unit disk, provided that the Green's function for the Dirichlet problem on is known.
9
IDENTIFICATION OF A CORE
4. Smoothness, Di erentiability and Uniqueness. In this section we consider certain properties of the operator F as de ned in (29). We investigate smoothness properties of F (a) (') and calculate the Frechet derivative F 0 (a). The domain for F is de ned as the set dom F = fa 2 H 1 (T) : 0 a(') 1 for all ' 2 Tg:
(38)
Frequently we shall consider the interior (dom F ) of dom F in H 1 (T), given by (dom F ) = fa 2 H 1 (T) : 0 < a(') < 1 for all ' 2 Tg:
(39)
Proposition 4.1. Let F be as de ned in (29) and suppose that dom F and (dom F ) are given as in (38) and (39) respectively. Then
F (a) 2 H l (T) for every a 2 dom F and l < 12 . If a 2 (dom F ) , then F (a) 2 C 1 (T): Moreover, for every a 2 (dom F ) , there exists a neighbourhood U (a) of a in H 1 (T) such that
F : U (a) H 1 (T) ! H l (T) is Lipschitz continuous for every l 2 R, with Lipschitz constant given by k = j1 ,2 j c~, where c~ is some constant depending on U (a) and l. Proof. From (14) and (29) follows
(40) (41)
2
l
Z 2 2 X a( )jnj+2 e,in d , 1 2 2 l (1 + n ) H (T) = 2 jnj + 2 0 n2Z 2l X j1 , 2 j2 ((1jn+j +n2))2 : n2Z
kF (a)k2 l
2
n) l 2l,2 for all l 0 and jnj 1, the series (40) converges if Since ((1+ jnj+2)2 2 jnj 2l , 2 < ,1, i.e. if l < 12 . Let us now suppose that a 2 (dom F ) . Since a 2 H 1 (T) C (T), it follows that sup'2T ja(')j = kak1 < 1. We have
Z
2
jnj+2
a( ) e,in d jnjk F (a) b(n) = jnjk 1 2, 2 jnj + 2 0 k (42) jnjjn+j 2 j1 , 2 j kakj1nj+2 ! 0 as jnj ! 1 since kak1 < 1, and hence F (a) 2 C 1 (T), due to (17). The boundedness of the embedding H 1 (T) ,! C (T) implies that we can nd a neighbourhood U (a) of a 2 (dom F ) in H 1 (T) such that m := supf kbk1 : b 2 U (a) g < 1. We consider F : U (a) H 1 (T) ! C 1 (T) H l (T) ,
10
W. RING
with l 2 R. Then kF (a1 ) , F (a2 )kH l (T) =
X
n2Z
2
Z
2
0
a1 ( )jnj+2 , a2 ( )jnj+2 e,in d jnj + 2
2 1 2
j+1 (1 + n2 )l Z 2 ja ( ) , a ( )j jnX k j n j +1,k d 2 a ( ) a ( ) 1 2 1 2 2 0 n2Z (jnj + 2) k=0 X 1 j1 , 2 j (1 + n2 )l m2(jnj+1) 2 ka1 , a2 k1 k ka1 , a2 kH 1 (T);
j1 2, 2 j (43)
(1 + n2 )l 1 , 2
X
12
n2Z
where
k := j1 , 2 j
X
(1 + n2 )l m2(jnj+1)
n2Z
21
c1 < 1
with c1 such that kf k1 c1 kf kH 1 (T) for all f 2 H 1 (T). The convergence of the series in (43) is easily seen by the quotient criterion. Hence we have proved that F : U (a) H 1 (T) ! H l (T) is locally Lipschitz continuous for every l 2 R. Remark 4.2. If (25) has a solution a 2 (dom F ) , then Proposition 4.1 and (25) imply that although each of the boundary data z2 and (Dz2 ) is only in H , 12 (T) (and, in general, will not be in any space H l (T) with l > , 21 ) their di erence is in C 1 (T). We also point out that in Proposition 4.1 we have also proved the Lipschitz continuity of F with respect to the C (T)-norm on (dom F ) . Proposition 4.3. F : (dom F ) H 1 (T) ! H l (T), with l 2 R, is Frechet di erentiable at every point a 2 (dom F ) with derivative given by [F 0 (a) h](') = 1 2, 2
2
Z
0
a( ) P (a( ); ' , ) h( ) d
(44)
for all h 2 H 1 (T). The Fourier transform of F 0 (a) h is given by
[F 0 (a) h]b(n) = 1 , 2
2
Z
2
0
a( )jnj+1 h( ) e,in d
(45)
for h 2 H 1 (T) and n 2 N . Proof. First of all we show that the Fourier transform of (44) is in fact given by (45). Since for a 2 (dom F ) , 2 P a( ); ' , (46) 1 , kak1 < 1; we can exchange summation and integral in (44) and obtain
F 0 (a)h = 1 2, 2
XZ
n2Z 0
2
a( )jnj+1 h( ) e,in d ein'
11
IDENTIFICATION OF A CORE
which proves that (45) is correct. To prove Frechet di erentiability it is sucient (c.f. Berger [1, p 68, Th 2.1.13]) to ensure that 1. F 0 (a) : H 1 (T) ! H l (T) is a bounded linear Operator. 2. F 0 depends locally Lipschitz continuously on a, i.e. there exists an L > 0 (depending on a) such that
0
F (a) , F 0 (b) (47) L(H 1 (T);H l(T)) Lka , bkH 1 (T) for all b in some neighbourhood of a. (This is even more than required in [1].) 3. F 0 (a) is the Gateaux derivative of F at a, i.e. lim kF (a + #h) , F (a) , # F 0 (a) hkH l (T) = 0; #!0
where # 2 [0; 1) for every h 2 H l (T). We start with 1. and show that F 0 as de ned in (44) is a bounded linear operator from H 1 (T) into H l (T) for some arbitrary but xed l 2 R. Using (45) we nd
kF 0 (a) hkH l (T) =
X
(1 + n2 )l 1 , 2
2
n2Z
X
(1 , 2 ) c
,X
2
Z
0
2 12 j n j +1 , in a( ) h( ) e d
jnj+1) jh^ (n)j2 (1 + n2 )l kak2( 1
12
n2Z
1
jh^ (n)j2 2 = c khkL2(T)
n2Z
where c is chosen large enough such that (1 , 2 ) (1 + n2 ) 2l kakj1nj+1 c for all n 2 Z. (Recall that a 2 (dom F ) and therefore kak1 < 1.) Hence F 0 (a) : H 1 (T) L2 (T) ! H l (T) is a bounded linear operator. We note that this is a stronger result than required in 1. and that F 0 (a) can be uniquely extended to a bounded linear operator from L2 (T) into H l (T). For the proof of 2. we choose a neighbourhood V (a) of a in H 1 (T) such that kbk1 m < 1 for all b 2 V (a). Then we have kF 0 (a) h , F 0 (b) hkH l (T) =
X
(1 + n2 )l 1 , 2
n2Z
(1 , 2 )
2
X
2 ,
Z
0
a( ) , b( )
jnj
X
k=0
a( )k b( )jnj,k h( ) e,in
d
2 12
1
(1 + n2 )l (jnj + 1)2 m2jnj jh^ (n)j2 2 ka , bk1:
n2Z
If we choose constants c1 , c2 such that ka , bk1 c1 ka , bkH 1(T) for all b 2 V (a) and for all n 2 Z, we obtain
(1 + n2 )l,1 (jnj + 1)2 m2jnj c22 X
kF 0(a) h , F 0 (b) hkH l (T) (1 , 2 ) c1 c2 ka , bkH 1 (T) = L ka , bkH 1 (T) khkH 1 (T);
(1 + n2 ) jh^ (n)j2
n2Z
12
12
W. RING
with
L = (1 , 2 ) c1 c2 and consequently
(48)
F 0 (a) , F 0 (b) L(H 1 (T);H l(T)) L ka , bkH 1 (T)
for all b 2 V (a). It remains to prove that F 0 is the Gateaux derivative of F . We nd
kF (a + #h) , F (a) , # F 0 (a) hkH l (T) X 2 = (1 + n2 )l 1 , 2 n2Z Z
2
a( ) + #h( ) jnj+2 , a( )jnj+2 jnj+1 h( ) e,in d , # a ( ) jnj + 2
2 ,
0
= 1 2, 2 X
n2Z
# (1 , 2 )
2
Z
(1 + n2 )l X
0
j+2 , i 1 h jnX jnj+2#k h( )k a( )jnj+2,k e,in d jnj + 2 k=2 k
(1 + n2 )l
n2Z
2 12
1
a + p#jhj
jnj+2 2 1 jnj + 2
12
2 12
;
k where we p used #k,1 # 2 for all k 2 and # 2 (0; 1). The last series above converges if ka + # jhj k1 < 1 by the quotient criterion, hence the whole expression tends to 0 if # ! 0. We now address the problem of uniqueness of a solution a 2 dom F of F (a) = z for z 2 rg F . In Section 1. we proved that every solution a 2 C (T) with 0 a(') 1 is also a solution of the equivalent inverse potential problem (3). Hence, to prove injectivity of the operator F , it is sucient to prove unique solvability of (3). The question of uniqueness for (3) is a well investigated topic in the literature. It was rst shown by Novikov [14] that ~ in (2) is uniquely determined by the potential uj@ in the class of star-shaped open sets with smooth boundary. Note that a = f(r; ') : r < a(')g is star shaped with respect to the origin. A proof of a unique-solvability result for (2) in Rd (d 2) for star-shaped regions can be found in [9, p 55, Th 3.1.1]. The following theorem is an immediate consequence of this result. Theorem 4.4 (Uniqueness). Suppose that a; b 2 dom F , satisfying a(') > 0 or b(') > 0 for all ' 2 T and suppose that
F (a) = F (b) holds in H l (T) with l < 21 . Then
a=b in H 1 (T). Next we shall prove the injectivity of F 0 .
13
IDENTIFICATION OF A CORE
Theorem 4.5. The operator F 0 (a) : H 1 (T) ! H l (T) with l 2 R is injective for all a 2 (dom F ) . Proof. Suppose that f 2 ker F 0 (a) H 1 (T), i.e. (c.f. (45) ) Z
0
2
a(')jnj+1 f (') ein' d' = 0 for all n 2 Z:
(49)
We shall rst prove that (49) implies 2
Z
0
,
a(') f (') h a('); ' d' = 0
(50)
for all functions h which are harmonic on some neighbourhood V of a = (r; ') : r a(') . For this purpose x h harmonic on V and let hc : V ! R be a harmonic conjugate of h, i.e. suppose hc is chosen such that the function H (x + iy) = h(x; y) + ihc(x; y) is analytic on V . Here we have changed for the moment from the usual polar coordinates (r; ') to Cartesian coordinates (x; y). For a proof of the existence of a harmonic conjugate to a given harmonic function see Kellogg [10, p 345, Th III]. H can be approximated by polynomials
Qn (z ) =
mX (n) k=0
ck;n z k ;
with z = x + iy, ck;n 2 C and m(n) 2 N , uniformly on a by Runge's theorem (c.f. [17, p 271, Th 13.7]). Then
We set and
(n) m X , Re Qn (z ) , H (z ) = 12 ck;n z k + ck;n z k , h(z ): k=0
(51)
d0;n = 21 (c0;n + c0;n ) dk;n
(
1 2 ck;n 1 ck;n 2
for k 1 : for k ,1
With z = r ein' we can write (51) as Re (Qn , H ) =
mX (n) k=,m(n)
dk;n rjkj eik' , h(r; '):
(52)
Since jRe(z )j jz j for every complex number z , (52) implies that h can be approximated by a sequence of `harmonic polynomials'
qn (r; ') =
mX (n) k=,m(n)
dk;n rjkj eik'
14
W. RING
uniformly on a . Hence we obtain Z
2
0
,
a(') f (') h a('); ' d'
Z
2
0
(53)
,
,
a(') f (') h a('); ' , qn a('); ' d'
+
mX (n)
2
Z
k=,m(n)
dk;n
0
a(')jnj+1 f (') ein' d' :
The second term on the right-hand side of (53) is 0 due to (49) and the rst one tends to 0 since pn ! h uniformly on @ a . Hence the left-hand side of (53) must be 0 and (50) is proved. We consider now the Dirichlet problem (
w = 0 on
a , w @ a = w a('); ' = f (') :
(54)
since a 2 (dom F ) , there exists > 0 satisfying a(') for all ' 2 T. This together with f 2 H 1 (T) C (T) implies that f : @ a ! C in (54) is a continuous function on @ a . From this and from geometric properties of the boundary @ a (c.f. [6, p 336, Th 1 and p 340, Prop 4]) the existence of a classical solution w for (54) follows, that is to say w = 0 on a and lim w(x) = f (z ) x!z ,
for x 2 a and z 2 @ a. We de ne wk (r; ') = w ,(1 , k1 )r; ' for k 2 N . From the properties of a classical solution it follows ,that wk a('); ' ! f (') pointwise on T and by the maximum modulus principle w a('); ' is bounded by jf (')j. Hence, by ,k the dominated convergence theorem, wk a('); ' ! f (') in L1 (T). Moreover, wk is harmonic on Vk = (1 + k1 ) a . Employing (50) and the last arguments we obtain 0=
2
Z
0
,
a(') f (') wk a('); ' d' !
2
Z
0
a(') jf (')j2 d' 2 kf k2L2(T);
where 0 < a(') for all ' 2 T. Hence we nd f = 0 in H 1 (T). R Remark 4.6. The fact that 02 a(') f (') a(')jnj ein' d' = 0 for all n 2 Z implies that f = 0 in H 1 (T) is equivalent to the assertion that the functions fa(')jnj ein' gn2Z are dense in H 1 (T) (and hence in L2 (T)). It would be interesting to prove this consequence of Theorem 4.5 directly and not by solving boundary value problems which seems somewhat unnatural in this context. 5. Ill-Posedness, Least-Squares Formulation, and Regularization. In this section we consider the problem F (a) = z; (55) within the concept of least-squares methods. Solutions of (55) do not depend continuously on the data z , not even if we replace the parameter space H 1 (T) by L2(T) as the following example shows. We consider the sequence of functions fan gn2N dom F de ned by an (') = 12 + 2 cos n' = 12 + (ein' + e,in')
15
IDENTIFICATION OF A CORE
with 0 < < 14 and the constant function
r
a(') = 14 + 2 2: It is easy to see that r
2 kan , ak2L2(T) = 21 , 14 + 2 2 + 2 2 = const > 0; and hence an 6! a for n ! 1 in L2 (T). On the other hand we nd Z 2 h , , 1 1 + (ein + e,in )2 , 1 , 2 2i d = 0 1 2 b F (an ) , F (a) (0) = 2 2 2 4
0
for all n 1, and (recall that F (a) is constant) ,
Z
F (an ) , F (a) b(k) = 1 2, 2
2
0
1 1 + (ein + e,in )jkj+2 e,ik d jk j + 2 2
for n 1. The term 12 + (ein + e,in ) jkj+2 can be written in the form P all k 6= 0 and imn with coecients cm;n 2 C , cm;n 6= 0 only for nitely many m 2 Z. m2Zcm;n e , ik 2 Hence the orthogonality of the functions e in L (T) implies that F (an ) b(k) = 0 for all k 2 Z with 1 jkj < n. For jkj n we have
,
,
F (a
Z 2 , 1 1 + 2 jkj+2 d = 1 , 2 qjkj+2 1 2 k 2 jk j + 2 0 jk j + 2 2
n ) b(
)
where q = 21 + 2 < 1. We therefore obtain
kF (an ) , F (a)k2H l (T) =
X
,
,
k2Z
(1 , 2 )
X
jkjn
jkjn
jkj+2
(1 + k2 )l jqkj + 2 :
Using the quotient criterion it is easy to see that the series converges. Thus we have X
(1 + k2 )l F (an ) b(k) , F (a) b(k)
2 l qjkj+2 k2Z(1 + k ) jkj+2
P
jkj+2
(1 + k2 )l jqkj + 2 ! 0
for n ! 1 and consequently F (an ) ! F (a) in H l (T). Since the kth derivative Dk is a continuous operator from H k+1 (T) into H 1 (T) it follows that Dk F (an ) ! Dk F (a) in H 1 (T) and due to the compact embedding H 1 (T) ,! C (T) we have
Dk F (an ) ! Dk F (a) uniformly on T for all k 2 N . From this we can conclude that F (an ) ! F (a) even in the topology on C 1(T).
16
W. RING
The lack of continuous dependence of the interface a on the observation z in (55) is the reason why we consider the inverse problem of determining a from z as `illposed' (in the sense of Hadamard). An ecient and well investigated tool to deal with ill-posed problems is Tikhonov's regularization method. (c.f. Engl, Kunisch and Neubauer [8], Seidman and Vogel [18], the work [12].) The method consists in approximating solutions of (55) by minimizers of the Tikhonov functional n
o
min kF (a) , z k2H , 12 (T) + kD(a , a )k2L2 (T) ; a2dom F
(56)
where D : H 1 (T) ! L2 (T) and is some (small) regularization parameter. The function a in the regularization term plays the r^ole of an `origin' with respect to which the semidistance kD( , a )kL2 (T) is measured. It should be chosen from some a-priori information close to the actual solution of (55). In the following we assume that a 2 dom F H 1 (T) and z 2 H , 12 (T). Although we could use any H l (T)-space with l 2 (,1; 12 ) as output-space (c.f. Proposition 4.1), we consider only the case where l = , 21 . This is motivated by the fact that the data on the right-hand side of (55) have the form z = z2 , (Dz1 ) , 22 , where the functions z2 and (Dz1 ) are in H , 12 (T), only their di erence is smoother, e.g. in C 1 (T) if a 2 (dom F ) . If (Dz1 ) and z2 result from real data, we cannot expect the di erence z2 , (Dz1 ) to be reasonable accurate in the C 1 (T)-topology or in any H l (T)-norm with large l but only in the topology which is `natural' to both terms 12 , of the di erence, which is H (T1). Moreover, we assume that z is given by 1some inaccurate measurement z 2 H , 2 (T) of `unperturbed' data z0 2 rg F H , 2 (T), with kz0 , z kH , 21 (T) ; (57) where we suppose that a0 2 (dom F ) exists, satisfying F (a0 ) = z0 . By Theorem 4.4, a0 is then uniquely determined by z0 . In [12] conditions were obtained which ensure solvability, stability with respect to noisy data, and convergence of solutions of (56) | with a certain convergence rate | towards a solution of (55) if the noise level and the parameter go to 0. Applying the results [12] leads to the following 1 Theorem 5.1. Assume that F : dom F H 1 (T) ! H , 2 (T) is given as in (29) and (38) respectively, and a 2 H 1 (T), z 2 H , 12 (T) and > 0. Then there exists a solution a = a(1 ; z ) 2 dom F of (56). Solutions of (56) are stable in the sense that for zn ! z in H , 2 (T), the sequence an = a( ; zn ) has a H 1 (T)-convergent subsequence, and the limit of every weakly convergent subsequence of fang is solution of n
o
min kF (a) , z k2H , 12 (T) + kD(a , a )k2L2 (T) : a2dom F If z0 = F (a0 ) and kz0 , z kH , 12 (T) and if we choose = () then a( (); z ) ! a0 in H 1 (T). Moreover, if there exists ! 2 H 12 (T) such that
,
D2 (a0 , a ) (
or equivalently
X
Z 2 , , 1 2 ) = 2 a0 ( ) P a0 ( ); ' , 0
n2 (a0 , a )b(n) ein = (1 , 2 )
n2Z
X
!(') d'
!^ (n) a0 ( )jnj+1 ein :
n2Z
(58) (59)
17
IDENTIFICATION OF A CORE
and if k!kH 12 (T) is bounded by a certain constant M depending on a0 , then
p ka( (); z ) , a0 kH 1 (T) = O( ) and kF a( (); z ) , z0 kH , 12 (T) = O():
Remark 5.2. A detailed description of the dependence of the constant M on the solution a0 is ,given inthe Appendix. Condition (58) is in fact the explicit form of D D(a0 , a ) = F 0 (a0 ) !. We demonstrate that (58) includes certain smoothness , requirements on the di erence a , a0 . From (59) we see that the range of F 0 (a0 ) consists of functions whichPcan be expressed as series in terms of a0 ( )jnj+1 ein with coecients cn satisfying n2Zjnj jcn j2 < 1. Using the constant function 1 ,1 2 we , nd that a0 is always in the range of F 0 (a0 ) . It is easily, seen that the series 0 on the right-hand side of (59) converges uniformly, therefore F (a0 ) g 2 C (T) for all g 2 H 12 (T). Thus, if (58) is ful lled, we have D2 (a0 , a ) 2 C (T) and hence a0 , a 2 C 2 (T). This, however, is only the minimal smoothness requirement as the following example shows. Suppose that a0 = q < 1 is constant. Then X F 0 (a0 ) g = (1 , 2 ) g^(n) qjnj+1 ein 2 C 1 (T)
,
n2Z
for all g 2 H 21 (T). Hence in this case we have rg F 0 (a0 ) C 1 (T). 6. Numerical Results. As a starting point for the numerical implementation we use the system of equations (33). We introduce the set SN of real-valued trigonometric polynomials of degree less or equal than d = N2,1 given by SN = f N =
d X
j =,d
,
cNj eij' : cN0 2 R; cNj 2 C ; cNj = cN,j for j 6= 0
with odd N 2 N . We use the abbreviations 'Nk := 2N k and fkN := f N ('Nk ) for f N 2 SN and k = 0; ; N , 1. Given f N 2 SN there is a one-to-one correspondence between the Fourier-coecients fcN0 ; ; cNd g and the interpolation points ff0N ; ; fNN,1 g via the discrete Fourier transform (c.f. [5, pp 60, 61]), hence we can use ff0N ; ; fNN,1g as coordinate-vector for f N 2 SN instead of fcN,d; ; cNd g. In the following it is assumed that F acts on SN , i.e. we consider F (aN ) and F (aN ) b(n) for aN 2 SN . Since the exact evaluation of the integral in (33) for aN 2 SN and large values of jnj is very complicated, we used a numerical quadrature by means of the rectangular rule for the computation of the integral F (aN ) b(n) ' 1 , 2
N
NX ,1 (aN )jnj+2 k e,in'Nk : j n j + 2 k=0
(60)
Note that the fact that aN (0) = aN (2) implies that the rectangular rule and the trapezoid rule give the same numerical quadrature formula (60). We also mention that the dimension N of SN also de nes the mesh-width and therefore also the accuracy of the numerical integration (60). It follows from the de nition that the operator F depends linearly on the di erence = 1 , 2 , hence (56) can be transformed into h
i
min 2 kFe (a) , ze k2H , 12 (T) + e kD(a , a )k2L2 (T) ; a2dom F
18
W. RING
where Fe = 1 F , ze = 1 z and e = 2 are independent of . Therefore, if 6= 0, the optimization problem (56) depends on only as a scaling factor. This justi es to x
1 = 2 and 2 = 1 for all numerical examples. In order to obtain numerical data for the right-hand side of (33) from some known interface function a, we also used the numerical integration (60). We set NX d,1 (aNd )jnj+2 ,in'Nk d =: z^Nd (n); k e d k=0 jnj + 2 for all jnj Nd2,1 , and z^Nd (n) = 0 for jnj > Nd2,1 . For this (forward) calculation large values for Nd (usually Nd = 3003) were used to obtain reasonable accurate `real' data, while for the calculation of faNk gNk=0 from the data z^Nd , mostly N = 91 was chosen. The accuracy of the data was then arti cially destroyed by adding a random noise zNd , where
z^(n) ' 1N, 2
zNd ('Nk d ) = Rand(0; 1) , 12 h
i
with Rand(0; 1) a random number in [0; 1] and the noise level is chosen as one of the values in f0:05; 0:01; 0:005; 0:001; 0:0005g. All interface functions under consideration were centered around r = 0:8. This gave data z with kz k1 0:32. Hence the noise level varied between 0.075% and 7.5% with respect to the kk1 -norm. Moreover, the random function zNd was the same for all interface functions. A di erent possibility to construct data z from a known interface a would be to solve the forward problem u = fa by means of the fundamental solution F (x; ) = , 21 log jx , jRand to evaluate both Dirichlet and Neumann boundary data of the solution u(x) = F (x; ) fa() d at the points xk = (cos 'Nk d ; sin 'Nk d ). We could then use the discrete Fourier transform to obtain the Fourier transform of the data z = z2 , (Dz1 ) , 22 . With these approximations for the operator F and the data z , we can enter (56) and calculate the interface a as the minimizer of the approximated Tikhonov functional
J N (aN ) = (61)
jnj M2,1
+
, 2 (1 + n2 ), 12 1 N
X
X
jnj M2,1
NX ,1 (aN )jnj+2 ,in'Nk k e k=0 jnj + 2
2
, z^(n)
jnj2 a^N (n) , a^ (n) 2 :
The rst thing to come to mind is perhaps to choose M = N in (61). With this choice it was observed in the numerical computations that large values of M increased the computation time considerably. We were therefore bound to x M somewhere below 100. However, it turned out that (for xed M ) the value of N (the number of interpolation points) could be increased beyond M without increasing the computational error and without loosing too much computation time. This may be due to
19
IDENTIFICATION OF A CORE
the fact that [F (a)]b(n) decreases exponentially as jnj ! 1 (c.f. (42) with k = 0), and hence the in uence of the terms for large values of jnj to the minimum of (61) (the stabilized problem!) is small. For this reason we used `more variables than equations', in fact we xed N = 91 and M = 45 for all computations. Moreover for 0 < n M2,1 we have F N (aN ) b(,n) = F N (aN ) b(n), z^N (,n) = z^N (n) and (aN , a )b(,n) = (aN , a )b(n); hence each of the terms
,1 (aN )jnj+2 1 , 2 NX ,in'Nk , z^(n) 2 k e N k=0 jnj + 2
and N ^ (
a n) , a^ (n) 2
occurs twice in (61). If we identify aN 2 S N with its coordinate vector (aN0 ; aNN ,1 ) 2 RN and if we ignore the constraints 0 aN 1 we obtain for (61)
min
aN 2RN
+2 (62)
M2,1 X
n=1
,1 jaN j2 1 , 2 NX k , z^(0) 2 N k=0 2
, 2 (1 + n2 ), 12 1 N
M2,1
+2
X
n=1
NX ,1 jaN jn+2 ,in'Nk k e k=0 n + 2
jnj2 a^N (n) , a^ (n) 2
2
, z^(n)
:
This optimization problem was implemented using the numerical package Matlab and the optimization routine leastsq provided there, which is based on the Levenberg-Marquard algorithm for nonlinear least-squares problems. The Fourier transform (aN , a )b(n) occurring in the regularization term was computed by means of Matlab's fast Fourier transform routine fft. Formulation (62) also allows to calculate an analytic gradient which was used in the numerical computations. As a rst example we consider the interface given by the function
ac (') = 0:8 + 0:1 cos ';
(63)
which we named `Cosine'. The a-priori guess a waschosen 0. In Figure 1 the true ,1 of (62) as small interface ac as a solid line, and the solution-vector aNk (z; ) Nk=0 crosses, are plotted vs the angle '. The noise level was chosen to be 0 and the regularization parameter = 10,11 was (experimentally) determined such that the error
kaN (z; ) , ackL2 (T) '
X (
jnjN
21
aN (z; ))b(n) , a^c (n) 2
is minimal. This fact is indicated in the plot by the word `optimal'. Also the absolute errors with respect to the L2 (T)- and the H 1 (T)-norm are given. Figure 2 shows the result for the same interface function ac with perturbed data z N + zN , where = 0:05 (i.e. 7.5% noise). As in the previous example the regularization parameter
20
W. RING true interface (-) and solution of the regularized problem (+) 0.95
0.9 +++++++ 0.85
r
0.8
0.75
0.7
0.65
+++++ ++ ++ ++ ++ + ++ ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ++ + ++ ++ + guess=0 ++ + +++ ++ ++++++++
number of points=91 number of equations=45 regularization parameter alpha=1e-11, optimal H(-0.5)-noise=0, L2-noise=0 L2-error=2.907e-05, H1-error=9.476e-05
0.6 0
1
2
3
4
5
6
phi
Fig. 1. Cosine without Noise true interface (-) and solution of the regularized problem (+) 0.95
0.9 ++++++ 0.85
r
0.8
0.75
0.7
0.65
++ +++ ++ ++ ++ ++ + ++ ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ++ + ++ ++ ++ + guess=0 + ++ ++++++++++++
number of points=91 number of equations=45 regularization parameter alpha=0.002, optimal H(-0.5)-noise=0.001021, L2-noise=0.0143 L2-error=0.001338, H1-error=0.003781
0.6 0
1
2
3
4
5
6
phi
Fig. 2. Cosine with Noise ( = 0:05)
is optimal. Note that the optimal regularization parameter increases for growing noiselevel. It is also worth mentioning that the di erence between the L2 (T)-error and the H 1 (T)-error is very small which can be attributed to the fact that the representation of ac as a Fourier series has only two modes (ac 2 S3 ). An error in the higher modes is therefore not likely to occur since in (62) they are penalized by the regularization term with comparably large . The four plots in Figure 3 show numerical experiments with the interface function (
0:8 for ' 2 [0; , 1 ] [ [ + 1 ; 2] ah(') = 0:8 , 0:27 exp , 2 1 2 for ' 2 ( , 1 ; + 1 ) : 1,4(',)2 2 2
(64)
21
IDENTIFICATION OF A CORE
We called this function `Hole' and we investigated the behaviour of solutions of (62) for di erent noise-levels = 0; = 0:0005; = 0:005 and = 0:05. According to Theorem 5.1 the regularization parameter was chosen of order ; i.e. we choose = 10,8; = 10,7 and = 10,6 for the values 6= 0, where the value = 10,6 for = 0:05 is optimal. For this interface the optimal parameter is considerably smaller than for Cosine (with the same noise level), which is due to the fact that higher frequencies are needed to represent ah , hence it is not good to suppress them by a strong regularization term. true interface (-) and solution of the regularized problem (+)
true interface (-) and solution of the regularized problem (+)
0.95
0.95
0.9
0.9
0.85
0.85
0.8 ++++++++++++++++++++++++++++++++++++++++
+++++++++++++++++++++++++++++++++++++++
r
+
0.75
+
0.65
0.6 0
2
+
0.7
0.65
3
4
5
6
0.6 0
1
2
3
true interface (-) and solution of the regularized problem (+)
0.9
0.9
0.85
0.85
0.8 +++++++++++++++++++++++++++++++++++++++
++++++++++++++++++++++++++++++++++++++ +
r
+
0.6 0
2
++ ++ +++++++++++++ ++++++++++++++++++ ++ + + +
+ +
0.75
3
0.7
0.65
4
5
6
0.6 0
+
+ +
+
+ guess=0 + + + ++ number of points=91 number of equations=45 regularization parameter alpha=9e-07 H(-0.5)-noise=0.0001021, L2-noise=0.00143 L2-error=0.002288, H1-error=0.03717
1
++++ + +
0.8 ++++++++++++++++++++++++++++++++++
+
+
+
6
true interface (-) and solution of the regularized problem (+) 0.95
+
5
Noise ( = 0:0005)
0.95
0.65
4
phi
no Noise
0.7
+
guess=0 + + + + ++ number of points=91 number of equations=45 regularization parameter alpha=9e-08, optimal H(-0.5)-noise=1.021e-05, L2-noise=0.000143 L2-error=0.001425, H1-error=0.02805
phi
0.75
+
+
0.75
+
guess=0 + + + + ++ number of points=91 number of equations=45 regularization parameter alpha=1e-11 H(-0.5)-noise=0, L2-noise=0 L2-error=0.0004521, H1-error=0.01477 1
++++++++++++++++++++++++++++++++++++++ + +
+
+
+
0.7
0.8 ++++++++++++++++++++++++++++++++++++++++
+
+
+ + guess=0 + + ++ number of points=91 number of equations=45 regularization parameter alpha=9e-06 H(-0.5)-noise=0.001021, L2-noise=0.0143 L2-error=0.005702, H1-error=0.05715
1
2
phi
3
4
5
6
phi
Noise ( = 0:005)
Noise ( = 0:05)
Fig. 3. Hole with Increasing Noise and
In Table 1 and Table 2 we listed the H , 21 (T)-noise , the regularization parameter 2, the H 21 (T)- and the L2 (T)-error and ~ of the computed interface, and the rates and ~ for the preceding examples Cosine and Hole. The noise-levels in the H , 12 (T)-norm correspond to values for 2 f0p :05; 0:01; 0:005; 0:001; 0:0005g. Note that there seems to be a convergence rate for Cosine, while for Hole we can not observe this convergence rate neither in H 1 (T) nor in, L2 (T).Most likely this can be attributed to a lack of containment of D Dah in rg F 0 (ah ) (c.f. Theorem 5.1
22
W. RING
and Appendix, condition (R4) ).
H 21 -noise Reg. parameter H 1 -error 2 = 0.001021 0.0002043 0.0001021 0.00002043 0.00001021
0.002 0.0004 0.0002 0.00004 0.00002
0.003781 0.001556 0.0009800 0.0003477 0.0002492
0.0140 0.0119 0.0094 0.0059 0.0061
L2 -error ~
0.001338 0.0003865 0.0002247 0.00006623 0.00004551
~2 =
0.0018 0.00073 0.00049 0.00021 0.00020
Table 1
Error Rates for Cosine
H , 12 -noise Reg. parameter H 1 -error 2 = 0.001021 0.0002043 0.0001021 0.00002043 0.00001021
0.000009 0.0000018 0.0000009 0.0000018 0.00000009
0.05715 0.04042 0.03717 0.03089 0.02804
L2-error ~ ~2 =
3.199 7.997 13.532 46.705 77.007
0.005702 0.002681 0.002288 0.001685 0.001425
0.0318 0.0352 0.0513 0.1390 0.1989
Table 2
Error Rates for Hole
true interface (-) and solution of the regularized problem (+) 0.95
0.9
+ +++ + ++ + +++ + + + +++ ++ + +++ + 0.85 +++ +++ + ++ + + ++ +++ + +++ + 0.8 + +
r
+++ +++ ++ +++ + + ++ ++ +++ ++ + +++ + +++ ++ + +++ +++ + ++ ++
0.75
0.7
0.65
0.6 0
guess=0 number of points=91 number of equations=45 regularization parameter alpha=1e-11, optimal H(-0.5)-noise=0, L2-noise=0 L2-error=0.0007609, H1-error=0.0248 1
2
3
4
5
6
phi
Fig. 4. Zigzag without Noise
For the plots in Figures 4{7 we considered the piecewise linear interface function
23
IDENTIFICATION OF A CORE true interface (-) and solution of the regularized problem (+) 0.95
0.9
r
++++++ ++ ++ + ++++ + ++ + +++ + + 0.85 ++ ++ + ++ + ++ ++++++ + + ++ 0.8 ++ +
++ ++ ++++ +++ ++ + + +++ ++++ + ++++ +++ + + + + + + + + +++++
+
0.75
0.7
0.65
guess=0 number of points=91 number of equations=45 regularization parameter alpha=7e-05, optimal H(-0.5)-noise=0.001021, L2-noise=0.0143 L2-error=0.004186, H1-error=0.04132
0.6 0
1
2
3
4
5
6
phi
Fig. 5. Zigzag with Noise ( = 0:05)
`Zigzag', de ned by
8 4 >
