CUBO A Mathematical Journal Vol.10, N o ¯ 02, (15–30). July 2008 Dynamical Inverse Problem for the Equation utt − ∆u − ∇ ln ρ · ∇u = 0 (the BC Method) M.I. Belishev Saint-Petersburg Department of the Steklov Mathematical Institute (POMI), 27 Fontanka, St. Petersburg 191023, Russia, email: belishev@pdmi.ras.ru ABSTRACT A dynamical system of the form utt − ∆u − ∇ ln ρ · ∇u = 0, in R n + × (0, T ) u|t=0 = ut|t=0 = 0, in R n + uxn = f on ∂R n + × [0, T ], is considered, where R n + := { x = {x1, . . . , xn}| xn > 0 } ; ρ = ρ(x) is a smooth positive function (density) such that ρ, 1 ρ are bounded in R n +; f is a (Neumann) boundary control of the class L2(∂R n + × [0, T ]); u = u f (x, t) is a solution (wave). With the system one associates a response operator RT : f 7→ uf |∂Rn + ×[0,T ]. A dynamical inverse problem is to determine the density from the given response operator. Fix an open subset σ ⊂ ∂Rn+; let L2(σ×[0, T ]) be the subspace of controls supported on σ. A partial response operator RT σ acts in this subspace by the rule RT σ f = uf |σ×[0,T ]; let R2T σ be the operator corresponding to the same system considered on the doubled time interval [0, 2T ]. Denote BT σ := { x ∈ Rn+| {x 1, . . . , xn−1, 0} ∈ σ, 0 < xn < T } and assume ρ|σ to be known. We show that R 2T σ determines ρ|BT σ and propose an efficient 16 M.I. Belishev CUBO 10, 2 (2008) procedure recovering the density. The procedure is available for constructing numerical algorithms. The instrument for solving the problem is the boundary control method which is an approach to inverse problems based on their relations with control theory (Belishev, 1986). Our presentation is elementary and can serve as introduction to the BC method. RESUMEN Consideramos el sistema dinámico utt − ∆u − ∇ ln ρ · ∇u = 0, en R n + × (0, T ) u|t=0 = ut|t=0 = 0, en R n + uxn = f sobre ∂R n + × [0, T ], donde R n + := { x = {x1, . . . , xn}| xn > 0 } ; ρ = ρ(x) es una función positiva suave (densidad) tal que ρ, 1 ρ son limitada en R n +; f es un control en la frontera (Neumann) de clase L2(∂R n + × [0, T ]); u = u f (x, t) es la solución (Onda). Con el sistema asociamos un operador respuesta RT : f 7→ uf |∂Rn + ×[0,T ]. Un problema dinámico inverso consiste en determinar la densidad desde el operador respuesta. Fije un subconjunto abierto σ ⊂ ∂Rn+; sea L2(σ × [0, T ]) el subespacio de los controles soportados em σ. Un operador respuesta parcial RT σ actua en este subespacio mediante la regla RT σ f = uf |σ×[0,T ]; sea R 2T σ el operador correspondiente al mismo sistema consid- erado en el intervalo de tiempo [0, 2T ]. Denote BT σ := {x ∈ Rn+| {x 1, . . . , xn−1, 0} ∈ σ, 0 < xn < T } y suponga que ρ|σ es conocido. Nosotros mostramos que R 2T σ determina ρ| B T σ y es propuesto un procedimento eficiente de recuperar la densidad. El proced- imiento es encontrado por construción de algoritmos númericos. El instrumento de resolver el problema es el método de control en la frontera este es un abordage para problemas inversos basado en sus relaciones con teoria de control (Belishev, 1986). Nuestra presentación es elemental y puede servir como introducción al método BC. Key words and phrases: Dynamical inverse problem, response operator, determination of den- sity, boundary control method. Math. Subj. Class.: 35Bxx, 35R30 CUBO 10, 2 (2008) Dynamical Inverse Problem for the Equation ... 17 1 About the paper The problem under consideration comes from geophysics. We deal with a dynamical system of the form utt − ∆u − ∇ ln ρ · ∇u = 0 in R n + × (0, T ) (1) u|t=0 = ut|t=0 = 0 in R n + (2) uxn = f on Γ × [0, T ], (3) where R n + := { x = {x1, . . . , xn}| xn > 0 } simulates the Earth, Γ := ∂Rn+ is the Earth surface, ρ = ρ(x) is a smooth function (density) satisfying 0 < ρ∗ ≤ ρ(·) ≤ ρ ∗ , f is a (Neumann) boundary control of the class L2(Γ×[0, T ]), u = u f (x, t) is a solution. The solution describes a wave initiated by the control and propagating into the Earth from the surface. With the system one associates a response operator RT : f 7→ uf |∂Rn + ×[0,T ]. From the physical viewpoint, f is a force applied at the surface, whereas uf |∂Rn + ×[0,T ] is a displacement measured at the same surface. Thus, R T is an ”input 7→ output” map representing the measurements, which the external observer implements at the surface. The dynamical inverse problem, which the paper deals with, is to determine the density from the given response operator. 2 Results Begin with certain of the notations. With a point x ∈ Rn+ we associate a pair (γ, τ ) : γ := {x1, . . . , xn−1, 0} ∈ Γ, τ := xn ≥ 0 of its semigeodesic coordinates and write x = x(γ, τ ). Fix ξ > 0, let σ ⊂ Γ be an open subset at the surface; the set Bξ σ := {x ∈ Rn+| x = x(γ, τ ), γ ∈ σ, τ ∈ (0, ξ)} is called a tube with a bottom σ and a top σξ = {x(γ, τ )| γ ∈ σ, τ = ξ}. Also, introduce a subdomain Ω ξ σ := {x ∈ Rn+| dist (x, σ) < ξ} 1 (contoured with cdef on Fig 1). The tube (shadowed on Fig 1) is the part of the subdomain illuminated with rays emanating from σ in normal direction to the surface. Figure 1: The tube Bξ σ and the subdomain Ωξ σ 1dist is the standard Euclidean distance in Rn + 18 M.I. Belishev CUBO 10, 2 (2008) Let σ ⊂ Γ and T > 0 be fixed. Consider system (1)–(3) with the final moment t = 2T and introduce a partial response operator R2T σ acting in the (sub)space L2(σ × [0, 2T ]) of controls supported on σ by the rule R2T σ f := uf |σ×[0,2T ] 2 . As is well-known, the waves in system (1)–(3) propagate with the unit velocity. By this, the response operator depends on the density locally: R2T σ is determined by the behavior of ρ in the subdomain ΩT σ only 3 . Such a locality motivates the setup of the inverse problem: given operator R2T σ to recover ρ|ΩT σ . However, since the substitution ρ → cρ with a constant c > 0 does not change system (1)–(2), the unique determination of density is impossible. Avoiding such a nonuniqueness, it is natural to assume the boundary values of ρ to be known. The main result of the paper is Theorem 1 Let T > 0 be fixed, the operator R2T σ given, and the function ρ|σ known. Then these data uniquely determine ρ|BT σ . The proof is constructive: we propose an efficient procedure recovering the density in the tube. Moreover, the procedure is provided with an additional option that is visualization of waves: given f we recover uf | B T σ . 3 Motivation and comments There are two reasons to deal namely with the version (1) of the general wave equation with variable coefficients. First, the interest is motivated by possible applications in geophysics (see, e.g., [6]). The second reason is the following. The instrument for solving the problem is the boundary control method (BCm), which is an approach to inverse problems based on their relations with control theory (Belishev, 1986). In comparison with another versions (see [1] – [4]), the variant of the BCm available for equation (1) is the simplest one. As such, it has good chances for numerical realization. Our presentation is elementary: along with the paper [3], this one can serve as an introduction to the BCm. The plan of the paper is as follows: • in section 4, the basic notions and objects (spaces, operators etc) are introduced • in section 5, we present the so-called amplitude formula (AF) which solves the inverse prob- lem: it recovers ρ|BT σ via R2T σ • sections 6–10 are devoted to the derivation of the AF • in section 11, a certain additional option of the BCm is described: the AF enables one to recover the solutions uf in the tube BT σ that is what we call a wave visualization • section 12 contains the concluding remarks; also, the extension of our results is shortly discussed. 2this operator represents the measurements implemented at the part of the Earth surface 3in other words, R2Tσ does not depend on ρ|Rn + \ΩT σ CUBO 10, 2 (2008) Dynamical Inverse Problem for the Equation ... 19 Simplifying the notations, we accept the convention: unless the otherwise is specified, the subset σ ⊂ Γ is assumed fixed and we write R2T instead of R2T σ , Ω ξ instead of Ω ξ σ , Bξ instead of Bξ σ , etc. Also, without loss of generality, we assume σ to be bounded and ∂σ smooth. 4 Dynamical system With system (1)–(3) one associates (i) an outer space FT := {f ∈ L2, ρ0 (Γ × [0, T ])| supp f ⊂ σ × [0, T ]} of controls acting from σ with the inner product (f, g) F T = ∫ σ×[0,T ] f (γ, t) g(γ, t) ρ0(γ)dΓdt , where ρ0 := ρ|Γ, dΓ is the Euclidean surface element on ∂R n +. In F T we single out a family of subspaces F T, ξ := { f ∈ FT | supp f ⊂ σ × [T − ξ, T ] } , ξ ∈ [0, T ] consisting of the delayed controls (T − ξ is the value of delay, ξ is the action time; FT, 0 = {0}, FT, T = FT ); (ii) an inner space of states (waves) HT := L2, ρ(Ω T ) with the product (y, w) H T = ∫ ΩT y(x) w(x) ρ(x)dx and the family of its subspaces 4 H ξ := { y ∈ HT | supp y ⊂ Ωξ } , ξ ∈ [0, T ] . Since the waves described by a hyperbolic system (1)–(3) propagate (from σ into Rn+) with the speed 1, the inclusion supp uf (·, t) ⊂ Ωt holds for all t ∈ [0, T ], whereas ΩT is the subdomain filled with waves at the final moment t = T . Correspondingly, we consider the waves as time depended elements of the space HT ; (iii) a control operator W T : FT → FT , W T f := uf (·, T ). This operator is continuous [7] and injective for any T > 0 [1]. By the above-mentioned hyperbolicity, for f ∈ FT, ξ one has supp uf (·, T ) ⊂ Ωξ that yields the embedding W T FT, ξ ⊂ Hξ; (iv) a response operator RT :FT → FT , RT f := uf |σ×[0,T ], which is also a continuous map 5 ; (v) a connecting operator CT : FT → FT , CT := (W T )∗W T . For f, g ∈ FT , one has ( uf (·, T ), ug(·, T ) ) H T = (W T f, W T g) H T = (CT f, g) F T , (4) 4recall that Ωξ = {x ∈ Rn + | dist (x, σ) < ξ} 5moreover, RT is a compact operator [7] 20 M.I. Belishev CUBO 10, 2 (2008) so that CT connects the metrics of the outer and inner spaces. By injectivity of W T , the operator CT is also injective. One of central points of the BCm is an explicit relation between the response and connecting operators. Denote by ST : FT → F2T the operator extending controls from σ × [0, T ] to σ × [0, 2T ] as odd (with respect to t = T ) functions of time; let J 2T : F2T → F2T be the integration: (J 2T f )(·, t) := ∫ t 0 f (·, s) ds; R2T : F2T → F2T the response operator of system (1)–(3) with the final moment t = 2T . Lemma 1 The relation CT = 1 2 (ST )∗J 2T R2T ST (5) holds. Proof Choose f, g ∈ FT and denote f− := S T f . Assume f , g to be such that uf− , ug are the classical solutions. Blagoveschenskii’s function b(s, t) := ∫ R n + uf− (·, s) ug(·, t) ρdx s, t ∈ [0, 2T ] × [0, T ] is well defined and satisfies btt(s, t) − bss(s, t) = ∫ Ω [ u f − (·, s) u g tt (·, t) − u f − ss (·, s) u g (·, t) ] ρdx = ∫ Ω [ uf− (·, s) divρ∇ug(·, t) − divρ∇uf− (·, s) ug(·, t) ] dx = ∫ Γ [ −uf− (·, s) u g xn (·, t) + u f − xn (·, s) ug(·, t) ] ρ0dΓ = − ∫ Γ [ (R2T f−)(·, s) g(·, t) − f−(·, s) (R T g)(·, t) ] ρ0dΓ =: F (s, t) (in the second equality we have used equation (1) in the form ρutt = div(ρ∇u). Finding b by the D’Alembert formula (with regard to the initial conditions b(·, 0) = bt(·, 0) = 0), putting t = T , and taking into account the oddness of f−, we get b(T, T ) = 1 2 ∫ T 0 dt ∫ 2T −t t F (s, t) ds = − 1 2 ∫ T 0 dt ∫ 2T −t t ds ∫ Γ (R 2T f−)(·, s) g(·, t) ρ0dΓ := ∫ Γ×[0,T ] 1 2 {∫ t 0 (R 2T f−)(·, s) ds − ∫ 2T −t 0 (R 2T f−)(·, s) ds } g(·, t) ρ0dΓdt that can be easily transformed to b(T, T ) = ( 1 2 (ST )∗J 2T R2T ST f, g ) F T . (6) CUBO 10, 2 (2008) Dynamical Inverse Problem for the Equation ... 21 On the other hand, we have b(T, T ) = ( u f − (·, T ), u g (·, T ) ) H T = ( u f (·, T ), u g (·, T ) ) H T = 〈see (4)〉 = (C T f, g) F T . (7) Comparing (6) with (7) and taking into account the density (in FT ) of f, g used, we arrive at (5). � 5 Amplitude formula Here we present a relation that determines the density from the response operator. Recall that x(γ, τ ) denotes the point in Rn+ with the semigeodesic coordinates γ and τ , the subset σ ∋ γ is fixed, ρ0 = ρ|Γ. Fix ξ ∈ (0, T ); let f̌ ξ := {f ξ k } ∞ k=1 ⊂ F T, ξ be a linearly independent complete system 6 of delayed controls acting from σ and such that the corresponding solutions uf ξ k are classical. With the system we associate its Gram matrix {Gik} ∞ i,k=1 : Gik := (C T f ξ i , f ξ k ) F T = ∫ σ×[T −ξ,T ] (CT f ξ i )(γ, t) f ξ k (γ, t) ρ0(γ)dΓdt (8) and a sequence of numbers {βi} ∞ i=1 : βi := − (κ T , f ξ i ) F T = − ∫ σ×[T −ξ,T ] (T − t) f ξ i (γ, t) ρ0(γ)dΓdt , (9) where κ T = κ T (γ, t) := T − t. For any integer N ≥ 1, a linear algebraic system N∑ k=1 Gik α N k = βi, i = 1, . . . , N (10) is uniquely solvable (w.r.t. αN1 , . . . , α N N ) by injectivity of the operator CT . Lemma 2 For a fixed (γ, ξ) ∈ σ × (0, T ), the relation ρ (x(γ, ξ)) = ρ0(γ) {[ ∂ ∂t lim N→∞ N∑ k=1 α N k ( C T f ξ k )] (γ, t) ∣∣∣∣ t=T −ξ+0 t=T −ξ−0 }2 (11) holds. Relation (11) is a relevant version of the so-called amplitude formula (AF), which is one of the main tools for solving inverse problems by the BCm: see [1], [2], [3]. The assertion of Theorem 1 easily follows from the AF. Indeed, if ρ0 is known and R 2T is given, then one can find CT by (5) and, choosing a system f̌ ξ, determine the r.h.s. of (11). As (γ, ξ) runs over σ × (0, T ), the points x(γ, ξ) exhaust the tube BT . Hence, ρ|BT is determined by R 2T . In sections 6–10 we derive the AF. 6the completeness means that closFT span f̌ ξ = FT, ξ 22 M.I. Belishev CUBO 10, 2 (2008) 6 Propagation of jumps The well-known and typical for hyperbolic problems fact is that discontinuity of a control f in system (1)–(3) implies discontinuity of the corresponding wave uf , the latter discontinuity (singu- larity) propagating along the space-time rays and being supported on the characteristic surfaces. Below we recall certain details. Let θ0(t) := 1 2 [1 + sign t], t ∈ R be the Heavyside function; the sequence {θk(·)} k=∞ k=−∞ : dθk dt = θk−1 is usually referred to as a smoothness scale 7 . Choose a smooth control f ∈ FT and fix ξ ∈ (0, T ); by fξ(γ, t) := f (γ, t)θ0 (t − (T − ξ)) we denote its cut-off function on σ × [T − ξ, T ], so that fξ is an element of the subspace F T, ξ . In the generic case, the control fξ has a jump at the cross-section σ × {t = T − ξ} (see ab on Fig 2), the value (amplitude) of the jump being equal to fξ( · , t)| t=T −ξ+0 t=T −ξ−0 = f ( · , T − ξ) − 0 = f ( · , T − ξ). t=0 t=T t=T-x Figure 2: Jumps in system (1)–(3) Jumps of a Neumann control induce jumps of a wave velocity. Namely, the velocity u fξ t turns out to be discontinuous; its jump is supported on the characteristic surface {(x, t)| t = (T − ξ)+ dist (x, σ)}. This surface consists of the plane part abde and two conic parts bcd and aef (see Fig 2; the arrows pick up the space-time rays, which the discontinuity propagates along). The jump on the conic parts is weaker than the one on the plane part and plays no role in the further considerations. The jump on the plane part is described as follows. In a (space-time) neighborhood of abde, the solution is sought in the form ufξ = Ap + wp , Ap(x, t) := p∑ k=1 Ak(γ, τ ) θk (t − (T − ξ) − τ ) , (12) where x in the l.h.s. is x(γ, τ ); Ap is an ansatz, which is a function of the class C p−1 loc ; Ak are the so- called amplitude functions; wp ∈ C p loc is a smoother reminder. Substituting such a representation in (1)–(3), one derives a recurrent system of ODE’s 8 for the amplitude functions. As result, one 7θk with negative k’s are understood in the sense of distributions 8the so-called transport equations CUBO 10, 2 (2008) Dynamical Inverse Problem for the Equation ... 23 arrives at the representation (for p = 1) of the form ufξ (x, t) = − ( ρ (x(γ, τ )) ρ0(γ, τ ) ) − 1 2 f (γ, T − ξ) θ1 (t − (T − ξ) − τ ) + w(x, t) θ2 (t − (T − ξ) − τ ) , (13) where x = x(γ, τ ) in the l.h.s. and w is a smooth function. For details of this technique see, e.g., [5]. At the final moment t = T , the wave ufξ (· , T ) is supported in the subdomain Ωξ contoured with cdef on Fig 2. The surface cdef coincides with the forward front of the wave. By (13), in a neighborhood of the plane part ed of the front, the representation u fξ (x, T ) = (W T fξ)(x) = − ( ρ (x(γ, τ )) ρ0(γ, τ ) ) − 1 2 f (γ, T − ξ) θ1 (ξ − τ ) + w(x) θ2(ξ − τ ) (14) holds with a smooth w. Correspondingly, the velocity of the wave has a jump at ed: u fξ t (x(γ, τ ), T ) ∣∣∣∣ τ =ξ+0 τ =ξ−0 = 0 − u fξ t (x(γ, ξ − 0)) = 〈see (13)〉 = − ( ρ (x(γ, ξ)) ρ0(γ, ξ) ) − 1 2 f (γ, T − ξ) , γ ∈ σ . (15) So, up to the factor −( ρ ρ0 ) − 1 2 , the shape of the velocity jump reproduces the shape of the control jump. 7 Dual system A dynamical system vtt − ∆v − ∇ ln ρ · ∇v = 0, in R n + × (0, T ) (16) v|t=0 = 0, vt|t=0 = y, in R n + (17) vxn = 0 on Γ × [0, T ], (18) is said to be dual to system (1)–(3); its solution v = vy(x, t) describes a wave initiated by a velocity perturbation y and propagating (in the reversed time) into Rn+. The term ’dual’ is motivated by the following relation between solutions of these systems. Lemma 3 For any square summable f and y compactly supported in Γ×[0, T ] and Rn+ respectively, the equality ∫ R n + u f (· , T ) y ρdx = ∫ Γ×[0,T ] f v y ρ0dΓdt (19) is valid. 24 M.I. Belishev CUBO 10, 2 (2008) Proof Let f and y be compactly supported and such that the solutions uf and vy are classical. The relations 0 = ∫ R n + ×[0,T ] [ u f tt − ∆uf − ∇ρ · ∇uf ] vy ρdx dt = ∫ R n + ×[0,T ] [ ρu f tt − div ρ∇u f ] v y ρdx dt = ∫ R n + [ u f t vy − uf v y t ] ∣∣∣∣ t=T t=0 ρdx + ∫ T 0 dt ∫ Γ [ u f xn vy − uf v y xn ] ρ0dΓ + ∫ R n + ×[0,T ] u f [v y tt − ∆v y − ∇ρ · ∇v y ] ρdx dt = 〈 see (2), (17) 〉 = − ∫ R n + uf (·, T ) y ρdx + ∫ Γ×[0,T ] f vy ρ0dΓdt imply (19). By the density of the chosen f ’s and y’s in the corresponding L2- spaces, the passage to the limit in the proper sense leads to the assertion of the lemma. � Taking in (17) y ∈ HT , introduce an observation operator OT : HT → FT , OT y := vy|σ×[0,T ]. For f ∈ FT , relation (19) can be written in the form (W T f, y) H T = (f, OT y) F T , which yields OT = ( W T ) ∗ (20) and clarifies the duality of the systems. 8 Jumps in dual system Here we consider the dual system provided with the specific Cauchy data (17): the velocity per- turbation y is assumed to be discontinuous. Recall the notations: Ωξ = {x ∈ Rn+| dist (x, σ) < ξ}, Bξ = {x ∈ Rn+| x = x(γ, τ ), γ ∈ σ, τ ∈ (0, ξ)}, ξ > 0. Fix ξ ∈ (0, T ) and take a function y ∈ C∞(ΩT ) ⊂ HT ; denote by yξ := { y in Ωξ 0 in ΩT \ Ωξ its cut-off function onto the subdomain Ω ξ . In the generic case, the function yξ has a jump at a surface ∂Ωξ ∩ Rn+ (see cdef on Fig 3). CUBO 10, 2 (2008) Dynamical Inverse Problem for the Equation ... 25 t=0 t=T t=T-x Figure 3: Jumps in the dual system (16)-(18) Return to the dual system and replace y by yξ in (17). The well-known fact is that discontin- uous data produce discontinuous waves. Namely, the jump of the data at cdef implies the jump of the wave velocity, the latter one being supported on the characteristic surfaces. In particular, these surfaces contain the plane parts abde and ednm consisting of the space-time rays (see the arrows) emanating from the set σξ × {t = T } (see ed). The amplitude of the jumps at these parts can be found by standard geometrical optics technique, i.e., by the use of the relevant analog of representation (12) for the solution vyξ . Omitting the details, the result is as follows. For a point (x0, t0) ∈ abde such that x0 = x(γ, τ ), γ ∈ σ, 0 < τ < ξ and t0 = T − (τ − ξ), one has v yξ t (x0, t) ∣∣∣∣ t=t0+0 t=t0−0 = 1 2 ( ρ (x(γ, ξ)) ρ (x(γ, τ )) ) 1 2 y (x(γ, ξ)) . (21) The meaning of this formula is quite transparent: the jump of data at the point x(γ, ξ) initiates the jump of velocity, which propagates (in the reversed time) along the ray {(x(γ, τ ), t) | t = T − (τ − ξ), τ ∈ [0, ξ]}, its amplitude being proportional to the jump of data up to the factor 1 2 ( ρ ρ ) 1 2 . 9 Further, at the moment t = T − ξ the jump reaches the boundary Γ and is reflected back into R n + (see the rays constituting the part ablk). As result, a trace of the velocity v yξ t on σ × [0, T ] turns out to be discontinuous, its jump being supported at the cross-section σ × {t = T − ξ} (see ab). The amplitude of this jump can be found from (21): the equality v yξ t (γ, t) ∣∣∣∣ t=T −ξ+0 t=T −ξ−0 = ( ρ (x(γ, ξ)) ρ (x(γ, 0)) ) 1 2 y (x(γ, ξ)) 9the part ednm also supports the jump moving into Rn + in opposite direction. This jump plays no role in the further considerations. 26 M.I. Belishev CUBO 10, 2 (2008) holds 10 . By the definition of the observation operator, this equality can be written as ( ∂ ∂t OT yξ ) (γ, t) ∣∣∣∣ t=T −ξ+0 t=T −ξ−0 = ( ρ (x(γ, ξ)) ρ0(γ) ) 1 2 y (x(γ, ξ)) . (22) At last, putting y = 1 (so that 1ξ is an indicator of the subdomain Ω ξ ), we obtain the important auxiliary relation ( ρ (x(γ, ξ)) ρ0(γ) ) 1 2 = ( ∂ ∂t OT 1ξ ) (γ, t) ∣∣∣∣ t=T −ξ+0 t=T −ξ−0 , (γ, ξ) ∈ σ × (0, T ) . (23) Completing the derivation of (11) in the next two sections, we show how to express the r.h.s. of (23) through the inverse data (operator R2T ). 9 Controllability. Wave basis Fix ξ ∈ (0, T ). In control theory, the set Uξ := Ran W T = {uf (·, T )| f ∈ FT, ξ} of waves produced by controls acting from σ (the action time is ξ) is said to be reachable (at the moment t = ξ). Since the wave propagation speed is equal to 1, the waves constituting U ξ are supported in the metric neighborhood Ω ξ of σ; as result, the embedding Uξ ⊂ Hξ holds. The property of system (1)–(3), which plays the key role in solving inverse problems, is that for any σ and ξ this embedding is dense: clos Uξ = Hξ. Control theory interprets this fact as a local approximate boundary controllability 11 of the system. Roughly speaking, it means that the reachable set is rich enough: any function supported in the subdomain Ω ξ filled with waves, can be approximated (in L2- metric) by a wave uf (·, T ) ∈ Uξ with the properly chosen control f ∈ FT, ξ. The proof of this property relays on the fundamental Holmgren–John–Tataru uniqueness theorem (see [1] for detail). Let f̌ ξ = {f ξ k }∞ k=1 ⊂ F T, ξ be a linearly independent complete system of controls acting from σ; denote by ǔξ = {u ξ k }∞ k=1 ⊂ H ξ, u ξ k := uf ξ k (·, T ) = W T f ξ k a system of the corresponding waves. By the controllability, the latter system is also complete: clos span ǔξ = Hξ. With a slight abuse of terms, we call ǔξ a wave basis. Denote by P ξ N the orthogonal projection in HT onto span {u ξ k }N k=1. Since the waves form a complete system, one has s-limN→∞ P ξ N = P ξ, where P ξ projects in HT onto Hξ, i.e., cuts off functions supported in Ω T onto Ω ξ . Recall that 1ξ denotes the indicator of Ω ξ . Representing 1ξ = P ξ 1T = lim N→∞ P ξ N 1T , P ξ N 1T = N∑ k=1 αN k u ξ k , (24) 10this equality can be also derived from (14), (15) by the use of duality relation (19). The factor 1 2 is doubled owing to the contribution of the reflected rays 11this motivates the name ”boundary control method” CUBO 10, 2 (2008) Dynamical Inverse Problem for the Equation ... 27 one can find the coefficients α ξ N as follows. Multiplying the last equality in (24) by u ξ i , one gets N∑ k=1 Gik α N k = βi, i = 1, . . . , N , (25) where Gik := (u ξ i , u ξ k ) H T = (W T f ξ i , W T f ξ k ) H T = 〈see (4)〉 = (CT f ξ i , f ξ k ) F T = ∫ σ×[T −ξ,T ] (CT f ξ i )(γ, t) f ξ k (γ, t) ρ0(γ)dΓdt (26) and βi := (P ξ N 1T , u ξ i ) H T = (1T , P ξ N u ξ i ) H T = (1T , u ξ i ) H T = ∫ ΩT uf ξ i (·, T ) ρdx = 〈see (2)〉 = ∫ ΩT ρdx ∫ T 0 (T − t) u f ξ i tt (·, T ) dt = ∫ T 0 dt (T − t) ∫ ΩT divρ∇uf ξ i (·, T ) dx = − ∫ Γ×[0,T ] (T − t) u f ξ i xn (·, T ) ρ0dΓdt = 〈see (3)〉 = − ∫ σ×[T −ξ,T ] (T − t) f ξ i (γ, t) ρ0dΓdt = − (κ T , f ξ i ) F T (27) with κ T (γ, t) := T − t. In the course of integration by parts, the integral over ∂ΩT ∩ Rn+ vanishes because the wave uf ξ i (·, T ) is supported into the smaller subdomain Ωξ. Also, note that the linear independence of the system f̌ ξ and the injectivity of the of the operator W T provide the linear independence of the system ǔξ in Hξ. By the latter, the Gram matrix G is invertible and, hence, system (25) is uniquely solvable w.r.t. αN1 , . . . , α N N . 10 Completing the proof. Applying OT to the both sides of (24), with regard to u ξ k = W T f ξ k , we have O T 1ξ = lim N→∞ N∑ k=1 α N k O T W T f ξ k = 〈see (20)〉 = lim N→∞ N∑ k=1 α N k C T f ξ k . (28) At last, substituting (28) into (23), we arrive at (11). The amplitude formula is derived and Theorem 1 is proven. � The obtained results can be presented in the form of a procedure recovering the density. If the values of ρ0|σ are known and the response operator R 2T is given (as result of measurements on the part σ of the Earth surface), the external observer can determine the density ρ in the tube BT as follows: Step 1 Find the operator CT by (5). Fix ξ ∈ (0, T ) and choose a complete system of controls f̌ ξ ⊂ FT, ξ. Find the Gram matrix {Gik} and the numbers βi by (26), (27). 28 M.I. Belishev CUBO 10, 2 (2008) Step 2 Solve system (25) and find αN1 , . . . , α N N for N = 1, 2, . . . . Fixing a γ ∈ σ, determine ρ (x(γ, ξ)) by the AF (11). Step 3 Varying (γ, ξ) ∈ σ × (0, T ), recover ρ| B T . 11 Visualization of waves The amplitude formula enables one to recover not only the density but the waves themselves. Choose a control f ∈ FT providing the wave uf to be a classical solution of (1)–(3). For a fixed (γ, ξ) ∈ σ × (0, T ), the representation u f (x(γ, ξ), T ) = ( ρ (x(γ, ξ)) ρ0(γ) ) − 1 2 {[ ∂ ∂t lim N→∞ N∑ k=1 η N k C T f ξ k ] (γ, t) ∣∣∣∣ t=T −ξ+0 t=T −ξ−0 } (29) holds, where the coefficients ηN1 , . . . , η N N satisfy the system N∑ k=1 Gik η N k = θi, i = 1, . . . , N with θi := (C T f, f ξ i ) F T . Indeed, putting y = uf (·, T ) in (22), we have ( ρ (x(γ, ξ)) ρ0(γ) ) 1 2 uf (x(γ, ξ), T ) = [ ∂ ∂t OT ( uf (·, T ) ) ξ ] (γ, t) ∣∣∣∣ t=T −ξ+0 t=T −ξ−0 . (30) Using the wave basis by analogy with (24), we can represent ( uf (·, T ) ) ξ = P ξuf (·, T ) = lim N→∞ P ξ N uf (·, T ) , P ξ N uf (·, T ) = N∑ k=1 ηN k u ξ k (31) and find the coefficients ηN1 , . . . , η N N from the Gram system N∑ k=1 Gikη N k = θi , i = 1, . . . , N with the r.h.s. θi := ( uf (·, T ), u ξ i ) H T = 〈see (4)〉 = (CT f, f ξ i ) F T . Then, substituting OT ( uf (·, T ) ) ξ = 〈see (31)〉 = lim N→∞ N∑ k=1 ηN k OT u ξ k = lim N→∞ N∑ k=1 ηN k CT f ξ k CUBO 10, 2 (2008) Dynamical Inverse Problem for the Equation ... 29 in (30), we arrive at (29). Thus, the external observer with the knowledge of the response operator R2T can make the waves visible in the tube BT under the part σ of the Earth surface. This is what we call a visualization (see [1]–[3]). 12 Comments • As regards the numerical realization of the procedure outlined in section 10, the most prob- lematic point is Step 2, which consists of solving system (25) for large N ’s. Namely, since CT is a compact operator, the condition number of the matrix {Gik} grows as N → ∞, so that (25) turns out to be an ill-posed system. The second problem is that the passage to the limit and the differentiation w.r.t. to t in (11) do not commute. Actually, the difficulties of this type are unavoidable: they reflect the well-known strong ill-posedness of multidimensional inverse problems. However, certain affirmative results in numerical testing do exist and show that elaboration of workable BC-algorithms is not a hopeless endeavor [2], [6]. • For numerical realization, an important issue is the proper choice of the system f̌ ξ. For many reasons, the controls producing the waves with sharp forward front are preferable. In the problem with Neumann boundary controls (2), it looks reasonable to simulate a complete system by the set {f ξ qp } M,N p=1, q=1: f ξ qp (γ, t) = δ′ ε ( t − [ (T − ξ) + p ξ M + 1 ]) fq(γ) , where δ′ ε (t) is a relevant regularization of the first derivative of the Dirac function δ(t), {fq} ∞ q=1 is an orthonormal basis in L2(σ), M and N are large integers. Also, note that in applications one deals as a rule with a certain prescribed set of ’standard’ controls. Therefore, it is important to develop the algorithms with controls simulating the real sources. • The amplitude formulae (11) and (29) can be easily extended to the case of a curved boundary Γ. Namely, assume that σ ⊂ Γ and T > 0 are such that the field of the normal rays, which form the tube BT , is regular. Then, the only correction required is to replace ρ0 by ρ0 J , where J = J(γ, τ ) is the Jacobian of the passage from the semideodesic coordinates in BT σ to the Cartesian coordinates (see [3]). 13 Acknowledgements I would like to thank Prof. C. Cuevas for kind invitation to write this paper for the Journal. I’m grateful to S.V.Belisheva and I.V.Kubyshkin for assistance in computer graphics. The work is supported by the RFBR grants No. 08-01-00511 and the Project NSh-8336.2006.1. Received: November 2007. Revised: January 2008. 30 M.I. Belishev CUBO 10, 2 (2008) References [1] M.I. Belishev, Boundary control in reconstruction of manifolds and metrics (the BC-method). Invers Problems., 13 (1997), No 5, R1–R45. [2] M.I. Belishev, V.Yu. Gotlib, Dynamical variant of the BC-method: theory and numerical testing. Journal of Inverse and Ill-Posed Problems, 7, no 3: 221–240, 1999. [3] M.I. Belishev, How to see waves under the Earth surface (the BC-method for geophysicists). Ill-Posed and Inverse Problems, S.I.Kabanikhin and V.G.Romanov (Eds). VSP, 55–72, 2002. [4] M.I. Belishev, Recent progress in the boundary control method. Invers Problems., 23 (2007), No 5, R1–R67. [5] M. Ikawa. Hyperbolic PDEs and Wave Phenomena, Translations of Mathematical Mono- graphs, v. 189 AMS; Providence. Rhode Island, 1997. [6] S.I. Kabanikhin, M.A. Shishlenin, A.D. Satybaev, Direct Methods of Solving Inverse Hyperbolic Problems. Utrecht, The Netherlands, VSP, 2004. [7] I. Lasiecka, J-L. Lions, R. Triggiani, Non homogeneous boundary value problems for second order hyperbolic operators. J. Math. Pures Appl, v. 65 (1986), no 3, 142–192. N2