IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY IN A ...
c 2007 Society for Industrial and Applied Mathematics
SIAM J. CONTROL OPTIM. Vol. 46, No. 2, pp. 694–713
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY IN A HEAT CONDUCTION PROCESS∗ SEMION GUTMAN† AND JUNHONG HA‡ Abstract. We study the identification and identifiability problems for heat conduction in a nonhomogeneous rod. The identifiability results are established for two different sets of observations. Given a sequence of distributed type observations, the identifiability is proved for conductivities in a piecewise smooth class of functions. In the case of observations taken at finitely many points the identifiability is established for piecewise constant conductivities. Such conductivities can be uniquely identified using the proposed marching algorithm. Key words. identification, identifiability, piecewise constant conductivity AMS subject classifications. 35R30, 93B30 DOI. 10.1137/060657364
1. Introduction. Consider the heat conduction in a nonhomogeneous insulated rod of a unit length, with the ends kept at zero temperature at all times. Our main interest is in the identification and identifiability of the discontinuous conductivity (thermal diffusivity) coefficient a(x), 0 ≤ x ≤ 1. The identification problem consists of finding a conductivity a(x) in an admissible set K for which the temperature u(x, t) fits given observations in a prescribed sense. Under a wide range of conditions one can establish the continuity of the objective function J(a) representing the best fit to the observations. Then the existence of the best fit to data conductivity follows if the admissible set K is compact in the appropriate topology. However, such an approach usually does not guarantee the uniqueness of the found conductivity a(x). Establishing such a uniqueness is referred to as the identifiability problem. If the conductivity is identifiable and one can design an algorithm for its reconstruction, then we say that a is constructively identifiable. From physical considerations the conductivity coefficients a(x) are assumed to be in (1.1)
Aad = {a ∈ L∞ (0, 1) : 0 < ν ≤ a(x) ≤ μ}.
The temperature u(a) = u(x, t; a) inside the rod satisfies (1.2)
ut − (a(x)ux )x = 0, Q = (0, 1) × (0, T ), u(0, t) = u(1, t) = 0, t ∈ (0, T ), u(x, 0) = g(x), x ∈ (0, 1),
where g ∈ L2 (0, 1). In general, the solution of (1.2) is understood in the weak sense. According to [11] for any a ∈ Aad there exists a unique weak solution u(a) ∈ L2 (0, 1; H01 (0, 1)) ∩C([0, 1]; L2 (0, 1)), and so the map a → u(a) is well defined. Moreover, this map is continuous from Aad equipped with the L2 (0, 1) topology into ∗ Received by the editors April 15, 2006; accepted for publication (in revised form) January 25, 2007; published electronically May 7, 2007. http://www.siam.org/journals/sicon/46-2/65736.html † Department of Mathematics, The University of Oklahoma, Norman, Oklahoma 73019 (sgutman@ ou.edu). ‡ School of Liberal Arts, Korea University of Technology and Education, Cheonan 330-708, South Korea ([email protected]).
694
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY
695
C([0, T ]; L2 (0, 1)). In fact, in [11] these results are established for multidimensional parabolic problems. The identification (parameter estimation) problem for (1.2) is as follows: Find a conductivity a ∈ Aad such that the solution u(a) of (1.2) fits a given observation z of the heat conduction process. For example, given z ∈ L2 (0, 1) one defines (1.3)
J(a) = u(x, T ; a) − z(x)L2 (0,1) .
Then the parameter estimation problem for (1.2) is reduced to the minimization of the objective function J over the admissible set Aad or its subset Kad : Find a ∈ Kad ⊂ Aad such that (1.4)
J(¯ a) = inf{J(a) : a ∈ Kad }.
The above-mentioned properties of the solutions u(a) imply that the objective function J(a) is continuous on Aad ∩ L2 (0, 1). Therefore the identification problem (1.4) has a solution if Aad is compact in L2 (0, 1). One such choice is Kad = {a ∈ Aad ∩ H 1 (0, 1) : aH 1 ≤ constant}; see [8]. However, a ∈ H 1 (0, 1) implies that the conductivity is continuous. Therefore this choice of Kad is not suitable for the study of the identification problems with discontinuous coefficients. To overcome this difficulty we have shown in [5] (in a multidimensional case) that one can take for Kad the set of functions in Aad which have a uniformly bounded variation. Such a set Kad is compact in L2 (0, 1), and the existence of solutions to the identification problem (1.4) follows. See [5] for additional details and numerical experiments for 2D parameter identification problems. A variety of identification problems is studied in [1] under very general assumptions on the problem’s parameters. The identifiability questions for partial differential equations are much more difficult, and there are just a few available results. Suppose that one is given an observation z(t) = u(p, t; a) of the heat conduction process (1.2) for t1 < t < t2 at some observation point 0 < p < 1. From the series solution for (1.2) and the uniqueness of the Dirichlet series expansion (see section 2), one can, in principle, recover all of the eigenvalues of the associated Sturm–Liouville problem. If one also knows the eigenvalues for the heat conduction process with the same coefficient a and different boundary conditions, then the classical results of Gelfand and Levitan [4] show that smooth coefficients a(x) can be uniquely identified from the knowledge of the two spectral sequences. Also, if the entire spectral function is known (i.e., the eigenvalues and the values of the derivatives of the normalized eigenfunctions at x = 0), then the conductivity is identifiable as well. However, such results have little practical value, since the observation data z(t) always contain some noise, and therefore one cannot hope to adequately identify more than just the few first eigenvalues of the problem. A different approach is taken in [6, 12, 13, 14]. These works show that one can identify a constant conductivity a in (1.2) from the measurement z(t) taken at one point p ∈ (0, 1). These works also discuss problems more general than (1.2), including problems with a broad range of boundary conditions, nonzero forcing functions, as well as elliptic and hyperbolic problems. In [7, 3] and references therein identifiability results are obtained for elliptic and parabolic equations with discontinuous parameters in a multidimensional setting. A typical assumption there is that one knows the normal derivative of the solution at the boundary of the region for every Dirichlet boundary input. The main result of this paper is contained in Theorem 4.6. This theorem describes and justifies the marching algorithm for the unique identification of piecewise constant conductivities from observations of (1.2) given at finitely many points pk ∈ (0, 1). We
696
SEMION GUTMAN AND JUNHONG HA
start by recalling some basic properties of (1.2) in section 2. Identifiability results for countably many distributed observations are given in section 3. Identifiability of piecewise constant conductivities a is discussed in section 4. Numerical results for the identifiability algorithms described in this paper require an extensive exposition, and they will be presented elsewhere. 2. Auxiliary results. In this section we collect some well-known results for the solutions u(x, t; a) of (1.2), as well as for its associated Sturm–Liouville problem. Since such results are scattered in the literature, some brief proof outlines are included as well. See [2, 9, 10, 11] for a detailed discussion. Definition 2.1. Function a(x) is said to belong to the class PS N if (i) a ∈ Aad = {a ∈ L∞ (0, 1) : 0 < ν ≤ a(x) ≤ μ} for some positive constants ν and μ; (ii) function a is piecewise smooth; that is, there exists a finite sequence of points 0 = x0 < x1 < · · · < xN −1 < xN = 1 such that both a(x) and a (x) are continuous on every open subinterval (xi , xi+1 ), i = 0, . . . , N − 1, and both can be continuously extended to the closed intervals [xi , xi+1 ], i = 0, . . . , N − 1. For definiteness, we assume that a and a are continuous from the right, i.e., a(x) = a(x+) and a (x) = a (x+) for all x ∈ [0, 1). Also let a(1) = a(1−). Definition 2.2. PS = ∪∞ N =1 PS N . Everywhere in the following the conductivities a are assumed to be in PS. If a ∈ PS N , then the regularity conditions on a and the uniqueness of the weak solutions imply that for any t > 0 the weak solution u(x, t; a) of (1.2) satisfies the equation in the classical sense on any subinterval (xi , xi+1 ), i = 0, . . . , N − 1. Also u satisfies the matching conditions for the continuity of the solution and its conormal derivative at xi ∈ (0, 1), i = 1, 2, . . . , N − 1: ut − (a(x)ux )x = 0, x = xi , t ∈ (0, T ), u(0, t) = u(1, t) = 0, t ∈ (0, T ), u(xi +, t) = u(xi −, t), a(xi +)ux (xi +, t) = a(xi −)ux (xi −, t), u(x, 0) = g(x), x ∈ (0, 1),
(2.1)
where g ∈ L2 (0, 1); see [11, 16]. Denote by ·, ·, · the norm and the inner product, respectively, in H = L2 (0, 1). Theorem 2.3. Let a ∈ PS. Then (i) the associated Sturm–Liouville problem (a(x)v(x) ) = −λv(x), x = xi , v(0) = v(1) = 0, v(xi +) = v(xi −), a(xi +)vx (xi +) = a(xi −)vx (xi −)
(2.2)
has infinitely many eigenvalues 0 < λ1 < λ2 < · · · → ∞. The eigenvalues {λk }∞ k=1 and the corresponding orthonormal set of eigenfunctions {vk }∞ k=1 satisfy (2.3) λk = inf
1 0
a(x)[v (x)]2 dx 1 : v ∈ H0 (0, 1), < v, vj = 0, j = 1, 2, . . . , k − 1 , 1 [v(x)]2 dx 0
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY
(2.4)
1
λk =
697
a(x)[vk (x)]2 dx.
0 2 The normalized eigenfunctions {vk }∞ k=1 form a basis in L [0, 1]. (ii) Each eigenvalue is simple. For each eigenvalue λk there exists a unique continuous, piecewise smooth normalized eigenfunction vk (x) such that vk (0+) > 0, and the function a(x)vk (x) is continuous on [0, 1]. (iii) Eigenvalues {λk }∞ k=1 satisfy the inequality
νπ 2 k 2 ≤ λk ≤ μπ 2 k 2 . (iv) The first eigenfunction v1 satisfies v1 (x) > 0 for any x ∈ (0, 1). (v) The first eigenfunction v1 has a unique point of maximum q ∈ (0, 1) : v1 (x) < v1 (q) for any x = q. (vi) For any fixed t > 0 the solution u of (2.1) is given by u(x, t; a) =
∞
g, vk e−λk t vk (x),
k=1
and the series converges uniformly and absolutely on [0, 1]. (vii) For any p ∈ (0, 1) the function z(t) = u(p, t; a),
t > 0,
is real analytic on (0, ∞). Proof. (i) The proof is standard; see, e.g., [10]. (ii) On any subinterval (xi , xi+1 ) the coefficient a(x) has a bounded continuous derivative. Therefore, on any such interval the initial value problem (a(x)v (x)) +λv = 0, v(xi ) = A, v (xi ) = B has a unique solution. Suppose that two eigenfunctions w1 (x) and w2 (x) correspond to the same eigenvalue λk . Then they both satisfy the condition w1 (0) = w2 (0) = 0. Therefore their Wronskian is equal to zero at x = 0. Consequently, the Wronskian is zero throughout the interval (x0 , x1 ), and the solutions are linearly dependent there. Thus w2 (x) = Cw1 (x) on (x0 , x1 ), w2 (x1 −) = Cw1 (x1 −), and w2 (x1 −) = Cw1 (x1 −). The linear matching conditions imply that w2 (x1 +) = Cw1 (x1 +) and w2 (x1 +) = Cw1 (x1 +). The uniqueness of solutions implies that w2 (x) = Cw1 (x) on (x1 , x2 ), etc. Thus w2 (x) = Cw1 (x) on (0, 1), and each eigenvalue λk is simple. In particular λ1 is a simple eigenvalue. The uniqueness and the matching conditions also imply that any solution of (a(x)v (x)) + λv = 0, v(0) = 0, v (0) = 0 must be identically equal to zero on the entire interval (0, 1). Thus no eigenfunction vk (x) satisfies vk (0) = 0. Assuming that the eigenfunction vk is normalized in L2 (0, 1), it leaves us with the choice of its sign for vk (0). Letting vk (0) > 0 makes the eigenfunction unique. (iii) The eigenvalues of (2.3) satisfy the min-max principle 1 2 a(x)[v (x)] dx 0 λk = min max : v ∈ Vk , 1 Vk [v(x)]2 dx 0 where Vk varies over all subspaces of H01 (0, 1) of finite dimension k; see [10]. Therefore (a) (b) a(x) ≤ b(x), x ∈ [0, 1] implies λk ≤ λk . Since the eigenvalues of (2.3) with a(x) = 1 2 2 are π k , the required inequality follows.
698
SEMION GUTMAN AND JUNHONG HA
(iv) Recall that v1 (x) is a continuous function on [0, 1]. Suppose that there exists p ∈ (0, 1) such that v1 (p) = 0. Let wl (x) = v1 (x) for 0 ≤ x < p and wl (x) = 0 for p ≤ x ≤ 1. Let wr (x) = v1 (x) − wl (x), x ∈ [0, 1]. Then wl , wr are continuous, and, moreover, wl , wr ∈ H01 (0, 1). Also
1
wl (x)wr (x)dx = 0
1
and
0
a(x)wl (x)wr (x)dx = 0.
0
Suppose that wl is not an eigenfunction for λ1 . Then
1
a(x)[wl (x)]2 dx
1
[wl (x)]2 dx.
> λ1
0
0
Since
1
a(x)[wr (x)]2 dx ≥ λ1
0
1
[wr (x)]2 dx, 0
we have λ1 =
1
1
0
0
a(x)[v1 (x)]2 dx = 1 2 dx [v (x)] 1 0
a(x)([wl (x)]2 + [wr (x)]2 )dx 1 ([wl (x)]2 + [wr (x)]2 )dx 0
1 >
0
(λ1 [wl (x)]2 + λ1 [wr (x)]2 )dx = λ1 . 1 ([wl (x)]2 + [wr (x)]2 )dx 0
This contradiction implies that wl (and wr ) must be an eigenfunction for λ1 . However, wl (x) = 0 for p ≤ x ≤ 1, and as in (ii) it implies that wl (x) = 0 for all x ∈ [0, 1], which is impossible. Since v1 (0) > 0 the conclusion is that v1 (x) > 0 for x ∈ (0, 1). (v) From part (i), any eigenfunction vk is continuous and satisfies (a(x)vk (x)) = −λk vk (x) for x = xi . Also the function a(x)vk (x) is continuous on [0, 1] because of the matching conditions at the points of discontinuity xi , i = 1, 2, . . . , N − 1 of a. The integration gives x a(x)vk (x) = a(p)vk (p) − λk vk (s)ds p
for any x, p ∈ (0, 1). Let p ∈ (0, 1) be a point of maximum of vk . If p = xi , then vk (p) = 0. If p = xi , then vk (xi −) ≥ 0 and vk (xi +) ≤ 0. Therefore limx→p a(x)vk (x) = 0 and vk (p+) = vk (p−) = 0 since a(x) ≥ ν > 0. In any case for such a point p we have (2.5)
a(x)vk (x)
= −λk
x
vk (s)ds,
x ∈ (0, 1).
p
Since v1 (x) > 0, a(x) > 0 on (0, 1), (2.5) implies that v1 (x) > 0 for any 0 ≤ x < p and v1 (x) < 0 for any p < x ≤ 1. Since the derivative of v1 is zero at any point of maximum, we have to conclude that such a maximum p is unique.
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY
699
(vi) We prove only the convergence part. Note that 1 2 νvk ≤ a(x)[vk (x)]2 dx = λk vk 2 = λk . 0
Thus
√ λk vk ≤ √ ν
and
x
|vk (x)| ≤
|vk (s)|ds
≤
vk
0
√ λk ≤ √ . ν
Bessel’s inequality implies that the sequence of Fourier coefficients g, vk is bounded. Therefore, denoting by C various constants and using the fact that the function √ s → se−σs is bounded on [0, ∞) for any σ > 0, one gets √ λk t λ k λk t λk t −λk t | g, vk e vk (x)| ≤ C √ e− 2 e− 2 ≤ Ce− 2 . ν From (iii) of this theorem λk ≥ νπ 2 k 2 . Thus ∞
| g, vk e−λk t vk (x)| ≤ C
k=1
∞
e−
νπ 2 k2 t 2
k=1
≤C
∞
e−
νπ 2 t 2
k
< ∞.
k=1
∞ −λk t0 From (vi), the series vk (p) (vii) Let t0 > 0 and p ∈ (0, k=1 g, vk e 1). ∞ −λk s converges absolutely. Therefore k=1 g, vk e vk (p) is analytic in the part of the complex plane {s ∈ C : Re s > t0 }, and the result follows. ∞ Series of the form k=1 Ck e−λk t are known as the Dirichlet series. The following lemma shows that the Dirichlet series representation of a function is unique. Lemma 2.4. Let μk > 0, k = 1, 2, . . . , be a strictly increasing sequence. Suppose ∞ that T1 ≥ 0 and k=1 |Ck | < ∞. If ∞
Ck e−μk t = 0
for all
t ∈ (T1 , T2 ),
k=1
then Ck = 0 for k = 1, 2, . . . . ∞ −μk z The result follows at once from the observation that the series k=1 Ck e converges uniformly in the Re z > 0 region of the complex plane, implying that it is an analytic function there. See Chapter 9 of [15] for additional results on Dirichlet series. Remark. According to Theorem 2.3(vi) for each fixed p ∈ (0, 1) the solution z(t) = u(p, t; a) of (2.1) is given by a Dirichlet series. However, Lemma 2.4 is not directly applicable since the coefficients Ck = g, vk vk (p) are only square summable. Nevertheless, the conclusion of Lemma 2.4 remains valid, since the exponents μk in the Dirichlet series are the eigenvalues λk which satisfy the growth condition stated in Theorem 2.3(iii). This allows one to conclude (Theorem 2.3(vii)) that the solution z(t) is a real analytic function on (0, ∞), and the uniqueness of such a representation follows. Thus it would be a mistake to simply refer to the standard results such as Lemma 2.4 for the uniqueness of the Dirichlet series representation to justify the paper’s conclusions.
700
SEMION GUTMAN AND JUNHONG HA
3. Identifiability by distributed measurements. Suppose that one is given some observations of the heat conduction process (2.1) with an unknown conductivity a(x) and that they coincide with the observations of the model process
(3.1)
m m m um t ∈ (0, T ), t − (a (x)ux )x = 0, x = xi , um (0, t) = um (1, t) = 0, t ∈ (0, T ), m m um (xm i +, t) = u (xi −, t), m m m m m m am (xm +)u (x x i i +, t) = a (xi −)ux (xi −, t), um (x, 0) = g(x), x ∈ (0, 1),
where g is the same as in (2.1). The conductivity a is said to be identifiable in some class of functions M if the coincidence of the measurements of the observed and the model processes implies that a = am , provided a, am ∈ M. 2 Theorem 3.1. Let {ψn }∞ n=1 be a complete orthonormal set in H = L (0, 1). Suppose that nonzero initial data g ∈ H and the observations zn (t) = u(x, t; a), ψn for n = 1, 2, . . . and 0 ≤ T1 < t < T2 of the heat conduction process (2.1) are given. Then the conductivity a(x) ∈ Aad is constructively identifiable in the class of piecewise smooth functions PS. Proof. To show the identifiability of a we give an algorithm for its reconstruction from the data zn (t), n = 1, 2, . . . , guaranteeing the uniqueness in each step. Using Theorem 2.3(vi) we have (3.2)
zn (t) =
∞
g, vk e−λk t vk , ψn
k=1
for each n = 1, 2, . . . and 0 ≤ T1 < t < T2 . Fix an n > 0. Since ψn ∈ H and {vk }∞ k=1 form a basis in H, the Bessel inequality implies that the sequence of the Fourier coefficients { vk , ψn }∞ k=1 is bounded. From Theorem 2.3(vi) one concludes that the above series converges absolutely. Note that some products g, vk vk , ψn may be equal to zero. The zero value products present a difficulty, since we would not know how to associate the sequence of exponents recovered from (3.2) with the (unknown) eigenvalues λk : Some eigenvalues may be missing from the sequence. Define (possibly empty) subsets Qn ⊂ N by Qn = {k ∈ N : g, vk vk , ψn = 0},
n = 1, 2, . . . .
For Qn = ∅ reindex (3.2) so that there would be no vanishing coefficients: (3.3)
zn (t) =
∞
Cl,n e−μl,n t ,
t ∈ (T1 , T2 ).
l=1
By Theorem 2.3(vii) the solutions zn (t) are real analytic. Therefore, since all of the coefficients Cl,n = 0, one can uniquely determine them and the sequences ∞ μl,n , l = 1, 2, . . . . Recall that {μl,n }∞ l=1 ⊂ {λk }k=1 for any n with a nonempty Qn so each μl,n ≥ λ1 > 0. For each n such that Qn = ∅ let γn = min{μl,n : l ∈ N}. Define γ = min{γn : Qn = ∅}
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY
and (3.4)
701
⎧ ⎨ C1,n if γn = γ, Qn = ∅, 0 if γn > γ, Qn = ∅, An = ⎩ 0 if Qn = ∅.
Let (3.5)
w(x) =
∞
An ψn (x).
n=1
We claim that w is a nonzero multiple of some eigenfunction vJ of (2.1). Indeed, let J be the smallest index for which g, vJ = 0. Such an index exists since g = 0, and the eigenfunctions form a basis in H. Now, since vJ = 0 and {ψn }∞ n=1 also form a basis in H, there exists n ∈ N such that vJ , ψn = 0. Thus g, vJ vJ , ψn = 0 and γ ≤ λJ . The choice of J implies that γ = λJ . By the definition An = g, vJ vJ , ψn for nonzero products. Therefore w(x) =
∞
g, vJ vJ , ψn ψn (x)
n=1
= g, vJ
∞
vJ , ψn ψn (x) = g, vJ vJ (x)
n=1
as claimed. Now we show that the set of points y ∈ (0, 1) where w (y) = 0 is finite. Assuming the opposite, and since w is a nonzero multiple of vJ , there exists a sequence yj ∈ (0, 1) such that vJ (yj ) = 0 and limj→∞ yj = c ∈ [0, 1]. The continuity of a(x)vJ (x) implies that vJ (c) = 0. From (a(x)vJ (x)) = −γvJ (x) one gets yj+1 0 = a(yj+1 )vJ (yj+1 ) − a(yj )vJ (yj ) = −γ (3.6) vJ (s)ds yj
and concludes that vJ cannot be strictly positive or strictly negative on (yj , yj+1 ). Let ζj ∈ (yj , yj+1 ) be such that vJ (ζj ) = 0. Then limj→∞ ζj = c ∈ [0, 1] and vJ (c) = 0. Now we have both vJ (c) = 0 and vJ (c) = 0. But then the uniqueness of the Cauchy problem for the second order linear equations, and the matching conditions (see the proof of Theorem 2.3(ii)) imply that vJ (x) = 0 for all x ∈ [0, 1], which is impossible. Let q be a point of maximum of w. Note that w may have several maxima, unless it is a multiple of v1 . In any case, (2.5) implies x (3.7) w(s)ds, x ∈ (0, 1). a(x)w (x) = −γ q
Then, outside of the finite sets {xi } and {yj }, the conductivity a(x) is uniquely determined from (3.7) by x w(s)ds q . a(x) = −γ w (x) Because a is assumed to be in PS, it can be uniquely extended to the entire interval [0, 1].
702
SEMION GUTMAN AND JUNHONG HA
√ Remark 1. In an application one can choose ψn = 2 sin πnx and the initial condition g(x) > 0 on (0, 1). Then g, v1 v1 , ψ1 = 0, since v1 (x) > 0 on (0, 1) by Theorem 2.3(iv). Thus J = 1 in this case, and w(x) = g, v1 v1 (x) in the above algorithm. Also, one can see from (3.6) that there is only one point y1 = q ∈ (0, 1) where v1 (q) = 0, and it is the point of maximum of v1 (x) on (0, 1). Indeed, if there were two such points, then by (3.6) v1 (x) would have to become negative between them, which would contradict v1 (x) > 0 on (0, 1). Remark 2. Since the system {ψn }∞ n=1 is complete, the conditions of Theorem 3.1 imply that for any t > 0 one knows u(x, t; a) almost everywhere on [0, 1]. Since u(x, t; a) is continuous in x and analytic in t, Theorem 3.1, in fact, assumes that the solution is known in [0, 1] × (0, T ) or (equivalently) in [0, 1] × (0, ∞). Thus, Theorem 3.1 can be stated under any of these conditions. However, a practical reconstruction of the conductivity a directly from the equation ut = (aux )x is extremely unstable. The algorithm presented in the above theorem does not reconstruct the entire solution u but just the first eigenfunction of the associated Sturm–Liouville problem. Its stability properties will be studied elsewhere. 4. Identifiability of piecewise constant conductivities from finitely many observations. The identifiability is sought within the following set. Definition 4.1. Let PC ⊂ PS be the class of piecewise constant conductivities, and let PC N = PC ∩ PS N . Functions a ∈ PC N have the form a(x) = ai for x ∈ [xi−1 , xi ), i = 1, 2, . . . , N . In this case the governing system (2.1) is
(4.1)
ut − ai uxx = 0, x ∈ (xi−1 , xi ), u(0, t) = u(1, t) = 0, t ∈ (0, T ), u(xi +, t) = u(xi −, t), ai+1 ux (xi +, t) = ai ux (xi −, t), u(x, 0) = g(x), x ∈ (0, 1),
t ∈ (0, T ),
where g ∈ L2 (0, 1) and i = 1, 2, . . . , N − 1. The associated Sturm–Liouville problem is
(4.2)
ai v (x) = −λv(x), x ∈ (xi−1 , xi ), v(0) = v(1) = 0, v(xi +) = v(xi −), ai+1 v (xi +) = ai v (xi −)
for i = 1, 2, . . . , N − 1. We are interested only in the first eigenfunction v1 of (4.2). Let λ1 be the first eigenvalue. Suppose that p∗ ∈ (xi−1 , xi ). Then v1 can be expressed on (xi−1 , xi ) as
π λ1 π ∗ (x − p ) + γ , A > 0, − < γ < . v1 (x) = A cos ai 2 2 The range for γ in the above representation follows from the fact that v1 (p∗ ) = A cos γ > 0 by Theorem 2.3(iv). The identifiability of piecewise constant conductivities is based on the following three lemmas. Lemma 4.2. Suppose that δ > 0. Assume Q1 , Q3 ≥ 0, Q2 > 0, and 0 < Q1 + Q3 < 2Q2 . Let π π π Γ = (A, ω, γ) : A > 0, 0 < ω < , − 0 and cos γ > 0. Therefore cos(ωδ − γ) sin γ Q1 = cos(ωδ) + sin(ωδ) = , cos γ cos γ Q2 Q3 cos(ωδ + γ) sin γ = . = cos(ωδ) − sin(ωδ) cos γ cos γ Q2
(4.3) (4.4)
Adding (4.3) and (4.4) yields cos ωδ =
Q1 + Q3 . 2Q2
Since 0 < ωδ < π2 and 0 < (Q1 +Q3 )/2Q2 < 1 the above equation is uniquely solvable. Now subtracting (4.4) from (4.3) yields tan γ =
Q1 − Q3 , 2Q2 sin ωδ
which is also uniquely solvable, since −π/2 < γ < π/2. Finally, we have A = Q2 / cos γ. Lemma 4.3. Suppose that δ > 0, 0 < p ≤ x1 < p + δ < 1, 0 < ω1 , ω2 < π/2δ. Let w(x), v(x), x ∈ [p, p + δ] be such that w(x) = A1 cos ω1 x + B1 sin ω1 x, v(x) = A2 cos ω2 x + B2 sin ω2 x. Suppose that v(x1 ) = w(x1 ), v (x1 ) > 0,
ω12 v (x1 ) = ω22 w (x1 ),
v(x1 ) > 0.
Then
(i) conditions v(p + δ) = w(p + δ), v (p + δ) ≥ 0, and ω1 ≤ ω2 imply ω1 = ω2 ; (ii) conditions v(p + δ) = w(p + δ), w (p + δ) ≥ 0, and ω1 ≥ ω2 imply ω1 = ω2 . Proof. Since v(x1 ) > 0, v (x1 ) > 0 we have v(x) = A sin[ω2 (x − x1 ) + γ],
0 0. For ω2 > ω1 and v (p + δ) ≥ 0 one has cos[ω2 (p + δ − x1 ) + γ] =
1 v (p + δ) ≥ 0. ω2
Therefore ω2 + ω1 π (p + δ − x1 ) + γ < ω2 (p + δ − x1 ) + γ ≤ 2 2 and
cos
ω2 + ω 1 (p + δ − x1 ) + γ > cos[ω2 (p + δ − x1 ) + γ] ≥ 0. 2
Thus v(p + δ) − w(p + δ) > 0, and the conclusion (i) of the lemma follows. The case ω2 < ω1 and w (p + δ) ≥ 0 is reduced to the already established one by interchanging ω1 with ω2 and w with v. Lemma 4.4. Let δ > 0, 0 < η ≤ 2δ, ω1 = ω2 , with 0 < ω1 δ, ω2 δ < π/2. Also let A, B > 0, 0 ≤ p < p + η ≤ 1, and w(x) = A cos[ω1 (x − p) + γ1 ], v(x) = B cos[ω2 (x − p − η) + γ2 ], with |γ1 |, |γ2 | < π/2. Then the system (4.5) (4.6)
w(q) = v(q), ω22 w (q) = ω12 v (q),
(4.7)
w(q) > 0,
v(q) > 0
admits at most one solution q on [p, p + η]. This unique solution q can be computed as follows: If γ1 ≥ 0, then
B 2 − A2 1 − γ1 . (4.8) arctan ω1 2 2 q =p+ ω1 A ω2 − B 2 ω12 If γ2 ≤ 0, then (4.9)
B 2 − A2 1 − γ2 . − arctan ω2 2 2 q =p+η+ ω2 A ω2 − B 2 ω12
705
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY
Otherwise, compute q1 and q2 according to (4.8) and (4.9) and discard the one that does not satisfy the conditions of the lemma. Proof. Let α > 0 and cos t c(t; α) = , t ∈ R. α sin t Vector function c(t, α) traverses the ellipse E(1, α) centered in the origin with the x semiaxis equal to 1 and the y semiaxis equal to α. This function can be rewritten as c(t; α) = P(α)M(t)e1 , where M(t) =
cos t sin t
− sin t cos t
,
P(α) =
1 0
0 α
,
e1 =
1 0
.
Note that M(t) is the counterclockwise rotation in R2 by the angle t, P(α) is the α-contraction (expansion) in R2 along the y axis, and e1 is the standard basis vector along the x axis. Furthermore (4.10)
c(t1 + t2 ; α) = P(α)M(t1 )M(t2 )e1 .
With this notation system (4.5)–(4.6) is Ac(ω1 (q − p) + γ1 ; 1/ω1 ) = Bc(ω2 (q − p − η) + γ2 ; 1/ω2 ) or (4.11)
P(ω1−1 )M[ω1 (q − p) + γ1 ]Ae1 = P(ω2−1 )M[ω2 (q − p − η) + γ2 ]Be1 .
If q is a solution of (4.5)–(4.6), then the vectors in the right and the left sides of (4.11) are identical. Thus they belong to the intersection of the ellipses E(A, A/ω1 ) and E(B, B/ω2 ), and this intersection is not empty. In general the ellipses intersect in four points: one in each quadrant. Suppose that q ∗ = q is another solution of (4.5)–(4.6) on [p, p+η]. We can assume that q ∗ = q + τ for some τ > 0, 0 < ω1 τ, ω2 τ < π. System (4.5)–(4.6) at x = q ∗ is (4.12) P(ω1−1 )M[ω1 (q + τ − p) + γ1 ]Ae1 = P(ω2−1 )M[ω2 (q + τ − p − η) + γ2 ]Be1 . Using (4.10) and P(α)P(α−1 ) = I the right side of (4.12) can be written as P(ω2−1 )M[ω2 (q + τ − p − η) + γ2 ]Be1 = P(ω2−1 )M[ω2 τ ]M[ω2 (q − p − η) + γ2 ]Be1 = P(ω2−1 )M[ω2 τ ]P(ω2 )P(ω2−1 )M[ω2 (q − p − η) + γ2 ]Be1 = P(ω2−1 )M[ω2 τ ]P(ω2 )P(ω1−1 )M[ω1 (q − p) + γ1 ]Ae1 . Similarly the left side of (4.12) can be written as P(ω1−1 )M[ω1 (q + τ − p) + γ1 ]Ae1 = P(ω1−1 )M[ω1 τ ]M[ω1 (q − p) + γ1 ]Ae1 = P(ω1−1 )M[ω1 τ ]P(ω1 )P(ω1−1 )M[ω1 (q − p) + γ1 ]Ae1 .
706
SEMION GUTMAN AND JUNHONG HA
Let v = P(ω1−1 )M[ω1 (q − p) + γ1 ]Ae1 and D = P(ω2−1 )M[ω2 τ ]P(ω2 ) − P(ω1−1 )M[ω1 τ ]P(ω1 ). Then (4.12) is Dv = 0. Since v = 0 we must have det(D) = 0. Note that 1 det(D) = 2ω1 ω2 − (ω12 + ω22 ) sin ω1 τ sin ω2 τ − 2ω1 ω2 cos ω1 τ cos ω2 τ ω ω 1 2 1 1 1 2 2 2ω1 ω2 − (ω1 + ω2 ) cos(ω1 − ω2 )τ + (ω1 − ω2 ) cos(ω1 + ω2 )τ . = ω1 ω2 2 2 Let us define f (τ ) on [0, π/ω1 ) ∩ [0, π/ω2 ) as 1 1 f (τ ) = 2ω1 ω2 + (ω1 − ω2 )2 cos(ω1 + ω2 )τ − (ω1 + ω2 )2 cos(ω1 − ω2 )τ. 2 2 Function f is smooth on (0, π/ω1 ) ∩ (0, π/ω2 ), and its first and second derivatives are 1 1 f (τ ) = − (ω1 − ω2 )2 (ω1 + ω2 ) sin(ω1 + ω2 )τ + (ω1 + ω2 )2 (ω1 − ω2 ) sin(ω1 − ω2 )τ, 2 2 f (τ ) = (ω1 − ω2 )2 (ω1 + ω2 )2 sin ω1 τ sin ω2 τ. Since f (0) = f (0) = 0, f (τ ) > 0, and f (τ ) > 0 on (0, π/ω1 )∩(0, π/ω2 ), we conclude that f (τ ) > 0 for all τ ∈ (0, π/ω1 ) ∩ (0, π/ω2 ). Thus det(D) = 0 if and only if τ = 0. This contradicts the assumption τ > 0. Therefore the solution q of (4.5)–(4.7) is unique on [p, p + η]. To obtain formulas (4.8) and (4.9) notice that the ellipses E(A, A/ω1 ) and E(B, B/ω2 ) are given by x2 + ω12 y 2 = A2
and x2 + ω22 y 2 = B 2 .
At the intersection points we have y2 =
B 2 − A2 ω22 − ω12
and x2 =
A2 ω22 − B 2 ω12 . ω22 − ω12
The polar angle of the intersection point in the first quadrant is B 2 − A2 . ζ = arctan 2 2 A ω − B 2 ω2 2
1
Since w(q) = v(q) > 0 the intersection points corresponding to the solution q are in either the first or the fourth quadrants, and 0 ≤ ζ < π/2. If w (p) ≤ 0, then γ1 ≥ 0. Therefore 0 ≤ γ1 ≤ ω1 (q − p) + γ1 . In this case the intersection point is in the first quadrant. Accordingly tan[ω1 (q − p) + γ1 ] = ω1 tan ζ. Thus
B 2 − A2 1 − γ1 . arctan ω1 2 2 q =p+ ω1 A ω2 − B 2 ω12 If v (p+η) ≥ 0, then γ2 ≤ 0 and ω2 (q−p−η)+γ2 ≤ γ2 ≤ 0. In this case the intersection point is in the fourth quadrant. Accordingly tan[ω2 (q − p − η) + γ2 ] = −ω2 tan ζ and one gets (4.9).
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY
707
Now we would like to define a class of piecewise constant conductivities with sufficiently separated points of discontinuity. Definition 4.5. By the definition of a ∈ PC there exists N ∈ N and a finite sequence 0 = x0 < x1 < · · · < xN −1 < xN = 1 such that a is a constant on each subinterval (xn−1 , xn ), n = 1, . . . , N . Let σ > 0. Define PC(σ) = {a ∈ PC : xn − xn−1 ≥ σ,
n = 1, 2, . . . , N }.
Note that a ∈ PC(σ) attains at most N = [[1/σ]] distinct values ai , 0 < ν ≤ ai ≤ μ. The following theorem is our main result. It describes and justifies the marching algorithm for the unique identification of piecewise constant conductivities in the class PC(σ). Theorem 4.6. Given σ > 0 let an integer M be such that 3 μ and M > 2 M≥ . σ ν Suppose that the initial data g(x) > 0, 0 < x < 1, and the observations zm (t) = u(pm , t; a), pm = m/M for m = 1, 2, . . . , M − 1 and 0 ≤ T1 < t < T2 of the heat conduction process (4.1) are given. Then the conductivity a ∈ Aad is constructively identifiable in the class of piecewise constant functions PC(σ). First, we present the marching algorithm for the unique identification of the conductivity a and then justify it. The algorithm marches from the left end x = 0 to a certain observation point pl−1 ∈ (0, 1) and identifies the values an and the discontinuity points xn of the conductivity a on [0, pl−1 ]. Then the algorithm marches from the right end point x = 1 to the left until it reaches the observation point pl+1 ∈ (0, 1) identifying the values and the discontinuity points of a on [pl+1 , 1]. Finally, the values of a and its discontinuity are identified on the interval [pl−1 , pl+1 ]. The overall goal of the algorithm is to determine the number N − 1 of the discontinuities of a on [0, 1], the discontinuity points xn , n = 1, 2, . . . , N − 1, and the values an of a on [xn−1 , xn ], n = 1, 2, . . . , N (x0 = 0, xN = 1). As a part of the process the algorithm determines certain functions Hn (x) defined on intervals [xn−1 , xn ], n = 1, 2, . . . , N . The resulting function H(x) defined on [0, 1] is a multiple of the first eigenfunction v1 . Marching Algorithm. (i) Represent the data zm (t) as (4.13)
zm (t) =
∞
ck,m e−λk t ,
m = 1, 2, . . . , M − 1,
0 ≤ T 1 < t < T2 ,
k=1
and use it to uniquely identify the first eigenvalue λ1 and the coefficients Gm = c1,m , m = 1, 2, . . . , M − 1. Let G0 = GM = 0. (ii) Find l, 0 < l < M , such that Gl = max{Gm : m = 1, 2, . . . , M − 1} and Gm < Gl for any 0 ≤ m < l. (iii) Let i = 1, m = 0. (iv) Use Lemma 4.2 to find Ai , ωi , and γi from the system ⎧ ⎨ Ai cos(ωi δ − γi ) = Gm , Ai cos γi = Gm+1 , (4.14) ⎩ Ai cos(ωi δ + γi ) = Gm+2 . Let Hi (x) = Ai cos(ωi (x − pm+1 ) + γi ).
708
SEMION GUTMAN AND JUNHONG HA
(v) If m + 3 ≥ l, then go to step (viii). If Hi (pm+3 ) = Gm+3 , or Hi (pm+3 ) = Gm+3 and Hi (pm+3 ) ≤ 0, then a has a discontinuity xi on interval [pm+2 , pm+3 ). Proceed to the next step (vi). If Hi (pm+3 ) = Gm+3 and Hi (pm+3 ) > 0, then let m := m + 1 and repeat this step (v). (vi) Use Lemma 4.2 to find Ai+1 , ωi+1 , and γi+1 from the system ⎧ ⎨ Ai+1 cos(ωi+1 δ − γi+1 ) = Gm+3 , Ai+1 cos γi+1 = Gm+4 , (4.15) ⎩ Ai+1 cos(ωi+1 δ + γi+1 ) = Gm+5 . Let Hi+1 (x) = Ai+1 cos(ωi+1 (x − pm+4 ) + γi+1 ). (vii) Use the formulas in Lemma 4.4 to find the unique discontinuity point xi ∈ [pm+2 , pm+3 ). The parameters and functions used in Lemma 4.4 are defined as follows. Let p = pm+2 , η = δ. To avoid confusion we are going to use the notation Ω1 , Ω2 , Γ1 , Γ2 for the corresponding parameters ω1 , ω2 , γ1 , γ2 required in Lemma 4.4. Let Ω1 = ωi , Ω2 = ωi+1 . For w(x) use function Hi (x) recentered at p = pm+2 ; i.e., rewrite Hi (x) in the form w(x) = Hi (x) = A cos(Ω1 (x − pm+2 ) + Γ1 ),
|Γ1 | < π/2.
For v(x) use function Hi+1 recentered at p + η = pm+3 ; i.e., v(x) = Hi+1 (x) = B cos(Ω2 (x − pm+3 ) + Γ2 ),
|Γ2 | < π/2.
Let i := i + 1, m := m + 3. If m < l, then return to step (v). If m ≥ l, then go to the next step (viii). (viii) Do steps (iii)–(vii) in the reverse direction of x, advancing from x = 1 to x = pl+1 . Identify the values and the discontinuity points of a on [pl+1 , 1], and determine the corresponding functions Hi (x). (ix) Using the notation introduced in (vii) let Hj (x) be the previously determined function H on interval [pl−2 , pl−1 ]. Recenter it at p = pl−1 ; i.e., w(x) = Hj (x) = A cos(Ω1 (x − pl−1 ) + Γ1 ). Let Hj+1 (x) be the previously determined function H on interval [pl+1 , pl+2 ]. Recenter it at pl+1 : v(x) = Hj+1 (x) = B cos(Ω2 (x−pl+1 )+Γ2 ). If Ω1 = Ω2 , then stop; otherwise, use Lemma 4.4 with η = 2δ and the above parameters to find the discontinuity xj ∈ [pl−1 , pl+1 ]. Stop. Proof. To prove Theorem 4.6 we need to justify the marching algorithm and to show the uniqueness of the identification in each step. (i) Using Theorem 2.3(vi) we get (4.16) zm (t) =
∞
gk e−λk t vk (pm ),
m = 1, 2, . . . , M − 1,
0 ≤ T1 < t < T2 ,
k=1
where gk = g, vk for k = 1, 2, . . . . By Theorem 2.3(iv) v1 (x) > 0 on interval (0, 1). Since g is positive on (0, 1) we conclude that g1 v1 (pm ) > 0. According to Theorem 2.3(vii) each observation zm (t) is a real analytic function. Thus one can uniquely determine the nonzero coefficients in (4.16) and the corresponding exponents. In particular, one determines the first eigenvalue λ1 and the values of (4.17)
Gm = g1 v1 (pm ) > 0,
pm = m/M, m = 1, 2, . . . , M − 1.
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY
709
Because of the zero boundary conditions we can let G0 = GM = 0. The crucial point M −1 is that the numbers {Gm }m=1 are not arbitrary but are the values (up to a nonzero multiplicative constant g1 ) of the still-undetermined eigenfunction v1 . (ii) Let index l be defined as in (ii) of the marching algorithm. By Theorem 2.3(v) there exists a unique point q ∗ of maximum of v1 on (0, 1). Note that q ∗ ∈ (pl−1 , pl+1 ). Thus Gl+1 ≤ Gl and Gm < Gl for m > l + 1. Also v1 (pm −) > 0 for m = 1, 2, . . . , l − 1 and v1 (pm −) < 0 for m = l + 1, l + 2, . . . , M − 1. (iii) Start at the left end point p0 = 0 and work on interval [0, x1 ], where x1 is the first discontinuity point of a. (iv) Let δ = 1/M . Since σ ≥ 3δ and a ∈ PC(σ) we conclude that [0, p2 ) ⊂ [0, x1 ) and a = a1 on [0, x1 ). To apply Lemma 4.2 we just need to check the conditions for Q1 , Q2 , Q3 required there. We have Q1 = G0 = 0, Q2 = G1 = g1 v1 (p1 ) > 0 , Q3 = G2 = g1 v1 (p2 ) > 0. Let λ1 ω1 = . a1 By Theorem 2.3(iii) 0 < λ1 ≤ μπ 2 . Since 0 < ν ≤ a1 we have μπ 2 1 ν π 0 < ω1 δ < = . ν 2 μ 2 This inequality and v1 (x) > 0 on (0, 1) imply that the first eigenfunction v1 of (4.2) can be represented on (0, x1 ) as v1 (x) = C1 cos(ω1 (x − p1 ) + γ1 ) for some (C1 , ω1 , γ1 ) ∈ Γ, where Γ was defined in Lemma 4.2. Also Q1 + Q3 = g1 C1 (cos(ω1 δ + γ1 ) + cos(ω1 δ − γ1 )) = 2g1 C1 cos(ω1 δ) cos γ1 < 2G1 = 2Q2 since 0 < ω1 δ < π/2; hence 0 < cos(ω1 δ) < 1. Now Lemma 4.2 guarantees a unique solution of the system ⎧ ⎨ g1 C1 cos(ω1 δ − γ1 ) = G0 , g1 C1 cos γ1 = G1 , (4.18) ⎩ g1 C1 cos(ω1 δ + γ1 ) = G2 . It also gives formulas for the computation of A1 = g1 C1 , γ1 , ω1 from the known values of G0 , G1 , and G2 . Thus one can determine a1 = λ1 /ω12 and obtain (4.19)
H1 (x) = g1 C1 cos(ω1 (x − p1 ) + γ1 ) = g1 v1 (x)
for x ∈ [0, x1 ). (v) Let ai be the value of a on the part of the interval [pm+1 , pm+2 ) adjacent to [pm+2 , pm+3 ). By construction this value and the associated function Hi (x) = g1 v1 (x) are already determined by the algorithm. If there is no discontinuity of a on interval [pm+2 , pm+3 ), then a has the same value ai on interval [pm+2 , pm+3 ) as well. Therefore Hi (x) = g1 v1 (x) on this interval, and we must have Gm+3 = Hi (pm+3 ) by (4.17). If one has Gm+3 = Hi (pm+3 ), then the implication is that there is a discontinuity of a on [pm+2 , pm+3 ), and one proceeds to step (vi). On the other hand, if Gm+3 = Hi (pm+3 ), then one cannot, in general, conclude that there is no discontinuity of a on [pm+2 , pm+3 ). However, since we have m + 3 < l
710
SEMION GUTMAN AND JUNHONG HA
then (ii) of the proof implies that v1 (pm+3 −) > 0. Then the assumption a = ai on [pm+2 , pm+3 ) implies Hi (pm+3 ) = g1 v1 (pm+3 −). Therefore the equality Gm+3 = Hi (pm+3 ) together with Hi (pm+3 ) ≤ 0 lead to a contradiction. The conclusion is that Gm+3 = Hi (pm+3 ) and Hi (pm+3 ) ≤ 0 imply a discontinuity of a on [pm+2 , pm+3 ), and one proceeds to step (vi). Finally, one uses Lemma 4.3 to conclude that m + 3 < l, Gm+3 = Hi (pm+3 ), and Hi (pm+3 ) > 0 imply that there is no discontinuity of a on [pm+2 , pm+3 ). Indeed, suppose that there is a discontinuity point xi of a on interval [pm+2 , pm+3 ). Then a = ai on [pm+1 , xi ) and a = ai+1 on [xi , pm+3 ]. We are going to use the notation ω1 , and ω2 used in Lemma 4.3. Let xi , Ω1 , Ω2 for the corresponding variables x1 , p = pm+2 , p + δ = pm+3 , Ω1 = λ1 /ai , Ω2 = λ1 /ai+1 , and w(x) = Hi (x) = g1 v1 (x) = A1 cos Ω1 x + B1 sin Ω1 x, x ∈ [p, xi ), v(x) = g1 v1 (x) = A2 cos Ω2 x + B2 sin Ω2 x, x ∈ [xi , p + δ]. Note that the condition Ω21 v (xi ) = Ω22 w (xi ) is just the matching condition (4.2) at x = xi . Since m + 3 < l, the maximum q ∗ of v1 satisfies q ∗ > pm+3 . Because w is a positive multiple of v1 , it implies w(xi ) > 0 and w (xi ) > 0. Therefore v(xi ) > 0 and v (xi ) > 0. Because v is a positive multiple of v1 , we have v (p+δ) > 0. The condition v(p + δ) = w(p + δ) means v(pm+3 ) = g1 v1 (pm+3 ) = Gm = Hi (pm+3 ) = w(pm+3 ). Suppose that Ω1 < Ω2 . We have v(p + δ) = w(p + δ) and v (p + δ) > 0. According to Lemma 4.3(i), this is impossible. Suppose that Ω1 > Ω2 . We have v(p+δ) = w(p+δ) and w (p+δ) = Hi (pm+3 ) > 0. According to Lemma 4.3(ii), this is also impossible. Thus the conclusion is that there is no point of discontinuity of a on [pm+2 , pm+3 ) in this case. By assigning m := m + 1 one advances to the next observation interval [pm+3 , pm+4 ) and repeats the analysis of (v). (vi) Since it is already determined that there is a discontinuity point on interval [pm+2 , pm+3 ), the assumption a ∈ PC(σ) implies that a is constant on [pm+3 , pm+5 ]. This value ai+1 of a can be uniquely determined from the system in (vi) similarly to the argument presented in (iv). Note that Hi+1 (x) = g1 v1 (x) on [pm+3 , pm+5 ]. (vii) One knows that the discontinuity xi ∈ [pm+2 , pm+3 ) as well as the values ai and ai+1 of a on the adjacent intervals [pm+1 , pm+2 ] and [pm+3 , pm+4 ] together with the corresponding functions Hi (x) and Hi+1 (x). According to Lemma 4.4 one can determine the unique location of the discontinuity xi by the formulas given there. (viii) The advance of the algorithm from x = 1 to x = pl+1 is justified by reducing it to (iii)–(vii) using the change of variables z = 1 − x. (ix) Lemma 4.4 is applicable with η = 2δ. Note that there can be only one discontinuity of a on [pl−1 , pl+1 ], since 2δ < σ. The values of a as well as the corresponding functions Hj and Hj+1 are already known on the adjacent intervals. The discontinuity of a exists on [pl−1 , pl+1 ) if ωj = ωj+1 . The marching algorithm of Theorem 4.6 requires measurements of the system at a possibly large number of observation points. Our next theorem shows that if a piecewise constant conductivity a is known to have just one point of discontinuity x1 , and its values a1 and a2 are known beforehand, then the discontinuity point x1 can be determined from just one measurement of the heat conduction process. Theorem 4.7. Let p ∈ (0, 1) be an observation point, g(x) > 0 on (0, 1), and the observation zp (t) = u(xp , t; a), t ∈ (T1 , T2 ), of the heat conduction process (4.1) be given. Suppose that the conductivity a ∈ Aad is piecewise constant and has only one (unknown) point of discontinuity x1 ∈ (0, 1). Given positive values a1 = a2 such that
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY
711
a(x) = a1 for 0 ≤ x < x1 and a(x) = a2 for x1 ≤ x < 1, the point of discontinuity x1 is constructively identifiable. Proof. Arguing as in the previous theorem, zp (t) =
∞
gk e−λk t vk (p),
0 ≤ T 1 < t < T2 ,
k=1
where gk = g, vk for k = 1, 2, . . . . Since g1 v1 (p) > 0, the uniqueness of the Dirichlet series representation implies that one can uniquely determine the first eigenvalue λ1 and the value of Gp = g1 v1 (p). Without loss of generality one can assume that a1 > a2 . In this case we show that the first eigenvalue λ1 is strictly increasing as a function of x1 ∈ [0, 1]. Indeed, suppose that 0 ≤ xa1 < xb1 ≤ 1; that is, a(x) =
a1 , 0 < x < xa1 a2 , xa1 < x < 1
and b(x) =
a1 , 0 < x < xb1 a2 , xb1 < x < 1.
By Theorem 2.3(i) 1 1 1 b(x)[v1,b (x)]2 dx a(x)[v1,b (x)]2 dx a(x)[v (x)]2 dx 0 0 b 0 λ1 = 1 > 1 ≥ inf = λa1 1 1 (0,1) 2 dx 2 dx 2 dx v∈H [v (x)] [v (x)] [v(x)] 0 0 0 1,b 0 1,b provided that the derivative v1,b (x) of the first eigenfunction v1,b (x) is not identically a b zero on (x1 , x1 ). But, from (b(x)v1,b (x)) = −λb1 v1,b (x), the assumption v1,b (x) = 0 a b a b on (x1 , x1 ) implies v1,b (x) = 0 on (x1 , x1 ), and this is impossible, since v1,b (x) > 0 on (0, 1). Thus there exists a unique conductivity of the type sought in the theorem for which its first eigenvalue is equal to λ1 ; i.e., a is identifiable. Now the unique discontinuity point x1 of a can be determined as follows. Let λ1 λ1 ω1 = , ω2 = . a1 a2
Then the first eigenfunction v1 is given by A sin ω1 x, 0 < x < x1 , v1 (x) = (4.20) B sin ω2 (1 − x), x1 < x < 1, for some A, B > 0. The matching conditions at x1 give A sin ω1 x1 = B sin ω2 (1 − x1 ) and
A B cos ω1 x1 = cos ω2 (1 − x1 ). ω1 ω2
Since v1 (x1 ) > 0 we have 0 < ω1 x1 < π and 0 < ω2 (1 − x1 ) < π. Therefore x1 satisfies 1 1 cot ω1 x = cot ω2 (1 − x). ω1 ω2 The existence and the uniqueness of the solution x1 of the above nonlinear equation follows from the monotonicity and the continuity of the cotangent functions. Practically, the value of x1 can be found by a numerical method.
712
SEMION GUTMAN AND JUNHONG HA
5. Conclusions. The prevalent approach to parameter identification (estimation) problems is to find such parameters from the best fit to data minimization. However, such an approach usually does not guarantee the uniqueness of the identified parameters. The identifiability problem consists of finding sufficient conditions assuring such a uniqueness, and there have been just a few results for the identifiability in distributed parameter systems. In this paper we have shown that in some cases a variable conductivity in a 1D heat conduction process can be uniquely identified from observations of this process. The identifiability has been established for two sets of observations. In one case it is assumed that the conductivity is piecewise smooth, and we are given a sequence of distributed observations of the form zn (t) = u(x, t; a), ψn for n = 1, 2, . . . on a finite 2 time interval, where functions {ψN }∞ n=1 form a basis in H = L (0, 1). An algorithm for the conductivity identification is proposed. Its numerical study will be reported elsewhere. In the second case it is assumed that the conductivity is piecewise constant with sufficiently separated points of discontinuity. The observations of the process are taken at equidistant points pm ∈ (0, 1). The total number of points needed for the unique conductivity identification can be computed from a priori known parameters of the process. A marching algorithm for the conductivity identification is presented and justified. In both cases the plant does not require a special external input for its identifiability; i.e., it is modeled by ut = (aux )x rather than by ut = (aux )x +f (x, t). It will be of interest to extend the developed methods to vibration and steady-state processes. Our current research shows that the methods described in this paper can be extended to identifiability problems for heat conduction processes admitting various boundary (e.g., periodic) inputs and to other cases. A numerical implementation shows that the marching algorithm achieves a perfect identification for observations with low noise levels. These results will be presented elsewhere. REFERENCES [1] N. U. Ahmed, Optimization and identification of systems governed by evolution equations on Banach space, Pitman Res. Notes Math. Ser. 184, Longman Scientific and Technical, Harlow, UK, 1988. [2] G. Birkhoff and G. Rota, Ordinary Differential Equations, Wiley & Sons, New York, 1989. [3] A. Elayyan and V. Isakov, On uniqueness of recovery of the discontinuous coefficient of a parabolic equation, SIAM J. Math. Anal., 28 (1997), pp. 49–59. [4] I. M. Gelfand and B. M. Levitan, On the determination of a differential equation from its spectral function, Amer. Math. Soc. Transl. Ser. 2, 2 (1955), pp. 253–304. [5] S. Gutman, Identification of discontinuous parameters in flow equations, SIAM J. Control Optim., 28 (1990), pp. 1049–1060. [6] S. Kitamura and S. Nakagiri, Identifiability of spatially-varying and constant parameters in distributed systems of parabolic type, SIAM J. Control Optim., 15 (1977), pp. 785–802. [7] R. Kohn and M. Vogelius, Determining conductivity by boundary measurements II: Interior results, Comm. Pure Appl. Math., 41 (1988), pp. 865–877. [8] C. Kravaris and J. H. Seinfeld, Identification of parameters in distributed parameter systems by regularization, SIAM J. Control Optim., 23 (1985), pp. 217–241. [9] O. A. Ladyˇ zenskaja, V. A. Solonnikov, and N. N. Uralceva, Linear and quasi-linear equations of parabolic type, Trans. Math. Monogr., 23, AMS, Providence, RI, 1968. ´e, Partial Differential Equations with Numerical Methods, Springer[10] S. Larsson and V. Thome Verlag, New York, 2003. [11] J. L. Lions, Optimal Control of Systems Governed by Partial Differential Equations, SpringerVerlag, New York, 1971.
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY
713
[12] S. Nakagiri, Review of Japanese work of the last 10 years on identifiability in distributed parameter systems, Inverse Problems, 9 (1993), pp. 143–191. [13] Y. Orlov and J. Bentman, Adaptive distributed parameter systems identification with enforceable identifiability conditions and reduced-order spatial differentiation, IEEE Trans. Automat. Control, 45 (2000), pp. 203–216. [14] A. Pierce, Unique identification of eigenvalues and coefficients in a parabolic problem, SIAM J. Control Optim., 17 (1979), pp. 494–499. [15] S. Saks and A. Zygmund, Analytic Functions, Monografie Matematyczne, Warsaw, 1965. [16] E. Zauderer, Partial Differential Equations of Applied Mathematics, Wiley & Sons, New York, 1983.
SIAM J. CONTROL OPTIM. Vol. 46, No. 2, pp. 694–713
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY IN A HEAT CONDUCTION PROCESS∗ SEMION GUTMAN† AND JUNHONG HA‡ Abstract. We study the identification and identifiability problems for heat conduction in a nonhomogeneous rod. The identifiability results are established for two different sets of observations. Given a sequence of distributed type observations, the identifiability is proved for conductivities in a piecewise smooth class of functions. In the case of observations taken at finitely many points the identifiability is established for piecewise constant conductivities. Such conductivities can be uniquely identified using the proposed marching algorithm. Key words. identification, identifiability, piecewise constant conductivity AMS subject classifications. 35R30, 93B30 DOI. 10.1137/060657364
1. Introduction. Consider the heat conduction in a nonhomogeneous insulated rod of a unit length, with the ends kept at zero temperature at all times. Our main interest is in the identification and identifiability of the discontinuous conductivity (thermal diffusivity) coefficient a(x), 0 ≤ x ≤ 1. The identification problem consists of finding a conductivity a(x) in an admissible set K for which the temperature u(x, t) fits given observations in a prescribed sense. Under a wide range of conditions one can establish the continuity of the objective function J(a) representing the best fit to the observations. Then the existence of the best fit to data conductivity follows if the admissible set K is compact in the appropriate topology. However, such an approach usually does not guarantee the uniqueness of the found conductivity a(x). Establishing such a uniqueness is referred to as the identifiability problem. If the conductivity is identifiable and one can design an algorithm for its reconstruction, then we say that a is constructively identifiable. From physical considerations the conductivity coefficients a(x) are assumed to be in (1.1)
Aad = {a ∈ L∞ (0, 1) : 0 < ν ≤ a(x) ≤ μ}.
The temperature u(a) = u(x, t; a) inside the rod satisfies (1.2)
ut − (a(x)ux )x = 0, Q = (0, 1) × (0, T ), u(0, t) = u(1, t) = 0, t ∈ (0, T ), u(x, 0) = g(x), x ∈ (0, 1),
where g ∈ L2 (0, 1). In general, the solution of (1.2) is understood in the weak sense. According to [11] for any a ∈ Aad there exists a unique weak solution u(a) ∈ L2 (0, 1; H01 (0, 1)) ∩C([0, 1]; L2 (0, 1)), and so the map a → u(a) is well defined. Moreover, this map is continuous from Aad equipped with the L2 (0, 1) topology into ∗ Received by the editors April 15, 2006; accepted for publication (in revised form) January 25, 2007; published electronically May 7, 2007. http://www.siam.org/journals/sicon/46-2/65736.html † Department of Mathematics, The University of Oklahoma, Norman, Oklahoma 73019 (sgutman@ ou.edu). ‡ School of Liberal Arts, Korea University of Technology and Education, Cheonan 330-708, South Korea ([email protected]).
694
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY
695
C([0, T ]; L2 (0, 1)). In fact, in [11] these results are established for multidimensional parabolic problems. The identification (parameter estimation) problem for (1.2) is as follows: Find a conductivity a ∈ Aad such that the solution u(a) of (1.2) fits a given observation z of the heat conduction process. For example, given z ∈ L2 (0, 1) one defines (1.3)
J(a) = u(x, T ; a) − z(x)L2 (0,1) .
Then the parameter estimation problem for (1.2) is reduced to the minimization of the objective function J over the admissible set Aad or its subset Kad : Find a ∈ Kad ⊂ Aad such that (1.4)
J(¯ a) = inf{J(a) : a ∈ Kad }.
The above-mentioned properties of the solutions u(a) imply that the objective function J(a) is continuous on Aad ∩ L2 (0, 1). Therefore the identification problem (1.4) has a solution if Aad is compact in L2 (0, 1). One such choice is Kad = {a ∈ Aad ∩ H 1 (0, 1) : aH 1 ≤ constant}; see [8]. However, a ∈ H 1 (0, 1) implies that the conductivity is continuous. Therefore this choice of Kad is not suitable for the study of the identification problems with discontinuous coefficients. To overcome this difficulty we have shown in [5] (in a multidimensional case) that one can take for Kad the set of functions in Aad which have a uniformly bounded variation. Such a set Kad is compact in L2 (0, 1), and the existence of solutions to the identification problem (1.4) follows. See [5] for additional details and numerical experiments for 2D parameter identification problems. A variety of identification problems is studied in [1] under very general assumptions on the problem’s parameters. The identifiability questions for partial differential equations are much more difficult, and there are just a few available results. Suppose that one is given an observation z(t) = u(p, t; a) of the heat conduction process (1.2) for t1 < t < t2 at some observation point 0 < p < 1. From the series solution for (1.2) and the uniqueness of the Dirichlet series expansion (see section 2), one can, in principle, recover all of the eigenvalues of the associated Sturm–Liouville problem. If one also knows the eigenvalues for the heat conduction process with the same coefficient a and different boundary conditions, then the classical results of Gelfand and Levitan [4] show that smooth coefficients a(x) can be uniquely identified from the knowledge of the two spectral sequences. Also, if the entire spectral function is known (i.e., the eigenvalues and the values of the derivatives of the normalized eigenfunctions at x = 0), then the conductivity is identifiable as well. However, such results have little practical value, since the observation data z(t) always contain some noise, and therefore one cannot hope to adequately identify more than just the few first eigenvalues of the problem. A different approach is taken in [6, 12, 13, 14]. These works show that one can identify a constant conductivity a in (1.2) from the measurement z(t) taken at one point p ∈ (0, 1). These works also discuss problems more general than (1.2), including problems with a broad range of boundary conditions, nonzero forcing functions, as well as elliptic and hyperbolic problems. In [7, 3] and references therein identifiability results are obtained for elliptic and parabolic equations with discontinuous parameters in a multidimensional setting. A typical assumption there is that one knows the normal derivative of the solution at the boundary of the region for every Dirichlet boundary input. The main result of this paper is contained in Theorem 4.6. This theorem describes and justifies the marching algorithm for the unique identification of piecewise constant conductivities from observations of (1.2) given at finitely many points pk ∈ (0, 1). We
696
SEMION GUTMAN AND JUNHONG HA
start by recalling some basic properties of (1.2) in section 2. Identifiability results for countably many distributed observations are given in section 3. Identifiability of piecewise constant conductivities a is discussed in section 4. Numerical results for the identifiability algorithms described in this paper require an extensive exposition, and they will be presented elsewhere. 2. Auxiliary results. In this section we collect some well-known results for the solutions u(x, t; a) of (1.2), as well as for its associated Sturm–Liouville problem. Since such results are scattered in the literature, some brief proof outlines are included as well. See [2, 9, 10, 11] for a detailed discussion. Definition 2.1. Function a(x) is said to belong to the class PS N if (i) a ∈ Aad = {a ∈ L∞ (0, 1) : 0 < ν ≤ a(x) ≤ μ} for some positive constants ν and μ; (ii) function a is piecewise smooth; that is, there exists a finite sequence of points 0 = x0 < x1 < · · · < xN −1 < xN = 1 such that both a(x) and a (x) are continuous on every open subinterval (xi , xi+1 ), i = 0, . . . , N − 1, and both can be continuously extended to the closed intervals [xi , xi+1 ], i = 0, . . . , N − 1. For definiteness, we assume that a and a are continuous from the right, i.e., a(x) = a(x+) and a (x) = a (x+) for all x ∈ [0, 1). Also let a(1) = a(1−). Definition 2.2. PS = ∪∞ N =1 PS N . Everywhere in the following the conductivities a are assumed to be in PS. If a ∈ PS N , then the regularity conditions on a and the uniqueness of the weak solutions imply that for any t > 0 the weak solution u(x, t; a) of (1.2) satisfies the equation in the classical sense on any subinterval (xi , xi+1 ), i = 0, . . . , N − 1. Also u satisfies the matching conditions for the continuity of the solution and its conormal derivative at xi ∈ (0, 1), i = 1, 2, . . . , N − 1: ut − (a(x)ux )x = 0, x = xi , t ∈ (0, T ), u(0, t) = u(1, t) = 0, t ∈ (0, T ), u(xi +, t) = u(xi −, t), a(xi +)ux (xi +, t) = a(xi −)ux (xi −, t), u(x, 0) = g(x), x ∈ (0, 1),
(2.1)
where g ∈ L2 (0, 1); see [11, 16]. Denote by ·, ·, · the norm and the inner product, respectively, in H = L2 (0, 1). Theorem 2.3. Let a ∈ PS. Then (i) the associated Sturm–Liouville problem (a(x)v(x) ) = −λv(x), x = xi , v(0) = v(1) = 0, v(xi +) = v(xi −), a(xi +)vx (xi +) = a(xi −)vx (xi −)
(2.2)
has infinitely many eigenvalues 0 < λ1 < λ2 < · · · → ∞. The eigenvalues {λk }∞ k=1 and the corresponding orthonormal set of eigenfunctions {vk }∞ k=1 satisfy (2.3) λk = inf
1 0
a(x)[v (x)]2 dx 1 : v ∈ H0 (0, 1), < v, vj = 0, j = 1, 2, . . . , k − 1 , 1 [v(x)]2 dx 0
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY
(2.4)
1
λk =
697
a(x)[vk (x)]2 dx.
0 2 The normalized eigenfunctions {vk }∞ k=1 form a basis in L [0, 1]. (ii) Each eigenvalue is simple. For each eigenvalue λk there exists a unique continuous, piecewise smooth normalized eigenfunction vk (x) such that vk (0+) > 0, and the function a(x)vk (x) is continuous on [0, 1]. (iii) Eigenvalues {λk }∞ k=1 satisfy the inequality
νπ 2 k 2 ≤ λk ≤ μπ 2 k 2 . (iv) The first eigenfunction v1 satisfies v1 (x) > 0 for any x ∈ (0, 1). (v) The first eigenfunction v1 has a unique point of maximum q ∈ (0, 1) : v1 (x) < v1 (q) for any x = q. (vi) For any fixed t > 0 the solution u of (2.1) is given by u(x, t; a) =
∞
g, vk e−λk t vk (x),
k=1
and the series converges uniformly and absolutely on [0, 1]. (vii) For any p ∈ (0, 1) the function z(t) = u(p, t; a),
t > 0,
is real analytic on (0, ∞). Proof. (i) The proof is standard; see, e.g., [10]. (ii) On any subinterval (xi , xi+1 ) the coefficient a(x) has a bounded continuous derivative. Therefore, on any such interval the initial value problem (a(x)v (x)) +λv = 0, v(xi ) = A, v (xi ) = B has a unique solution. Suppose that two eigenfunctions w1 (x) and w2 (x) correspond to the same eigenvalue λk . Then they both satisfy the condition w1 (0) = w2 (0) = 0. Therefore their Wronskian is equal to zero at x = 0. Consequently, the Wronskian is zero throughout the interval (x0 , x1 ), and the solutions are linearly dependent there. Thus w2 (x) = Cw1 (x) on (x0 , x1 ), w2 (x1 −) = Cw1 (x1 −), and w2 (x1 −) = Cw1 (x1 −). The linear matching conditions imply that w2 (x1 +) = Cw1 (x1 +) and w2 (x1 +) = Cw1 (x1 +). The uniqueness of solutions implies that w2 (x) = Cw1 (x) on (x1 , x2 ), etc. Thus w2 (x) = Cw1 (x) on (0, 1), and each eigenvalue λk is simple. In particular λ1 is a simple eigenvalue. The uniqueness and the matching conditions also imply that any solution of (a(x)v (x)) + λv = 0, v(0) = 0, v (0) = 0 must be identically equal to zero on the entire interval (0, 1). Thus no eigenfunction vk (x) satisfies vk (0) = 0. Assuming that the eigenfunction vk is normalized in L2 (0, 1), it leaves us with the choice of its sign for vk (0). Letting vk (0) > 0 makes the eigenfunction unique. (iii) The eigenvalues of (2.3) satisfy the min-max principle 1 2 a(x)[v (x)] dx 0 λk = min max : v ∈ Vk , 1 Vk [v(x)]2 dx 0 where Vk varies over all subspaces of H01 (0, 1) of finite dimension k; see [10]. Therefore (a) (b) a(x) ≤ b(x), x ∈ [0, 1] implies λk ≤ λk . Since the eigenvalues of (2.3) with a(x) = 1 2 2 are π k , the required inequality follows.
698
SEMION GUTMAN AND JUNHONG HA
(iv) Recall that v1 (x) is a continuous function on [0, 1]. Suppose that there exists p ∈ (0, 1) such that v1 (p) = 0. Let wl (x) = v1 (x) for 0 ≤ x < p and wl (x) = 0 for p ≤ x ≤ 1. Let wr (x) = v1 (x) − wl (x), x ∈ [0, 1]. Then wl , wr are continuous, and, moreover, wl , wr ∈ H01 (0, 1). Also
1
wl (x)wr (x)dx = 0
1
and
0
a(x)wl (x)wr (x)dx = 0.
0
Suppose that wl is not an eigenfunction for λ1 . Then
1
a(x)[wl (x)]2 dx
1
[wl (x)]2 dx.
> λ1
0
0
Since
1
a(x)[wr (x)]2 dx ≥ λ1
0
1
[wr (x)]2 dx, 0
we have λ1 =
1
1
0
0
a(x)[v1 (x)]2 dx = 1 2 dx [v (x)] 1 0
a(x)([wl (x)]2 + [wr (x)]2 )dx 1 ([wl (x)]2 + [wr (x)]2 )dx 0
1 >
0
(λ1 [wl (x)]2 + λ1 [wr (x)]2 )dx = λ1 . 1 ([wl (x)]2 + [wr (x)]2 )dx 0
This contradiction implies that wl (and wr ) must be an eigenfunction for λ1 . However, wl (x) = 0 for p ≤ x ≤ 1, and as in (ii) it implies that wl (x) = 0 for all x ∈ [0, 1], which is impossible. Since v1 (0) > 0 the conclusion is that v1 (x) > 0 for x ∈ (0, 1). (v) From part (i), any eigenfunction vk is continuous and satisfies (a(x)vk (x)) = −λk vk (x) for x = xi . Also the function a(x)vk (x) is continuous on [0, 1] because of the matching conditions at the points of discontinuity xi , i = 1, 2, . . . , N − 1 of a. The integration gives x a(x)vk (x) = a(p)vk (p) − λk vk (s)ds p
for any x, p ∈ (0, 1). Let p ∈ (0, 1) be a point of maximum of vk . If p = xi , then vk (p) = 0. If p = xi , then vk (xi −) ≥ 0 and vk (xi +) ≤ 0. Therefore limx→p a(x)vk (x) = 0 and vk (p+) = vk (p−) = 0 since a(x) ≥ ν > 0. In any case for such a point p we have (2.5)
a(x)vk (x)
= −λk
x
vk (s)ds,
x ∈ (0, 1).
p
Since v1 (x) > 0, a(x) > 0 on (0, 1), (2.5) implies that v1 (x) > 0 for any 0 ≤ x < p and v1 (x) < 0 for any p < x ≤ 1. Since the derivative of v1 is zero at any point of maximum, we have to conclude that such a maximum p is unique.
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY
699
(vi) We prove only the convergence part. Note that 1 2 νvk ≤ a(x)[vk (x)]2 dx = λk vk 2 = λk . 0
Thus
√ λk vk ≤ √ ν
and
x
|vk (x)| ≤
|vk (s)|ds
≤
vk
0
√ λk ≤ √ . ν
Bessel’s inequality implies that the sequence of Fourier coefficients g, vk is bounded. Therefore, denoting by C various constants and using the fact that the function √ s → se−σs is bounded on [0, ∞) for any σ > 0, one gets √ λk t λ k λk t λk t −λk t | g, vk e vk (x)| ≤ C √ e− 2 e− 2 ≤ Ce− 2 . ν From (iii) of this theorem λk ≥ νπ 2 k 2 . Thus ∞
| g, vk e−λk t vk (x)| ≤ C
k=1
∞
e−
νπ 2 k2 t 2
k=1
≤C
∞
e−
νπ 2 t 2
k
< ∞.
k=1
∞ −λk t0 From (vi), the series vk (p) (vii) Let t0 > 0 and p ∈ (0, k=1 g, vk e 1). ∞ −λk s converges absolutely. Therefore k=1 g, vk e vk (p) is analytic in the part of the complex plane {s ∈ C : Re s > t0 }, and the result follows. ∞ Series of the form k=1 Ck e−λk t are known as the Dirichlet series. The following lemma shows that the Dirichlet series representation of a function is unique. Lemma 2.4. Let μk > 0, k = 1, 2, . . . , be a strictly increasing sequence. Suppose ∞ that T1 ≥ 0 and k=1 |Ck | < ∞. If ∞
Ck e−μk t = 0
for all
t ∈ (T1 , T2 ),
k=1
then Ck = 0 for k = 1, 2, . . . . ∞ −μk z The result follows at once from the observation that the series k=1 Ck e converges uniformly in the Re z > 0 region of the complex plane, implying that it is an analytic function there. See Chapter 9 of [15] for additional results on Dirichlet series. Remark. According to Theorem 2.3(vi) for each fixed p ∈ (0, 1) the solution z(t) = u(p, t; a) of (2.1) is given by a Dirichlet series. However, Lemma 2.4 is not directly applicable since the coefficients Ck = g, vk vk (p) are only square summable. Nevertheless, the conclusion of Lemma 2.4 remains valid, since the exponents μk in the Dirichlet series are the eigenvalues λk which satisfy the growth condition stated in Theorem 2.3(iii). This allows one to conclude (Theorem 2.3(vii)) that the solution z(t) is a real analytic function on (0, ∞), and the uniqueness of such a representation follows. Thus it would be a mistake to simply refer to the standard results such as Lemma 2.4 for the uniqueness of the Dirichlet series representation to justify the paper’s conclusions.
700
SEMION GUTMAN AND JUNHONG HA
3. Identifiability by distributed measurements. Suppose that one is given some observations of the heat conduction process (2.1) with an unknown conductivity a(x) and that they coincide with the observations of the model process
(3.1)
m m m um t ∈ (0, T ), t − (a (x)ux )x = 0, x = xi , um (0, t) = um (1, t) = 0, t ∈ (0, T ), m m um (xm i +, t) = u (xi −, t), m m m m m m am (xm +)u (x x i i +, t) = a (xi −)ux (xi −, t), um (x, 0) = g(x), x ∈ (0, 1),
where g is the same as in (2.1). The conductivity a is said to be identifiable in some class of functions M if the coincidence of the measurements of the observed and the model processes implies that a = am , provided a, am ∈ M. 2 Theorem 3.1. Let {ψn }∞ n=1 be a complete orthonormal set in H = L (0, 1). Suppose that nonzero initial data g ∈ H and the observations zn (t) = u(x, t; a), ψn for n = 1, 2, . . . and 0 ≤ T1 < t < T2 of the heat conduction process (2.1) are given. Then the conductivity a(x) ∈ Aad is constructively identifiable in the class of piecewise smooth functions PS. Proof. To show the identifiability of a we give an algorithm for its reconstruction from the data zn (t), n = 1, 2, . . . , guaranteeing the uniqueness in each step. Using Theorem 2.3(vi) we have (3.2)
zn (t) =
∞
g, vk e−λk t vk , ψn
k=1
for each n = 1, 2, . . . and 0 ≤ T1 < t < T2 . Fix an n > 0. Since ψn ∈ H and {vk }∞ k=1 form a basis in H, the Bessel inequality implies that the sequence of the Fourier coefficients { vk , ψn }∞ k=1 is bounded. From Theorem 2.3(vi) one concludes that the above series converges absolutely. Note that some products g, vk vk , ψn may be equal to zero. The zero value products present a difficulty, since we would not know how to associate the sequence of exponents recovered from (3.2) with the (unknown) eigenvalues λk : Some eigenvalues may be missing from the sequence. Define (possibly empty) subsets Qn ⊂ N by Qn = {k ∈ N : g, vk vk , ψn = 0},
n = 1, 2, . . . .
For Qn = ∅ reindex (3.2) so that there would be no vanishing coefficients: (3.3)
zn (t) =
∞
Cl,n e−μl,n t ,
t ∈ (T1 , T2 ).
l=1
By Theorem 2.3(vii) the solutions zn (t) are real analytic. Therefore, since all of the coefficients Cl,n = 0, one can uniquely determine them and the sequences ∞ μl,n , l = 1, 2, . . . . Recall that {μl,n }∞ l=1 ⊂ {λk }k=1 for any n with a nonempty Qn so each μl,n ≥ λ1 > 0. For each n such that Qn = ∅ let γn = min{μl,n : l ∈ N}. Define γ = min{γn : Qn = ∅}
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY
and (3.4)
701
⎧ ⎨ C1,n if γn = γ, Qn = ∅, 0 if γn > γ, Qn = ∅, An = ⎩ 0 if Qn = ∅.
Let (3.5)
w(x) =
∞
An ψn (x).
n=1
We claim that w is a nonzero multiple of some eigenfunction vJ of (2.1). Indeed, let J be the smallest index for which g, vJ = 0. Such an index exists since g = 0, and the eigenfunctions form a basis in H. Now, since vJ = 0 and {ψn }∞ n=1 also form a basis in H, there exists n ∈ N such that vJ , ψn = 0. Thus g, vJ vJ , ψn = 0 and γ ≤ λJ . The choice of J implies that γ = λJ . By the definition An = g, vJ vJ , ψn for nonzero products. Therefore w(x) =
∞
g, vJ vJ , ψn ψn (x)
n=1
= g, vJ
∞
vJ , ψn ψn (x) = g, vJ vJ (x)
n=1
as claimed. Now we show that the set of points y ∈ (0, 1) where w (y) = 0 is finite. Assuming the opposite, and since w is a nonzero multiple of vJ , there exists a sequence yj ∈ (0, 1) such that vJ (yj ) = 0 and limj→∞ yj = c ∈ [0, 1]. The continuity of a(x)vJ (x) implies that vJ (c) = 0. From (a(x)vJ (x)) = −γvJ (x) one gets yj+1 0 = a(yj+1 )vJ (yj+1 ) − a(yj )vJ (yj ) = −γ (3.6) vJ (s)ds yj
and concludes that vJ cannot be strictly positive or strictly negative on (yj , yj+1 ). Let ζj ∈ (yj , yj+1 ) be such that vJ (ζj ) = 0. Then limj→∞ ζj = c ∈ [0, 1] and vJ (c) = 0. Now we have both vJ (c) = 0 and vJ (c) = 0. But then the uniqueness of the Cauchy problem for the second order linear equations, and the matching conditions (see the proof of Theorem 2.3(ii)) imply that vJ (x) = 0 for all x ∈ [0, 1], which is impossible. Let q be a point of maximum of w. Note that w may have several maxima, unless it is a multiple of v1 . In any case, (2.5) implies x (3.7) w(s)ds, x ∈ (0, 1). a(x)w (x) = −γ q
Then, outside of the finite sets {xi } and {yj }, the conductivity a(x) is uniquely determined from (3.7) by x w(s)ds q . a(x) = −γ w (x) Because a is assumed to be in PS, it can be uniquely extended to the entire interval [0, 1].
702
SEMION GUTMAN AND JUNHONG HA
√ Remark 1. In an application one can choose ψn = 2 sin πnx and the initial condition g(x) > 0 on (0, 1). Then g, v1 v1 , ψ1 = 0, since v1 (x) > 0 on (0, 1) by Theorem 2.3(iv). Thus J = 1 in this case, and w(x) = g, v1 v1 (x) in the above algorithm. Also, one can see from (3.6) that there is only one point y1 = q ∈ (0, 1) where v1 (q) = 0, and it is the point of maximum of v1 (x) on (0, 1). Indeed, if there were two such points, then by (3.6) v1 (x) would have to become negative between them, which would contradict v1 (x) > 0 on (0, 1). Remark 2. Since the system {ψn }∞ n=1 is complete, the conditions of Theorem 3.1 imply that for any t > 0 one knows u(x, t; a) almost everywhere on [0, 1]. Since u(x, t; a) is continuous in x and analytic in t, Theorem 3.1, in fact, assumes that the solution is known in [0, 1] × (0, T ) or (equivalently) in [0, 1] × (0, ∞). Thus, Theorem 3.1 can be stated under any of these conditions. However, a practical reconstruction of the conductivity a directly from the equation ut = (aux )x is extremely unstable. The algorithm presented in the above theorem does not reconstruct the entire solution u but just the first eigenfunction of the associated Sturm–Liouville problem. Its stability properties will be studied elsewhere. 4. Identifiability of piecewise constant conductivities from finitely many observations. The identifiability is sought within the following set. Definition 4.1. Let PC ⊂ PS be the class of piecewise constant conductivities, and let PC N = PC ∩ PS N . Functions a ∈ PC N have the form a(x) = ai for x ∈ [xi−1 , xi ), i = 1, 2, . . . , N . In this case the governing system (2.1) is
(4.1)
ut − ai uxx = 0, x ∈ (xi−1 , xi ), u(0, t) = u(1, t) = 0, t ∈ (0, T ), u(xi +, t) = u(xi −, t), ai+1 ux (xi +, t) = ai ux (xi −, t), u(x, 0) = g(x), x ∈ (0, 1),
t ∈ (0, T ),
where g ∈ L2 (0, 1) and i = 1, 2, . . . , N − 1. The associated Sturm–Liouville problem is
(4.2)
ai v (x) = −λv(x), x ∈ (xi−1 , xi ), v(0) = v(1) = 0, v(xi +) = v(xi −), ai+1 v (xi +) = ai v (xi −)
for i = 1, 2, . . . , N − 1. We are interested only in the first eigenfunction v1 of (4.2). Let λ1 be the first eigenvalue. Suppose that p∗ ∈ (xi−1 , xi ). Then v1 can be expressed on (xi−1 , xi ) as
π λ1 π ∗ (x − p ) + γ , A > 0, − < γ < . v1 (x) = A cos ai 2 2 The range for γ in the above representation follows from the fact that v1 (p∗ ) = A cos γ > 0 by Theorem 2.3(iv). The identifiability of piecewise constant conductivities is based on the following three lemmas. Lemma 4.2. Suppose that δ > 0. Assume Q1 , Q3 ≥ 0, Q2 > 0, and 0 < Q1 + Q3 < 2Q2 . Let π π π Γ = (A, ω, γ) : A > 0, 0 < ω < , − 0 and cos γ > 0. Therefore cos(ωδ − γ) sin γ Q1 = cos(ωδ) + sin(ωδ) = , cos γ cos γ Q2 Q3 cos(ωδ + γ) sin γ = . = cos(ωδ) − sin(ωδ) cos γ cos γ Q2
(4.3) (4.4)
Adding (4.3) and (4.4) yields cos ωδ =
Q1 + Q3 . 2Q2
Since 0 < ωδ < π2 and 0 < (Q1 +Q3 )/2Q2 < 1 the above equation is uniquely solvable. Now subtracting (4.4) from (4.3) yields tan γ =
Q1 − Q3 , 2Q2 sin ωδ
which is also uniquely solvable, since −π/2 < γ < π/2. Finally, we have A = Q2 / cos γ. Lemma 4.3. Suppose that δ > 0, 0 < p ≤ x1 < p + δ < 1, 0 < ω1 , ω2 < π/2δ. Let w(x), v(x), x ∈ [p, p + δ] be such that w(x) = A1 cos ω1 x + B1 sin ω1 x, v(x) = A2 cos ω2 x + B2 sin ω2 x. Suppose that v(x1 ) = w(x1 ), v (x1 ) > 0,
ω12 v (x1 ) = ω22 w (x1 ),
v(x1 ) > 0.
Then
(i) conditions v(p + δ) = w(p + δ), v (p + δ) ≥ 0, and ω1 ≤ ω2 imply ω1 = ω2 ; (ii) conditions v(p + δ) = w(p + δ), w (p + δ) ≥ 0, and ω1 ≥ ω2 imply ω1 = ω2 . Proof. Since v(x1 ) > 0, v (x1 ) > 0 we have v(x) = A sin[ω2 (x − x1 ) + γ],
0 0. For ω2 > ω1 and v (p + δ) ≥ 0 one has cos[ω2 (p + δ − x1 ) + γ] =
1 v (p + δ) ≥ 0. ω2
Therefore ω2 + ω1 π (p + δ − x1 ) + γ < ω2 (p + δ − x1 ) + γ ≤ 2 2 and
cos
ω2 + ω 1 (p + δ − x1 ) + γ > cos[ω2 (p + δ − x1 ) + γ] ≥ 0. 2
Thus v(p + δ) − w(p + δ) > 0, and the conclusion (i) of the lemma follows. The case ω2 < ω1 and w (p + δ) ≥ 0 is reduced to the already established one by interchanging ω1 with ω2 and w with v. Lemma 4.4. Let δ > 0, 0 < η ≤ 2δ, ω1 = ω2 , with 0 < ω1 δ, ω2 δ < π/2. Also let A, B > 0, 0 ≤ p < p + η ≤ 1, and w(x) = A cos[ω1 (x − p) + γ1 ], v(x) = B cos[ω2 (x − p − η) + γ2 ], with |γ1 |, |γ2 | < π/2. Then the system (4.5) (4.6)
w(q) = v(q), ω22 w (q) = ω12 v (q),
(4.7)
w(q) > 0,
v(q) > 0
admits at most one solution q on [p, p + η]. This unique solution q can be computed as follows: If γ1 ≥ 0, then
B 2 − A2 1 − γ1 . (4.8) arctan ω1 2 2 q =p+ ω1 A ω2 − B 2 ω12 If γ2 ≤ 0, then (4.9)
B 2 − A2 1 − γ2 . − arctan ω2 2 2 q =p+η+ ω2 A ω2 − B 2 ω12
705
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY
Otherwise, compute q1 and q2 according to (4.8) and (4.9) and discard the one that does not satisfy the conditions of the lemma. Proof. Let α > 0 and cos t c(t; α) = , t ∈ R. α sin t Vector function c(t, α) traverses the ellipse E(1, α) centered in the origin with the x semiaxis equal to 1 and the y semiaxis equal to α. This function can be rewritten as c(t; α) = P(α)M(t)e1 , where M(t) =
cos t sin t
− sin t cos t
,
P(α) =
1 0
0 α
,
e1 =
1 0
.
Note that M(t) is the counterclockwise rotation in R2 by the angle t, P(α) is the α-contraction (expansion) in R2 along the y axis, and e1 is the standard basis vector along the x axis. Furthermore (4.10)
c(t1 + t2 ; α) = P(α)M(t1 )M(t2 )e1 .
With this notation system (4.5)–(4.6) is Ac(ω1 (q − p) + γ1 ; 1/ω1 ) = Bc(ω2 (q − p − η) + γ2 ; 1/ω2 ) or (4.11)
P(ω1−1 )M[ω1 (q − p) + γ1 ]Ae1 = P(ω2−1 )M[ω2 (q − p − η) + γ2 ]Be1 .
If q is a solution of (4.5)–(4.6), then the vectors in the right and the left sides of (4.11) are identical. Thus they belong to the intersection of the ellipses E(A, A/ω1 ) and E(B, B/ω2 ), and this intersection is not empty. In general the ellipses intersect in four points: one in each quadrant. Suppose that q ∗ = q is another solution of (4.5)–(4.6) on [p, p+η]. We can assume that q ∗ = q + τ for some τ > 0, 0 < ω1 τ, ω2 τ < π. System (4.5)–(4.6) at x = q ∗ is (4.12) P(ω1−1 )M[ω1 (q + τ − p) + γ1 ]Ae1 = P(ω2−1 )M[ω2 (q + τ − p − η) + γ2 ]Be1 . Using (4.10) and P(α)P(α−1 ) = I the right side of (4.12) can be written as P(ω2−1 )M[ω2 (q + τ − p − η) + γ2 ]Be1 = P(ω2−1 )M[ω2 τ ]M[ω2 (q − p − η) + γ2 ]Be1 = P(ω2−1 )M[ω2 τ ]P(ω2 )P(ω2−1 )M[ω2 (q − p − η) + γ2 ]Be1 = P(ω2−1 )M[ω2 τ ]P(ω2 )P(ω1−1 )M[ω1 (q − p) + γ1 ]Ae1 . Similarly the left side of (4.12) can be written as P(ω1−1 )M[ω1 (q + τ − p) + γ1 ]Ae1 = P(ω1−1 )M[ω1 τ ]M[ω1 (q − p) + γ1 ]Ae1 = P(ω1−1 )M[ω1 τ ]P(ω1 )P(ω1−1 )M[ω1 (q − p) + γ1 ]Ae1 .
706
SEMION GUTMAN AND JUNHONG HA
Let v = P(ω1−1 )M[ω1 (q − p) + γ1 ]Ae1 and D = P(ω2−1 )M[ω2 τ ]P(ω2 ) − P(ω1−1 )M[ω1 τ ]P(ω1 ). Then (4.12) is Dv = 0. Since v = 0 we must have det(D) = 0. Note that 1 det(D) = 2ω1 ω2 − (ω12 + ω22 ) sin ω1 τ sin ω2 τ − 2ω1 ω2 cos ω1 τ cos ω2 τ ω ω 1 2 1 1 1 2 2 2ω1 ω2 − (ω1 + ω2 ) cos(ω1 − ω2 )τ + (ω1 − ω2 ) cos(ω1 + ω2 )τ . = ω1 ω2 2 2 Let us define f (τ ) on [0, π/ω1 ) ∩ [0, π/ω2 ) as 1 1 f (τ ) = 2ω1 ω2 + (ω1 − ω2 )2 cos(ω1 + ω2 )τ − (ω1 + ω2 )2 cos(ω1 − ω2 )τ. 2 2 Function f is smooth on (0, π/ω1 ) ∩ (0, π/ω2 ), and its first and second derivatives are 1 1 f (τ ) = − (ω1 − ω2 )2 (ω1 + ω2 ) sin(ω1 + ω2 )τ + (ω1 + ω2 )2 (ω1 − ω2 ) sin(ω1 − ω2 )τ, 2 2 f (τ ) = (ω1 − ω2 )2 (ω1 + ω2 )2 sin ω1 τ sin ω2 τ. Since f (0) = f (0) = 0, f (τ ) > 0, and f (τ ) > 0 on (0, π/ω1 )∩(0, π/ω2 ), we conclude that f (τ ) > 0 for all τ ∈ (0, π/ω1 ) ∩ (0, π/ω2 ). Thus det(D) = 0 if and only if τ = 0. This contradicts the assumption τ > 0. Therefore the solution q of (4.5)–(4.7) is unique on [p, p + η]. To obtain formulas (4.8) and (4.9) notice that the ellipses E(A, A/ω1 ) and E(B, B/ω2 ) are given by x2 + ω12 y 2 = A2
and x2 + ω22 y 2 = B 2 .
At the intersection points we have y2 =
B 2 − A2 ω22 − ω12
and x2 =
A2 ω22 − B 2 ω12 . ω22 − ω12
The polar angle of the intersection point in the first quadrant is B 2 − A2 . ζ = arctan 2 2 A ω − B 2 ω2 2
1
Since w(q) = v(q) > 0 the intersection points corresponding to the solution q are in either the first or the fourth quadrants, and 0 ≤ ζ < π/2. If w (p) ≤ 0, then γ1 ≥ 0. Therefore 0 ≤ γ1 ≤ ω1 (q − p) + γ1 . In this case the intersection point is in the first quadrant. Accordingly tan[ω1 (q − p) + γ1 ] = ω1 tan ζ. Thus
B 2 − A2 1 − γ1 . arctan ω1 2 2 q =p+ ω1 A ω2 − B 2 ω12 If v (p+η) ≥ 0, then γ2 ≤ 0 and ω2 (q−p−η)+γ2 ≤ γ2 ≤ 0. In this case the intersection point is in the fourth quadrant. Accordingly tan[ω2 (q − p − η) + γ2 ] = −ω2 tan ζ and one gets (4.9).
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY
707
Now we would like to define a class of piecewise constant conductivities with sufficiently separated points of discontinuity. Definition 4.5. By the definition of a ∈ PC there exists N ∈ N and a finite sequence 0 = x0 < x1 < · · · < xN −1 < xN = 1 such that a is a constant on each subinterval (xn−1 , xn ), n = 1, . . . , N . Let σ > 0. Define PC(σ) = {a ∈ PC : xn − xn−1 ≥ σ,
n = 1, 2, . . . , N }.
Note that a ∈ PC(σ) attains at most N = [[1/σ]] distinct values ai , 0 < ν ≤ ai ≤ μ. The following theorem is our main result. It describes and justifies the marching algorithm for the unique identification of piecewise constant conductivities in the class PC(σ). Theorem 4.6. Given σ > 0 let an integer M be such that 3 μ and M > 2 M≥ . σ ν Suppose that the initial data g(x) > 0, 0 < x < 1, and the observations zm (t) = u(pm , t; a), pm = m/M for m = 1, 2, . . . , M − 1 and 0 ≤ T1 < t < T2 of the heat conduction process (4.1) are given. Then the conductivity a ∈ Aad is constructively identifiable in the class of piecewise constant functions PC(σ). First, we present the marching algorithm for the unique identification of the conductivity a and then justify it. The algorithm marches from the left end x = 0 to a certain observation point pl−1 ∈ (0, 1) and identifies the values an and the discontinuity points xn of the conductivity a on [0, pl−1 ]. Then the algorithm marches from the right end point x = 1 to the left until it reaches the observation point pl+1 ∈ (0, 1) identifying the values and the discontinuity points of a on [pl+1 , 1]. Finally, the values of a and its discontinuity are identified on the interval [pl−1 , pl+1 ]. The overall goal of the algorithm is to determine the number N − 1 of the discontinuities of a on [0, 1], the discontinuity points xn , n = 1, 2, . . . , N − 1, and the values an of a on [xn−1 , xn ], n = 1, 2, . . . , N (x0 = 0, xN = 1). As a part of the process the algorithm determines certain functions Hn (x) defined on intervals [xn−1 , xn ], n = 1, 2, . . . , N . The resulting function H(x) defined on [0, 1] is a multiple of the first eigenfunction v1 . Marching Algorithm. (i) Represent the data zm (t) as (4.13)
zm (t) =
∞
ck,m e−λk t ,
m = 1, 2, . . . , M − 1,
0 ≤ T 1 < t < T2 ,
k=1
and use it to uniquely identify the first eigenvalue λ1 and the coefficients Gm = c1,m , m = 1, 2, . . . , M − 1. Let G0 = GM = 0. (ii) Find l, 0 < l < M , such that Gl = max{Gm : m = 1, 2, . . . , M − 1} and Gm < Gl for any 0 ≤ m < l. (iii) Let i = 1, m = 0. (iv) Use Lemma 4.2 to find Ai , ωi , and γi from the system ⎧ ⎨ Ai cos(ωi δ − γi ) = Gm , Ai cos γi = Gm+1 , (4.14) ⎩ Ai cos(ωi δ + γi ) = Gm+2 . Let Hi (x) = Ai cos(ωi (x − pm+1 ) + γi ).
708
SEMION GUTMAN AND JUNHONG HA
(v) If m + 3 ≥ l, then go to step (viii). If Hi (pm+3 ) = Gm+3 , or Hi (pm+3 ) = Gm+3 and Hi (pm+3 ) ≤ 0, then a has a discontinuity xi on interval [pm+2 , pm+3 ). Proceed to the next step (vi). If Hi (pm+3 ) = Gm+3 and Hi (pm+3 ) > 0, then let m := m + 1 and repeat this step (v). (vi) Use Lemma 4.2 to find Ai+1 , ωi+1 , and γi+1 from the system ⎧ ⎨ Ai+1 cos(ωi+1 δ − γi+1 ) = Gm+3 , Ai+1 cos γi+1 = Gm+4 , (4.15) ⎩ Ai+1 cos(ωi+1 δ + γi+1 ) = Gm+5 . Let Hi+1 (x) = Ai+1 cos(ωi+1 (x − pm+4 ) + γi+1 ). (vii) Use the formulas in Lemma 4.4 to find the unique discontinuity point xi ∈ [pm+2 , pm+3 ). The parameters and functions used in Lemma 4.4 are defined as follows. Let p = pm+2 , η = δ. To avoid confusion we are going to use the notation Ω1 , Ω2 , Γ1 , Γ2 for the corresponding parameters ω1 , ω2 , γ1 , γ2 required in Lemma 4.4. Let Ω1 = ωi , Ω2 = ωi+1 . For w(x) use function Hi (x) recentered at p = pm+2 ; i.e., rewrite Hi (x) in the form w(x) = Hi (x) = A cos(Ω1 (x − pm+2 ) + Γ1 ),
|Γ1 | < π/2.
For v(x) use function Hi+1 recentered at p + η = pm+3 ; i.e., v(x) = Hi+1 (x) = B cos(Ω2 (x − pm+3 ) + Γ2 ),
|Γ2 | < π/2.
Let i := i + 1, m := m + 3. If m < l, then return to step (v). If m ≥ l, then go to the next step (viii). (viii) Do steps (iii)–(vii) in the reverse direction of x, advancing from x = 1 to x = pl+1 . Identify the values and the discontinuity points of a on [pl+1 , 1], and determine the corresponding functions Hi (x). (ix) Using the notation introduced in (vii) let Hj (x) be the previously determined function H on interval [pl−2 , pl−1 ]. Recenter it at p = pl−1 ; i.e., w(x) = Hj (x) = A cos(Ω1 (x − pl−1 ) + Γ1 ). Let Hj+1 (x) be the previously determined function H on interval [pl+1 , pl+2 ]. Recenter it at pl+1 : v(x) = Hj+1 (x) = B cos(Ω2 (x−pl+1 )+Γ2 ). If Ω1 = Ω2 , then stop; otherwise, use Lemma 4.4 with η = 2δ and the above parameters to find the discontinuity xj ∈ [pl−1 , pl+1 ]. Stop. Proof. To prove Theorem 4.6 we need to justify the marching algorithm and to show the uniqueness of the identification in each step. (i) Using Theorem 2.3(vi) we get (4.16) zm (t) =
∞
gk e−λk t vk (pm ),
m = 1, 2, . . . , M − 1,
0 ≤ T1 < t < T2 ,
k=1
where gk = g, vk for k = 1, 2, . . . . By Theorem 2.3(iv) v1 (x) > 0 on interval (0, 1). Since g is positive on (0, 1) we conclude that g1 v1 (pm ) > 0. According to Theorem 2.3(vii) each observation zm (t) is a real analytic function. Thus one can uniquely determine the nonzero coefficients in (4.16) and the corresponding exponents. In particular, one determines the first eigenvalue λ1 and the values of (4.17)
Gm = g1 v1 (pm ) > 0,
pm = m/M, m = 1, 2, . . . , M − 1.
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY
709
Because of the zero boundary conditions we can let G0 = GM = 0. The crucial point M −1 is that the numbers {Gm }m=1 are not arbitrary but are the values (up to a nonzero multiplicative constant g1 ) of the still-undetermined eigenfunction v1 . (ii) Let index l be defined as in (ii) of the marching algorithm. By Theorem 2.3(v) there exists a unique point q ∗ of maximum of v1 on (0, 1). Note that q ∗ ∈ (pl−1 , pl+1 ). Thus Gl+1 ≤ Gl and Gm < Gl for m > l + 1. Also v1 (pm −) > 0 for m = 1, 2, . . . , l − 1 and v1 (pm −) < 0 for m = l + 1, l + 2, . . . , M − 1. (iii) Start at the left end point p0 = 0 and work on interval [0, x1 ], where x1 is the first discontinuity point of a. (iv) Let δ = 1/M . Since σ ≥ 3δ and a ∈ PC(σ) we conclude that [0, p2 ) ⊂ [0, x1 ) and a = a1 on [0, x1 ). To apply Lemma 4.2 we just need to check the conditions for Q1 , Q2 , Q3 required there. We have Q1 = G0 = 0, Q2 = G1 = g1 v1 (p1 ) > 0 , Q3 = G2 = g1 v1 (p2 ) > 0. Let λ1 ω1 = . a1 By Theorem 2.3(iii) 0 < λ1 ≤ μπ 2 . Since 0 < ν ≤ a1 we have μπ 2 1 ν π 0 < ω1 δ < = . ν 2 μ 2 This inequality and v1 (x) > 0 on (0, 1) imply that the first eigenfunction v1 of (4.2) can be represented on (0, x1 ) as v1 (x) = C1 cos(ω1 (x − p1 ) + γ1 ) for some (C1 , ω1 , γ1 ) ∈ Γ, where Γ was defined in Lemma 4.2. Also Q1 + Q3 = g1 C1 (cos(ω1 δ + γ1 ) + cos(ω1 δ − γ1 )) = 2g1 C1 cos(ω1 δ) cos γ1 < 2G1 = 2Q2 since 0 < ω1 δ < π/2; hence 0 < cos(ω1 δ) < 1. Now Lemma 4.2 guarantees a unique solution of the system ⎧ ⎨ g1 C1 cos(ω1 δ − γ1 ) = G0 , g1 C1 cos γ1 = G1 , (4.18) ⎩ g1 C1 cos(ω1 δ + γ1 ) = G2 . It also gives formulas for the computation of A1 = g1 C1 , γ1 , ω1 from the known values of G0 , G1 , and G2 . Thus one can determine a1 = λ1 /ω12 and obtain (4.19)
H1 (x) = g1 C1 cos(ω1 (x − p1 ) + γ1 ) = g1 v1 (x)
for x ∈ [0, x1 ). (v) Let ai be the value of a on the part of the interval [pm+1 , pm+2 ) adjacent to [pm+2 , pm+3 ). By construction this value and the associated function Hi (x) = g1 v1 (x) are already determined by the algorithm. If there is no discontinuity of a on interval [pm+2 , pm+3 ), then a has the same value ai on interval [pm+2 , pm+3 ) as well. Therefore Hi (x) = g1 v1 (x) on this interval, and we must have Gm+3 = Hi (pm+3 ) by (4.17). If one has Gm+3 = Hi (pm+3 ), then the implication is that there is a discontinuity of a on [pm+2 , pm+3 ), and one proceeds to step (vi). On the other hand, if Gm+3 = Hi (pm+3 ), then one cannot, in general, conclude that there is no discontinuity of a on [pm+2 , pm+3 ). However, since we have m + 3 < l
710
SEMION GUTMAN AND JUNHONG HA
then (ii) of the proof implies that v1 (pm+3 −) > 0. Then the assumption a = ai on [pm+2 , pm+3 ) implies Hi (pm+3 ) = g1 v1 (pm+3 −). Therefore the equality Gm+3 = Hi (pm+3 ) together with Hi (pm+3 ) ≤ 0 lead to a contradiction. The conclusion is that Gm+3 = Hi (pm+3 ) and Hi (pm+3 ) ≤ 0 imply a discontinuity of a on [pm+2 , pm+3 ), and one proceeds to step (vi). Finally, one uses Lemma 4.3 to conclude that m + 3 < l, Gm+3 = Hi (pm+3 ), and Hi (pm+3 ) > 0 imply that there is no discontinuity of a on [pm+2 , pm+3 ). Indeed, suppose that there is a discontinuity point xi of a on interval [pm+2 , pm+3 ). Then a = ai on [pm+1 , xi ) and a = ai+1 on [xi , pm+3 ]. We are going to use the notation ω1 , and ω2 used in Lemma 4.3. Let xi , Ω1 , Ω2 for the corresponding variables x1 , p = pm+2 , p + δ = pm+3 , Ω1 = λ1 /ai , Ω2 = λ1 /ai+1 , and w(x) = Hi (x) = g1 v1 (x) = A1 cos Ω1 x + B1 sin Ω1 x, x ∈ [p, xi ), v(x) = g1 v1 (x) = A2 cos Ω2 x + B2 sin Ω2 x, x ∈ [xi , p + δ]. Note that the condition Ω21 v (xi ) = Ω22 w (xi ) is just the matching condition (4.2) at x = xi . Since m + 3 < l, the maximum q ∗ of v1 satisfies q ∗ > pm+3 . Because w is a positive multiple of v1 , it implies w(xi ) > 0 and w (xi ) > 0. Therefore v(xi ) > 0 and v (xi ) > 0. Because v is a positive multiple of v1 , we have v (p+δ) > 0. The condition v(p + δ) = w(p + δ) means v(pm+3 ) = g1 v1 (pm+3 ) = Gm = Hi (pm+3 ) = w(pm+3 ). Suppose that Ω1 < Ω2 . We have v(p + δ) = w(p + δ) and v (p + δ) > 0. According to Lemma 4.3(i), this is impossible. Suppose that Ω1 > Ω2 . We have v(p+δ) = w(p+δ) and w (p+δ) = Hi (pm+3 ) > 0. According to Lemma 4.3(ii), this is also impossible. Thus the conclusion is that there is no point of discontinuity of a on [pm+2 , pm+3 ) in this case. By assigning m := m + 1 one advances to the next observation interval [pm+3 , pm+4 ) and repeats the analysis of (v). (vi) Since it is already determined that there is a discontinuity point on interval [pm+2 , pm+3 ), the assumption a ∈ PC(σ) implies that a is constant on [pm+3 , pm+5 ]. This value ai+1 of a can be uniquely determined from the system in (vi) similarly to the argument presented in (iv). Note that Hi+1 (x) = g1 v1 (x) on [pm+3 , pm+5 ]. (vii) One knows that the discontinuity xi ∈ [pm+2 , pm+3 ) as well as the values ai and ai+1 of a on the adjacent intervals [pm+1 , pm+2 ] and [pm+3 , pm+4 ] together with the corresponding functions Hi (x) and Hi+1 (x). According to Lemma 4.4 one can determine the unique location of the discontinuity xi by the formulas given there. (viii) The advance of the algorithm from x = 1 to x = pl+1 is justified by reducing it to (iii)–(vii) using the change of variables z = 1 − x. (ix) Lemma 4.4 is applicable with η = 2δ. Note that there can be only one discontinuity of a on [pl−1 , pl+1 ], since 2δ < σ. The values of a as well as the corresponding functions Hj and Hj+1 are already known on the adjacent intervals. The discontinuity of a exists on [pl−1 , pl+1 ) if ωj = ωj+1 . The marching algorithm of Theorem 4.6 requires measurements of the system at a possibly large number of observation points. Our next theorem shows that if a piecewise constant conductivity a is known to have just one point of discontinuity x1 , and its values a1 and a2 are known beforehand, then the discontinuity point x1 can be determined from just one measurement of the heat conduction process. Theorem 4.7. Let p ∈ (0, 1) be an observation point, g(x) > 0 on (0, 1), and the observation zp (t) = u(xp , t; a), t ∈ (T1 , T2 ), of the heat conduction process (4.1) be given. Suppose that the conductivity a ∈ Aad is piecewise constant and has only one (unknown) point of discontinuity x1 ∈ (0, 1). Given positive values a1 = a2 such that
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY
711
a(x) = a1 for 0 ≤ x < x1 and a(x) = a2 for x1 ≤ x < 1, the point of discontinuity x1 is constructively identifiable. Proof. Arguing as in the previous theorem, zp (t) =
∞
gk e−λk t vk (p),
0 ≤ T 1 < t < T2 ,
k=1
where gk = g, vk for k = 1, 2, . . . . Since g1 v1 (p) > 0, the uniqueness of the Dirichlet series representation implies that one can uniquely determine the first eigenvalue λ1 and the value of Gp = g1 v1 (p). Without loss of generality one can assume that a1 > a2 . In this case we show that the first eigenvalue λ1 is strictly increasing as a function of x1 ∈ [0, 1]. Indeed, suppose that 0 ≤ xa1 < xb1 ≤ 1; that is, a(x) =
a1 , 0 < x < xa1 a2 , xa1 < x < 1
and b(x) =
a1 , 0 < x < xb1 a2 , xb1 < x < 1.
By Theorem 2.3(i) 1 1 1 b(x)[v1,b (x)]2 dx a(x)[v1,b (x)]2 dx a(x)[v (x)]2 dx 0 0 b 0 λ1 = 1 > 1 ≥ inf = λa1 1 1 (0,1) 2 dx 2 dx 2 dx v∈H [v (x)] [v (x)] [v(x)] 0 0 0 1,b 0 1,b provided that the derivative v1,b (x) of the first eigenfunction v1,b (x) is not identically a b zero on (x1 , x1 ). But, from (b(x)v1,b (x)) = −λb1 v1,b (x), the assumption v1,b (x) = 0 a b a b on (x1 , x1 ) implies v1,b (x) = 0 on (x1 , x1 ), and this is impossible, since v1,b (x) > 0 on (0, 1). Thus there exists a unique conductivity of the type sought in the theorem for which its first eigenvalue is equal to λ1 ; i.e., a is identifiable. Now the unique discontinuity point x1 of a can be determined as follows. Let λ1 λ1 ω1 = , ω2 = . a1 a2
Then the first eigenfunction v1 is given by A sin ω1 x, 0 < x < x1 , v1 (x) = (4.20) B sin ω2 (1 − x), x1 < x < 1, for some A, B > 0. The matching conditions at x1 give A sin ω1 x1 = B sin ω2 (1 − x1 ) and
A B cos ω1 x1 = cos ω2 (1 − x1 ). ω1 ω2
Since v1 (x1 ) > 0 we have 0 < ω1 x1 < π and 0 < ω2 (1 − x1 ) < π. Therefore x1 satisfies 1 1 cot ω1 x = cot ω2 (1 − x). ω1 ω2 The existence and the uniqueness of the solution x1 of the above nonlinear equation follows from the monotonicity and the continuity of the cotangent functions. Practically, the value of x1 can be found by a numerical method.
712
SEMION GUTMAN AND JUNHONG HA
5. Conclusions. The prevalent approach to parameter identification (estimation) problems is to find such parameters from the best fit to data minimization. However, such an approach usually does not guarantee the uniqueness of the identified parameters. The identifiability problem consists of finding sufficient conditions assuring such a uniqueness, and there have been just a few results for the identifiability in distributed parameter systems. In this paper we have shown that in some cases a variable conductivity in a 1D heat conduction process can be uniquely identified from observations of this process. The identifiability has been established for two sets of observations. In one case it is assumed that the conductivity is piecewise smooth, and we are given a sequence of distributed observations of the form zn (t) = u(x, t; a), ψn for n = 1, 2, . . . on a finite 2 time interval, where functions {ψN }∞ n=1 form a basis in H = L (0, 1). An algorithm for the conductivity identification is proposed. Its numerical study will be reported elsewhere. In the second case it is assumed that the conductivity is piecewise constant with sufficiently separated points of discontinuity. The observations of the process are taken at equidistant points pm ∈ (0, 1). The total number of points needed for the unique conductivity identification can be computed from a priori known parameters of the process. A marching algorithm for the conductivity identification is presented and justified. In both cases the plant does not require a special external input for its identifiability; i.e., it is modeled by ut = (aux )x rather than by ut = (aux )x +f (x, t). It will be of interest to extend the developed methods to vibration and steady-state processes. Our current research shows that the methods described in this paper can be extended to identifiability problems for heat conduction processes admitting various boundary (e.g., periodic) inputs and to other cases. A numerical implementation shows that the marching algorithm achieves a perfect identification for observations with low noise levels. These results will be presented elsewhere. REFERENCES [1] N. U. Ahmed, Optimization and identification of systems governed by evolution equations on Banach space, Pitman Res. Notes Math. Ser. 184, Longman Scientific and Technical, Harlow, UK, 1988. [2] G. Birkhoff and G. Rota, Ordinary Differential Equations, Wiley & Sons, New York, 1989. [3] A. Elayyan and V. Isakov, On uniqueness of recovery of the discontinuous coefficient of a parabolic equation, SIAM J. Math. Anal., 28 (1997), pp. 49–59. [4] I. M. Gelfand and B. M. Levitan, On the determination of a differential equation from its spectral function, Amer. Math. Soc. Transl. Ser. 2, 2 (1955), pp. 253–304. [5] S. Gutman, Identification of discontinuous parameters in flow equations, SIAM J. Control Optim., 28 (1990), pp. 1049–1060. [6] S. Kitamura and S. Nakagiri, Identifiability of spatially-varying and constant parameters in distributed systems of parabolic type, SIAM J. Control Optim., 15 (1977), pp. 785–802. [7] R. Kohn and M. Vogelius, Determining conductivity by boundary measurements II: Interior results, Comm. Pure Appl. Math., 41 (1988), pp. 865–877. [8] C. Kravaris and J. H. Seinfeld, Identification of parameters in distributed parameter systems by regularization, SIAM J. Control Optim., 23 (1985), pp. 217–241. [9] O. A. Ladyˇ zenskaja, V. A. Solonnikov, and N. N. Uralceva, Linear and quasi-linear equations of parabolic type, Trans. Math. Monogr., 23, AMS, Providence, RI, 1968. ´e, Partial Differential Equations with Numerical Methods, Springer[10] S. Larsson and V. Thome Verlag, New York, 2003. [11] J. L. Lions, Optimal Control of Systems Governed by Partial Differential Equations, SpringerVerlag, New York, 1971.
IDENTIFIABILITY OF PIECEWISE CONSTANT CONDUCTIVITY
713
[12] S. Nakagiri, Review of Japanese work of the last 10 years on identifiability in distributed parameter systems, Inverse Problems, 9 (1993), pp. 143–191. [13] Y. Orlov and J. Bentman, Adaptive distributed parameter systems identification with enforceable identifiability conditions and reduced-order spatial differentiation, IEEE Trans. Automat. Control, 45 (2000), pp. 203–216. [14] A. Pierce, Unique identification of eigenvalues and coefficients in a parabolic problem, SIAM J. Control Optim., 17 (1979), pp. 494–499. [15] S. Saks and A. Zygmund, Analytic Functions, Monografie Matematyczne, Warsaw, 1965. [16] E. Zauderer, Partial Differential Equations of Applied Mathematics, Wiley & Sons, New York, 1983.
