FRÉCHET DERIVATIVES FOR SOME BILINEAR INVERSE PROBLEMS

THOMAS DIERKES ^{∗}, OLIVER DORN ^{†}, FRANK NATTERER ^{‡}, VICTOR PALAMODOV ^{§}, AND HELMUT SIELSCHOTT ^{¶}

Abstract. In many inverse problems a functional of u is given by measurements where u solves a partial differential equation of the type L(p)u + Su = q. Here, q is a known source term and L(p), S are operators with p as unknown parameter of the inverse problem. For the numerical reconstruction of p often the heuristically derived Fréchet derivative R′ of the mapping R : p → ’measurement functional of u’ is used.

We show for three problems — a transport problem in optical tomography, an elliptic equation governing near infrared tomography, and the wave equation in moving media — that R′ is the derivative in the strict sense. Our method is applicable in more general problems than established methods for similar inverse problems.

AMS subject classifications. 35R30, 35R25, 65M32, 65N21, 92C55, 92F05, 35L20, 35J25

1. Introduction. This paper is motivated by a problem which often comes up in the numerical solution of inverse problems for partial differential equations, namely the differentiability of the operator R that maps the (sought for) coefficients of the partial differential equation to the (observable) quantity, e.g. the values of the normal derivative of the solution of the partial differential equation at the boundary. For many numerical methods — dual field methods [6], output least squares methods, Newton type methods, Kaczmarz’s method [10, 12, 13, 24, 25, 30] — the actual computation of the derivative R′ of R is of primary importance. It is fairly easy to derive expressions for R′ without actually proving that R′ is really a derivative, and most numerically oriented papers use these heuristically found derivatives with considerable success. In this paper we show that these heuristically obtained derivatives are in fact derivatives in the strict sense.

Our inverse problems are of the following bilinear type. Let L(p), S be linear differential operators in a domain X that may be a spatial domain Ω ⊆ IR^{n}, a domain Ω × (0, T) or even Ω × S^{n−1} × (0, T). Consider the problem

L(p)u + Su = q in X , u = f on ∂Ω. (1.1)

Here, L(p) is depending linearly on p, the sought-for parameter. We assume (1.1) to be uniquely solvable for u for any given q, f, p. The inverse problem calls for the determination of p from

R(p) = Mu

where M is the (linear) measurement operator.

Heuristically we can derive an expression for R′(p) as follows. Let u = u(p) be the solution of (1.1) and assume h to be a small perturbation of p. Then,

L(p + h)u(p + h) + Su(p + h) = q in X , u(p + h) = f on ∂Ω . (1.2)

Neglecting second order terms we obtain u(p + h) = u(p) + v where v satisfies

L(p)v + Sv = −L(h)u(p) in X , v = 0 on ∂Ω . (1.3)

∗Department of Computer Science, University College London, Gower Street, London WC1E 6BT, UK, t.dierkes@cs.ucl.ac.uk

†Department of Computer Science, University of British Columbia, 201-2366 Main Mall, Vancouver, B.C., V6T1Z4, Canada, dorn@cs.ubc.ca

‡Institut für Numerische und instrumentelle Mathematik, Universität Münster, Germany, frank.natterer@math.uni-muenster.de

§Department of Mathematics, Tel Aviv University, Israel, palamodo@math.tau.ac.il

¶Institut für Numerische und instrumentelle Mathematik, Universität Münster, Germany, helmut@sielschott.de

With v the — uniquely determined — solution to (1.3) it is natural to assume that the operator R^{0}(p) defined by

R^{0}(p)h = Mv (1.4)

is the derivative R′(p) of R at p. The purpose of this paper is to prove that for three selected problems this is actually the case.

In principle it is clear how to prove that (1.4) is in fact the derivative: We have to show that, in suitable function spaces,

∥R(p + h) − R(p) − R^{0}(p)h∥ = o(∥h∥)

as h → 0 for an appropriate norm ∥ · ∥. Subtracting (1.1) from (1.2) we obtain

u(p + h) − u(p) = −(L(p) + S)^{−1}L(h)u(p + h)
= −(L(p) + S)^{−1}L(h)u(p) − (L(p) + S)^{−1}L(h)(u(p + h) − u(p))
= v − (L(p) + S)^{−1}L(h)(u(p + h) − u(p)) .

It follows that

R(p + h) − R(p) − R^{0}(p)h = −M(L(p) + S)^{−1}L(h)(u(p + h) − u(p)) .

It remains to be shown that the right hand side is o(∥h∥).

This is fairly easy if S is an elliptic operator, L(p) of order not exceeding that of S, and M the normal derivative on ∂Ω and could be done along the lines of [8]. In § 3 we carry out the details for the problem of NIR tomography. For the case of optical tomography in time domain S is the transport operator. This case is treated in § 2. The much more difficult case of a hyperbolic operator S is treated in § 4 on ultrasound tomography in moving media. In all three cases we need very explicit a-priori estimates for the operators S the proofs of which are rather cumbersome.

Even though the use of derivatives is quite common for the numerical solution of inverse problems, the pertinent literature is surprisingly sparse. Many papers deal with the onedimensional case [3, 20], with L(p) being a lower order operator [4, 5], or with problems in which only the boundary of an object has to be identified [15, 17, 18, 26]. Fréchetdifferentiability for an inverse conductivity problem is proven in [11]. More recently, an inverse scattering problem for Maxwell’s equations in an infinite domain is investigated in [27, 28] and Fréchet-differentiability is shown by using an integral equation formulation. A differential equation satisfied by the derivative is then deduced from this result. A similar problem is treated in [1] directly in a differential equation framework for the case where near-field measurements are taken. We believe that the method followed in our paper is powerful enough to deal with quite general inverse problems. We will present three practically important examples in the following sections.

2. A transport problem in optical tomography.

2.1. The linear transport equation. We consider the time-dependent linear transport equation

\frac{∂u}{∂t} + θ · ∇u(x, θ, t) + a(x)u(x, θ, t) − b(x) ∫_{S^{n−1}} η(θ · θ′)u(x, θ′, t)dθ′ = q(x, θ, t) (2.1)

in Ω × S^{n−1} × [0, T] with initial condition

u(x, θ, 0) = 0 in Ω × S^{n−1} (2.2) and boundary condition

u(x, θ, t) = 0 on Γ_{−}. (2.3)

Here and in the following, Ω is a bounded domain in IR^{n}, n ≥ 2, with C^{2}-smooth boundary ∂Ω, ν(x) denotes the outward unit normal to ∂Ω at the point x ∈ ∂Ω, and Γ_{±} is defined as

Γ_{±} := {(x, θ, t) ∈ ∂Ω × S^{n−1} × [0, T], ±ν(x) · θ > 0} .

The function u(x, θ, t) is the density function of particles (photons) which travel in Ω at time t through the point x in the direction θ. Analogously, q is the density function of sources. Both u and q are assumed to be nonnegative and integrable functions in the domain X := Ω̄ × S^{n−1} × [0, T]. The velocity c of the particles is normalized to c = 1cms^{−1} and has been neglected in the formulation of (2.1).

The scattering function η(θ · θ′) describes the probability for a particle entering a scattering process with the direction of propagation θ′ to leave this process with the direction θ. It is assumed to be a real, nonnegative function and is normalized to

∫_{S^{n−1}} η(θ · θ′)dθ = 1, (2.4)

which expresses particle conservation in pure scattering events. The dot-product in the argument depends only on the cosine of the scattering angle cos ϑ := θ ·θ′ and, in particular, is independent of the location of the scattering event x. From now on we assume that the scattering function η is given and fixed, and that it satisfies (2.4).

The parameter a(x) is the total cross section or attenuation, and b(x) is the scattering cross section. The difference μ := a − b has the physical meaning of an absorption cross section which models the loss of particles from the system due to interactions with the background material. The parameters a, b and μ are assumed to be real, nonnegative and bounded functions of the position x.

In the application of optical tomography, the photons propagating in the medium are physically created by laser sources. These can be modeled mathematically either as an inhomogeneous incoming boundary condition, or as a nonzero source term q. Since an inhomogeneous boundary condition can be reformulated as an equivalent source distribution which is concentrated at the boundary [7], we have chosen to use the homogeneous boundary condition (2.3), and all photons in the medium are created by the source q. The modifications which are necessary for treating the alternative situation of a nonzero incoming boundary condition and a vanishing source term are straightforward, such that we will not discuss them here in particular. The initial condition (2.2) indicates that there are no photons moving inside of Ω at the starting time of our experiment. The inverse problem can now be formulated as follows.

Inverse problem. Assume that we measure for ℓ different source distributions q_{j}, j = 1, …, ℓ, some given functionals of the outgoing densities u_{j}|_{Γ_{+}} at the boundary ∂Ω. Provided with this information and knowing the scattering function η, how can we reconstruct both coefficients a and b of (2.1) inside of Ω simultaneously?

More information and references about this and related inverse problems are given in Arridge [2], Dierkes [10], and Dorn [12, 13].

2.2. Continuity. We say that a pair of parameters (a, b) is admissible if

a, b ∈ L^{∞}(Ω) with a > b > 0.

It is shown in [7] that, for admissible parameters (a, b) and with (2.4), the solution of the problem (2.1)-(2.3) is uniquely determined by specifying the initial condition (2.3), the incoming radiation at the boundary (2.2), and the source q. Furthermore, it can be shown that for admissible (a, b) and for sources q ∈ L^{1}(Ω × S^{n−1} × [0, T]) a solution of (2.1)-(2.4) exists which is in L^{1}(Ω × S^{n−1} × [0, T]). For details we refer to [9, 19]. In the present paper, however, the following result will be used.

Lemma 2.1. Let η be positive with (2.4), a(x) and b(x) be admissible, and assume that for the nonnegative source density q we have

Q̂ := ∫_{X} q(x, θ, t) dxdθdt < ∞. (2.5)

Then there exists a nonnegative solution of (2.1)-(2.3) that satisfies for any t ≥ 0 the inequality

∫_{S^{n−1}} ∫_{Ω} u(x, θ, t) dxdθ ≤ Q̂ exp(βt) (2.6)

with β := sup |b|. This solution is unique.

Proof. We prove the lemma by constructing the Neumann series solution to (2.1)-(2.3) as it is given for example in [7]. We have

u(x, θ, t) = Q(x, θ, t) + Ku(x, θ, t)

with

Q(x, θ, t) := ∫_{0}^{t} q(x − (t − s)θ, θ, s)e^{− ∫_{s}^{t} a(x−(t−τ)θ)dτ} ds, (2.7)

and

Ku(x, θ, t) := ∫_{0}^{t} ∫_{Ω} ∫_{S^{n−1}} δ(x′ − x + θ(t − t′)) (2.8)

× exp[ − ∫_{t′}^{t} a(x − θ(t − τ))dτ]b(x′)η(θ · θ′)u(x′, θ′, t′) dθ′dx′dt′

where q, a and b are assumed to be extended by zero outside of X. The Neumann series therefore has the form

u_{0}(x, θ, t) = Q(x, θ, t),

u_{j}(x, θ, t) = Ku_{j−1}(x, θ, t).

Integrating (2.7) against the density dxdθ and taking into account (2.5), we get

∫_{S^{n−1}} ∫_{Ω} u_{0}(x, θ, t) dxdθ ≤ Q̂.

Iterating this procedure and using (2.4), (2.5) and (2.8) yields

∫_{S^{n−1}} ∫_{Ω} u_{j}(x, θ, t) dxdθ ≤ Q̂\frac{β^{j}t^{j}}{j!}

for all integers j ≥ 0. This implies

∫_{S^{n−1}} ∫_{Ω} u(x, θ, t) dxdθ = ∑_{j=0}^{∞} ∫_{S^{n−1}} ∫_{Ω} u_{j}(x, θ, t) dxdθ ≤ Q̂ exp(βt),

which is just the estimate (2.6) of the lemma. The fact that u(x, θ, t) is nonnegative follows from the fact that every term in the series is nonnegative.

From now on we will assume that the conditions of Lemma 2.1 are satisfied. We introduce the following function spaces.

U := { u ∈ L^{1}(X) , (∂_{t} + θ · ∇)u ∈ L^{1}(X),
u(x, θ, 0) = 0 on Ω × S^{n−1} , u(x, θ, t) = 0 on Γ_{−} },
P := L^{∞}(Ω) × L^{∞}(Ω),  Z := L^{1}(∂Ω × [0, T]),

where U, P are equipped with the norms

∥u∥_{U} := ∥u∥_{1} + ∥(∂_{t} + θ · ∇)u∥_{1}, ∥v∥_{1} := ∫_{X} |v|dxdtdθ,

∥(a, b)∥_{P} := ∥a∥_{∞} + ∥b∥_{∞}.

Here P denotes the space of parameters (a, b), Z the space of data or measurements g, and U the space of states u. Moreover, we denote the space of source functions, equipped with the L^{1}-norm, by Y.

We mention that, in particular, we require in the definition of U that along almost every ray (x(s), t(s)) with x(s) = x_{0} + sθ, t(s) = t_{0} + s, s ≥ 0, (x_{0}, t_{0}) ∈ Ω × [0, T], the solution u(x(s), θ, t(s)) as well as its derivative (∂_{t} + θ · ∇)u(x(s), θ, t(s)) are integrable. Together with Fubini’s theorem this will ensure that, along almost every ray, u(x(s), θ, t(s)) is absolutely continuous. We will make use of this fact later when defining boundary values for u.

Notice that one reason for the choice of using L^{1}-norms for describing the states u and sources q in our problem is their physical significance in linear transport theory. L^{1}-norms are considered to be the ‘physically relevant’ norms here since the L^{1}-norm of u at a fixed time t_{0} describes the total number of particles N(t_{0}) moving at this time inside of Ω,

N(t_{0}) = ∫_{Ω} ∫_{S^{n−1}} u(x, θ, t_{0}) dxdθ.

An analogous interpretation holds for the sources q ∈ Y and for the measurements g ∈ Z which will be defined below. Notice that, according to Lemma 2.1, these quantities are in fact nonnegative since q is nonnegative. There is no obvious physical interpretation for e.g. the corresponding L^{2} norms.

The following continuity result will be needed in the proof of Fréchet differentiability. Lemma 2.2. Let (a, b) be admissible and q ≥ 0 with (2.5). For solutions u of (2.1)-(2.3)

we have the two estimates

∥u∥_{1} ≤ T ∥q∥_{1}, ∥u∥_{U} ≤ C_{1} ∥q∥_{1}, (2.9) where the constant C_{1} > 0 depends on a, b and T.

Proof. To start with we mention that

∫_{S^{n−1}} b(x) ( ∫_{S^{n−1}} η(θ · θ′)u(x, θ′, t)dθ′ ) dθ (2.10)

= b(x) ∫_{S^{n−1}} u(x, θ′, t) ( ∫_{S^{n−1}} η(θ · θ′)dθ) dθ′ = b(x) ∫_{S^{n−1}} u(x, θ′, t) dθ′.

Therefore, integration of (2.1) over Ω × S^{n−1} yields

∫_{S^{n−1}} ∫_{Ω} (∂_{t} + θ · ∇)u dxdθ = ∫_{S^{n−1}} ∫_{Ω} (−a + b)udxdθ + ∫_{S^{n−1}} ∫_{Ω} qdxdθ. (2.11)

The first term on the right hand side is non-positive, hence

∫_{S^{n−1}} ∫_{Ω} (∂_{t} + θ · ∇)u dxdθ ≤ ∫_{S^{n−1}} ∫_{Ω} q dxdθ. (2.12)

Take an arbitrary point (x, t) ∈ Ω × [0, T], and a vector θ ∈ S^{n−1} and consider the ray (x′(s), t′(s)) with x′(s) = x − sθ, t′(s) = t − s, s ≥ 0. Let (x_{0}, t_{0}) = (x′(s_{0}), t′(s_{0})) be the first point where the ray meets the boundary of the space-time domain Ω̄ × [0, T]. We have either t_{0} = 0 or x_{0} ∈ ∂Ω. We put s_{0} = t − t_{0}. Integration along the ray yields

u(x, θ, t) − u(x_{0}, θ, t_{0}) = ∫_{0}^{s_{0}}(∂_{t} + θ · ∇)u(x′, θ, t′) ds. (2.13)

Notice that the values on the left hand side of (2.13) are well-defined for almost every ray since (∂_{t} + θ · ∇)u(x′, θ, t′) ∈ L^{1} along the ray. In the first case we have u(x_{0}, θ, 0) = 0 by (2.2). In the second case we have either ν(x_{0}) · θ < 0 or ν(x_{0}) · θ = 0. In the case of the strict inequality, the function u vanishes by (2.3) in (x_{0}, θ, t_{0}). The equation ν(x_{0}) · θ = 0 means that the ray is tangent to the boundary of Ω. The set of the corresponding points (x, θ) is of zero measure and can be neglected when we integrate (2.13):

∫_{S^{n−1}} ∫_{Ω} udxdθ = ∫_{S^{n−1}} ∫_{Ω} ∫_{0}^{s_{0}}(∂_{t} + θ · ∇)u dsdxdθ ≤ ∫_{S^{n−1}} ∫_{Ω} ∫_{0}^{s_{0}}q dsdxdθ.

By changing variables x′ = x − sθ, t′ = t − s we get the inequality

∫_{S^{n−1}} ∫_{Ω} u(x, θ, t)dxdθ ≤ ∫_{S^{n−1}} ∫_{Ω} ∫_{0}^{t} q(x, θ, s) dsdxdθ.

Integration with respect to the time variable yields

∫_{X} u(x, θ, t) dxdθdt ≤ T ∫_{X} q(x, θ, t) dxdθdt,

and therefore the first inequality of (2.9). Integrating now (2.11) we get furthermore

∫_{X} |(∂_{t} + θ · ∇)u| dxdθdt ≤ ∫_{X} (a − b)|u| dxdθdt + ∫_{X} q dxdθdt

≤ ∥a − b∥_{∞} ∥u∥_{1} + ∥q∥_{1} ≤ (T ∥a − b∥_{∞} + 1)∥q∥_{1}.

It follows finally

∥u∥_{U} = ∥u∥_{1} + ∥ (∂_{t} + θ · ∇)u ∥_{1} ≤ (T ∥a − b∥_{∞} + T + 1)∥q∥_{1},

which concludes the proof.

For functions u ∈ U we define the measurement operator M by

M : U ⟶ Z, Mu(x, t) := ∫_{S^{n−1}} ν(x) · θ u(x, θ, t) dθ. (2.14)

The expression (2.14) describes typical measurements in the application of optical tomography, see for example [16].

Lemma 2.3. If the boundary ∂Ω is C^{2}-smooth, the operator M : U → Z is continuous and we have

∥Mu∥_{Z} ≤ ∥u∥_{U}. (2.15)

Proof. First note that the integration in (2.14) is actually performed only on those directions θ with ν(x) · θ > 0. All other directions do not contribute to the integration, because of the boundary condition u|_{Γ_{−}} = 0. For (y, θ, t) ∈ Γ_{+} we define

τ(y, −θ) := sup{ t , y − τ′θ ∈ Ω for 0 ≤ τ′ < t }.

The quantity τ(y, −θ) is just the travel time of a particle entering the domain Ω at a point y ∈ ∂Ω with direction −θ (and having normalized velocity) until the first arrival at the boundary of Ω. Consider the mapping

(2.16) φ : ∂Ω × S^{n−1} × IR × IR → IR^{n} × S^{n−1} × IR, φ(y, θ, t′, s) = (y + sθ, θ, t′ + s)

such that for a point (x, θ, t) ∈ Ω×S^{n−1} ×[0, T] we have (x, t) = (y, t′)+s(θ, 1) with y ∈ ∂Ω. Writing ũ = u ∘ φ it follows

u(x, θ, t) = ũ(y, θ, t′, s).

By the chain rule, the derivative of ũ with respect to the variable s is just

\frac{∂}{∂s}ũ(y, θ, t′, s) = (∂_{t} + θ · ∇) u(x, θ, t)

The pull back of the volume element dxdθdt is equal to

φ^{∗}(dxdθdt) = |J(φ)| dydθdt′ds,

where J(φ) is the determinant of the Jacobian of (2.16). It can be verified easily that |J(φ)| = |(ν · θ)^{−1}| on Γ. We integrate the function ũ from s = −τ(y, −θ) to s = 0. Doing this we get for each direction θ such that τ(y, −θ) > 0 the equality

|ũ(y, θ, t′, 0) − ũ(y, θ, t′, −τ(y, −θ))| = | ∫_{−τ(y,−θ)}^{0} \frac{∂ũ}{∂s} ds |

We have ũ(y, θ, t′, −τ(y, −θ)) = 0 because of the boundary condition u|_{Γ_{−}} = 0. Therefore

|ũ(y, θ, t′, 0)| ≤ ∫_{−τ(y,−θ)}^{0} |\frac{∂ũ}{∂s}| ds.

Integrate this inequality against the density ν · θ dydθdt′ over Γ_{+}:

∫_{Γ_{+}} |ν · θ ũ(y, θ, t′, 0)| dydθ dt′ ≤ ∫_{Γ_{+}} ∫_{−τ(y,−θ)}^{0} |\frac{∂ũ}{∂s}| ds ν · θ dydθdt′.

Returning to the coordinates (x, θ, t) we get finally

∫_{Γ_{+}} (ν · θ) u(y, θ, t)dydθdt ≤ ∫_{X} |(∂_{t} + θ · ∇)u| dxdθdt ≤ ∥u∥_{U} < ∞.

This is the estimate (2.15) of the lemma.

2.3. The Fréchet derivative. In the following we will use consequently the short notation p = (a, b) and h = (h_{a}, h_{b}) for elements of the parameter space P.

For functions u ∈ U and admissible p ∈ P we define the transport operator

Su = ∂_{t}u + θ · ∇u

and the scattering operator

L(p)u = au − b ∫_{S^{n−1}} ηu dθ′.

Equation (2.1) gets the form

L(p)u + Su = q. (2.17)

Let us assume now that we are given physical data g ∈ Z which correspond to the true parameter distribution p̃. The measurements which correspond to our guess p are given by Mu(p) where M is the measurement operator introduced in (2.14) and u(p) solves (2.17) with parameters p. For these functions we define the (nonlinear) residual operator R by

R : P ⟶ Z, (R(p))(x, t) := (Mu(p))(x, t) − g(x, t).

When solving our inverse problem we want to find a parameter distribution p̂ ∈ P such that the corresponding residual R(p̂) is minimized in some given norm. More information about a possible strategy for solving the inverse problem can be found in [13]. In the present paper we will concentrate on studying the properties of the residual operator R, in particular its Fréchet differentiability.

For a given p in P, we denote by U_{0}(p) a sufficiently small neighbourhood of p such that all points p + h ∈ U_{0}(p) are also admissible. For a point h = (h_{a}, h_{b}) with p + h ∈ U_{0}(p), we consider now the problem

\begin{matrix}\frac{∂w}{∂t} + θ · ∇w(x, θ, t) + a(x)w(x, θ, t) − b(x) ∫_{S^{n−1}} η(θ · θ′)w(x, θ′, t)dθ′ \\ = −h_{a}u + h_{b} ∫_{S^{n−1}} η(θ · θ′)u(x, θ′, t)dθ′ \end{matrix} (2.18)

in Ω × S^{n−1} × [0, T] with initial condition

w(x, θ, 0) = 0 in Ω × S^{n−1}

and boundary condition

w(x, θ, t) = 0 on Γ_{−},

where u is the solution of (2.1)-(2.3). Notice that the right hand side of (2.18) admits the physical interpretation of a ‘scattering source’ which is caused by the parameter perturbation h. With these functions, we define the (linear) operator R^{0}(p) by

R^{0}(p) : P ⟶ Z

(R^{0}(p)h)(x, t) = ∫_{S^{n−1}} ν · θ w(x, θ, t) dθ.

In the remainder of this section we will prove the following theorem.

Theorem 2.4. Let Ω be bounded with ∂Ω ∈ C^{2}. Let furthermore p = (a, b) ∈ P be admissible. For a given nonnegative source q ∈ Y, we assume that q ≥ 0 with (2.5). Then the mapping R is Fréchet differentiable at p, and the Fréchet derivative R′_{p} of R in p is given by R^{0}(p).

2.4. Proof of Fréchet differentiability. In this section we will prove Theorem 2.4. We consider the combined mapping

\begin{aligned}R : & P → & U → & Z \\ & p → & u(p) → & Mu(p) − g, \end{aligned}

with an arbitrary function g ∈ Z. Let p, h ∈ P with h so small that p + h is in U_{0}(p). Let u, v, w be solutions of

(2.19)  L(p)u + Su = q,
(2.20)  L(p + h)v + Sv = q,
(2.21)  L(p)w + Sw = −L(h)u.

Notice that (2.21) ist just (2.18). Since u ∈ U and p, p + h ∈ P are admissible, Lemma 2.1 can be applied to each of the problems (2.19)–(2.21), where, however, the solutions w of (2.21) are now allowed to assume negative values as well.

Lemma 2.5. We have the estimate

∥v − u∥_{U} ≤ C_{2} ∥h∥_{P} , (2.22)

with a constant C_{2} > 0 which depends on p, T and q.

Proof. Taking the difference between (2.20) and (2.19) yields

S(v − u) + L(p)(v − u) = −L(h)v.

Using (2.9), and a calculation analogous to (2.10), we get successively

∥v − u∥_{U} ≤ C_{1} ∥L(h)v∥_{1} ≤ 2C_{1}∥h∥_{P} ∥v∥_{1} ≤ 2C_{1}∥h∥_{P} T ∥q∥_{1}.

With C_{2} := 2T C_{1}∥q∥_{1} the statement of the lemma follows.

Lemma 2.6. Let z := v − u − w and let ∂Ω be C^{2}. Then the following estimate holds

∥Mz∥_{Z} ≤ C_{3} ∥h∥_{P}^{2} , (2.23) with a constant C_{3} > 0 which depends on p, T and q.

Proof. Taking the difference ‘(2.20)−(2.19)−(2.21)’ yields

Sz + L(p)z = −L(h)(v − u).

With (2.9), (2.22), and a calculation analogous to (2.10), we get successively

∥z∥_{U} ≤ C_{1} ∥L(h)(v − u)∥_{1} ≤ 2 C_{1} ∥h∥_{P} ∥v − u∥_{1}
≤ 2 C_{1} ∥h∥_{P} C_{2}∥h∥_{P} = C_{3}∥h∥_{P}^{2} ,

with the constant C_{3} := 2C_{1}C_{2}. Applying Lemma 2.3 we arrive at (2.23).

The statement of Theorem 2.4 follows now directly from the equation R(p + h) − R(p) − R^{0}(p)h = Mz and Lemma 2.6. This completes the proof of the theorem.

3. Near infrared tomography by an elliptic equation. Near infrared (NIR) or optical tomography has come to mean to identify regions of both higher scattering and absorption by transilluminating highly diffusive media, such as human tissue, with laser light in the wavelength range between 700 and 1000nm. Considering photons propagating through tissue as particles leads to the time-dependent linear transport equation as mathematical model, see § 2. A general closed solution formula to this integro-differential equation in higher dimensions is not known. Therefore, the linear transport equation is often reduced by using some approximation. Not only from a numerical point of view these reduction steps are interesting (permitting faster algorithms), but also theoretically within our context (transferring the differentiability result from § 2 to the approximation considered in this section).

In this section we will deal with the diffusion approximation of the linear transport equation. A more detailed derivation of this approximation from the linear transport equation, in time domain as well as in frequency domain, can be found in [16, 12, 2]. Here, we will first present this diffusion approximation in frequency domain, and will then prove the Fréchet differentiability of the corresponding forward operator for this reduced problem.

3.1. Governing equation. Let Ω ⊂ IR^{n} be a bounded, simply connected C^{∞}-domain. Let D : Ω → IR be a L^{∞}(Ω) function with a strictly positive lower bound, and let μ : Ω → IR be a strictly positive L^{∞}(Ω) function. Then, we consider the following form of the diffusion approximation in frequency domain,

−∇ · (D(x)∇u) + (μ(x) + i k)u = 0 in Ω, u|_{∂Ω} = f. (3.1)

The boundary distribution f ∈ H^{1/2}(∂Ω) in (3.1) is the spatial component of a timeharmonic modulated signal at the boundary. Accordingly, the solution u ∈ H^{1}(Ω) to (3.1) approximately gives the spatial component (integrated over all directions θ ∈ S^{n−1}) of the photon density function u(x, t, θ) of (2.1) in frequency domain. The real number k is the wave number k = ω/c of the modulation frequency ω.

The following theorem is needed for our differentiability proof.

Theorem 3.1. Let k ≥ 0 be a real number and f ∈ H^{1/2}(∂Ω). Let (D, μ) be a pair of two real-valued L^{∞}(Ω) functions that satisfies, for some m_{D}, M_{D}, M_{μ} > 0,

0 < m_{D} ≤ D(x) ≤ M_{D}, 0 < μ(x) ≤ M_{μ}.

Then, there exists a unique solution u ∈ H^{1}(Ω) for the Dirichlet problem (3.1). Moreover, there is a constant C_{Ω} > 0 only depending on ∂Ω such that

∥u∥_{H^{1}(Ω)} ≤ C_{Ω} (β C_{α,k} + 1) ∥f∥_{H^{1/2}(∂Ω)} (3.2) where C_{α,k} = 4α + 2k, α = max{ M_{μ}, M_{D} } and β = (1 + diam^{2} Ω)/m_{D}.

Proof. Define the differential operator A ≡ −∇ · (D(x)∇) + (μ(x) + i k), and let v = E_{1}f ∈ H^{1}(Ω) where E_{1} is a linear and bounded extension operator from H^{1/2}(∂Ω) to H^{1}(Ω) such that E_{1}f|_{∂Ω} = f (see e.g. [21] for existence of such extension operators). Consider the Dirichlet problem

Aw = h in Ω, w|_{∂Ω} = 0 (3.3)

where h = −Av. Multiplying the equation in (3.3) by an arbitrary \overline{z} ∈ H_{0}^{1}(Ω) and integrating over Ω yields the weak formulation

a(w, z) = b(z), z ∈ H_{0}^{1}(Ω)

for (3.3) where

a(w, z) = ∫_{Ω} (D(x) ∇w · \overline{∇z} + (μ(x) + i k) w\overline{z}) dx

and

b(z) = − ∫_{Ω} (D(x) ∇v · \overline{∇z} + (μ(x) + i k) v\overline{z}) dx.

Since Ω is open and bounded, Poincaré ’s inequality ∥z∥_{L^{2}} ≤ (diam Ω) ∥∇z∥_{L^{2}} holds for all z ∈ H_{0}^{1}(Ω), and the estimate

Re a(z, z) ≥ m_{D}∥∇z∥_{L^{2}}^{2} ≥ \frac{1}{β} ∥z∥_{H^{1}}^{2}, z ∈ H_{0}^{1}(Ω)

is true with β^{−1} = m_{D}/(1 + diam^{2} Ω). Furthermore, since we can estimate

|a(w, z)| ≤ M_{D} ∫_{Ω} |∇w| |∇z| dx + |M_{μ} + ik| ∫_{Ω} |w| |z| dx

≤ (2α + k) ∫_{Ω} (|∇w| + |w|)(|∇z| + |z|) dx
≤ 2 (2α + k) ∥w∥_{H^{1}} ∥z∥_{H^{1}},  w, z ∈ H^{1}(Ω)

with α = max{ M_{μ}, M_{D} }, it follows from these two estimates that there exists a unique solution w ∈ H_{0}^{1}(Ω) for the Dirichlet problem (3.3) by the theorem of Lax–Milgram. Therefore, the sum u := w + v ∈ H^{1}(Ω) solves the Dirichlet problem (3.1) in the weak sense. Additionally, the Lax–Milgram theorem yields

∥w∥_{H^{1}} ≤ β ∥h∥_{H^{−1}} ≤ β C_{α,k} ∥v∥_{H^{1}}

with C_{α,k} = 2(2α + k). Since v = E_{1}f where E_{1} is the bounded linear operator from H^{1/2}(∂Ω) to H^{1}(Ω) we have the estimate

∥v∥_{H^{1}(Ω)} ≤ C_{Ω} ∥f∥_{H^{1/2}(∂Ω)}

with the constant C_{Ω} = ∥E_{1}∥ > 0. Combining the last two inequalities we get (3.2). The uniqueness of the solution u ∈ H^{1}(Ω) follows because the difference between two different solutions must be zero by (3.2), and the theorem is proved.

3.2. Dirichlet-to-Neumann map Λ_{D,μ}. Define the Dirichlet-to-Neumann (DtN) map corresponding to D, μ by

Λ_{D,μ} : H^{1/2}(∂Ω) → H^{−1/2}(∂Ω), (3.4)

⟨Λ_{D,μ}f, g⟩ = ∫_{Ω} (D(x) ∇u · \overline{∇v} + (μ(x) + i k) u\overline{v}) dx

where v is some H^{1}(Ω) function with trace g on the boundary and u ∈ H^{1}(Ω) is a solution to the Dirichlet problem (3.1). Since this mapping is characterized by

Λ_{D,μ} : f ↦ (D ∂_{ν}u)|_{∂Ω}

it models boundary measurements of the outward photon flux which is given by the derivative in the direction of the outward normal of ∂Ω.

It is common for the application which we are addressing in this section to consider the whole DtN map rather than the forward operator corresponding to a single source position. Therefore, we slightly generalize the meaning of the operator R(p) which was introduced in § 1. It is now related to the whole DtN map, and is for this section defined by

R(p) = Λ_{D,μ}

where p = (D, μ) denotes the pair of unknown coefficients.

Clearly, the Dirichlet-to-Neumann map defined in (3.4) shows a bilinear structure with respect to the coefficients D, μ on the one hand and the solution u on the other. Thus, we will consider D as a parameter independent from μ in the following and prove the differentiability with respect to D and μ.

3.3. Continuity of Λ_{D,μ}. In order to show that a map is differentiable continuity must be settled first. This will be done in this section for the Dirichlet-to-Neumann map Λ_{D,μ} with respect to D, μ ∈ L^{∞}(Ω). By ∥ · ∥_{s,t} we denote the usual linear operator norm of a mapping from H^{s} to H^{t}. We will use the norm ∥ · ∥_{∞} for a pair p = (D, μ) which is given by ∥p∥_{∞} := ∥D∥_{L^{∞}} + ∥μ∥_{L^{∞}}. The next theorem can be considered as an extension of a result in [31, Thm. 0.2].

Theorem 3.2. Let k ≥ 0 be a real number. Suppose that p_{1} = (D_{1}, μ_{1}), p_{2} = (D_{2}, μ_{2}) are two pairs of real-valued L^{∞}(Ω) functions that satisfy, for some m_{D}, M_{D}, M_{μ} > 0

0 < m_{D} ≤ D_{1,2}(x) ≤ M_{D}, 0 < μ_{1,2}(x) ≤ M_{μ}.

Then, there exists a constant C_{Ω,α,β,k} > 0, independent of p_{1} and p_{2}, such that

∥R(p_{2}) − R(p_{1})∥_{1/2,−1/2} ≤ C_{Ω,α,β,k} ∥p_{2} − p_{1}∥_{∞} (3.5)

where α and β are the numbers given in Theorem 3.1 and R(p_{j}) = Λ_{D_{j}, μ_{j}}, j = 1, 2, denote the DtN maps with wave number k corresponding to each function pair for short.

Proof. Let f ∈ H^{1/2}(∂Ω) be arbitrary. Consider the difference z := u_{2} − u_{1} ∈ H_{0}^{1}(Ω) between two solutions to the Dirichlet problem (3.1) with respect to p_{1} = (D_{1}, μ_{1}) and p_{2} = (D_{2}, μ_{2}), respectively. This difference z is a (weak) solution to the inhomogeneous Dirichlet problem

−∇ · (D_{2}(x)∇z) + (μ_{2}(x) + i k)z = ∇ · ((D_{2} − D_{1})(x)∇u_{1}) − (μ_{2} − μ_{1})(x) u_{1} in Ω

with homogeneous boundary condition z|_{∂Ω} = 0.

Applying the theorem of Lax–Milgram to the weak formulation of the differential equation we see that the inequality

∥z∥_{H^{1}} = ∥u_{2} − u_{1}∥_{H^{1}} ≤ 2 β ∥p_{2} − p_{1}∥_{∞} ∥u_{1}∥_{H^{1}}

holds for the difference z. By using Theorem 3.1 we arrive at the estimate

∥u_{2} − u_{1}∥_{H^{1}} ≤ 2 β C_{Ω} (1 + β C_{α,k}) ∥p_{2} − p_{1}∥_{∞} ∥f∥_{H^{1/2}} (3.6)

with C_{α,k} = 4α + 2k and C_{Ω} > 0. Likewise, we have for the solution u_{2} ∈ H^{1}(Ω) the estimate

∥u_{2}∥_{H^{1}} ≤ C_{Ω} (1 + β C_{α,k}) ∥f∥_{H^{1/2}} (3.7)

by Theorem 3.1.

Now let g ∈ H^{1/2}(∂Ω) be an arbitrary boundary density with unit norm, ∥g∥_{H^{1/2}} = 1, and let v be an arbitrary H^{1}(Ω) function with trace v|_{∂Ω} = g. Then we have

sup_{∥g∥_{1/2}=1} |⟨(R(p_{2}) − R(p_{1}))f, g⟩|

= sup_{∥g∥_{1/2}=1} | ∫_{Ω} ((D_{2} − D_{1})(x)∇u_{2} · \overline{∇v} + D_{1}(x)∇(u_{2} − u_{1}) · \overline{∇v}

+ (μ_{2} − μ_{1})(x) u_{2}\overline{v} + (μ_{1}(x) + i k) (u_{2} − u_{1})\overline{v}) dx |
≤ 2 ∥p_{2} − p_{1}∥_{∞} ∥u_{2}∥_{H^{1}} sup_{∥g∥_{1/2}=1} ∥v∥_{H^{1}}

+ C_{α,k} ∥u_{2} − u_{1}∥_{H^{1}} sup_{∥g∥_{1/2}=1} ∥v∥_{H^{1}} (3.8)

(3.9)  ≤ 2 C_{Ω} (1 + β C_{α,k})^{2} ∥p_{2} − p_{1}∥_{∞} ∥f∥_{H^{1/2}}

where we have used (3.6) and (3.7) for the step from (3.8) to (3.9).

Thus, we can set C_{Ω,α,β,k} = 2 C_{Ω}(1 + βC_{α,k})^{2} which is greater than zero. Now, the inequality (3.5) follows from (3.9), and the theorem is proved.

3.4. Fréchet derivative of Λ_{D,μ}. Denote by L(H^{s}, H^{t}) the space of linear operators from H^{s} to H^{t}. Let h = (h_{1}, h_{2}) be an arbitrary pair of L^{∞}(Ω) functions and let p = (D, μ) be a pair of strictly positive L^{∞}(Ω) functions where the function D is bounded from below by a positive constant. With these settings, define a linear map

R^{0}(p) : L^{∞}(Ω) × L^{∞}(Ω) → L(H^{1/2}(∂Ω), H^{−1/2}(∂Ω))

by

⟨(R^{0}(p)h)f, g⟩ = ∫_{Ω} (D(x)∇w · \overline{∇v} + (μ(x) + i k) w\overline{v} + h_{1}(x)∇u · \overline{∇v} + h_{2}(x) u\overline{v}) dx

where v is some H^{1}(Ω) function with trace g on the boundary as in section 3.2, and w ∈ H^{1}(Ω) solves the inhomogeneous Dirichlet problem

−∇ · (D(x)∇w) + (μ(x) + i k)w = ∇ · (h_{1}(x)∇u) − h_{2}(x)u in Ω, w|_{∂Ω} = 0. (3.10)

The function u ∈ H^{1}(Ω) is a solution to the Dirichlet problem (3.1) with boundary values u|_{∂Ω} = f.

The next theorem gives the main result of this section.

Theorem 3.3. With the notations and assumptions of Theorem 3.2, the following is true. Independent of the difference h = p_{2} − p_{1}, there exists a constant C′_{Ω,α,β,k} > 0 such that

∥R(p_{2}) − R(p_{1}) − R^{0}(p_{1})h∥_{1/2, −1/2} ≤ C′_{Ω,α,β,k} ∥h∥_{∞}^{2} (3.11)

where R(p_{j}) = Λ_{D_{j}, μ_{j}}, j = 1, 2.

In particular, R^{0}(p_{1}) is the Fréchet derivative at p_{1} = (D_{1}, μ_{1}) of Λ_{D,μ} with respect to D, μ ∈ L^{∞}(Ω).

Proof. Let f ∈ H^{1/2}(∂Ω) be arbitrary, and let u_{1,2} ∈ H^{1}(Ω) be the solutions to the Dirichlet problem (3.1) with D(x) = D_{1,2}(x) and μ(x) = μ_{1,2}(x), respectively. Furthermore, let w ∈ H_{0}^{1}(Ω) be the solution to (3.10) with u = u_{1}. Considering the weak formulation of the differential equation which the difference u_{2} − u_{1} − w ∈ H^{1}(Ω) solves, we conclude as in the proof of Theorem 3.2 by the theorem of Lax–Milgram that the inequality

∥u_{2} − u_{1} − w∥_{H^{1}} ≤ 2 β ∥h∥_{∞} ∥u_{2} − u_{1}∥_{H^{1}} ≤ (2 β ∥h∥_{∞})^{2} ∥u_{1}∥_{H^{1}} (3.12)

holds. Thus, with v being any H^{1}(Ω) function with trace v|_{∂Ω} = g, we can estimate

sup_{∥g∥_{1/2}=1} |⟨(R(p_{2}) − R(p_{1}) − R^{0}(p_{1})h)f, g⟩|

= sup_{∥g∥_{1/2}=1} | ∫_{Ω} (D_{2}∇u_{2} · \overline{∇v} − D_{1}∇u_{1} · \overline{∇v} − D_{1}∇w · \overline{∇v} − h_{1}∇u_{1} · \overline{∇v}

+ (μ_{2} + i k) u_{2}\overline{v} − (μ_{1} + i k) u_{1}\overline{v} − (μ_{1} + i k) w\overline{v} − h_{2} u_{1}\overline{v}) dx |

= sup_{∥g∥_{1/2}=1} | ∫_{Ω} (h_{1}∇(u_{2} − u_{1}) · \overline{∇v} + h_{2} (u_{2} − u_{1})\overline{v}

+ D_{1} ∇(u_{2} − u_{1} − w) · \overline{∇v} + (μ_{1} + i k) (u_{2} − u_{1} − w)\overline{v}) dx |
(3.13)  ≤ (2 ∥h∥_{∞} ∥u_{2} − u_{1}∥_{H^{1}} + C_{α,k} ∥u_{2} − u_{1} − w∥_{H^{1}}) sup_{∥g∥_{1/2}=1} ∥v∥_{H^{1}}
(3.14)  ≤ 4 β (1 + β C_{α,k}) ∥h∥_{∞}^{2} ∥u_{1}∥_{H^{1}}
(3.15)  ≤ 4 β C_{Ω} (1 + β C_{α,k})^{2} ∥h∥_{∞}^{2} ∥f∥_{H^{1/2}}

where the constant C_{α,k} = 4α + 2k > 0 in (3.13) is given by Theorem 3.1. The step from (3.13) to (3.14) makes use of the estimate (3.12), and (3.15) arises from applying Theorem 3.1 to the solution u_{1} ∈ H^{1}(Ω), yielding the additional factor C_{Ω} (1 + β C_{α,k}) > 0.

With C′_{Ω,α,β,k} = 4 β C_{Ω}(1+β C_{α,k})^{2} and (3.15), we finally arrive at the inequality (3.11). This completes our proof of Theorem 3.3.

4. An inverse problem for the wave equation in moving media.

4.1. The wave equation in moving media. In a medium moving with velocity a(x), under assumptions of which some will be mentioned below, the acoustic wave equation is given by

u_{tt}(x, t) + 2a(x) · ∇u_{t}(x, t) − c^{2}(x)Δu(x, t) = q(x, t). (4.1)

As usual, c(x) denotes the speed of sound, u(x, t) is the acoustic pressure and q(x, t) is a source term. Equation (4.1) is hyperbolic. In the case a(x) = 0 it is the usual wave equation.

We consider (4.1) in a bounded domain Ω ⊂ IR^{n}, n ∈ {2, 3}, with homogeneous initial values and homogeneous Dirichlet boundary conditions. Equation (4.1) is deduced in [30] from the conservation laws of mass and moment. We mention that the standard reference for this derivation, [22], contains an error that is explained in [30]. The assumptions for the validity of (4.1) are

(i) div a = 0, (ii) ∥a∥_{∞} ≪ ∥c∥_{∞}, and (iii) the density of the medium is constant.

The above assumptions are nearly satisfied in acoustic pyrometry. This is a temperature measurement technique for burning gas. It is used in power stations where coal or waste are burned, see [29]. They match even better in ocean acoustic tomography (see [23]) as water is incompressible and the quotient ∥a∥_{∞}/∥c∥_{∞} is smaller.

We need to know that the homogeneous initial boundary value problem of (4.1) has a unique solution. This problem is usually considered for vanishing a(x). For the proof of Fréchet differentiability we also need some information about the regularity of the solution. The following theorems state the results on the forward problem. These results are well known for the wave equation in a medium at rest. We consider a moving medium and we show that the energy inequalities are valid uniformly in a neighbourhood of the parameter c.

Denote by C_{+}^{k}(Ω̄), k ∈ IN, the subset of strictly positive functions in C^{k}(Ω̄) and with ν the exterior normal on ∂Ω. Here, C^{k}(Ω̄) denotes the set of functions from C^{k}(Ω) that are continuous in Ω̄ and have continuous derivatives in Ω̄. We define X = Ω×]0, T [ and Γ = ∂Ω×[0, T]. We use the concept of generalized solutions as described in [14, sec. 7.2.1b].

In the following, we study the problem

u_{tt}(x, t) + 2a(x) · ∇u_{t}(x, t) − c^{2}(x)Δu(x, t) = q(x, t) in X,
u = 0  on Γ,
u = 0,  u_{t} = 0  on Ω × {0}.

(4.2)

This initial boundary value problem usually is not considered in text books as the term involving ∇u_{t} appears. Problem (4.2) is hyperbolic and existence of a generalized solution can be proved in the standard fashion, see [14, sec. 7.2]. However, we will need “energy estimates” of the type ∥u∥_{H^{m+1}(X)} ≤ C∥q∥_{H^{m}(X)} that are valid uniformly in a neighbourhood of the parameter c. In Theorem 4.1, we give an explicit upper bound for C in the case m = 1. In Theorem 4.3 we study the case m > 1.

Theorem 4.1. Let ∂Ω ∈ C^{1}, c ∈ C_{+}^{1}(Ω̄) with 0 < c_{0} ≤ c(x), c_{0} < 1, and |div c^{2}(x)| ≤ c_{1} in Ω̄. Furthermore, let a ∈ (C^{1}(Ω̄))^{n}, with ν · a|_{Γ} = 0 and |div a(x)| ≤ a_{1} in Ω̄, q ∈ L^{2}(X), and let u ∈ H^{1}(X) be a generalized solution to (4.2).

Then, for C_{0} := exp (c_{0}^{−2}(2 + 2a_{1} + c_{1})T ) we have

∥u∥_{H^{1}(X)}^{2} ≤ C_{0}∥q∥_{L^{2}(X)}^{2}. (4.3)

Proof. For short, we use the notation ∥u(·, t)∥_{0} = ∥u(·, t)∥_{L^{2}(Ω)}. Multiplying the differential equation by u_{t} and integrating over Ω×]0, τ[ with τ ∈]0, T] we obtain

∫_{Ω} ∫_{0}^{τ} (u_{t}u_{tt} + u_{t}a · ∇u_{t} − c^{2}u_{t}Δu) dxdt = ∫_{Ω} ∫_{0}^{τ} qu_{t} dxdt . (4.4)

16 T. Dierkes, O. Dorn, F. Natterer, V. Palamodov, H. Sielschott

We consider each of the terms on the left hand side separately. To begin with,

∫_{Ω} ∫_{0}^{τ} u_{t}u_{tt} dxdt = \frac{1}{2} ∫_{Ω} ∫_{0}^{τ} \frac{∂}{∂t}(u_{t})^{2} dxdt = \frac{1}{2} ∫_{Ω} u_{t}^{2}(x, τ) dx . (4.5)

Similarly, by Green’s formula and using that a(x) is tangent to ∂Ω,

2 ∫_{Ω} ∫_{0}^{τ} (a(x) · ∇u_{t})u_{t} dxdt = − ∫_{Ω} ∫_{0}^{τ} u_{t}^{2} div a(x) dxdt. (4.6)

Finally, since u = 0 on ∂Ω,

∫_{Ω} ∫_{0}^{τ} c^{2}u_{t}Δu dxdt = − ∫_{Ω} ∫_{0}^{τ} ∇u · ∇(c^{2}u_{t}) dxdt
= − ∫_{Ω} ∫_{0}^{τ} u_{t}∇u · ∇c^{2} dxdt − ∫_{Ω} ∫_{0}^{τ} c^{2}∇u · ∇u_{t} dxdt.

In the latter integral we can do an integration by parts with respect to t, yielding

∫_{Ω} ∫_{0}^{τ} c^{2}u_{t}Δu dxdt = − ∫_{Ω} ∫_{0}^{τ} u_{t}∇u · ∇c^{2} dxdt + ∫_{Ω} ∫_{0}^{τ} c^{2}∇u_{t} · ∇u dxdt

− ∫_{Ω} |∇u|^{2}(x, τ)c^{2} dx .

Adding the last two formulas and dividing by 2 we obtain

∫_{Ω} ∫_{0}^{τ} c^{2}u_{t}Δu dxdt = − ∫_{Ω} ∫_{0}^{τ} u_{t}∇u · ∇c^{2} dxdt − \frac{1}{2} ∫_{Ω} |∇u|^{2}(x, τ)c^{2} dx . (4.7)

We have (with Cauchy-Schwarz)

∫_{Ω} u^{2}(x, τ) dx = 2 ∫_{Ω} ∫_{0}^{τ} u(x, t)u_{t}(x, t) dxdt ≤ 2 ∫_{0}^{τ} ∥u∥_{0} ∥u_{t}∥_{0} dt . (4.8)

Inserting (4.5), (4.6), (4.7) in (4.4) yields

\frac{1}{2} ∫_{Ω} u_{t}^{2}(x, T) dx + \frac{1}{2} ∫_{Ω} |∇u|^{2}(x, T)c^{2} dx

= ∫_{Ω} ∫_{0}^{τ} qu_{t} dxdt + ∫_{Ω} ∫_{0}^{τ} u_{t}^{2} div a(x) dxdt − ∫_{Ω} ∫_{0}^{τ} u_{t}∇u · ∇c^{2} dxdt .

Multiplying by 2, using Cauchy-Schwarz and adding (4.8) we obtain

∫_{Ω} u^{2}(x, τ) dx + ∫_{Ω} u_{t}^{2}(x, τ) dx + c_{0}^{2} ∫_{Ω} |∇u|^{2}(x, τ) dx

≤ 2 ∫_{0}^{τ} ∥q∥_{0}∥u_{t}∥_{0} dt + 2a_{1} ∫_{0}^{τ} ∥u_{t}∥_{0}^{2} dt + 2c_{1} ∫_{0}^{τ} ∥u_{t}∥_{0}∥∇u∥_{0} dt + 2 ∫_{0}^{τ} ∥u∥_{0}∥u_{t}∥_{0} dt .

With U(τ) := ∫_{Ω} u^{2}(x, τ) dx + ∫_{Ω} u_{t}^{2}(x, τ) dx + ∫_{Ω} |∇u|^{2}(x, τ) dx the sum on the left is greater

than c_{0}^{2}U(τ). Using 2αβ ≤ α^{2} + β^{2} on the right hand side yields

U(τ) ≤ c_{0}^{−2} ∫_{0}^{τ} ∥q∥_{0}^{2} dt + c_{0}^{−2} ∫_{0}^{τ} (∥u∥_{0}^{2} + (2 + 2a_{1} + c_{1})∥u_{t}∥_{0}^{2} + c_{1}∥∇u∥_{0}^{2}) dt

≤ c_{0}^{−2} ∫_{0}^{τ} ∥q∥_{0}^{2} dt + \frac{2 + 2a_{1} + c_{1}}{c_{0}^{2}} ∫_{0}^{τ} (∥u∥_{0}^{2} + ∥u_{t}∥_{0}^{2} + ∥∇u∥_{0}^{2}) dt.

Using the Lemma of Gronwall ([14, Ineq. B.2.k]) we get

U(τ) ≤ exp (c_{0}^{−2}(2 + 2a_{1} + c_{1})τ)∥q(·, τ)∥_{0}^{2}.

Integrating this inequality over [0, T] we arrive at (4.3).

Theorem 4.2. Let the assumptions of Theorem 4.1 be satisfied.

Then, there exists a generalized solution u ∈ H^{1}(X) to (4.2). This solution is unique. Proof. We do not prove this result here. It can be deduced from Theorem 4.1 with a

Galerkin method. For the wave equation this is shown in [14, sec. 7.2.2]. For more details see [30].

Theorem 4.3. Let the assumptions of Theorem 4.1 be satisfied. In addition, we assume for some k ≥ 1: ∂Ω ∈ C^{k}, c ∈ C_{+}^{k}(Ω̄), a ∈ (C^{k}(Ω̄))^{n}, q ∈ H^{k}(X), and \frac{∂^{j}q}{∂t^{j}} = 0 for j = 0, …, k. Then, the generalized solution of (4.2) satisfies u ∈ H^{k+1}(X) and there is a constant C > 0 independent of q with

∥u∥_{H^{k+1}(X)} ≤ C_{k}∥q∥_{H^{k}(X)}. (4.9)

For each closed and bounded neighbourhood of c in C_{+}^{k}(Ω) we find a constant C_{k} such that (4.9) is true in this neighbourhood.

Proof. We consider the case k = 1 and formally differentiate the problem (4.2) with respect to t. We denote (formally) u_{t} by z:

z_{tt}(x, t) + 2a(x) · ∇z_{t}(x, t) − c^{2}(x)Δz(x, t) = q_{t}(x, t) in X,
z = 0  on Γ,
z = 0,  z_{t} = 0  on Ω × {0}.

(4.10)

The previous theorem states that a generalized solution z ∈ H^{1}(X) to (4.10) exists. Integrating z with respect to t we see z = u_{t} and we have

u_{tt} ∈ L^{2}(X) and u_{tx_{i}} ∈ L^{2}(X) (4.11)

for i = 1, …, n and from Theorem 4.1

∥u_{t}∥_{H^{1}(X)} ≤ C_{0}∥q_{t}∥_{L^{2}(X)}. (4.12)

Next, note that u is for each t ∈ [0, T] a solution to the elliptic boundary value problem

\begin{aligned}−Δu(·, t) & = & q_{e}(·, t) & in Ω, \\ u & = & 0 & on ∂Ω \end{aligned}

with

q_{e}(x, t) = c^{−2}(x)(q(x, t) − u_{tt}(x, t) − 2a(x) · ∇u_{t}(x, t)).

For this situation and q_{e}(·, t) ∈ L^{2}(Ω) we know from Theorem 4 in [14, sec. 6.3] that u(·, t) ∈ H^{2}(Ω) and that there is a number C_{e} > 0 with

∥u(·, t)∥_{H^{2}(Ω)} ≤ C_{e}(∥q_{e}(·, t)∥_{L^{2}(Ω)} + ∥u(·, t)∥_{L^{2}(Ω)}).

Together with (4.11) we conclude u ∈ H^{2}(X). Using Theorem 4.1 and then (4.12) we get

∥u(·, t)∥_{H^{2}(Ω)} ≤ C_{e}(c_{0}^{−2}(x)(q(x, t) − u_{tt}(x, t) − 2a(x) · ∇u_{t}(x, t)) + C_{0}∥q∥_{L^{2}(X)})
≤ C_{e}c_{0}^{−2}C_{0}∥q∥_{H^{1}(X)}.

This proves (4.9) for k = 1 with C = C_{e}c_{0}^{−2}C_{0}.

The constant C depends on the upper and lower bounds of c and its derivatives. For any closed and bounded neighbourhood of c in C_{+}^{1}(Ω̄), by increasing C we make sure that (4.9) is true in this neighbourhood.

For k > 1 the result follows by repeating the above steps, thus, by an induction approach.

4.2. Fréchet differentiability of the residual operator. To state the main result, we put p(x) = c^{2}(x). The space of parameters is given by P = C_{+}^{2}(Ω̄) with ∥·∥_{P} = ∥·∥_{H^{2}(Ω)} and the space of boundary data is Z = H^{1/2}(Γ). We denote by ν the exterior normal on ∂Ω.

Theorem 4.4. Let T > 0, Ω ⊂ IR^{n} a bounded domain with n ∈ {2, 3}, ∂Ω ∈ C^{4}, a ∈ (C^{2}(Ω̄))^{n} with ν·a|_{Γ} = 0, q ∈ H^{2}(X), and g ∈ Z. For p ∈ P we denote by u(p) ∈ H^{2}(X) the solution to (4.2). The residual operator

R : P ⟶ Z, p ↦ u(p)|_{Γ} − g

is Fréchet-differentiable for any p ∈ P. The derivative is given by R′(p) = R^{0}(p), with R^{0}(p) defined as

R^{0}(p) : P → Z,
R^{0}(p)h = w|_{Γ}.

Here, w ∈ H^{1}(X) is the solution to

w_{tt} + 2a · ∇w_{t} − pΔw = −hΔu(p) in X,
w = 0  on Γ,
w = 0, w_{t} = 0  in Ω × {0}.

Remark: We formulate the inverse problem. In acoustic pyrometry, we have a good guess of the velocity of the medium from vector tomography, see [29]. There is no (physical) need for improving this guess as we know ∥a∥_{∞} ≪ ∥c∥_{∞}. We also assume to know the source q(x, t), so the unknown parameter is p. By p̃ we denote the true value of p. The measurement g on ∂Ω is given by

g(x, t) = u(p̃)(x, t) on Γ.

The inverse problem asks for a function p ∈ P that minimizes ∥R(p)∥_{Z}.

Proof of Theorem 4.4. Choose p ∈ P fixed. We have u := u(p) ∈ H^{3}(X). Let h ∈ P be given with p(x)+h(x) > 0 a.e. in Ω̄. Let w be chosen as in the theorem. From hΔu ∈ H^{1}(X) we get w ∈ H^{2}(X).

Now consider the solution v ∈ H^{3}(X) to the problem

v_{tt} + 2a · ∇v_{t} − (p + h)Δv = q

with homogeneous initial values and boundary conditions. According to Theorem 4.3 there exists a real number C_{1} > 0 with

∥v∥_{H^{3}(X)} ≤ C_{1}∥q∥_{H^{2}(X)} . (4.13)

As this inequality is true in any closed neighbourhood of p (by enlarging C_{1}), the constant C_{1} is independent of h. The function v − u ∈ H^{3}(X) is in X a solution to the wave equation

(v − u)_{tt} + 2a · ∇(v − u)_{t} − pΔ(v − u) = hΔv

with homogeneous initial values and Dirichlet boundary data. As we have hΔv ∈ H^{1}(X) we know from Theorem 4.3 that there is a real number C_{2} > 0 with

∥v − u∥_{H^{2}(X)} ≤ C_{2} ∥hΔv∥_{H^{1}(X)}
≤ C_{2} ∥h∥_{P} ∥Δv∥_{H^{1}(X)}
≤ C_{2} ∥h∥_{P} ∥v∥_{H^{3}(X)}
(4.14)  ≤ C_{1}C_{2} ∥h∥_{P} ∥q∥_{H^{2}(X)} ,

where we have used (4.13). The function v − u − w ∈ H^{2}(X) solves in X the wave equation

(v − u − w)_{tt} + 2a · ∇(v − u − w)_{t} − pΔ(v − u − w) = hΔ(v − u)

with homogeneous initial values and Dirichlet boundary data. Because of hΔ(v−u) ∈ L^{2}(X) we know from Theorem 4.3 that there is a real number C_{3} > 0 with

∥v − u − w∥_{H^{1}(X)} ≤ C_{3} ∥hΔ(v − u)∥_{L^{2}(X)}
≤ C_{3} ∥h∥_{P} ∥Δ(v − u)∥_{L^{2}(X)}
≤ C_{3} ∥h∥_{P} ∥v − u∥_{H^{2}(X)}
≤ C_{1}C_{2}C_{3} ∥q∥_{H^{2}(X)} ∥h∥_{P}^{2} .

For the final inequality we used (4.14). We have thus proved Fréchet-differentiability:

R(p + h) − Rp − R^{0}(p)h = (v − u − w)|_{Γ} = o(∥h∥_{P})

as the trace theorem for n ∈ {2, 3} gives us a real number C_{4} > 0 with

\begin{aligned}∥(v − u − w)|_{Γ}∥_{Z} & ≤ & C_{4}∥v − u − w∥_{H^{1}(X)} \\ & ≤ & C_{1}C_{2}C_{3}C_{4} ∥q∥_{H^{2}(X)} ∥h∥_{P}^{2} . \end{aligned}

This result is not achievable with the implicit function approach of [8], as solutions to the wave equation are not smooth enough to satisfy the requirements of Theorem 2.1 in [8].

To see this, note that for q ∈ H^{k}(X) in Theorem 4.3 we only have u ∈ H^{k+1}(X). The function w solves the wave equation with the source term q = −hΔu and so from u ∈ H^{k+1}(X) we only get w ∈ H^{k}(X). For the contraction argument that proves Theorem 2.1 in [8] it would be necessary to have w ∈ H^{k+1}(X). This is satisfied in elliptic problems, but in this hyperbolic problem one degree of smoothness would be lost at each step of the contraction. To apply the implicit function approach we would have to assume p = c^{2}, h, and q to be infinitely often differentiable.

The advantage of our proof is the finite number of such steps: We only use the inequality (4.9) twice.

REFERENCES [1] H. Ammari and G. Bao, Analysis of the scattering map of a linearized inverse medium problem for

electromagnetic waves, Inverse Problems, 17 (2001), pp. 219–234.

[2] S. R. Arridge, Optical tomography in medical imaging, Inverse Problems, 15 (1999), pp. R41–R93.

[3] A. Bamberger, G. Chavent, and P. Lailly, About the stability of the inverse problem in 1–d wave equations — application to the interpretation of seismic profiles, Appl. Math. Optim., 5 (1979), pp. 1–47.

[4] G. Bao and W. W. Symes, On the sensitivity of hyperbolic equations to the coefficient, Comm. Partial Differential Equations, 21 (1996), pp. 395–422.

[5] G. Bao and W. W. Symes, Regularity of an inverse problem in wave propagation, in Inverse Problems of Wave Propagation and Diffraction, G. Chavent and P. C. Sabatier (eds.), Springer, 1997.

[6] H. D. Bui, Inverse Problems in the Mechanics of Materials, CRC Press, 1994.

[7] K. M. Case and P. F. Zweifel, Linear Transport Theory, Plenum Press, New York, 1967.

[8] T. J. Connolly and J. N. Wall, On Fréchet differentiability of some nonlinear operators occuring in inverse problems: an implicit function theorem approach, Inverse Problems, 6 (1990), pp. 949–966.

[9] R. Dautray and J. L. Lions, Evolution Problems II, Vol. 6 of Mathematical Analysis and Numerical Methods for Science and Technology, Springer, Berlin, 1993.

[10] T. Dierkes, Rekonstruktionsverfahren zur optischen Tomographie, Ph.D. thesis, Preprints Angewandte Mathematik und Informatik, 16/00-N, Münster, 2000.

Available at www.math.uni-muenster.de/num/Preprints/2000/dierkes .

[11] D. Dobson, Convergence of a reconstruction method for the inverse conductivity problem, SIAM J. Appl. Math., 52 (1992), pp. 442–458.

[12] O. Dorn, Das inverse Transportproblem in der Lasertomographie, Ph.D. thesis, Preprints Angewandte Mathematik und Informatik, 7/97-N, Münster, 1997. Available at

www.math.uni-muenster.de/num/Preprints/1997/dorn

[13] O. Dorn, A transport-backtransport method for optical tomography, Inverse Problems, 14 (1998), pp. 1107–1130.

[14] L. C. Evans, Partial Differential Equations, Graduate Studies in Mathematics, Vol. 19, AMS, 1998.

[15] F. Hettlich, Frechet derivatives in inverse obstacle scattering, Inverse Problems, 11 (1995), pp. 371– 382.

[16] J. P. Kaltenbach and M. Kaschke, Frequency and Time-Domain Modelling of Light Transport in Random Media, in Medical Optical Tomography: Functional Imaging and Monitoring, G. Müller (ed.), SPIE Optical Engineering Press, Vol. IS11 (1993), pp. 65–86.

[17] A. Kirsch, The domain derivative and two applications in inverse scattering theory, Inverse Problems, 9 (1993), pp. 81–96.

[18] R. Kress and L. Päivärinta, On the far field operator in obstacle scattering, SIAM J. Appl. Math., 59 (1999), pp. 1413–26.

[19] E. W. Larsen, Solution of Neutron Transport Problems in L_{1}, Comm. Pure Appl. Math., XXVIII (1975), pp. 729–746.

[20] R. M. Lewis and W. W. Symes, On the relation between the velocity coefficient and boundary value for solutions of the one-dimensional wave equation, Inverse Problems, 7 (1991), pp. 597–631.

[21] J. L. Lions and E. Magenes, Non-homogeneous Boundary Value Problems and Applications, Vol. I, Springer, New York, 1972.

[22] P. M. Morse and K. U. Ingard, Theoretical Acoustics, Princeton University Press, 1986.

[23] W. Munk and P. Worchester, Ocean acoustic tomography, Cambridge University Press, 1995.

[24] F. Natterer, Numerical Solution of Bilinear Inverse Problems, Preprints Angewandte Mathematik und Informatik, 19/96-N, Münster, 1996.

Available at www.math.uni-muenster.de/num/Preprints/1997/natterer 3 .

[25] F. Natterer and F. Wübbeling, A propagation-backpropagation method for ultrasound tomography, Inverse Problems, 11 (1995), pp. 1225–1232.

[26] R. Potthast, Fréchet differentiability of boundary integral operators, Inverse Problems, 10 (1994), pp. 431–447.

[27] R. Potthast, Electromagnetic scattering from an orthotropic medium, J. Int. Eq. Appl., 11 (2) (1999), pp. 197–215.

[28] R. Potthast, Integral equation methods in electromagnetic scattering from anisotropic media, Mathematical Methods in the Applied Sciences, 23 (2000), pp. 1145–1159.

[29] H. Sielschott, Measurement of horizontal flow in a large scale furnace using acoustic vector tomography, Flow Meas. Instrum.,8 (1997), pp. 191–197.

[30] H. Sielschott, Rückpropagationsverfahren für die Wellengleichung in bewegtem Medium, Ph.D. thesis, Preprints Angewandte Mathematik und Informatik, 15/00-N, Münster, 2000.

Available at www.math.uni-muenster.de/num/Preprints/2000/sielsch .

[31] J. Sylvester and G. Uhlmann, Inverse Boundary Value Problems at the Boundary — Continuous Dependence, Comm. Pure Appl. Math., 41 (1988), pp. 197–219.