error analysis of an equation error method for the ... - CiteSeerX
ERROR ANALYSIS OF AN EQUATION ERROR METHOD FOR THE IDENTIFICATION OF THE DIFFUSION COEFFICIENT IN A QUASILINEAR PARABOLIC DIFFERENTIAL EQUATION MARTIN HANKEy AND OTMAR SCHERZERz
Abstract. We consider the inverse problem of reconstructing the di usion coecient in a quasilinear parabolic di erential equation in divergence form from measurements of the solution at a nite number of points in the interior of the domain. An equation error method is developed which transforms the inverse problem into a system of linear operator equations for the di usion coecient, and which can be solved by the conjugate gradient method in a very ecient and stable manner. A detailed error analysis relates the required number of measurements with their accuracy. Numerical results illustrate the performance of the method. Key words. Inverse problems, quasilinearparabolicequation, parameter identi cation,ill-posed problems, regularization methods, conjugate gradient method AMS subject classi cations. 35R30, 65M30, 65M15
1. Introduction. We consider the quasilinear parabolic initial-boundary value
problem (1.1)
?
ut = a(u)ux x ; u(0; x) = !0 ; u(t; 0) = f(t) ; u(t; 1) = g(t) ;
with t 2 I = [0; T] and x 2 = [0; 1]. Here, a 2 H 2 is a positive function, !0 a constant and f and g are given boundary data. Given a, the direct problem consists in solving the di erential equation for u. Here we are concerned with the following inverse problem: Given u for certain abscissa 0 = x0 < : : : < xm = 1 and all t 2 [0; T], determine the di usion coecient a(u). We refer to Dummel and Handrock-Meyer [6] concerning the identi ability of an analytic di usion coecient from the given data. Inverse problems of similar sort arise in a variety of di erent contexts. This particular one occured in an industrial application where an ingot casting manufacturer was asking for the thermal conductivity of the form sand he is using in the manufacturing process. Due to di erent heat transport phenomena, the thermal conductivity of the sand is likely to change in the presence of very high temperatures. In this application the available data for the inverse problem is the temperature evolution measured by a few thermocouples placed at increasing distances from the heat source. The standard approach for the solution of this inverse problem is the output least squares method; cf., e.g., Nabokov [13] for an application of output least squares to this particular problem. This method requires a subroutine to solve the forward problem, e.g., a Crank-Nicolson type method with adaptive time steps that can handle the nonlinearity of the di erential equation; since a few tens or even hundreds of forward problems have to be solved during the optimization of the di usion coecient, the output least squares method has a tendency to become very expensive. Equation error methods, on the other hand, use the given data as input for the di erential equation (1.1), in order to solve it for a. This idea goes back at least to This research was partly supported by the Christian{Doppler{Society, Austria. y Fachbereich Mathematik, Universit at Kaiserslautern D{67653 Kaiserslautern, and
Institut fur Praktische Mathematik, Universitat Karlsruhe, D{76128 Karlsruhe, Germany. z Institut f ur Industriemathematik, Johannes-Kepler-Universitat Linz, A{4040 Linz, Austria. 1
Richter [16] who employed an equation error scheme to identify the di usion coecient in the stationary case. The approach we take, however, is di erent from Richter's: Integrating (1.1) from xi?1 to xi yields (1.2)
?
?
a u(t; xi) ux (t; xi) ? a u(t; xi?1) ux(t; xi?1) =
Z xi
ut(t; ) d ;
xi?1
i = 1; : : :; m. We choose those m linear equations for a as the basic problem to solve. In Sects. 2 and 3 we investigate a proper setting for this system of equations, before we set up a numerical algorithm in Sects. 4 and 6, and present an error analysis under the assumption that the exact coecient ay has additional smoothness, i.e., ay 2 C 3; . Under the same assumption we also suggest an optimal spacing h of the thermocouples, depending on the accuracy of the measured data: if is an upper bound in L2(I ) for the data error in u( ; xi) for each i = 0; : : :; m, then the optimal spacing is given by h 1=4. On the basis of these results we present in Sect. 5 a heuristic argument which indicates that asymptotically the error in H 2 of the computed reconstruction of the di usion coecient will decay roughly like 1=4 as ! 0. A numerical example is given at the end of the paper. To conclude this introduction we mention that the inverse problem of reconstructing a(u) in (1.1) has already been considered in somewhat di erent settings in the literature, e.g., in [4, 5, 8, 9]; in those works data are either given in the whole domain I , or are restricted entirely to the boundary, like additional measurements of the
ux, for example. Equation error methods for quasilinear parabolic problems similar to (1.1) seem to have only been employed by Hinestroza and Murio [8] and Cannon and DuChateau [3]; this latter reference is also the one which is closest to the present work. Aside of these two papers equation error methods have mainly been investigated in the context of linear elliptic problems, cf., e.g., [2, 10, 16, 21] and references therein, to name just a few works. 2. The forward problem. Let f and g in (1.1) be continuously di erentiable functions in I with (2.1)
f(0) = g(0) = !0 and f 0 (0) = g0(0) = 0 :
Denote by
U := ff(t) : t 2 I g [ fg(t) : t 2 I g
(2.2)
the interval of all boundary values f and g, and assume that a 2 H 2 (U ). Then the boundary value problem (1.1) has a classical solution u which is continuously di erentiable with respect to t and two times continuously di erentiable with respect to x. Furthermore, u satis es a maximum principle, and hence, all values u(t; x) belong to U for t 2 I and x 2 . Details are provided in Appendix A. Given this solution u we de ne linear operators Ti : C(U ) ! C(I ) via Ti : a 7! yi ;
(2.3) where (2.4)
?
i = 1; : : :; m ;
?
yi = a u( ; xi) ux ( ; xi) ? a u( ; xi?1) ux( ; xi?1) : 2
Furthermore, we set ? (2.5) T : C(U ) ?! C(I ) m ; T : a 7! (T1 a; : : :; Tm a) : As we will see, Ti and T are continuous operators between appropriate Hilbert spaces. Proposition 2.1. Given? u; v2 L2(I ) with u(t) 2 U for almost every t 2 I we de ne the operator S : a 7! a (u() v(). Then S is a continuous linear operator from C(U ) to L2(I ). Proof. For almost every t 2 I and a; ~a 2 C(U ) we conclude that ja~(u(t)) ? a(u(t))j ka~ ? ak C (U ) ; and hence,
kS(~a ? a)k 22 =
Z T? 0
~a(u(t)) ? a(u(t)) 2 v2 (t) dt ka~ ? ak 2C (U ) kvk 2L2 (I ) :
This shows the continuity of S : C(U ) ! L2 (I ) with kS k C (U )!L2 (I ) kvk L2 (I ). As we have mentioned before the values of the solution u of (1.1) belong to U . Hence, we immediately obtain from Proposition 2.1 and the Sobolev embedding theorems (cf., e.g., Adams [1]) the following corollary: Corollary 2.2. Let u be the solution of the di erential equation (1.1) with di usion coecient ay 2 H 2 (U ) and ux be its partial derivative with respect to the space variable. Then ? the operator T de ned in (2.3)-(2.5) is continuous as an operator from H 2(U ) to L2(I ) m . As? a continuous operator between Hilbert spaces T admits an adjoint operator T : L2 (I ) m ! H 2 (U ). We give an analytic representation of this operator in Appendix B. Note that in the formulation of Corollary 2.2 we have chosen H 2 (U ) for the domain of T because this is the Hilbert space in which the di usion coecient ay , which we want to reconstruct, lives. If we were merely interested in continuity properties of T we could of course relax the assumptions on the domain of T. In practice we do not have the exact solution u of (1.1) nor its partial derivative ux . We therefore have to study the impact of perturbations in u and v on the operator S de ned in Proposition 2.1 Theorem 2.3. Assume that v 2 C(I ) but instead of u and v we are given perturbed functions u~; v~ 2 L2(I ) with ku~ ? uk L2(I ) 1 and kv~ ? vk L2 (I ) 2 . Assume further that the values u(t) and u~(t) belong to U for almost every t 2? I . De ne S as in Proposition 2.1, and consider the perturbed operator S~ : a 7! a (~u() v~(). Then
?
kS~ ? S k H 2 (U )!L2 (I ) c kvk C (I ) 1 + 2 ;
where c is a constant depending only on the embedding operator from H 2 (U ) into C 1 (U ). Proof. For almost any xed t 2 I and any a 2 H 2 (U ) C 1 (U ) we have
? ? Z u~(t) 0 a u~(t) ? a u(t) = a (w) dw ka0k C(U ) u~(t) ? u(t) ;
and hence,
u(t)
Z T? Z T ? ? 2 2 0 2 2 a u~(t) ? a u(t) v (t) dt ka k C (U ) kvk C (I ) u~(t) ? u(t) 2 dt 0
0
= ka0k 2C (U ) kvk 2C (I ) ku~ ? uk 2L2(I ) : 3
It follows that ~ ? Sak L2 (I ) = ka(~u)~v ? a(u)vk L2 (I ) kSa
?
k a(~u) ? a(u) vk L2 (I ) + ka(~u)(~v ? v)k L2 (I ) kak C 1(U ) kvk C (I ) 1 + kak C (U ) 2 : This implies the theorem with a certain constant c which only depends on the embedding operator from H 2(U ) into C 1 (U ). Note that if S in Theorem 2.3 denotes any one of the two terms on the right-hand side of (2.4) then the assumptions of the theorem are ful lled because the solution u of (1.1) is continuously di erentiable with respect to x with values in U . Consequently, given only perturbed data u~ for the solution u (where the values of u~ are restricted to U ), we obtain an approximate operator T~ such that (2.6)
?
kT~ ? T k H 2 (U )!(L2 (I ))m O m(1 + 2 ) ;
if we can guarantee bounds 1 and 2 for u~i and some v~i , i = 0; : : :; m, such that (2.7) ku~i ? u( ; xi)k L2 (I ) 1 ; kv~i ? ux ( ; xi)k L2 (I ) 2 : 3. The inverse problem. As already described in the introduction the exact di usion coecient ay satis es the linear system (1.2) which we can now rewrite as (3.1) Tay = y ;
?
where T : H 2(U ) !? L2(I) m is the operator (2.5) from the previous section and y is an element from L2(I ) m whose ith component is the function
(3.2)
yi (t) =
Z xi
ut (t; ) d 2 L2 (I ) ;
xi?1
i = 1; : : :; m :
Given inexact data u~i we obtain an approximate operator T~ instead of T as indicated in the previous section, but in addition, we can only compute an approximate righthand side y~ for (3.1), e.g., by using a quadrature formula for the evaluation of the integral in (3.2). This will give us another error to deal with, i.e., we assume that (3.3) ky~i ? yi k L2 (I ) 3 ; i = 1; : : :; m : The inverse problem now consists in solving the perturbed problem ~ = y~ (3.4) Ta to nd an approximation a of ay. Recall from (2.5) that the range of T consists only of m-tuples of continuous functions, even continuous di erentiable functions if ay 2 ?C 2; , say (see Appendix A). In other words, the range of T embeds compactly into L2 (I ) m which shows that the inverse problem is ill-posed: in general, (3.4) need not have a solution, and even if it does, this solution might be arbitrarily far o (in H 2) from the true coecient ay. Therefore one has to apply regularization methods to this problem, a variety of which are discussed in [7, 12], for example. For reasons that we shall explain in Sect. 6 4
the method of conjugate gradients with early termination is very appealing for this particular problem. By this we refer to the conjugate gradient iteration applied to the normal equation system of (3.4), i.e., ~ = T~ y~ ; T~ Ta and regularization is only incorporated by terminating the iteration with an appropriate stopping rule, cf. (3.5) below. Note that the analysis of the conjugate gradient method as given in [7, 12] only deals with the case of errors in y but not with errors in T; however, rewriting the perturbed problem as ~ y = Tay + (T~ ? T)ay = y~ + (y ? y~) + (T~ ? T)ay =: y^ ; Ta ~ y = y^, we conclude from (2.6) and (3.3) that ay is the exact solution of the problem Ta where
ky^ ? y~k (L2 (I ))m m3 + c kayk H 2 (U ) m(1 + 2 ) : Here, c again only depends on the Sobolev embedding constants. An appropriate stopping criterion for the conjugate gradient method is the discrepancy principle: for a given > 1 we terminate the conjugate gradient method with the kth iterate ak , as soon as ~ k k (L2(I ))m ky^ ? y~k (L2 (I ))m (3.5) ky~ ? Ta for the rst time. It is shown in [7, 12] (a more sophisticated investigation of the perturbed operator case is contained in the original analysis by Nemirovskii [14]) that this yields a regularization method; i.e., as the data errors in T~ and y~ go to zero the resulting reconstructions of the di usion coecient tend to a solution of (3.1). Of course, in practice we only have rough estimates for j , j = 1; 2; 3 (see the following section), and we therefore terminate the iteration as soon as ~ k k (L2 (I ))m with := m maxf1 ; 2; 3 g : ky~ ? Ta (3.6)
4. Numerical approximation and error analysis. We will now investigate in more detail the magnitude of the data errors 1 ; 2, and 3 , under the assumption that ay 2 C 3; (U ), i.e., ay is three times di erentiable and (ay)000 is Holder continuous with some exponent 2 (0; 1). Moreover, we shall assume that f and g have Holder continuous second derivatives with f 00 (0) = g00 (0) = 0. As we indicate in Appendix A 2 4 @ @ this implies that @x4 u and @t2 u are continuous in I . We can therefore make use of the seminorms Z @4
2 1=2 ; kuk L2(I ;C 4 ( )) := @x u(t; : ) 4 1 I kuk H 2(I ;C ( )) :=
Z @2
2 u(t; : )
21 1=2 : @t I
In the sequel c always denotes a generic positive constant independent of the problem under consideration. 5
Of course 2 and 3 depend on 1 itself, but also on the particular way in which we numerically approximate ux ( ; xi) and the right-hand side functions yi of (3.2). This shall be described next. For simplicity we restrict our analysis to the case that the thermocouples form an equidistant mesh on with meshsize h, i.e., that xi ?xi?1 = h = 1=m for i = 1; : : :; m. Our numerical scheme is based on cubic interpolation of the given data u~i (t) u(t; xi) for xed t 2 I . It is well-known (cf., e.g., Stoer and Bulirsch [18, Sect. 2.4]) that for xed t 2 I there is a unique interpolating cubic spline s(t; ) satisfying 8 s(t; xi ) = u(t; xi) ; i = 0; : : :; m ; > > < 1 ??11u(t; x ) + 18u(t; x ) ? 9u(t; x ) + 2u(t; x ) ; (4.1) > sx (t; x0) = 6h 0 1 2 3 > ? 1 : sx(t; xm) = 11u(t; xm) ? 18u(t; xm?1) + 9u(t; xm?2) ? 2u(t; xm?3) ; 6h and a corresponding cubic spline s~(t; ) interpolating the data u~i (t), i = 0; : : :; m. For the splines in (4.1) the following error estimate holds, cf. Swartz and Varga [19],
@4 u(t; )
h4 ; js(t; ) ? u(t; )j c @x 4 1 (4.2) 4
@ u(t; )
h3 ; jsx (t; ) ? ux(t; )j c @x 4 1 valid for all 2 I . It remains to investigate the e ect of data errors in u~. To this end we rst observe that s~ ? s interpolates u~ ? u, and recall that the moments Mi (i.e., the second derivatives of (~s ? s)(t; ) at the mesh points xi ) are obtained from solving the linear system 0 1 1=2 0 1 0 M0 1 0 d0 1 CC BB ... CC BB ... CC B 1=4 1 1=4 B CC BB . CC BB . CC B ... B (4.3) 1=4 1 C BB .. CC = BB .. CC ; B B . . . . . . 1=4 C CA B@ .. CA B@ .. CA B @ . . M d 0 1=2 1 m m where the di are certain di erence quotient approximations for the second derivatives of (~u ? u)(t; ); moreover, from their explicit de nition which can be found in [18] one easily obtains the upper bound j +1 X ju~k (t) ? u(t; xk)j ; jdj j 3 h2 k=j ?1 valid for j = 1; : : :; m ? 1, while for j = 0 and j = m, 3 m ju~ (t) ? u(t; x )j X X k k : jd0j 9 ju~k (t) ?h2u(t; xk )j ; jdmj 9 2 h k=0 k=m?3 Now, using the quasi-local property of interpolating cubic splines (cf., e.g., [20, Sect. 14.7]) we obtain for some c > 0 1 1 X X ?j j j (4.4) jMij c 2 jdi?j j c~ 2?jj j ju~i?j (t) ?h2u(t; xi?j )j ; jj j=0 jj j=0 6
where the terms in the series are understood to be zero when i ? j 2= f0; : : :; mg. From this and formula (2.4.2.2) in [18] it follows readily that
j(~s ? s)x (t; xi)j c
(4.5)
1 X
jj j=0
2?jj j ju~i?j (t) ?hu(t; xi?j )j :
Now, since u~i is a function in L2(I ) we can integrate (4.2) and (4.5) to obtain ks~x ( ; xi) ? ux ( ; xi)k L2(I )
ks~x ( ; xi) ? sx ( ; xi)k L2(I ) + ksx ( ; xi) ? ux( ; xi)k L2(I ) c
1 X
jj j=0
ku~ ? u( ; xi?j )k L2 (I ) + h3 kuk 2 4 2?jj j i?j L (I ;C ( )) h
c h1 + h3 kuk L2(I ;C 4 ( )) ;
recall that c is used throughout as a generic positive constant. Thus, if we take v~i = s~x ( ; xi) in Section 2 then we can take 2 = c h3 kuk L2(I ;C 4 ( )) + h1 in (2.7). This expression becomes minimal when h 11=4, in which case
2 13=4 ; 1 ! 0 : (4.6) Next we turn our attention to a numerical approximation of the right-hand side functions yi of (3.2), i = 1; : : :; m. The idea is to rewrite Z xi d (4.7) u(t; ) d ; yi (t) = dt Yi (t) with Yi (t) = xi?1 and to replace the integrand u by the spline approximation s~ computed above, i.e., we de ne Z xi d ~ ~ (4.8) s~(t; ) d : y~i (t) = dt Yi (t) with Yi (t) = xi?1 In the same way as before we obtain from (4.2) and (4.4) that
jY~i(t) ? Yi (t)j
Z xi
js(t; ) ? u(t; )j d +
xi?1
4
Z xi
js~(t; ) ? s(t; )j d
xi?1
1 X
@ u(t; ) + h c h5 @x 2?jj jju~i?j (t) ? u(t; xi?j )j 4 1 jj j=0 and it follows that for h 11=4, ? (4.9) kY~i ? Yi k L2(I ) c h5 kuk L2(I ;C 4 ( )) + h1 c15=4 : 7
Recall that utt is continuous: therefore we have @2 Z xi @2 2
dt2 Yi(t) = x dt2 u(t; ) d h dt@ 2 u(t; ) 1 i?1 for almost every t 2 I , and hence, (4.10) Yi 2 H 2(I ) with kYi k H 2 (I ) h kuk H 2(I ;C ( )) : Now, extending Y~i by zero for negative t, we can use backward di erences to de ne ~ ~ y~i (t) := Yi(t) ? Yi (t ? ) ; a.e. t 2 I : The precise value of > 0 will be speci ed later. Using a similar analysis as is exempli ed in [7, Sect. 1.1] for the maximum norm we obtain ky~i ? yi k L2 (I ) c 1 kY~i ? Yi k L2 (I ) + kYi k H 2 (I ) : Consequently, for 15=8=h1=2 we obtain from (4.9) and (4.10) the bound
ky~i ? yi k L2 (I ) ch1=215=8 ;
so that, for h 11=4, we can choose 3 in (3.3) to be
(4.11) 3 13=4 ; 1 ! 0 : Concerning problem (3.1) we can now summarize our ndings as follows. Theorem 4.1. Let f; g 2 C 2; (I ) for some 2 (0; 1) satisfy (2.1) and f 00 (0) = 00 g (0) = 0, and assume that the di usion coecient ay belongs to C 3; (U ). If the measurements u~i at the individual thermocouples xi have values in U and satisfy ku~i ? u( ; xi)k L2 (I ) 1 ; i = 0; : : :; m ; then the optimal spacing between the thermocouples, h 11=4 , yields a perturbed linear system (3.4) instead of (3.1) with
p
p
kT~ ? T k (L2(I ))m c 1 ;
ky~ ? yk (L2 (I ))m c 1 ;
where c depends only on the exact temperature distribution u.
5. Towards error bounds for the computed reconstructions. Running ~ = y~ as described in Section 3 the conjugate gradient method on the linear system Ta we obtain an approximation a~ of the true di usion coecient ay. In this section we give a heuristic argument to support the claim that asymptotically ka~ ? ay k H 2 (U ) 11=4 h might hold provided that ay 2 W 4;1(U ). In fact, if we know that ay 2 R(T ) then the conjugate gradient theory in [7, 12, (5.1)
14] predicts an error bound (5.2)
p ka~ ? ayk H 2 (U ) = O( ) ; 8
where = m maxf1 ; 2 ; 3g is the same as in (3.6). Essentially, the condition ay 2 R(T ) or R(T~ ) requires ay to belong to a suitable Sobolev space, namely ay 2 W 4;1(U ), see Appendix B. Note that W 4;1(U ) C 3 (U ) (cf., e.g., [1, Thm. 5.4]) but it is not clear whether the third derivative of a function ay 2 R(T ) will also be Holder continuous. In fact, we cannot guarantee that an arbitrary function ay 2 C 3; will belong to R(T~ ), nor vice versa. Nevertheless, on the basis of the previous arguments one might guess that the gap between these two linear spaces is pretty small. In view of Theorem 4.1 this indicates that, essentially,
p
1 ; so that (5.2) coincides with (5.1) as desired. 6. Numerical implementation. In Section 4 we have described how we can ~ = y~ of (3.4). obtain v~i and y~i to set up the linear system Ta For our numerical computations we choose the discrete values a(u) at the integers u 2 U \ Z as unknowns, and stack them into a vector a. For each time t = tj for which data are available, and for each i = 1; : : :; m, the identity (1.2) yields one equation for these unknowns where, for simplicity, we round the arguments u(tj ; xi) and u(tj ; xi?1) of a to the nearest integers. In this way each such equation couples at most two of the entries of vector a, and we obtain a linear system (6.1) Aa = y for a with an extremely sparse matrix A which has at most two nonzeros (i.e., values of the approximate spatial derivatives v~i of u) per row. The vector y in (6.1) contains the respective values yi (tj ) from the right-hand side of (1.2). To take the required smoothness a 2 H 2 (U ) into account we use nite di erences to approximate
kak 2H 2 (U ) aT B a with
0 1 ?2 1 1 BB ?2 5 ?4 1 CC BB 1 ?4 6 ?4 1 CC 1 B CC + B = (u) B BB 1 ?4 6 ?4 1 C C @ 1 ?4 5 ?2 A 3
1 ?2 1
0 1 ?1 1 BB ?1 2 ?1 CC B CC ? 1 2 1 1B B C + u I : u B ?1 2 ?1 C B@ C ?1 2 ?1 A ?1 1
Here, the role of u is crucial: A rescaling of the solution u of the di erential equation (1.1) by a factor yields a temperature distribution u^ := u corresponding to a di usion coecient ^a(^u) := a(u); rescalings do not a ect the solution a of the linear system (6.1) but they do a ect the H 2 norm of the corresponding di usion coecients. It would be a desirable feature of the numerical inversion scheme to yield reconstructions of a that are independent of such rescalings. This is achieved, for example, by choosing u as the ratio of the scale of the a-axis over the scale of the u-axis. For instance, in the example of the following section, the di usion coecient is of the order of one whereas the temperature ranges between 20 and 1100 [ C]; we therefore take u := 1=1100 in this example. 9
Computing a Cholesky decomposition B = LLT of B we obtain a `standard form discrete ill-posed problem' (cf., e.g., [7, Chapter 9]) (6.2)
AL?1 b = y ;
a = L?1 b ;
for which we have to nd a regularized solution b. Although A and L are quite large they have at most three nonzeros per row. Therefore, iterative regularization methods like, for example, the conjugate gradient method are the most appropriate choice. Note that each iteration requires multiplications with A and AT , and the solution of one linear system with L and LT , respectively. Therefore the total amount of work per iteration is bounded by the number of nonzeros in A and L and so, because of what we have said before, each iteration takes only O(n) operations, where n is the size of vector a (in our example, n = 1100). Since only a few conjugate gradient iterations are necessary to obtain a good approximate solution a (see below) this yields an extremely cheap algorithm for the reconstruction of the di usion coecient. 7. Numerical results. Consider the following experiment: At time t = 0 hot metal is poured into the form sand at x = ?5 [cm] and temperature measurements are taken at m + 1 equidistant thermocouples located at xi = 20i=m [cm], i = 0; : : :; m. The temperature measurements at x0 = 0 and xm = 20 are considered as Dirichlet data for the di erential equation (1.1). The initial temperature is !0 = 20 [C]. To simulate this experiment we have chosen a phantom di usion coecient ay after consulting the collaborating engineers who also provided realistic Dirchlet data f and g for (1.1). The di erential equation was solved up to time T = 6 104 [s] with the subroutine ode15s from Matlab's ode suite, cf. Shampine and Reichelt [17]; to this end we used a semidiscrete approach with 121 piecewise linear nite elements for the discretization of the spatial domain. To study the e ect of measurement errors, the data have been perturbed by additive noise of up to 1%. A crucial part of the algorithm is the numerical di erentiation of the data u( ; xi) with respect to time. Rather than taking nite di erences we have interpolated u( ; xi) by cubic splines on relatively coarse grids, and then di erentiated the splines. For the linear system (6.1) only the equations (1.2) corresponding to i = 1; : : :; m and the time steps t = tj chosen by ode15s are considered. Still, not all of these data were actually used. For some reason the stability of our inversion scheme seems to improve when we throw away data corresponding to the initial time phase of the experiment. (We believe that the steep gradients ux in this area are badly captured by the interpolating splines; possibly, a nonequidistant distribution of the thermocouples might be more appropriate for this purpose.) As a consequence only the data within the shaded region in Figure 7.1 have been used for the reconstruction. They correspond to 22 time steps and thus lead to a linear system with 22 m equations for 1100 unknowns. We have run the algorithm with m = 5; 8, and 15, and hence, all linear systems have been strongly underdetermined. This is a problem that is taken care of by the regularization via early termination of the conjugate gradient iteration. Figure 7.2 shows the iteration history in the case of 0:1% additive noise. Note that the right-hand side plot shows the errors only in L2 , since the phantom ay itself is not H 2 , see Figure 7.3 below. While the residual is monotonically decreasing (as it should) the error history exhibits the typical semiconvergence phenomenon in illposed problems: the iterates seem to converge in the beginning before they start to 10
1000 800 600 400 200 0 0
1
2
3
4
5
6e4
Fig. 7.1. Temperature histories at nine thermocouples (m = 8) Relative Residuals
Relative Errors
1
1
m=5 m=8 m = 15
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0 0
5
10
m=5 m=8 m = 15
0.8
15
0 0
20
5
10
15
20
Fig. 7.2. Iteration history: residuals and iteration errors vs. iteration count (0:1% noise)
diverge eventually. In most cases the best approximation is obtained after just a few iterations. Reconstructions for several m and various noise levels are shown in Figure 7.3. We observe that the algorithm is able to reconstruct the phantom qualitatively correct provided that m is appropriately coupled to the noise level, . In addition to these gures Table 7.1 shows the number of iterations and the corresponding (L2 ) errors of the optimal approximations. A note on these numbers is in order: since the temperature data within the shaded region of Figure 7.1 only range between 20 and 940 [ C] it seems to be fair just to measure the errors in this interval (20; 940) rather than in the whole interval shown in Figure 7.3. The boldfaced entries in the table correspond to the optimal number of thermocouples depending on the amount of noise; the more noise the less thermocouples should be taken. Moreover, it can be seen that a large number of thermocouples (i.e., m = 15) is only useful with no or very little noise; with 1% noise the reconstruction corresponding to m = 15 is completely wrong, and is therefore omitted in Figure 7.3. It should be stressed that the numbers for the amount of noise must not be 11
no noise 0.65 0.6 0.55
0.1% noise 0.65
ay m=5 m=8 m = 15
0.6 0.55
0.5
0.5
0.45
0.45
0.4
0.4
0.35 0
200
400
600
800
0.35 0
1000
ay m=5 m=8 m = 15
200
400
0.3% noise 0.65 0.6 0.55
ay m=5 m=8 m = 15
0.6 0.55 0.5
0.45
0.45
0.4
0.4
200
400
600
800
1000
800
1000
1% noise 0.65
0.5
0.35 0
600
800
0.35 0
1000
ay m=5 m=8
200
400
600
Fig. 7.3. Reconstructions for various noise levels Table 7.1
Relative errors and iteration counts
m=5 m=8 m = 15
no noise 0.0223 6 its.
0:1% noise 0.0254 6 its.
24 its. 0.0071 38 its.
7 its. 0.0123 49 its.
0.0046
0.0118
0:3% noise
0.0347
5 its. 0.0436 3 its. 0.2662 4 its.
1% noise
0.0461
5 its. 0.2165 3 its. 0.7062 1 it.
Table 7.2
Discrepancy k y ? Aay k 2 = k y k 2
m=5 m=8 m = 15
no noise 0.1468 0.0583 0.0815
0:1% noise 0.1549 0.0759 0.1750 12
0:3% noise 0.1806 0.1853 0.4690
1% noise 0.2750 0.4891 1.5001
confused with the discrepancy of the linear system Aa = y of (6.1), i.e., the relative t ky ? Aay k 2= kyk 2 given the exact di usion coecient ay . To emphasize this point we have included those discrepancies in Table 7.2. As can be seen from this table, even with no noise added upon the temperature data the linear system (6.1) su ers from a discretization error between 5 and 15% of the actual right-hand side. We nally remark that with more than 0:3% noise the reconstructions do not match perfectly with the given phantom, but still it is possible to see the local minimum of ay and to detect the sharp bend near u = 600. On the other hand, those reconstructions are extremely cheap to compute: The whole algorithm including the setup of the linear system and seven conjugate gradient iterations takes only between 0:3 (m = 5) and 0:5 M ops (m = 15) in Matlab, with roughly 75% of this work amount being spent on solving the linear system. For comparison, the solution of the parabolic di erential equation (1.1) for the simulation of data (which is precisely the amount of work for solving the forward problem, i.e., the building block of all output least squares type methods for reconstructing the di usion coecient) took already more than 4 M ops with 121 nite elements, and still 2 M ops with 64 nite elements. We conclude this section with a comment on matrix A. As it turns out, A itself is not terribly ill-conditioned, i.e., the condition number of A remained between 400 and 700, independent of m and of the actual noise level. However, when m gets larger (in our experiments, when m = 8 or 15) then A becomes rank de cient. As we have indicated already in Sect. 3 the ill-posedness of the linear system is mostly due to the imbedding operator into H 2 , corresponding to the matrix L?1 in the discrete setting (6.2); in fact, cond (AL?1 ) 107, and hence, (6.2) is a severely ill-conditioned problem. Appendix A. We brie y recollect some of the results from the book by Ladyzenskaja, Solonnikov, and Ural'ceva [11] concerning the solution u of the forward problem (1.1). Note that a solution of (1.1) satis es a maximum principle (cf. Protter and Weinberger [15, Sect. 3.2]), i.e., the values of u over I belong to the set U de ned in (2.2). We set
?
'(t; x) := f(t) + x g(t) ? f(t) ; and transform (1.1) to
?
vt = a(v + ')vx x + a0 (v + ')(vx + 'x )'x ? 't ; v(0; x) = 0 ; v(t; 0) = 0 ; v(t; 1) = 0 ;
(7.1)
where u = v + '. From Thm. 6.1 in Chapter V of [11] it follows that (7.1) has a solution from the class C2+ of functions u with continuous derivatives ut, ux and uxx in I , which in addition have certain Holder continuity properties if one of the two arguments (t; x) is kept xed. Assume now that a 2 C 2; (U ) and that the Dirichlet data f; g 2 C 2; (I ) satisfy the compatibility condition (2.1), and f 00 (0) = g00 (0). Then, using the fact that w = vt is a solution of the di erential equation
?
wt = a(u)wx x + a0 (u)ux wx + 2a0(u)ux 'xt + a0 (u)utvxx + a00 (u)utu2x ? 'tt ; w(0; x) = 0 ; w(t; 0) = 0 ; w(t; 1) = 0 ; 13
we conclude from Thm. 5.2 in Chapter IV of [11] that vt also belongs to the class C2+ , and hence, utt is continuous in I . With similar arguments (di erentiating (7.1) with respect to x) we obtain that vx and vxx belong to C2+ provided that a 2 C 3; , and hence, @x@ 44 u is continuous in I . ? Appendix B. Here we derive the adjoint operator T : L2(I) m ! H 2 (U ) of the operator T of (2.5). To this end we consider once again the operator ? S : H 2 (U ) ?! L2 (I ) ; S : a 7! a (u() v() ; with u; v 2 L2(I ) and u has values in U almost everywhere. We also introduce the operator R : L2 (U ) ?! H 2 (U ) ; R : b 7! w ; where w 2 H 2(U ) is the weak solution of the di erential equation w ? w00 + w0000 = b with natural boundary conditions. In other words, we have hR a; biL2(U ) = ha; RbiH 2 (U ) = ha; biL2 (U ) for all a 2 H 2(U ) and all b 2 L2(U ), i.e., the adjoint R of R is the embedding operator : H 2(U ) ,! L2(U ). Note, however, that R extends in the usual way to a continuous operator R : H ?2(U ) ! H 2 (U ). Assume now for the moment that u 2 C 1(I ) and consider a time interval I1 I where u0 6= 0. Denote by U1 the interval of function values u(I1). Multiplying Sa by ' 2 L2(I ) we obtain Z Z Z (Sa)(t)'(t) dt = a(u(t))v(t)'(t) dt = a(u) v(t(u))'(t(u)) j u0(t(u)) j du ; I1 I1 U1 and we note that 8 v(t(u))'(t(u)) < u 2 U1 ; b1 (u) := : j u0(t(u)) j (7.2) 0 u 62 U1 ; belongs to L1(U ), since Z v(t(u))'(t(u)) Z Z Z (7.3) jb1(u)j du = jb1(u)j du = u0 (t(u)) du = jv(t)'(t)j dt : U1 U1 U I1 We proceed in a similar way for all further subintervals I2; I3; : : : of I in which u is monotone, thus obtaining further L1 functions b2; b3; : : :. In the remaining intervals D1 ; D2; ::: of I we have that u is a constant, i.e., ujDj = !j and
Z
Dj
(Sa)(t)'(t) dt =
Z
Dj
a(!j )v(t)'(t) dt = a(!j ) j
This yields
Z I
(Sa)(t)'(t) dt =
1 Z X i=1 Ii
(Sa)(t)'(t) dt + 14
Z
with j := v(t)'(t) dt : Dj
1Z X j =1 Dj
(Sa)(t)'(t) dt
= =
1Z X
a(u)bi (u) du +
i=1 U
Z
U
a(u)
1 X i=1
bi(u) du +
1 X j =1
j a(!j )
1 X j =1
j a(!j ) :
P b (u) converges in L1(U ), since From (7.3) then follows that the function b(u) := 1 i=1 i Z
U
jb(u)j du
1 Z X
i=1 U
jbi (u)j du =
1Z X
i=1 Ii
jv(t)'(t)j dt =
Z
I
jv(t)'(t)j dt
kvk L2 (I ) k'k L2 (I ) : P Since H s(U ) LP1 (U ) for s > 1=2 we have L1(U ) H ?2(U ). Similarly, j j j j < 1 and hence d := j j (!j ) converges as a distribution in L1(U ); here, (!j ) denotes the {distribution at !j . It therefore follows that R(b + d) 2 H 2 (U ), and hence,
hSa; 'iL2 (I ) =
Z
I
(Sa)(t)'(t) dt =
Z
U
a(u)(b(u) + d(u)) du = ha; R(b + d)iH 2 (U ) :
Furthermore, since b + d 2 L1(U ) it follows from the di erential equation for w that w = R(b+d) 2 W 4;1(U ) which embeds continuously into C 3(U), cf., e.g., [1, Thm. 5.4]. Thus we have shown that (7.4) S ' = R(b + d) 2 W 4;1(U ) ; with b and d as de ned above. The general case of u 2 L2(I ) follows by continuity of the mapping u 7! S according to Theorem 2.3. We denote by Si , i = 0; : : :; m, the operator S which is obtained when u in (7.2) is the trace u( ; xi) of the solution u of (1.1) along x = xi , and likewise v in (7.2) is the trace of the partial derivative ux of u. Then it follows that Ti = Si ? Si?1 is the operator de ned in (2.3), (2.4), and Ti = Si ? Si?1 is its adjoint. The adjoint ? 2 operator T : L (I ) m ! H 2(U ) of T in (2.5) is then given as T ('1 ; : : :; 'm ) = T1 '1 + : : : + Tm 'm ; 'i 2 L2(I ) ; and it follows from (7.4) that R(T ) W 4;1(U ). Acknowledgements. We like to thank Martin Bruhl for pointing out a gap in an earlier version of the derivation of T in Appendix B. In addition, we are grateful to an anonymous referee for bringing reference [3] to our attention. REFERENCES [1] R. A. Adams, Sobolev Spaces, Academic Press, New York, San Francisco, Melbourne, 1975. [2] H. T. Banks and K. Kunisch, Estimation Techniques for Distributed Parameter Systems, Birkhauser Verlag, Boston, Basel, Berlin, 1989. [3] J. R. Cannon and P. DuChateau, Design of an experiment for the determination of an unknown coecient in a nonlinear conduction-di usion equation, Internat. J. Engrg. Sci., 25 (1987), pp. 1067-1078. 15
[4] G. Chavent and P. Lemonnier, Identi cation de la non-linearite d'une equation parabolique quasilineaire, Appl. Math. Optim., 1 (1974), pp. 121-162. [5] P. DuChateau, An inverse problem for the hydraulic properties of porous media, SIAM J. Math. Anal., 28 (1997), pp. 611-632. [6] S. Dummel and S. Handrock-Meyer, Uniqueness of the solution of an inverse problem for a quasilinear parabolic equation in divergence form, Z. Anal. Anwendungen, 7 (1988), pp. 241-246. [7] H. W. Engl, M. Hanke, and A. Neubauer, Regularization of Inverse Problems, Kluwer, Dordrecht, 1996. [8] D. Hinestroza and D. Murio, Identi cation of transmissivity coecients by molli cation techniques. Part I: One-dimensional elliptic and parabolic problems, Comput. Math. Appl., 25 (1993), pp. 59-79. [9] T. Karkkainen, A linearization technique and error estimates for distributed parameter identi cation in quasilinear problems, Numer. Funct. Anal. Optim., 17 (1996), pp. 345364. [10] T. Karkkainen, An equation error method to recover di usion from the distributed observation, Inverse Problems, 13 (1997), pp. 1033-1051. [11] O. A. Ladyzenskaja, V. A. Solonnikov, and N. N. Ural'ceva, Linear and Quasilinear Equations of Parabolic Type, Amer. Math. Soc., Providence, RI, 1968. [12] A. K. Louis, Inverse und schlecht gestellte Probleme, Teubner, Stuttgart, 1989. [13] R. Nabokov, An inverse problem for porous medium equation, in Parameter Identi cation and Inverse Problems in Hydrology, Geology and Ecology, J. Gottlieb and P. DuChateau, eds., Kluwer, Dordrecht, Boston, London, 1996, pp. 155-163. [14] A. S. Nemirovskii, The regularizing properties of the adjoint gradient method in ill-posed problems, USSR Comput. Math. Math. Phys., 26 (1986), pp. 7-16. [15] M. H. Protter and H. F. Weinberger, Maximum Principles in Di erential Equations, Prentice-Hall, Englewood Cli s, New Jersey, 1967. [16] G. R. Richter, Numerical identi cation of a spatially varying di usion coecient, Math. Comp., 36 (1981), pp. 375-386. [17] L. F. Shampine and M. W. Reichelt, The Matlab ode suite, SIAM J. Sci. Comput., 18 (1997), pp. 1-22. [18] J. Stoer and R. Bulirsch, Introduction to Numerical Analysis, Springer-Verlag, New York, Heidelberg, Berlin, 1980. [19] B. K. Swartz and R. S. Varga, Error bounds for spline and L-spline interpolation, J. Approx. Theory, 6 (1972), pp. 6-49. [20] E. E. Tyrtyshnikov, A Brief Introduction to Numerical Analysis, Birkhauser, Boston, 1997. [21] G. M. Vainikko, Inverse problem of ground water ltration: identi ability, discretization and regularization, in Inverse Problems in Di usion Processes, H. W. Engl and W. Rundell, eds., SIAM, Philadelphia, 1995, pp. 90-107.
16
Abstract. We consider the inverse problem of reconstructing the di usion coecient in a quasilinear parabolic di erential equation in divergence form from measurements of the solution at a nite number of points in the interior of the domain. An equation error method is developed which transforms the inverse problem into a system of linear operator equations for the di usion coecient, and which can be solved by the conjugate gradient method in a very ecient and stable manner. A detailed error analysis relates the required number of measurements with their accuracy. Numerical results illustrate the performance of the method. Key words. Inverse problems, quasilinearparabolicequation, parameter identi cation,ill-posed problems, regularization methods, conjugate gradient method AMS subject classi cations. 35R30, 65M30, 65M15
1. Introduction. We consider the quasilinear parabolic initial-boundary value
problem (1.1)
?
ut = a(u)ux x ; u(0; x) = !0 ; u(t; 0) = f(t) ; u(t; 1) = g(t) ;
with t 2 I = [0; T] and x 2 = [0; 1]. Here, a 2 H 2 is a positive function, !0 a constant and f and g are given boundary data. Given a, the direct problem consists in solving the di erential equation for u. Here we are concerned with the following inverse problem: Given u for certain abscissa 0 = x0 < : : : < xm = 1 and all t 2 [0; T], determine the di usion coecient a(u). We refer to Dummel and Handrock-Meyer [6] concerning the identi ability of an analytic di usion coecient from the given data. Inverse problems of similar sort arise in a variety of di erent contexts. This particular one occured in an industrial application where an ingot casting manufacturer was asking for the thermal conductivity of the form sand he is using in the manufacturing process. Due to di erent heat transport phenomena, the thermal conductivity of the sand is likely to change in the presence of very high temperatures. In this application the available data for the inverse problem is the temperature evolution measured by a few thermocouples placed at increasing distances from the heat source. The standard approach for the solution of this inverse problem is the output least squares method; cf., e.g., Nabokov [13] for an application of output least squares to this particular problem. This method requires a subroutine to solve the forward problem, e.g., a Crank-Nicolson type method with adaptive time steps that can handle the nonlinearity of the di erential equation; since a few tens or even hundreds of forward problems have to be solved during the optimization of the di usion coecient, the output least squares method has a tendency to become very expensive. Equation error methods, on the other hand, use the given data as input for the di erential equation (1.1), in order to solve it for a. This idea goes back at least to This research was partly supported by the Christian{Doppler{Society, Austria. y Fachbereich Mathematik, Universit at Kaiserslautern D{67653 Kaiserslautern, and
Institut fur Praktische Mathematik, Universitat Karlsruhe, D{76128 Karlsruhe, Germany. z Institut f ur Industriemathematik, Johannes-Kepler-Universitat Linz, A{4040 Linz, Austria. 1
Richter [16] who employed an equation error scheme to identify the di usion coecient in the stationary case. The approach we take, however, is di erent from Richter's: Integrating (1.1) from xi?1 to xi yields (1.2)
?
?
a u(t; xi) ux (t; xi) ? a u(t; xi?1) ux(t; xi?1) =
Z xi
ut(t; ) d ;
xi?1
i = 1; : : :; m. We choose those m linear equations for a as the basic problem to solve. In Sects. 2 and 3 we investigate a proper setting for this system of equations, before we set up a numerical algorithm in Sects. 4 and 6, and present an error analysis under the assumption that the exact coecient ay has additional smoothness, i.e., ay 2 C 3; . Under the same assumption we also suggest an optimal spacing h of the thermocouples, depending on the accuracy of the measured data: if is an upper bound in L2(I ) for the data error in u( ; xi) for each i = 0; : : :; m, then the optimal spacing is given by h 1=4. On the basis of these results we present in Sect. 5 a heuristic argument which indicates that asymptotically the error in H 2 of the computed reconstruction of the di usion coecient will decay roughly like 1=4 as ! 0. A numerical example is given at the end of the paper. To conclude this introduction we mention that the inverse problem of reconstructing a(u) in (1.1) has already been considered in somewhat di erent settings in the literature, e.g., in [4, 5, 8, 9]; in those works data are either given in the whole domain I , or are restricted entirely to the boundary, like additional measurements of the
ux, for example. Equation error methods for quasilinear parabolic problems similar to (1.1) seem to have only been employed by Hinestroza and Murio [8] and Cannon and DuChateau [3]; this latter reference is also the one which is closest to the present work. Aside of these two papers equation error methods have mainly been investigated in the context of linear elliptic problems, cf., e.g., [2, 10, 16, 21] and references therein, to name just a few works. 2. The forward problem. Let f and g in (1.1) be continuously di erentiable functions in I with (2.1)
f(0) = g(0) = !0 and f 0 (0) = g0(0) = 0 :
Denote by
U := ff(t) : t 2 I g [ fg(t) : t 2 I g
(2.2)
the interval of all boundary values f and g, and assume that a 2 H 2 (U ). Then the boundary value problem (1.1) has a classical solution u which is continuously di erentiable with respect to t and two times continuously di erentiable with respect to x. Furthermore, u satis es a maximum principle, and hence, all values u(t; x) belong to U for t 2 I and x 2 . Details are provided in Appendix A. Given this solution u we de ne linear operators Ti : C(U ) ! C(I ) via Ti : a 7! yi ;
(2.3) where (2.4)
?
i = 1; : : :; m ;
?
yi = a u( ; xi) ux ( ; xi) ? a u( ; xi?1) ux( ; xi?1) : 2
Furthermore, we set ? (2.5) T : C(U ) ?! C(I ) m ; T : a 7! (T1 a; : : :; Tm a) : As we will see, Ti and T are continuous operators between appropriate Hilbert spaces. Proposition 2.1. Given? u; v2 L2(I ) with u(t) 2 U for almost every t 2 I we de ne the operator S : a 7! a (u() v(). Then S is a continuous linear operator from C(U ) to L2(I ). Proof. For almost every t 2 I and a; ~a 2 C(U ) we conclude that ja~(u(t)) ? a(u(t))j ka~ ? ak C (U ) ; and hence,
kS(~a ? a)k 22 =
Z T? 0
~a(u(t)) ? a(u(t)) 2 v2 (t) dt ka~ ? ak 2C (U ) kvk 2L2 (I ) :
This shows the continuity of S : C(U ) ! L2 (I ) with kS k C (U )!L2 (I ) kvk L2 (I ). As we have mentioned before the values of the solution u of (1.1) belong to U . Hence, we immediately obtain from Proposition 2.1 and the Sobolev embedding theorems (cf., e.g., Adams [1]) the following corollary: Corollary 2.2. Let u be the solution of the di erential equation (1.1) with di usion coecient ay 2 H 2 (U ) and ux be its partial derivative with respect to the space variable. Then ? the operator T de ned in (2.3)-(2.5) is continuous as an operator from H 2(U ) to L2(I ) m . As? a continuous operator between Hilbert spaces T admits an adjoint operator T : L2 (I ) m ! H 2 (U ). We give an analytic representation of this operator in Appendix B. Note that in the formulation of Corollary 2.2 we have chosen H 2 (U ) for the domain of T because this is the Hilbert space in which the di usion coecient ay , which we want to reconstruct, lives. If we were merely interested in continuity properties of T we could of course relax the assumptions on the domain of T. In practice we do not have the exact solution u of (1.1) nor its partial derivative ux . We therefore have to study the impact of perturbations in u and v on the operator S de ned in Proposition 2.1 Theorem 2.3. Assume that v 2 C(I ) but instead of u and v we are given perturbed functions u~; v~ 2 L2(I ) with ku~ ? uk L2(I ) 1 and kv~ ? vk L2 (I ) 2 . Assume further that the values u(t) and u~(t) belong to U for almost every t 2? I . De ne S as in Proposition 2.1, and consider the perturbed operator S~ : a 7! a (~u() v~(). Then
?
kS~ ? S k H 2 (U )!L2 (I ) c kvk C (I ) 1 + 2 ;
where c is a constant depending only on the embedding operator from H 2 (U ) into C 1 (U ). Proof. For almost any xed t 2 I and any a 2 H 2 (U ) C 1 (U ) we have
? ? Z u~(t) 0 a u~(t) ? a u(t) = a (w) dw ka0k C(U ) u~(t) ? u(t) ;
and hence,
u(t)
Z T? Z T ? ? 2 2 0 2 2 a u~(t) ? a u(t) v (t) dt ka k C (U ) kvk C (I ) u~(t) ? u(t) 2 dt 0
0
= ka0k 2C (U ) kvk 2C (I ) ku~ ? uk 2L2(I ) : 3
It follows that ~ ? Sak L2 (I ) = ka(~u)~v ? a(u)vk L2 (I ) kSa
?
k a(~u) ? a(u) vk L2 (I ) + ka(~u)(~v ? v)k L2 (I ) kak C 1(U ) kvk C (I ) 1 + kak C (U ) 2 : This implies the theorem with a certain constant c which only depends on the embedding operator from H 2(U ) into C 1 (U ). Note that if S in Theorem 2.3 denotes any one of the two terms on the right-hand side of (2.4) then the assumptions of the theorem are ful lled because the solution u of (1.1) is continuously di erentiable with respect to x with values in U . Consequently, given only perturbed data u~ for the solution u (where the values of u~ are restricted to U ), we obtain an approximate operator T~ such that (2.6)
?
kT~ ? T k H 2 (U )!(L2 (I ))m O m(1 + 2 ) ;
if we can guarantee bounds 1 and 2 for u~i and some v~i , i = 0; : : :; m, such that (2.7) ku~i ? u( ; xi)k L2 (I ) 1 ; kv~i ? ux ( ; xi)k L2 (I ) 2 : 3. The inverse problem. As already described in the introduction the exact di usion coecient ay satis es the linear system (1.2) which we can now rewrite as (3.1) Tay = y ;
?
where T : H 2(U ) !? L2(I) m is the operator (2.5) from the previous section and y is an element from L2(I ) m whose ith component is the function
(3.2)
yi (t) =
Z xi
ut (t; ) d 2 L2 (I ) ;
xi?1
i = 1; : : :; m :
Given inexact data u~i we obtain an approximate operator T~ instead of T as indicated in the previous section, but in addition, we can only compute an approximate righthand side y~ for (3.1), e.g., by using a quadrature formula for the evaluation of the integral in (3.2). This will give us another error to deal with, i.e., we assume that (3.3) ky~i ? yi k L2 (I ) 3 ; i = 1; : : :; m : The inverse problem now consists in solving the perturbed problem ~ = y~ (3.4) Ta to nd an approximation a of ay. Recall from (2.5) that the range of T consists only of m-tuples of continuous functions, even continuous di erentiable functions if ay 2 ?C 2; , say (see Appendix A). In other words, the range of T embeds compactly into L2 (I ) m which shows that the inverse problem is ill-posed: in general, (3.4) need not have a solution, and even if it does, this solution might be arbitrarily far o (in H 2) from the true coecient ay. Therefore one has to apply regularization methods to this problem, a variety of which are discussed in [7, 12], for example. For reasons that we shall explain in Sect. 6 4
the method of conjugate gradients with early termination is very appealing for this particular problem. By this we refer to the conjugate gradient iteration applied to the normal equation system of (3.4), i.e., ~ = T~ y~ ; T~ Ta and regularization is only incorporated by terminating the iteration with an appropriate stopping rule, cf. (3.5) below. Note that the analysis of the conjugate gradient method as given in [7, 12] only deals with the case of errors in y but not with errors in T; however, rewriting the perturbed problem as ~ y = Tay + (T~ ? T)ay = y~ + (y ? y~) + (T~ ? T)ay =: y^ ; Ta ~ y = y^, we conclude from (2.6) and (3.3) that ay is the exact solution of the problem Ta where
ky^ ? y~k (L2 (I ))m m3 + c kayk H 2 (U ) m(1 + 2 ) : Here, c again only depends on the Sobolev embedding constants. An appropriate stopping criterion for the conjugate gradient method is the discrepancy principle: for a given > 1 we terminate the conjugate gradient method with the kth iterate ak , as soon as ~ k k (L2(I ))m ky^ ? y~k (L2 (I ))m (3.5) ky~ ? Ta for the rst time. It is shown in [7, 12] (a more sophisticated investigation of the perturbed operator case is contained in the original analysis by Nemirovskii [14]) that this yields a regularization method; i.e., as the data errors in T~ and y~ go to zero the resulting reconstructions of the di usion coecient tend to a solution of (3.1). Of course, in practice we only have rough estimates for j , j = 1; 2; 3 (see the following section), and we therefore terminate the iteration as soon as ~ k k (L2 (I ))m with := m maxf1 ; 2; 3 g : ky~ ? Ta (3.6)
4. Numerical approximation and error analysis. We will now investigate in more detail the magnitude of the data errors 1 ; 2, and 3 , under the assumption that ay 2 C 3; (U ), i.e., ay is three times di erentiable and (ay)000 is Holder continuous with some exponent 2 (0; 1). Moreover, we shall assume that f and g have Holder continuous second derivatives with f 00 (0) = g00 (0) = 0. As we indicate in Appendix A 2 4 @ @ this implies that @x4 u and @t2 u are continuous in I . We can therefore make use of the seminorms Z @4
2 1=2 ; kuk L2(I ;C 4 ( )) := @x u(t; : ) 4 1 I kuk H 2(I ;C ( )) :=
Z @2
2 u(t; : )
21 1=2 : @t I
In the sequel c always denotes a generic positive constant independent of the problem under consideration. 5
Of course 2 and 3 depend on 1 itself, but also on the particular way in which we numerically approximate ux ( ; xi) and the right-hand side functions yi of (3.2). This shall be described next. For simplicity we restrict our analysis to the case that the thermocouples form an equidistant mesh on with meshsize h, i.e., that xi ?xi?1 = h = 1=m for i = 1; : : :; m. Our numerical scheme is based on cubic interpolation of the given data u~i (t) u(t; xi) for xed t 2 I . It is well-known (cf., e.g., Stoer and Bulirsch [18, Sect. 2.4]) that for xed t 2 I there is a unique interpolating cubic spline s(t; ) satisfying 8 s(t; xi ) = u(t; xi) ; i = 0; : : :; m ; > > < 1 ??11u(t; x ) + 18u(t; x ) ? 9u(t; x ) + 2u(t; x ) ; (4.1) > sx (t; x0) = 6h 0 1 2 3 > ? 1 : sx(t; xm) = 11u(t; xm) ? 18u(t; xm?1) + 9u(t; xm?2) ? 2u(t; xm?3) ; 6h and a corresponding cubic spline s~(t; ) interpolating the data u~i (t), i = 0; : : :; m. For the splines in (4.1) the following error estimate holds, cf. Swartz and Varga [19],
@4 u(t; )
h4 ; js(t; ) ? u(t; )j c @x 4 1 (4.2) 4
@ u(t; )
h3 ; jsx (t; ) ? ux(t; )j c @x 4 1 valid for all 2 I . It remains to investigate the e ect of data errors in u~. To this end we rst observe that s~ ? s interpolates u~ ? u, and recall that the moments Mi (i.e., the second derivatives of (~s ? s)(t; ) at the mesh points xi ) are obtained from solving the linear system 0 1 1=2 0 1 0 M0 1 0 d0 1 CC BB ... CC BB ... CC B 1=4 1 1=4 B CC BB . CC BB . CC B ... B (4.3) 1=4 1 C BB .. CC = BB .. CC ; B B . . . . . . 1=4 C CA B@ .. CA B@ .. CA B @ . . M d 0 1=2 1 m m where the di are certain di erence quotient approximations for the second derivatives of (~u ? u)(t; ); moreover, from their explicit de nition which can be found in [18] one easily obtains the upper bound j +1 X ju~k (t) ? u(t; xk)j ; jdj j 3 h2 k=j ?1 valid for j = 1; : : :; m ? 1, while for j = 0 and j = m, 3 m ju~ (t) ? u(t; x )j X X k k : jd0j 9 ju~k (t) ?h2u(t; xk )j ; jdmj 9 2 h k=0 k=m?3 Now, using the quasi-local property of interpolating cubic splines (cf., e.g., [20, Sect. 14.7]) we obtain for some c > 0 1 1 X X ?j j j (4.4) jMij c 2 jdi?j j c~ 2?jj j ju~i?j (t) ?h2u(t; xi?j )j ; jj j=0 jj j=0 6
where the terms in the series are understood to be zero when i ? j 2= f0; : : :; mg. From this and formula (2.4.2.2) in [18] it follows readily that
j(~s ? s)x (t; xi)j c
(4.5)
1 X
jj j=0
2?jj j ju~i?j (t) ?hu(t; xi?j )j :
Now, since u~i is a function in L2(I ) we can integrate (4.2) and (4.5) to obtain ks~x ( ; xi) ? ux ( ; xi)k L2(I )
ks~x ( ; xi) ? sx ( ; xi)k L2(I ) + ksx ( ; xi) ? ux( ; xi)k L2(I ) c
1 X
jj j=0
ku~ ? u( ; xi?j )k L2 (I ) + h3 kuk 2 4 2?jj j i?j L (I ;C ( )) h
c h1 + h3 kuk L2(I ;C 4 ( )) ;
recall that c is used throughout as a generic positive constant. Thus, if we take v~i = s~x ( ; xi) in Section 2 then we can take 2 = c h3 kuk L2(I ;C 4 ( )) + h1 in (2.7). This expression becomes minimal when h 11=4, in which case
2 13=4 ; 1 ! 0 : (4.6) Next we turn our attention to a numerical approximation of the right-hand side functions yi of (3.2), i = 1; : : :; m. The idea is to rewrite Z xi d (4.7) u(t; ) d ; yi (t) = dt Yi (t) with Yi (t) = xi?1 and to replace the integrand u by the spline approximation s~ computed above, i.e., we de ne Z xi d ~ ~ (4.8) s~(t; ) d : y~i (t) = dt Yi (t) with Yi (t) = xi?1 In the same way as before we obtain from (4.2) and (4.4) that
jY~i(t) ? Yi (t)j
Z xi
js(t; ) ? u(t; )j d +
xi?1
4
Z xi
js~(t; ) ? s(t; )j d
xi?1
1 X
@ u(t; ) + h c h5 @x 2?jj jju~i?j (t) ? u(t; xi?j )j 4 1 jj j=0 and it follows that for h 11=4, ? (4.9) kY~i ? Yi k L2(I ) c h5 kuk L2(I ;C 4 ( )) + h1 c15=4 : 7
Recall that utt is continuous: therefore we have @2 Z xi @2 2
dt2 Yi(t) = x dt2 u(t; ) d h dt@ 2 u(t; ) 1 i?1 for almost every t 2 I , and hence, (4.10) Yi 2 H 2(I ) with kYi k H 2 (I ) h kuk H 2(I ;C ( )) : Now, extending Y~i by zero for negative t, we can use backward di erences to de ne ~ ~ y~i (t) := Yi(t) ? Yi (t ? ) ; a.e. t 2 I : The precise value of > 0 will be speci ed later. Using a similar analysis as is exempli ed in [7, Sect. 1.1] for the maximum norm we obtain ky~i ? yi k L2 (I ) c 1 kY~i ? Yi k L2 (I ) + kYi k H 2 (I ) : Consequently, for 15=8=h1=2 we obtain from (4.9) and (4.10) the bound
ky~i ? yi k L2 (I ) ch1=215=8 ;
so that, for h 11=4, we can choose 3 in (3.3) to be
(4.11) 3 13=4 ; 1 ! 0 : Concerning problem (3.1) we can now summarize our ndings as follows. Theorem 4.1. Let f; g 2 C 2; (I ) for some 2 (0; 1) satisfy (2.1) and f 00 (0) = 00 g (0) = 0, and assume that the di usion coecient ay belongs to C 3; (U ). If the measurements u~i at the individual thermocouples xi have values in U and satisfy ku~i ? u( ; xi)k L2 (I ) 1 ; i = 0; : : :; m ; then the optimal spacing between the thermocouples, h 11=4 , yields a perturbed linear system (3.4) instead of (3.1) with
p
p
kT~ ? T k (L2(I ))m c 1 ;
ky~ ? yk (L2 (I ))m c 1 ;
where c depends only on the exact temperature distribution u.
5. Towards error bounds for the computed reconstructions. Running ~ = y~ as described in Section 3 the conjugate gradient method on the linear system Ta we obtain an approximation a~ of the true di usion coecient ay. In this section we give a heuristic argument to support the claim that asymptotically ka~ ? ay k H 2 (U ) 11=4 h might hold provided that ay 2 W 4;1(U ). In fact, if we know that ay 2 R(T ) then the conjugate gradient theory in [7, 12, (5.1)
14] predicts an error bound (5.2)
p ka~ ? ayk H 2 (U ) = O( ) ; 8
where = m maxf1 ; 2 ; 3g is the same as in (3.6). Essentially, the condition ay 2 R(T ) or R(T~ ) requires ay to belong to a suitable Sobolev space, namely ay 2 W 4;1(U ), see Appendix B. Note that W 4;1(U ) C 3 (U ) (cf., e.g., [1, Thm. 5.4]) but it is not clear whether the third derivative of a function ay 2 R(T ) will also be Holder continuous. In fact, we cannot guarantee that an arbitrary function ay 2 C 3; will belong to R(T~ ), nor vice versa. Nevertheless, on the basis of the previous arguments one might guess that the gap between these two linear spaces is pretty small. In view of Theorem 4.1 this indicates that, essentially,
p
1 ; so that (5.2) coincides with (5.1) as desired. 6. Numerical implementation. In Section 4 we have described how we can ~ = y~ of (3.4). obtain v~i and y~i to set up the linear system Ta For our numerical computations we choose the discrete values a(u) at the integers u 2 U \ Z as unknowns, and stack them into a vector a. For each time t = tj for which data are available, and for each i = 1; : : :; m, the identity (1.2) yields one equation for these unknowns where, for simplicity, we round the arguments u(tj ; xi) and u(tj ; xi?1) of a to the nearest integers. In this way each such equation couples at most two of the entries of vector a, and we obtain a linear system (6.1) Aa = y for a with an extremely sparse matrix A which has at most two nonzeros (i.e., values of the approximate spatial derivatives v~i of u) per row. The vector y in (6.1) contains the respective values yi (tj ) from the right-hand side of (1.2). To take the required smoothness a 2 H 2 (U ) into account we use nite di erences to approximate
kak 2H 2 (U ) aT B a with
0 1 ?2 1 1 BB ?2 5 ?4 1 CC BB 1 ?4 6 ?4 1 CC 1 B CC + B = (u) B BB 1 ?4 6 ?4 1 C C @ 1 ?4 5 ?2 A 3
1 ?2 1
0 1 ?1 1 BB ?1 2 ?1 CC B CC ? 1 2 1 1B B C + u I : u B ?1 2 ?1 C B@ C ?1 2 ?1 A ?1 1
Here, the role of u is crucial: A rescaling of the solution u of the di erential equation (1.1) by a factor yields a temperature distribution u^ := u corresponding to a di usion coecient ^a(^u) := a(u); rescalings do not a ect the solution a of the linear system (6.1) but they do a ect the H 2 norm of the corresponding di usion coecients. It would be a desirable feature of the numerical inversion scheme to yield reconstructions of a that are independent of such rescalings. This is achieved, for example, by choosing u as the ratio of the scale of the a-axis over the scale of the u-axis. For instance, in the example of the following section, the di usion coecient is of the order of one whereas the temperature ranges between 20 and 1100 [ C]; we therefore take u := 1=1100 in this example. 9
Computing a Cholesky decomposition B = LLT of B we obtain a `standard form discrete ill-posed problem' (cf., e.g., [7, Chapter 9]) (6.2)
AL?1 b = y ;
a = L?1 b ;
for which we have to nd a regularized solution b. Although A and L are quite large they have at most three nonzeros per row. Therefore, iterative regularization methods like, for example, the conjugate gradient method are the most appropriate choice. Note that each iteration requires multiplications with A and AT , and the solution of one linear system with L and LT , respectively. Therefore the total amount of work per iteration is bounded by the number of nonzeros in A and L and so, because of what we have said before, each iteration takes only O(n) operations, where n is the size of vector a (in our example, n = 1100). Since only a few conjugate gradient iterations are necessary to obtain a good approximate solution a (see below) this yields an extremely cheap algorithm for the reconstruction of the di usion coecient. 7. Numerical results. Consider the following experiment: At time t = 0 hot metal is poured into the form sand at x = ?5 [cm] and temperature measurements are taken at m + 1 equidistant thermocouples located at xi = 20i=m [cm], i = 0; : : :; m. The temperature measurements at x0 = 0 and xm = 20 are considered as Dirichlet data for the di erential equation (1.1). The initial temperature is !0 = 20 [C]. To simulate this experiment we have chosen a phantom di usion coecient ay after consulting the collaborating engineers who also provided realistic Dirchlet data f and g for (1.1). The di erential equation was solved up to time T = 6 104 [s] with the subroutine ode15s from Matlab's ode suite, cf. Shampine and Reichelt [17]; to this end we used a semidiscrete approach with 121 piecewise linear nite elements for the discretization of the spatial domain. To study the e ect of measurement errors, the data have been perturbed by additive noise of up to 1%. A crucial part of the algorithm is the numerical di erentiation of the data u( ; xi) with respect to time. Rather than taking nite di erences we have interpolated u( ; xi) by cubic splines on relatively coarse grids, and then di erentiated the splines. For the linear system (6.1) only the equations (1.2) corresponding to i = 1; : : :; m and the time steps t = tj chosen by ode15s are considered. Still, not all of these data were actually used. For some reason the stability of our inversion scheme seems to improve when we throw away data corresponding to the initial time phase of the experiment. (We believe that the steep gradients ux in this area are badly captured by the interpolating splines; possibly, a nonequidistant distribution of the thermocouples might be more appropriate for this purpose.) As a consequence only the data within the shaded region in Figure 7.1 have been used for the reconstruction. They correspond to 22 time steps and thus lead to a linear system with 22 m equations for 1100 unknowns. We have run the algorithm with m = 5; 8, and 15, and hence, all linear systems have been strongly underdetermined. This is a problem that is taken care of by the regularization via early termination of the conjugate gradient iteration. Figure 7.2 shows the iteration history in the case of 0:1% additive noise. Note that the right-hand side plot shows the errors only in L2 , since the phantom ay itself is not H 2 , see Figure 7.3 below. While the residual is monotonically decreasing (as it should) the error history exhibits the typical semiconvergence phenomenon in illposed problems: the iterates seem to converge in the beginning before they start to 10
1000 800 600 400 200 0 0
1
2
3
4
5
6e4
Fig. 7.1. Temperature histories at nine thermocouples (m = 8) Relative Residuals
Relative Errors
1
1
m=5 m=8 m = 15
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0 0
5
10
m=5 m=8 m = 15
0.8
15
0 0
20
5
10
15
20
Fig. 7.2. Iteration history: residuals and iteration errors vs. iteration count (0:1% noise)
diverge eventually. In most cases the best approximation is obtained after just a few iterations. Reconstructions for several m and various noise levels are shown in Figure 7.3. We observe that the algorithm is able to reconstruct the phantom qualitatively correct provided that m is appropriately coupled to the noise level, . In addition to these gures Table 7.1 shows the number of iterations and the corresponding (L2 ) errors of the optimal approximations. A note on these numbers is in order: since the temperature data within the shaded region of Figure 7.1 only range between 20 and 940 [ C] it seems to be fair just to measure the errors in this interval (20; 940) rather than in the whole interval shown in Figure 7.3. The boldfaced entries in the table correspond to the optimal number of thermocouples depending on the amount of noise; the more noise the less thermocouples should be taken. Moreover, it can be seen that a large number of thermocouples (i.e., m = 15) is only useful with no or very little noise; with 1% noise the reconstruction corresponding to m = 15 is completely wrong, and is therefore omitted in Figure 7.3. It should be stressed that the numbers for the amount of noise must not be 11
no noise 0.65 0.6 0.55
0.1% noise 0.65
ay m=5 m=8 m = 15
0.6 0.55
0.5
0.5
0.45
0.45
0.4
0.4
0.35 0
200
400
600
800
0.35 0
1000
ay m=5 m=8 m = 15
200
400
0.3% noise 0.65 0.6 0.55
ay m=5 m=8 m = 15
0.6 0.55 0.5
0.45
0.45
0.4
0.4
200
400
600
800
1000
800
1000
1% noise 0.65
0.5
0.35 0
600
800
0.35 0
1000
ay m=5 m=8
200
400
600
Fig. 7.3. Reconstructions for various noise levels Table 7.1
Relative errors and iteration counts
m=5 m=8 m = 15
no noise 0.0223 6 its.
0:1% noise 0.0254 6 its.
24 its. 0.0071 38 its.
7 its. 0.0123 49 its.
0.0046
0.0118
0:3% noise
0.0347
5 its. 0.0436 3 its. 0.2662 4 its.
1% noise
0.0461
5 its. 0.2165 3 its. 0.7062 1 it.
Table 7.2
Discrepancy k y ? Aay k 2 = k y k 2
m=5 m=8 m = 15
no noise 0.1468 0.0583 0.0815
0:1% noise 0.1549 0.0759 0.1750 12
0:3% noise 0.1806 0.1853 0.4690
1% noise 0.2750 0.4891 1.5001
confused with the discrepancy of the linear system Aa = y of (6.1), i.e., the relative t ky ? Aay k 2= kyk 2 given the exact di usion coecient ay . To emphasize this point we have included those discrepancies in Table 7.2. As can be seen from this table, even with no noise added upon the temperature data the linear system (6.1) su ers from a discretization error between 5 and 15% of the actual right-hand side. We nally remark that with more than 0:3% noise the reconstructions do not match perfectly with the given phantom, but still it is possible to see the local minimum of ay and to detect the sharp bend near u = 600. On the other hand, those reconstructions are extremely cheap to compute: The whole algorithm including the setup of the linear system and seven conjugate gradient iterations takes only between 0:3 (m = 5) and 0:5 M ops (m = 15) in Matlab, with roughly 75% of this work amount being spent on solving the linear system. For comparison, the solution of the parabolic di erential equation (1.1) for the simulation of data (which is precisely the amount of work for solving the forward problem, i.e., the building block of all output least squares type methods for reconstructing the di usion coecient) took already more than 4 M ops with 121 nite elements, and still 2 M ops with 64 nite elements. We conclude this section with a comment on matrix A. As it turns out, A itself is not terribly ill-conditioned, i.e., the condition number of A remained between 400 and 700, independent of m and of the actual noise level. However, when m gets larger (in our experiments, when m = 8 or 15) then A becomes rank de cient. As we have indicated already in Sect. 3 the ill-posedness of the linear system is mostly due to the imbedding operator into H 2 , corresponding to the matrix L?1 in the discrete setting (6.2); in fact, cond (AL?1 ) 107, and hence, (6.2) is a severely ill-conditioned problem. Appendix A. We brie y recollect some of the results from the book by Ladyzenskaja, Solonnikov, and Ural'ceva [11] concerning the solution u of the forward problem (1.1). Note that a solution of (1.1) satis es a maximum principle (cf. Protter and Weinberger [15, Sect. 3.2]), i.e., the values of u over I belong to the set U de ned in (2.2). We set
?
'(t; x) := f(t) + x g(t) ? f(t) ; and transform (1.1) to
?
vt = a(v + ')vx x + a0 (v + ')(vx + 'x )'x ? 't ; v(0; x) = 0 ; v(t; 0) = 0 ; v(t; 1) = 0 ;
(7.1)
where u = v + '. From Thm. 6.1 in Chapter V of [11] it follows that (7.1) has a solution from the class C2+ of functions u with continuous derivatives ut, ux and uxx in I , which in addition have certain Holder continuity properties if one of the two arguments (t; x) is kept xed. Assume now that a 2 C 2; (U ) and that the Dirichlet data f; g 2 C 2; (I ) satisfy the compatibility condition (2.1), and f 00 (0) = g00 (0). Then, using the fact that w = vt is a solution of the di erential equation
?
wt = a(u)wx x + a0 (u)ux wx + 2a0(u)ux 'xt + a0 (u)utvxx + a00 (u)utu2x ? 'tt ; w(0; x) = 0 ; w(t; 0) = 0 ; w(t; 1) = 0 ; 13
we conclude from Thm. 5.2 in Chapter IV of [11] that vt also belongs to the class C2+ , and hence, utt is continuous in I . With similar arguments (di erentiating (7.1) with respect to x) we obtain that vx and vxx belong to C2+ provided that a 2 C 3; , and hence, @x@ 44 u is continuous in I . ? Appendix B. Here we derive the adjoint operator T : L2(I) m ! H 2 (U ) of the operator T of (2.5). To this end we consider once again the operator ? S : H 2 (U ) ?! L2 (I ) ; S : a 7! a (u() v() ; with u; v 2 L2(I ) and u has values in U almost everywhere. We also introduce the operator R : L2 (U ) ?! H 2 (U ) ; R : b 7! w ; where w 2 H 2(U ) is the weak solution of the di erential equation w ? w00 + w0000 = b with natural boundary conditions. In other words, we have hR a; biL2(U ) = ha; RbiH 2 (U ) = ha; biL2 (U ) for all a 2 H 2(U ) and all b 2 L2(U ), i.e., the adjoint R of R is the embedding operator : H 2(U ) ,! L2(U ). Note, however, that R extends in the usual way to a continuous operator R : H ?2(U ) ! H 2 (U ). Assume now for the moment that u 2 C 1(I ) and consider a time interval I1 I where u0 6= 0. Denote by U1 the interval of function values u(I1). Multiplying Sa by ' 2 L2(I ) we obtain Z Z Z (Sa)(t)'(t) dt = a(u(t))v(t)'(t) dt = a(u) v(t(u))'(t(u)) j u0(t(u)) j du ; I1 I1 U1 and we note that 8 v(t(u))'(t(u)) < u 2 U1 ; b1 (u) := : j u0(t(u)) j (7.2) 0 u 62 U1 ; belongs to L1(U ), since Z v(t(u))'(t(u)) Z Z Z (7.3) jb1(u)j du = jb1(u)j du = u0 (t(u)) du = jv(t)'(t)j dt : U1 U1 U I1 We proceed in a similar way for all further subintervals I2; I3; : : : of I in which u is monotone, thus obtaining further L1 functions b2; b3; : : :. In the remaining intervals D1 ; D2; ::: of I we have that u is a constant, i.e., ujDj = !j and
Z
Dj
(Sa)(t)'(t) dt =
Z
Dj
a(!j )v(t)'(t) dt = a(!j ) j
This yields
Z I
(Sa)(t)'(t) dt =
1 Z X i=1 Ii
(Sa)(t)'(t) dt + 14
Z
with j := v(t)'(t) dt : Dj
1Z X j =1 Dj
(Sa)(t)'(t) dt
= =
1Z X
a(u)bi (u) du +
i=1 U
Z
U
a(u)
1 X i=1
bi(u) du +
1 X j =1
j a(!j )
1 X j =1
j a(!j ) :
P b (u) converges in L1(U ), since From (7.3) then follows that the function b(u) := 1 i=1 i Z
U
jb(u)j du
1 Z X
i=1 U
jbi (u)j du =
1Z X
i=1 Ii
jv(t)'(t)j dt =
Z
I
jv(t)'(t)j dt
kvk L2 (I ) k'k L2 (I ) : P Since H s(U ) LP1 (U ) for s > 1=2 we have L1(U ) H ?2(U ). Similarly, j j j j < 1 and hence d := j j (!j ) converges as a distribution in L1(U ); here, (!j ) denotes the {distribution at !j . It therefore follows that R(b + d) 2 H 2 (U ), and hence,
hSa; 'iL2 (I ) =
Z
I
(Sa)(t)'(t) dt =
Z
U
a(u)(b(u) + d(u)) du = ha; R(b + d)iH 2 (U ) :
Furthermore, since b + d 2 L1(U ) it follows from the di erential equation for w that w = R(b+d) 2 W 4;1(U ) which embeds continuously into C 3(U), cf., e.g., [1, Thm. 5.4]. Thus we have shown that (7.4) S ' = R(b + d) 2 W 4;1(U ) ; with b and d as de ned above. The general case of u 2 L2(I ) follows by continuity of the mapping u 7! S according to Theorem 2.3. We denote by Si , i = 0; : : :; m, the operator S which is obtained when u in (7.2) is the trace u( ; xi) of the solution u of (1.1) along x = xi , and likewise v in (7.2) is the trace of the partial derivative ux of u. Then it follows that Ti = Si ? Si?1 is the operator de ned in (2.3), (2.4), and Ti = Si ? Si?1 is its adjoint. The adjoint ? 2 operator T : L (I ) m ! H 2(U ) of T in (2.5) is then given as T ('1 ; : : :; 'm ) = T1 '1 + : : : + Tm 'm ; 'i 2 L2(I ) ; and it follows from (7.4) that R(T ) W 4;1(U ). Acknowledgements. We like to thank Martin Bruhl for pointing out a gap in an earlier version of the derivation of T in Appendix B. In addition, we are grateful to an anonymous referee for bringing reference [3] to our attention. REFERENCES [1] R. A. Adams, Sobolev Spaces, Academic Press, New York, San Francisco, Melbourne, 1975. [2] H. T. Banks and K. Kunisch, Estimation Techniques for Distributed Parameter Systems, Birkhauser Verlag, Boston, Basel, Berlin, 1989. [3] J. R. Cannon and P. DuChateau, Design of an experiment for the determination of an unknown coecient in a nonlinear conduction-di usion equation, Internat. J. Engrg. Sci., 25 (1987), pp. 1067-1078. 15
[4] G. Chavent and P. Lemonnier, Identi cation de la non-linearite d'une equation parabolique quasilineaire, Appl. Math. Optim., 1 (1974), pp. 121-162. [5] P. DuChateau, An inverse problem for the hydraulic properties of porous media, SIAM J. Math. Anal., 28 (1997), pp. 611-632. [6] S. Dummel and S. Handrock-Meyer, Uniqueness of the solution of an inverse problem for a quasilinear parabolic equation in divergence form, Z. Anal. Anwendungen, 7 (1988), pp. 241-246. [7] H. W. Engl, M. Hanke, and A. Neubauer, Regularization of Inverse Problems, Kluwer, Dordrecht, 1996. [8] D. Hinestroza and D. Murio, Identi cation of transmissivity coecients by molli cation techniques. Part I: One-dimensional elliptic and parabolic problems, Comput. Math. Appl., 25 (1993), pp. 59-79. [9] T. Karkkainen, A linearization technique and error estimates for distributed parameter identi cation in quasilinear problems, Numer. Funct. Anal. Optim., 17 (1996), pp. 345364. [10] T. Karkkainen, An equation error method to recover di usion from the distributed observation, Inverse Problems, 13 (1997), pp. 1033-1051. [11] O. A. Ladyzenskaja, V. A. Solonnikov, and N. N. Ural'ceva, Linear and Quasilinear Equations of Parabolic Type, Amer. Math. Soc., Providence, RI, 1968. [12] A. K. Louis, Inverse und schlecht gestellte Probleme, Teubner, Stuttgart, 1989. [13] R. Nabokov, An inverse problem for porous medium equation, in Parameter Identi cation and Inverse Problems in Hydrology, Geology and Ecology, J. Gottlieb and P. DuChateau, eds., Kluwer, Dordrecht, Boston, London, 1996, pp. 155-163. [14] A. S. Nemirovskii, The regularizing properties of the adjoint gradient method in ill-posed problems, USSR Comput. Math. Math. Phys., 26 (1986), pp. 7-16. [15] M. H. Protter and H. F. Weinberger, Maximum Principles in Di erential Equations, Prentice-Hall, Englewood Cli s, New Jersey, 1967. [16] G. R. Richter, Numerical identi cation of a spatially varying di usion coecient, Math. Comp., 36 (1981), pp. 375-386. [17] L. F. Shampine and M. W. Reichelt, The Matlab ode suite, SIAM J. Sci. Comput., 18 (1997), pp. 1-22. [18] J. Stoer and R. Bulirsch, Introduction to Numerical Analysis, Springer-Verlag, New York, Heidelberg, Berlin, 1980. [19] B. K. Swartz and R. S. Varga, Error bounds for spline and L-spline interpolation, J. Approx. Theory, 6 (1972), pp. 6-49. [20] E. E. Tyrtyshnikov, A Brief Introduction to Numerical Analysis, Birkhauser, Boston, 1997. [21] G. M. Vainikko, Inverse problem of ground water ltration: identi ability, discretization and regularization, in Inverse Problems in Di usion Processes, H. W. Engl and W. Rundell, eds., SIAM, Philadelphia, 1995, pp. 90-107.
16
