The quasi-reversibility method to numerically solve an inverse source problem for hyperbolic equations
TThe quasi-reversibility method to numerically solve an inverse sourceproblem for hyperbolic equations
Thuy T. Le ∗ Loc H. Nguyen ∗ Thi-Phong Nguyen † William Powell ∗ Abstract
We propose a numerical method to solve an inverse source problem of computing the initial condi-tion of hyperbolic equations from the measurements of Cauchy data. This problem arises in thermo-and photo- acoustic tomography in a bounded cavity, in which the reflection of the wave makes thewidely-used approaches, such as the time reversal method, not applicable. In order to solve this in-verse source problem, we approximate the solution to the hyperbolic equation by its Fourier series withrespect to a special orthonormal basis of L . Then, we derive a coupled system of elliptic equationsfor the corresponding Fourier coefficients. We solve it by the quasi-reversibility method. The desiredinitial condition follows. We rigorously prove the convergence of the quasi-reversibility method as thenoise level tends to 0. Some numerical examples are provided. In addition, we numerically prove thatthe use of the special basic above is significant. Key words: inverse source problem, hyperbolic equation, quasi-reversibility method
AMS subject classification: 35R30, 65M32
We consider an inverse source problem arising from biomedical imaging based on the photo-acousticand thermo-acoustic effects, which are named as photo-acoustic tomography (PAT) and thermo-acoustictomograpy (TAT) respectively. In PAT, [1, 2], non-ionizing laser pulses are sent to a biological tissueunder inspection (for instance, woman’s breast in mamography). A part of the energy will be absorbedand converted into heat, causing a thermal expansion and a subsequence ultrasonic wave propagating inspace. The ultrasonic pressures on a surface around the tissue are measured. The experimental set upfor TAT, [3], is similar to PAT except the use of microwave other than laser pulses. Finding some initialinformation of the pressures from these measurements yields structure inside this tissue.Due to the important real-world applications, the inverse source problem PAT/TAT has been studiedintensively. There are several methods to solve them available. In the case when the waves propagate inthe free space, one can find explicit reconstruction formulas in [4, 5, 6, 7], the time reversal method [8, 9,10, 11, 12], the quasi-reversibility method [13] and the iterative methods [14, 15, 16]. The publicationsabove study PAT/TAT for simple models for non-damping and isotropic media. The reader can findpublications about PAT/TAT for more complicated model involving a damping term or attenuation term[17, 18, 19, 20, 21, 22, 23, 24, 25]. The time reversibility requires an approximation of the wave at a “final ∗ Department of Mathematics and Statistics, University of North Carolina at Charlotte, Charlotte, NC 28223, USA,[email protected], [email protected], [email protected] † Department of Mathematics, Purdue University, West Lafayette, IN, 47907, [email protected] (corresponding au-thor). a r X i v : . [ m a t h . NA ] J a n tage” in the whole medium. This “internal data” might be known assuming that the reflection of thewave on the measured surface is negligible. However, there are many circumstances where this assumptionis no longer true. For example, when the biological tissue under consideration is located inside glass, thewaves reflects as in a resonant cavity [26, 27, 28]. In this case, measuring or approximating final stageof the wave inside the tissue is impossible. We draw the reader’s attention to [29] for a method to solvePAT/TAT in a bounded cavity with wave reflection. Our contribution in this paper is to apply the quasi-reversibility method to solve the inverse source problem of PAT/TAT for damping and nonhomogeneousmedia. In this case, the model is a full hyperbolic equation in a bounded domain. By this, the reflectionof the waves at the measurement sites is allowed.The uniqueness and the stability for the inverse source problem for general hyperbolic equations inthe damping and inhomogeneous media can be proved by using a Carleman estimate. These importantresults are the extensions of the uniqueness and stability for the inverse problem for a simple hyperbolicequation in [13] in the non-damping case. However, since the arguments are very similar to the ones inthat paper using Carleman estimate for hyperbolic operators, for the brevity, we do not write the proofhere.Instead of using the direct optimal control method, to find the initial value of general hyperbolicequations, we derive a system of elliptic partial differential equations, which are considered as an approx-imation model for our method. Solution of this system is the vector consisting of Fourier coefficients ofthe wave with respect to a special basis. This system and the given Cauchy boundary data can be solvedby the quasi-reversibility method. The quasi-reversibility method was first proposed by Latt`es and Lions[30] in 1969. Since then it has been studied intensively [31, 32, 33, 34, 13, 35, 36, 37]. The application ofCarleman estimates for proofs of convergence of those minimizers was first proposed in [38] for Laplace’sequation. In particular, [39] is the first publication where it was proposed to use Carleman estimates toobtain Lipschitz stability of solutions of hyperbolic equations with lateral Cauchy data. We draw thereader’s attention to the paper [40] that represents a survey of the quasi-reversibility method. Using aCarleman estimate, we prove Lipschitz-like convergence rate of regularized solutions generated by thequasi-reversibility method to the true solution of that Cauchy problem.It seems, in theory, that our method of approximation works for any orthonormal basis of L . However,this observation is not true in the numerical sense. That means, the special basis we use to establishthe approximation model is crucial. The basis we use in this paper was first introduced by Klibanov in[48], called { Ψ n } n ≥ . It has a very important property that Ψ (cid:48) n is not identically 0 in an open intervalwhile other bases; for e.g., trigonometric functions and orthonormal polynomials, do not. In this paper,we prove numerically that choosing this basis is optimal for our method.As mentioned, we establish in this paper the Lipschitz convergence of the quasi-reversibility method.Our main contribution to this field is to relax a technical condition on the noise. In our previous works[41, 42, 43, 44, 45] and references therein, we assumed that the noise contained in the boundary datacan be “smoothly extended” as a function defined on the domain Ω. This condition implies that thenoise must be smooth. Motivated by the fact that this assumption is not always true, we employ aCarleman estimate involving the boundary integrals to obtain the new convergence without imposingthis “extension” condition of the noise function.The paper is organized as follows. In Section 2, we state the inverse problems under consideration andderive an approximation model whose solution directly yields their solutions. In Section 3, we introducesome auxiliary results and prove the Carleman estimate, which plays an important role in our analysis.In Section 4, we implement the quasi-reversibility method to solve the system of elliptic equations andprove the convergence of the solution as the noise level tends to 0. Section 5 is for the numerical studies.Section 6 is to provide some numerical results for 1 D − problem and to show the significance of the used2rthonormal basis. Finally, section 7 is for the concluding remarks. Let Ω be a smooth and bounded domain in R d , where d ≥ T bea positive number. Let 1 ≤ a ∈ C (Ω), and a, c ∈ L ∞ (Ω) be functions defined on Ω. Let b be a d -dimensional vector valued function in L ∞ (Ω , R d ) . Define the elliptic operator Lφ := ∆ φ + b · ∇ φ + c φ (2.1)for all φ ∈ C (Ω). Let p ∈ H (Ω) represent a source, we consider the problems of solving u ( x , t ) ∈ H (Ω × (0 , T )) generated by the source p ( x ) ∈ C (Ω) and subjected to either homogeneous Dirichlet orNeumann boundary condition, which are given by a ( x ) u tt ( x , t ) + a ( x ) u t ( x ) = Lu ( x , t ) ( x , t ) ∈ Ω × (0 , T ) u t ( x ,
0) = 0 x ∈ Ω ,u ( x ,
0) = p ( x ) x ∈ Ω ,u ( x , t ) = 0 ( x , t ) ∈ ∂ Ω × [0 , T ] , (2.2)or a ( x ) u tt ( x , t ) + a ( x ) u t ( x ) = Lu ( x , t ) ( x , t ) ∈ Ω × (0 , T ) u t ( x ,
0) = 0 x ∈ Ω ,u ( x ,
0) = p ( x ) x ∈ Ω ,∂ ν u ( x , t ) = 0 ( x , t ) ∈ ∂ Ω × [0 , T ] . (2.3)respectively. Our interest is to determine the source p ( x ), x ∈ Ω, from some boundary observation of thewave. More precisely, the inverse source problems are formulated as:
Problem 2.1.
Determine the function p from the measurement of f ( x , t ) = ∂ ν u ( x , t ) ( x , t ) ∈ ∂ Ω × [0 , T ] (2.4) where u is the solution of (2.2) . Problem 2.2.
Determine the function p from the measurement of f ( x , t ) = u ( x , t ) ( x , t ) ∈ ∂ Ω × [0 , T ] (2.5) where u is the solution of (2.3) . The unique solvability of problems (2.2) and (2.3) can be obtained by Garlerkin approximations andenergy estimates as in Chapter 7, Section 7.2 in [46]. We now focus on our approach for solving thesetwo inverse problems.Let { Ψ n } n ≥ be an orthonormal basis of L (0 , T ) . The function u ( x , t ) can be represented as: u ( x , t ) = ∞ (cid:88) n =1 u n ( x )Ψ n ( t ) , for all ( x , t ) ∈ Ω × [0 , T ] (2.6)where u n ( x ) = (cid:90) T u ( x , t )Ψ n ( t ) dt, n ≥ . (2.7)3onsider u N ( x , t ) := N (cid:88) n =1 u n ( x )Ψ n ( t ) for all ( x , t ) ∈ Ω × [0 , T ] , (2.8)for some cut-off number N . This number N is chosen numerically such that u N well-approximates thefunction u , see Section 5.1 for more details. We have, u Nt ( x , t ) = N (cid:88) n =1 u n ( x )Ψ (cid:48) n ( t ) , and u Ntt ( x , t ) = N (cid:88) n =1 u n ( x )Ψ (cid:48)(cid:48) n ( t ) , for all ( x , t ) ∈ Ω × [0 , T ] . (2.9)Plugging (2.8) and (2.9) into the first equation of problem (2.2), we get N (cid:88) n =1 a ( x ) u n ( x )Ψ (cid:48)(cid:48) n ( t ) + a ( x ) N (cid:88) n =1 u n ( x )Ψ (cid:48) n ( t ) = N (cid:88) n =1 Lu n ( x )Ψ n ( t ) (2.10)for all ( x , t ) ∈ Ω × [0 , T ]. Remark 2.1.
Equation (2.10) is actually an approximation model. We only solve Problem 2.1 andProblem 2.2 in this approximation context. Studying the behavior of (2.10) as N → ∞ is extremelychallenging and out of the scope of the paper. In case of interesting, the reader could follow the techniquesin [47] to investigate the accuracy of (2.11) as N tends to ∞ . Since { Ψ n } n ≥ is an orthonormal basis of L , multiplying Ψ m ( t ) to both sides of (2.10) and thenintegrating the resulting equation with respect to t , for each m ∈ { , . . . , N } , yields N (cid:88) n =1 s mn u n ( x ) = Lu m ( x ) , for all x ∈ Ω (2.11)where s mn ( x ) = (cid:90) T (cid:2) a ( x )Ψ (cid:48)(cid:48) n ( t ) + a ( x )Ψ (cid:48) n ( t ) (cid:3) Ψ m ( t ) dt. Furthermore, from (2.7) we have ∂ ν u n ( x ) = (cid:90) T ∂ ν u ( x , t )Ψ n ( t ) dt, ∀ x ∈ ∂ Ω . Therefore, the Cauchy data of u n , for all n = 1 , . . . , N on the boundary ∂ Ω are determined by:1. Regarding Problem 2.1 ∂ ν u n ( x ) = (cid:90) T f ( x , t )Ψ n ( t ) dt,u n ( x ) = 0 (2.12)2. Regarding Problem 2.2 u n ( x ) = (cid:90) T f ( x , t )Ψ n ( t ) dt,∂ ν u n ( x ) = 0 . (2.13)4 lgorithm 1 A numerical method to solve Problem 2.1 Choose the basis { Ψ } n ≥ and a cut-off number N . Solve the system (2.11) with Cauchy data (2.12) for the vector valued function u comp ( x ) = ( u comp1 , . . . , u comp N ) , x ∈ Ω . The function p comp ( x ) is given by p comp ( x ) = u comp ( x ,
0) = N (cid:88) n =1 u comp n ( x )Ψ n (0) , x ∈ Ω . Algorithm 2
A numerical method to solve Problem 2.2 Choose the basis { Ψ } n ≥ and a cut-off number N . Solve the system (2.11) with Cauchy data (2.13) for the vector valued function u comp ( x ) = ( u comp1 , . . . , u comp N ) , x ∈ Ω . The function p comp ( x ) is given by p comp ( x ) = u comp ( x ,
0) = N (cid:88) n =1 u comp n ( x )Ψ n (0) , x ∈ Ω . The system of elliptic partial differential equations (2.11) together with Cauchy data either (2.12) or (2.13)is our approximation model, see Remark 2.1. It allows to determine coefficients u n , for all n = 1 , . . . , N ,and then the approximation u N ( x , t ) of u ( x , t ). The source term will be given by u N ( x , . In summary,the numerical method for solving Problem 2.1 and Problem 2.1 is described in Algorithm 1 and Algorithm2 below respectively.
Remark 2.2.
In Step 1 of Algorithm 1 and Algorithm 2, we choose the basis { Ψ n } n ≥ taken from [48].The cut-off number N is chosen numerically. More details will be discussed in Section 5. In Step 2 ofthese algorithms, we apply the quasi-reversibility method to solve (2.11) and (2.12) and (2.11) and (2.13) .The analysis about about the quasi-reversibility method and its convergence as the noise in the given datatends to are discussed in Section 4. As mentioned in Remark 2.2 that solving (2.11) and (2.12) for Problem 2.1 and solving (2.11) and(2.13) for Problem 2.2 are interesting when the given data in (2.4) and (2.5) contain noise. We employ thequasi-reversibility method to do so. Let (cid:15) be a small positive number playing the role of the regularizationparameter. To solve (2.11) and (2.12) we minimize the following mismatch functional J ( u , . . . , u N ) = N (cid:88) m =1 (cid:90) Ω (cid:12)(cid:12) Lu m − N (cid:88) n =1 s mn u n (cid:12)(cid:12) d x + N (cid:88) m =1 (cid:90) ∂ Ω (cid:12)(cid:12) ∂ ν u m − (cid:90) T f ( x , t )Ψ n ( t ) dt (cid:12)(cid:12) dσ + (cid:15) (cid:107) u (cid:107) H (Ω) N (2.14)5here u = ( u , . . . , u N ) is in H (Ω) N satisfying u ( x ) = 0 for all x ∈ ∂ Ω. To solve (2.11) and (2.13), weminimize the following mismatch functional J ( u , . . . , u N ) = N (cid:88) m =1 (cid:90) Ω (cid:12)(cid:12) Lu m − N (cid:88) n =1 s mn u n (cid:12)(cid:12) d x + N (cid:88) m =1 (cid:90) ∂ Ω (cid:12)(cid:12) u m − (cid:90) T f ( x , t )Ψ n ( t ) dt (cid:12)(cid:12) dσ + (cid:15) (cid:107) u (cid:107) H (Ω) N (2.15)where u = ( u , . . . , u N ) is in H (Ω) N satisfying ∂ ν u ( x ) = 0 for all x ∈ ∂ Ω.In the following sections, we prove that J has a unique minimizer and that minimizer converges tothe true solution of (2.11) and (2.12) as the noise level tends to 0. The results for J can be obtained inthe same manner. We introduce some auxiliary results in the next section. Let X be a number in (0 , (cid:101) Ω := (cid:110) ˜ x = (˜ x , . . . , ˜ x d ) : 0 < ˜ x + X − d (cid:88) i =2 ˜ x i < (cid:111) . Then, there exists α ∈ (0 , /
2) such that the function ψ (˜ x ) := ˜ x + 12 X d (cid:88) i =2 ˜ x i + α < , ∀ ˜ x ∈ (cid:101) Ω . (3.1)We have the following lemma. Lemma 3.1.
There are two positive constants λ and β depending only on α such that for all λ > λ and β > β , we have λβX e λψ − β |∇ φ | + λ β ψ − β − e λψ − β | φ | ≤ − CλβX e λψ − β φ | ∆ φ | + Cψ β +2 e λψ − β | ∆ φ | + divΦ (3.2) for all function φ ∈ C ( (cid:101) Ω) where the vector Φ satisfies | Φ | ≤ Ce λψ − β (cid:16) λβX |∇ φ | + λ β X ψ − β − | φ | (cid:17) . (3.3)Lemma 3.1 is a direct consequence of [49, Chapter 4, §
1, Lemma 3] in which the function φ isindependent of the time variable. Let ˜ R be a positive number such that Ω ⊂ B (0 , ˜ R ), where B (0 , ˜ R )denotes a ball of the center at 0 and the radius ˜ R . Let p and q be two positive numbers such that p B (0 , ˜ R ) + q e ⊂ (cid:101) Ω , (3.4)where e is the unit direction vector of x axis. For any x ∈ Ω, we define ˜ x := p x + q e , then (3.4) yields˜ x ∈ (cid:101) Ω. By modifying constant C in Lemma 3.1 (using Cp instead of C ) we have that the Lemma 3.1holds true in domain Ω. From now on, we apply Lemma 3.1 for all function in the space H (Ω). Thefollowing result plays an important role in our analysis.6 roposition 3.1 (Carleman estimate) . There exist λ > , β > , both of which only depend on a and α such that for all λ > λ and β > β for all φ ∈ H (Ω) , we have (cid:90) Ω ψ β +2 e λψ − β | ∆ φ | ≥ (cid:90) Ω Ce λψ − β (cid:2) λβ |∇ φ | + λ β ψ − β − | φ | (cid:3) d x − C (cid:90) ∂ Ω e λψ − β (cid:2) λβ |∇ φ | + λ β ψ − β − | φ | (cid:3) dσ (3.5) where C is a generic constant depending only on a , d , Ω , X and α .Proof. Let λ and β be as in Lemma 3.1. Fix λ > λ and β > β . Using the inequality ab ≤ λβψ − β − a + ψ β +2 b / (2 λβ ), for all φ ∈ C (Ω) , we have − λβe λψ − β φ ∆ φ ≤ λ β ψ − β − e λψ − β | φ | + 12 ψ β +2 e λψ − β | ∆ φ | in Ω . (3.6)Combining (3.2) and (3.6), since ψ <
1, we have ψ β +2 e λψ − β | ∆ φ | ≥ Ce λψ − β (cid:104) λβ |∇ φ | + λ β ψ − β − | φ | − div(Φ) (cid:105) in Ω (3.7)where Φ is a vector satisfying (3.3). Integrate (3.7) over Ω and apply the divergence theorem. Recaling(3.3), we obtain (3.5).The Carleman estimate (3.5) plays a key role for us to estimate the error of the solution to the inverseproblem assuming that the given data contains noise. As mentioned in section 2, we only prove the convergence of the quasi-reversibility method to solve(2.11) with Cauchy boundary data (2.12). In this case, the objective functional J , now named as J , see(2.14), is written as J ( u , . . . , u N ) = N (cid:88) m =1 (cid:90) Ω | Lu m − N (cid:88) n =1 s mn u n | d x + N (cid:88) m =1 (cid:90) ∂ Ω (cid:12)(cid:12)(cid:12) ∂ ν u m − (cid:90) T f ( x , t )Ψ m ( t ) dt (cid:12)(cid:12)(cid:12) dσ + (cid:15) N (cid:88) m =1 (cid:107) u m (cid:107) H (Ω) (4.1)subject to the constraint u = . . . = u N = 0 on ∂ Ω. Let H = { f ∈ H (Ω) N : f | ∂ Ω = 0 } be a closed subspace of H (Ω) N . It is clear that J is strictly convex in H . We now prove that J is hasa unique minimizer in H . Proposition 4.1.
The functional J has a unique minimizer in H . roof. Let e = inf (cid:8) J ( u ) : u = ( u , . . . , u n ) ∈ H (cid:9) ≥ { u i } i ≥ be a sequence satisfying lim i →∞ J ( u i ) = e. Then { u i } i ≥ is bounded in H . In fact, by contradiction, assume that { u i } i ≥ is unbounded. Then,there exists a subsequence, still named as { u i } i ≥ , satisfying lim i →∞ (cid:107) u i (cid:107) H (Ω) N = ∞ . Hence, e = lim i →∞ J ( u i ) ≥ lim i →∞ (cid:15) (cid:107) u i (cid:107) H (Ω) N = ∞ , which is impossible. Due to the boundedness of { u i } i ≥ in H , there exists a subsequence of { u i } i ≥ ,still named as { u i } i ≥ , which weakly converges to a function u in H . That implies { ∂ ν u i } i ≥ weaklyconverges to ∂ ν u in H / ( ∂ Ω) N , and therefore, { ∂ ν u i } i ≥ strongly converges to ∂ ν u in L ( ∂ Ω) by thecompact imbedding of H / ( ∂ Ω) into L ( ∂ Ω). Furthermore, the fact that { u i } i ≥ weakly converges u in H implies (cid:16) L ( u i ) m − (cid:80) Nn =1 s mn ( u i ) n (cid:17) Nm =1 weakly converges to (cid:16) L ( u ) m − (cid:80) Nn =1 s mn ( u ) n (cid:17) Nm =1 in L (Ω) N . As a result, N (cid:88) m =1 (cid:90) Ω (cid:12)(cid:12)(cid:12) L ( u ) m − N (cid:88) n =1 s mn ( u ) n (cid:12)(cid:12)(cid:12) d x ≤ lim sup i →∞ N (cid:88) m =1 (cid:90) Ω (cid:12)(cid:12)(cid:12) L ( u i ) m − N (cid:88) n =1 s mn ( u i ) n (cid:12)(cid:12)(cid:12) d x . Therefore J ( u ) ≤ lim sup i →∞ J ( u i ) = e. Thus u is a minimizer of J. The uniqueness of u is deduced from the strict convexity of J . Definition 4.1.
The unique minimizer, denoted by u min = ( u min1 , . . . , u min N ) , of functional J is said tobe the regularized solution of the problem (2.11) – (2.12) . We now assume that the measured data contain noise, with a noise level δ . We next study theconvergence of u min as noise level δ tends to 0. Let denote f δ ( x , t ) the noisy data and f ∗ ( x , t ) thecorresponding noiseless data, ( x , t ) ∈ ∂ Ω × [0 , T ]. By noise level δ , we mean (cid:18)(cid:90) T (cid:90) ∂ Ω | f δ ( x , t ) − f ∗ ( x , t ) | dσdt (cid:19) / < δ. Since the truncation number N is a finite number, we can write N (cid:88) m =1 (cid:90) ∂ Ω (cid:12)(cid:12)(cid:12) (cid:90) T f δ ( x , t )Ψ m ( t ) dt − (cid:90) T f ∗ ( x , t )Ψ m ( t ) dt (cid:12)(cid:12)(cid:12) ≤ Cδ , (4.2)where C is a generic constant depending only on N, a , b and c . For each m ∈ { , . . . , N } , define f ∗ m ( x ) = (cid:90) T f ∗ ( x , t )Ψ m ( t ) dt and f δm ( x ) = (cid:90) T f δ ( x , t )Ψ m ( t ) dt, (4.3)for all x ∈ ∂ Ω . The following theorem guarantees the Lipschitz stability of the reconstructed methodwith respect to noise. 8 heorem 4.1.
Let u δ min ∈ H be the minimizer of the functional J δ ( u , . . . , u N ) = N (cid:88) m =1 (cid:104) (cid:90) Ω | Lu m − N (cid:88) n =1 s mn u n | d x + (cid:90) ∂ Ω | ∂ ν u m − f δm | dσ (cid:105) + (cid:15) N (cid:88) m =1 (cid:107) u m (cid:107) H (Ω) . (4.4) Assume that the system Lu m − N (cid:88) n =1 s mn u n = 0 in Ω ,∂ ν u m = f ∗ m on ∂ Ω ,u m = 0 on ∂ Ω m ∈ { , . . . , N } (4.5) has the unique solution u ∗ = ( u ∗ , . . . , u ∗ N ) ∈ H . Then, (cid:107) u δ min − u ∗ (cid:107) H (Ω) N ≤ C (cid:0) δ + (cid:15) (cid:107) u ∗ (cid:107) H (Ω) N (cid:1) (4.6) where C is a generic constant depending only on N, a , b and c . As a result, let u comp and u ∗ be thefunctions obtained by (2.8) with ( u , . . . , u N ) replaced by u δ min and u ∗ respectively and let p comp ( x ) and p ∗ ( x ) be u comp ( x , and u ∗ ( x , respectively, x ∈ Ω . We have (cid:107) p comp − p ∗ (cid:107) H (Ω) ≤ C (cid:0) δ + (cid:15) (cid:107) u ∗ (cid:107) H (Ω) N (cid:1) . (4.7) Proof.
Due to (2.8), we have p comp ( x ) = N (cid:88) n =1 u comp n ( x )Ψ n (0) and p ∗ ( x ) = N (cid:88) n =1 u ∗ n ( x )Ψ n (0) . x ∈ Ω. Hence, (4.6) implies (4.7). It is sufficient to prove (4.6). Since u δ min = ( u , . . . , u N ) is theminimizer of J δ , for all h = ( h , . . . , h N ) ∈ H , N (cid:88) m =1 (cid:68) Lu m − N (cid:88) n =1 s mn u n , Lh m − N (cid:88) n =1 s mn h n (cid:69) L (Ω) + N (cid:88) m =1 (cid:68) ∂ ν u m − f δm , ∂ ν h m (cid:69) L ( ∂ Ω) + (cid:15) N (cid:88) m =1 (cid:104) u m , h m (cid:105) H (Ω) = 0 . (4.8)Since u ∗ = ( u ∗ , . . . , u ∗ N ) is the true solution to (4.5), for all h = ( h , . . . , h N ) ∈ H , N (cid:88) m =1 (cid:68) Lu ∗ m − N (cid:88) n =1 s mn u ∗ n , Lh m − N (cid:88) n =1 s mn h n (cid:69) L (Ω) + N (cid:88) m =1 (cid:68) ∂ ν u ∗ m − f ∗ m , ∂ ν h m (cid:69) L ( ∂ Ω) + (cid:15) N (cid:88) m =1 (cid:104) u ∗ m , h m (cid:105) H (Ω) = (cid:15) N (cid:88) m =1 (cid:104) u ∗ m , h m (cid:105) H (Ω) . (4.9)9ence, by subtracting (4.8) from (4.9) and setting h = ( h , . . . , h N ) = u δ min − u ∗ , we have N (cid:88) m =1 (cid:13)(cid:13)(cid:13) Lh m − N (cid:88) n =1 s mn h n (cid:13)(cid:13)(cid:13) L (Ω) + N (cid:88) m =1 (cid:104) ∂ ν h m − ( f δm − f ∗ m ) , ∂ ν h m (cid:105) L ( ∂ Ω) + (cid:15) N (cid:88) m =1 (cid:107) h m (cid:107) H (Ω) = − (cid:15) N (cid:88) m =1 (cid:104) u ∗ m , h m (cid:105) H (Ω) . Equivalently, N (cid:88) m =1 (cid:13)(cid:13)(cid:13) ∆ h m + b · ∇ h m + ch m − N (cid:88) n =1 s mn h n (cid:13)(cid:13)(cid:13) L (Ω) + N (cid:88) m =1 (cid:107) ∂h m (cid:107) L ( ∂ Ω) + (cid:15) N (cid:88) m =1 (cid:107) h m (cid:107) H (Ω) = N (cid:88) m =1 (cid:104) f δm − f ∗ m , ∂ ν h m (cid:105) L ( ∂ Ω) − (cid:15) N (cid:88) m =1 (cid:104) u ∗ m , h m (cid:105) H (Ω) . Using the inequality | ab | ≤ a + b , we have N (cid:88) m =1 (cid:13)(cid:13)(cid:13) ∆ h m + b · ∇ h m + ch m − N (cid:88) n =1 s mn h n (cid:13)(cid:13)(cid:13) L (Ω) + 12 N (cid:88) m =1 (cid:107) ∂ ν h m (cid:107) L ( ∂ Ω) + (cid:15) N (cid:88) m =1 (cid:107) h m (cid:107) H (Ω) ≤ N (cid:88) m =1 (cid:107) f δm − f ∗ m (cid:107) L ( ∂ Ω) + (cid:15) N (cid:88) m =1 (cid:107) u ∗ m (cid:107) H (Ω) . (4.10)It follows from (4.2), (4.3) and (4.10) that N (cid:88) m =1 (cid:90) Ω (cid:12)(cid:12)(cid:12) ∆ h m + b · ∇ h m + ch m − N (cid:88) n =1 s mn h n (cid:12)(cid:12)(cid:12) d x + N (cid:88) m =1 (cid:107) ∂ ν h m (cid:107) L ( ∂ Ω) ≤ C (cid:0) δ + (cid:15) (cid:107) u ∗ (cid:107) H (Ω) N (cid:1) (4.11)for a constant C > . It follows from (4.11) that N (cid:88) m =1 (cid:107) ∂ ν h m (cid:107) L ( ∂ Ω) ≤ C (cid:0) δ + (cid:15) (cid:107) u ∗ (cid:107) H (Ω) N (cid:1) . (4.12)Since h m = 0 on ∂ Ω, 1 ≤ m ≤ N , the tangent derivative of h m on ∂ Ω is 0. Hence, by (4.12) N (cid:88) m =1 (cid:107)∇ h m (cid:107) L ( ∂ Ω) ≤ C (cid:0) δ + (cid:15) (cid:107) u ∗ (cid:107) H (Ω) N (cid:1) . (4.13)Recall λ , β as in Lemma 3.1 and the function ψ as in (3.1). Fix β = β . Applying the inequality10 a − b ) ≥ a / − b , we have for all λ > λ N (cid:88) m =1 (cid:90) Ω (cid:12)(cid:12)(cid:12) ∆ h m + b · ∇ h m + ch m − N (cid:88) n =1 s mn h n (cid:12)(cid:12)(cid:12) d x ≥ min x ∈ Ω (cid:8) e − λψ − β ψ − β − (cid:9) N (cid:88) m =1 (cid:90) Ω e λψ − β ψ β +2 (cid:12)(cid:12)(cid:12) ∆ h m + b · ∇ h m + ch m − N (cid:88) n =1 s mn h n (cid:12)(cid:12)(cid:12) d x ≥ min x ∈ Ω (cid:8) e − λψ − β ψ − β − (cid:9)(cid:104) N (cid:88) m =1 (cid:90) Ω e λψ − β ψ β +2 | ∆ h m | d x − N (cid:88) m =1 (cid:90) Ω e λψ − β ψ β +2 (cid:12)(cid:12)(cid:12) b · ∇ h m + ch m − N (cid:88) n =1 s mn h n (cid:12)(cid:12)(cid:12) d x (cid:105) . Thus, by (4.11),min x ∈ Ω (cid:8) e − λψ − β ψ − β − (cid:9)(cid:104) N (cid:88) m =1 (cid:90) Ω e λψ − β ψ β +2 | ∆ h m | d x − N (cid:88) m =1 (cid:90) Ω e λψ − β ψ β +2 (cid:12)(cid:12)(cid:12) b · ∇ h m + ch m − N (cid:88) n =1 s mn h n (cid:12)(cid:12)(cid:12) d x (cid:105) ≤ C (cid:0) δ + (cid:15) (cid:107) u ∗ (cid:107) H (Ω) N (cid:1) . Applying the Carleman estimate (3.5), we havemin x ∈ Ω (cid:8) e − λψ − β ψ − β − (cid:9)(cid:104) C N (cid:88) m =1 (cid:90) Ω e λψ − β [ λβ |∇ h m | + λ β ψ − β − | h m | ] d x − C (cid:90) ∂ Ω e λψ − β (cid:2) λβ |∇ h m | + λ β ψ − β − | h m | (cid:3) d ] σ − N (cid:88) m =1 (cid:90) Ω e λψ − β ψ β +2 (cid:12)(cid:12)(cid:12) b · ∇ h m + ch m − N (cid:88) n =1 s mn h n (cid:12)(cid:12)(cid:12) d x (cid:105) ≤ C (cid:0) δ + (cid:15) (cid:107) u ∗ (cid:107) H (Ω) N (cid:1) . Since β = β fixed, choosing λ sufficiently large,we have N (cid:88) m =1 (cid:90) Ω (cid:2) λ |∇ h m | + λ | h m | (cid:3) d x ≤ C (cid:90) ∂ Ω | λ |∇ h m | + λ | h m | dσ + C (cid:0) δ + (cid:15) (cid:107) u ∗ (cid:107) H (Ω) N (cid:1) (4.14)Since h = u δ − u ∗ ∈ H , h m = 0 on ∂ Ω . Hence, we obtain (4.6) by using (4.13) and (4.14).
Remark 4.1.
The conclusion of Theorem 4.1 is similar to some theorems about the quasi-reversibilitymethod we have developed, see e.g. [41, Theorem 5.1] and [42, Theorem 4.1]. The main difference isthat in Theorem 4.1, we relax a technical condition that there exists an error vector valued function E ,well-defined in the whole Ω , such that ∂ ν E = f δ − f ∗ and (cid:107)E(cid:107) H (Ω) = O ( δ ) . The analysis for the quasi-reversibility method to solve (2.11) and (2.13) is similar to the argumentsabove. We do not repeat the proof here. 11
Numerical studies
In this section, we set the dimension d = 2 and Ω = ( − R, R ) with R = 1. Define a grid of points onΩ as G = { ( x i = − R + ( i − h x , y j = − R + ( j − h x ) : 1 ≤ i, j ≤ N x } where h x = 2 R/ ( N x −
1) with N x = 81 . We set T = 2. On [0 , T ], we also define the uniform partition T = { t i = ( i − h t : 1 ≤ i ≤ N T } where h t = T / ( N T − N T = 201. To generate the simulated data, we solve (2.2)with a ( x, y ) = 1 + sin ( x + y ) , b ( x, y ) = (2 , ,c ( x, y ) = cos( x + y ) ,e ( x, y ) = 0 . x + y ) + sin( x + y )] ( x, y ) ∈ Ωby the finite difference method in the implicit scheme. Let u ∗ ( x , t ), x ∈ G and t ∈ T , be the obtainednumerical solution. We can extract the data B u ∗ ( x , t ) and F u ∗ ( x , t ) on ( ∂ Ω × [0 , T ]) ∩ ( G × T ). Thesefunctions serve as the data without noise. For δ >
0, the noisy data are given by B u ( x , t ) = B u ∗ ( x , t )(1 + δ rand( x , t )) and F u ( x , t ) = F u ∗ ( x , t )(1 + δ rand( x , t ))where rand is the function taking uniformly distributed random numbers in [ − , . We present in details the implementation of Steps 1 and 2 of Algorithm 1 to solve Problem 2.1 whilethe other steps can be implemented directly. The implementation for Problem 2.2 can be done in thesame manner.
Step 1 in Algorithm 1.
In our numerical studies, we employ the basis { Ψ n } n ≥ that was first introducedby Klibanov in [48]. For each n ≥ , we define Φ n ( t ) := t n − e t − T/ . The set { Φ n : n ≥ } is complete in L (0 , T ). We apply the Gram-Schmidt orthonormalization process on this set to obtain the orthonormalbasis { Ψ n } n ≥ of L (0 , T ) . Remark 5.1.
The basis { Ψ n } n ≥ was successfully used very often in our research group to solve a longlist of inverse problems including the nonlinear coefficient inverse problems for elliptic equations [50] andparabolic equations [51, 43, 44, 52], and ill-posed inverse source problems for elliptic equations [41] andparabolic equations [42], transport equations [45] and full transfer equations [53]. Another reason for usto employ this basis rather than the well-known basis of the Fourier series is that the first elements ofthis basis is a constant. Hence, when we plug (2.9) into (2.2) , the information of u ( x )Ψ (cid:48)(cid:48) ( t ) will be lost.As a result, the contribution of u ( x ) in (2.11) is less than the that of the corresponding u ( x ) obtainedby the basis { Ψ n } n ≥ . To choose N , we numerically compare u ∗ ( x ,
0) and N (cid:88) n =1 u ∗ n ( x )Ψ n (0) for x ∈ Ω where u ∗ is the truesolution to (2.2) and the source p ( x ) is given in Example 1 below. The number N is chosen such thatthe error (cid:12)(cid:12)(cid:12) u ∗ ( x , − N (cid:88) n =1 u ∗ n ( x )Ψ n (0) (cid:12)(cid:12)(cid:12)
12s small enough. We perform this procedure and choose N = 35, see Figure 1. This cut off number isused for all numerical examples in the paper. (a) N = 15 (b) N = 20 (c) N = 35 Figure 1: The function (cid:12)(cid:12)(cid:12) u ∗ ( x , − N (cid:88) n =1 u ∗ n ( x )Ψ n (0) (cid:12)(cid:12)(cid:12) , x ∈ Ω where the function u ∗ is the solution to (2.2)and the source is given in Example 2. It is evident that the larger N , the better approximation in (2.8)is. Step 2 in Algorithm 1.
In this step, we apply the quasi-reversibility method to solve the system(2.11) with Cauchy data (2.12). That means, we minimize the functional J defined in (2.14). The finitedifference version of J , still called J , is J ( u , . . . , u N ) = h x N (cid:88) m =1 N x − (cid:88) i,j =2 (cid:12)(cid:12)(cid:12) ∆ u m ( x i , y j ) + b ( x i , y j ) · ∇ u m ( x i , y j ) + c ( x i , y j ) u m ( x i , y j ) − N (cid:88) n =1 s mn u n ( x i , y j ) (cid:12)(cid:12)(cid:12) + h x N (cid:88) m =1 N x (cid:88) j =1 ( | − ∂ x u m ( x , y j ) − f m ( x , y j ) | + | ∂ x u m ( x N x , y j ) − f m ( x N x , y j ) | )+ h x N (cid:88) m =1 N x − (cid:88) i =2 ( | − ∂ y u m ( x i , y ) − f m ( x i , y ) | + | ∂ y u m ( x i , y N x ) − f m ( x i , y N x ) | )+ h x N (cid:88) m =1 N x (cid:88) j =1 ( | u m ( x , y j ) | + | u m ( x N x , y j ) | ) + h x N (cid:88) m =1 N x − (cid:88) i =2 ( | u m ( x i , y ) | + | u m ( x i , y N x ) | )+ (cid:15)h x N (cid:88) m =1 N x − (cid:88) i =2 | u m ( x i , y j ) | + |∇ u m ( x i , y j ) | + | ∆ u m ( x i , y j ) | . (5.1)Here, instead of imposing the constraint u m = 0 on ∂ Ω, we add additional term: h x (cid:80) Nm =1 (cid:80) N x j =1 ( | u m ( x , y j ) | + | u m ( x N x , y j ) | ) + h x (cid:80) Nm =1 (cid:80) N x − i =2 ( | u m ( x i , y ) | + | u m ( x i , y N x ) | ) to the right hand side of (5.1). Thistechnique significantly reduces the efforts in the implementation. We now identify the vector valuefunction ( u , . . . , u N ) with it “line up” version u i = u m ( x i , y j ) where i = ( i − N N x + ( j − N + m for 1 ≤ i, j ≤ N x and 1 ≤ m ≤ N . The data f is also line-up in the same manner. f i = f m ( x i , y j ) where i = ( i − N N x + ( j − N + m i ∈ { , N x } , 1 ≤ j ≤ N x or 1 ≤ i ≤ N x , j ∈ { , N x } . It is not hard to rewrite J in term of u as J ( u ) = |L u | + |N u − f | + |D u | + (cid:15) | ( | u | + | D x u | + | u | + | L u | ) (5.2)for some matrices L , N , D , D x , D y and L . The matrix L is such that( L u ) i = h x (cid:16) ∆ u m ( x i , y j ) + b ( x i , y j ) · ∇ u m ( x i , y j ) + c ( x i , y j ) u m ( x i , y j ) − N (cid:88) n =1 s mn u n ( x i , y j ) (cid:17) with i = ( i − N N x + ( j − N + m, ≤ i, j ≤ N x − ≤ m ≤ N. The matrix N and D are the matrices such that N u and D u respectively correspond to the Neumann and Dirichlet values of u m ( x i , y j ) where ( x i , y j ) is on ∂ Ω, 1 ≤ m ≤ N. The matrix D x , D y and L are such that D x u , D y u and L u correspond to ∂ x u m ( x i , y j ), ∂ y u m ( x i , y j ) and ∆ u m ( x i , y j ), 2 ≤ i, j ≤ N x −
1, 1 ≤ m ≤ N. Theexplicit forms of these matrices can be written similarly to [42, Section 5.1]. For the brevity, we do notrepeat the details here.Since u is the minimizer of J defined in (5.2), u satisfies( L T L + N T N + D T D + (cid:15) (Id + D T x D x + D T y D y + L T1 L )) u = N T f where the superscript T indicates the transpose of matrices. This linear system can be solve by anylinear algebra package. We employ the command “lsqlin” of MATLAB for this purpose. In all followingexamples, the regularization parameter (cid:15) is chosen to be (cid:15) = 10 − . Step 3 of Algorithm 1 . These steps can be implemented directly since they involve only explicitformulas.Again, the implementation for solving Problem 2.2 is similar to that for solving Problem 2.1. We donot repeat the full process for this case. For the brevity, we just provide some numerical results.
We now perform four numerical examples for both Problem 2.1 and Problem 2.2.
Example 1.
We consider the true source function given by p ∗ ( x, y ) = (cid:26) x − . + y < . , δ = 10% the maximal computedvalue of the source is 0.96115 (relative error 3.9%) while in the case δ = 100% , the maximal computedvalue of the source is 1.01114 (relative error 1.1%). Regarding to Problem 2.2, in the case δ = 10% themaximal computed value of the source is 0.99389 (relative error 0.6%) while in the case δ = 100% , themaximal computed value of the source is 0.98797 (relative error 1.2%). Example 2 . We consider a more complicated source function p ∗ ( x, y ) = x − . + y < . , {| x + 0 . | , | y − . |} < . , , where the support of the source function consists of a disk and a square, see Figure 3 (a).14 a) The true source function (b) The computed solution to Prob-lem 2.1 from data with 10% noise (c) The computed solution to Prob-lem 2.1 from data with 100% noise(d) The computed solution to Prob-lem 2.2 from data with 10% noise (e) The computed solution to Prob-lem 2.2 from data with 100% noise Figure 2: Example 1. The true source function and the reconstructions of source functions.The reconstructions of source p ∗ are displayed in Figure 3, which show the accurate reconstructionsof the square and the disk. The computed values of the source function are quite accurate. Regardingto Problem 2.1, in the case δ = 10% the maximal computed values of the source in the square andthe disk are 1.97386 (relative error 1.3%) and 0.9608 (relative error 3.9%) respectively while in the case δ = 100% , the corresponding maximal computed values of the source are 2.19941 (relative error 9.9%)and 1.157 (relative error 15.7%) . Regarding to Problem 2.2, in the case δ = 10%, the maximal computedvalues of the source in the square and the disk are 2.01843 (relative error 0.9%) and 0.9994 (relative error0.0%) while in the case δ = 100%, the corresponding maximal computed values of the source are 2.11776(relative error 5.9%) and 0.9733 (relative error 2.3%). We observe that when the noise level is 100%,the values of the source are well computed while and the reconstructed shapes of the inclusions start tobreak out. Example 3.
We next test the case where the support of the source has more complicated geometriesthan the one in Example 2, and the source has both positive and negative values. p ∗ ( x, t ) = { | x − . | , | y |} < . x − . + y ≥ . , − . x + 0 . + ( y − . ≤ . , . The support of the source involves a rectangle with a void and an ellipse, see Figure 4(a).The numerical solutions of Example 3 are displayed in Figure 4, which show the accurate reconstruc-tions of the rectangle with the void and the ellipse. The computed values of the source function are quiteaccurate. Regarding to Problem 2.1, in the case δ = 10% the maximal and minimal computed values of15 a) The true source function (b) The computed solution to Prob-lem 2.1 from data with 10% noise (c) The computed solution to Prob-lem 2.1 from data with 100% noise(d) The computed solution to Prob-lem 2.2 from data with 10% noise (e) The computed solution to Prob-lem 2.2 from data with 100% noise Figure 3: Example 2. The true source function and the reconstructions.the source is 3.22793 (relative error 7.6%) and − . δ = 100% , the maximal and minimal computed values of the source are 3.21003 (relative error 7.0%) and − . δ = 10%, the maximaland minimal computed values of the source are 3.04546 (relative error 1.5%) and − . δ = 100%, the maximal and minimal computed values of the sourceare 3.15653 (relative error 5.2%) and − . Example 4.
In this example, the true source function p ∗ is the characteristic function of the letter T .The numerical solutions of Example 4 are displayed in Figure 5, which show the accurate reconstruc-tions of the letter T . The computed values of the source function are quite accurate. Regarding toProblem 2.1, in the case δ = 10% the maximal computed value of the source are 0.97705 (relative error2.3%) while in the case δ = 100% , the maximal computed value of the source is 0.93572 (relative error6.4%). Regarding to Problem 2.2, in the case δ = 10%, the maximal computed value of the source is1.00401 (relative error 0.4%) and in the case δ = 100%, the maximal computed value of the source is0.94551 (relative error 5.4%). We observe that when the noise level is 100%, the reconstructed valuesand the T -shape of the source meet the expectation.16 a) The true source function (b) The computed solution to Prob-lem 2.1 from data with 10% noise (c) The computed solution to Prob-lem 2.1 from data with 100% noise(d) The computed solution to Prob-lem 2.2 from data with 10% noise (e) The computed solution to Prob-lem 2.2 from data with 100% noise Figure 4: Example 3. The true source function and the reconstructions. { Ψ n } n ≥ As mention in Remark 5.1, we choose the basis { Ψ n } n ≥ rather than the the “sin and cosine” basisof the well-known Fourier series. To numerically verifying that this choice is important, we comparesome 1D numerical solutions to Problem 2.1 obtained by our method with respect to two bases: (1) thetrigonometric Fourier expansion and (2) the basis { Ψ n } n ≥ . We will show that the numerical results incase (1) do not meet the expectation while the numerical results in case (2) do.For the simplicity, we drop the damping term au t in the governing equation. The governing model is a ( x ) u tt ( x, t ) = u xx ( x, t ) + b ( x ) u x ( x, t ) + c ( x ) u ( x, t ) ( x, t ) ∈ ( − , × [0 , T ] ,u t ( x, t ) = 0 x ∈ ( − , ,u ( x,
0) = p ( x ) x ∈ ( − , ,u ( x, t ) = 0 x ∈ {− , } . (6.1)Here, we choose T = 4. The aim of Problem 2.1 is to compute the initial condition p ( x ) from themeasurement of f ( ± , t ) = u x ( ± , t ) t ∈ [0 , T ] . (6.2)We now try the trigonometric Fourier expansion to solve Problem 2.1. In this section, we display thenumerical results with the cut-off number N = 35. We have tried the cases when N = 50 and N = 100but the quality of the computed sources does not improve. For each x ∈ ( − , u ( x, t )17 a) The true source function (b) The computed solution to Prob-lem 2.1 from data with 10% noise (c) The computed solution to Prob-lem 2.1 from data with 100% noise(d) The computed solution to Prob-lem 2.2 from data with 10% noise (e) The computed solution to Prob-lem 2.2 from data with 100% noise Figure 5: Example 4. The true source function and the reconstructions.by the N − partial sum of its Fourier series u ( x, t ) = u ( x ) + N (cid:88) n =1 u n ( x ) cos (cid:16) πntT (cid:17) + N (cid:88) n =1 v n ( x ) sin (cid:16) πntT (cid:17) (6.3)where u ( x ) = 1 T (cid:90) T u ( x, t ) dt,u n ( x ) = 2 T (cid:90) T u ( x, t ) cos (cid:16) πntT (cid:17) n ≥ ,v n ( x ) = 2 T (cid:90) T u ( x, t ) sin (cid:16) πntT (cid:17) n ≥ , for x ∈ ( − , . Differentiate (6.3) with respect to t , we have u tt ( x, t ) = − N (cid:88) n =1 (cid:16) πnT (cid:17) u n ( x ) cos (cid:16) πntT (cid:17) − N (cid:88) n =1 (cid:16) πnT (cid:17) v n ( x ) sin (cid:16) πntT (cid:17) (6.4)18lugging (6.3) and (6.4) into (6.1) gives − a ( x ) (cid:104) N (cid:88) n =1 (cid:16) πnT (cid:17) u n ( x ) cos (cid:16) πntT (cid:17) + N (cid:88) n =1 (cid:16) πnT (cid:17) v n ( x ) sin (cid:16) πntT (cid:17)(cid:105) = (cid:104) u (cid:48)(cid:48) ( x ) + N (cid:88) n =1 u (cid:48)(cid:48) n ( x ) cos (cid:16) πntT (cid:17) + N (cid:88) n =1 v (cid:48)(cid:48) n ( x ) sin (cid:16) πntT (cid:17)(cid:105) + b ( x ) (cid:104) u (cid:48) ( x ) + N (cid:88) n =1 a (cid:48) n ( x ) cos (cid:16) πntT (cid:17) + N (cid:88) n =1 v (cid:48) n ( x ) sin (cid:16) πntT (cid:17)(cid:105) + c ( x ) (cid:104) u ( x ) + N (cid:88) n =1 u n ( x ) cos (cid:16) πntT (cid:17) + N (cid:88) n =1 v n ( x ) sin (cid:16) πntT (cid:17)(cid:105) (6.5)for all x ∈ ( − , , t ∈ [0 , T ] . Hence, we have u (cid:48)(cid:48) ( x ) + b ( x ) u (cid:48) ( x ) + c ( x ) u ( x ) = 0 ,u (cid:48)(cid:48) n ( x ) + b ( x ) u (cid:48) n ( x ) + (cid:104)(cid:16) πnT (cid:17) a ( x ) + c ( x ) (cid:105) u n ( x ) = 0 n ≥ v (cid:48)(cid:48) n ( x ) + b ( x ) v (cid:48) n ( x ) + (cid:104)(cid:16) πnT (cid:17) a ( x ) + c ( x ) (cid:105) v n ( x ) = 0 n ≥ x ∈ ( − , . The boundary constraints of u , { u n } n ≥ and { v n } n ≥ are u ( ±
1) = u n ( ±
1) = v n ( ±
1) = 0 n ≥ , ( u ) x ( ±
1) = 1 T (cid:90) T f ( ± , t ) dt, ( u n ) x ( ±
1) = 2 T (cid:90) T f ( ± , t ) cos (cid:16) πtT (cid:17) dt n ≥ , ( v n ) x ( ±
1) = 2 T (cid:90) T f ( ± , t ) sin (cid:16) πtT (cid:17) dt n ≥ u , u n , v n for n ≥
1. Then, the function initial condition is given by u ( x,
0) where u ( x, t ) is given by (6.3). We skip presenting the implementation for this method. The implementation issimilar but simpler than the implementation for the 2D case.In the numerical tests, we set a ( x ) = 1+sin ( x ), b ( x ) = sin( πx ) and c ( x ) = cos(2 πx ) for x ∈ ( − , . In this section, we perform three (3) tests. The source functions p , p and p correspond to Test 1, Test2 and Test 3 are given below: p = (cid:26) exp( | x − . | / ( | x − . | − . )) | x − . | < . , .p = 1 − x ,p = sin( πx )We compute these source functions with noise level 5%. The numerical examples are displayed inFigure 6. It is evident that approximating u ( x, t ) with the basis { Ψ n } n ≥ provides much better solutionsto Problem 2.1 in comparison to approximating u ( x, t ) with the popular “sin and cosine” basis.19 a) Test 1 (b) Test 2 (c) Test 3(d) Test 1 (e) Test 2 (f) Test 3 Figure 6: The true and computed source functions. Row 1 are the results obtained by solving (6.6)–(6.7)and row 2 are the results by the 1D version of Algorithm 1. It is evident that solving Problem 2.1 byusing the “sin and cosine” basis is not succesful. In contrast, the reconstructions using the basis { Ψ n } n ≥ are quite accurate. 20 Concluding remarks
In this paper, we introduced a new approach to numerically compute the source function for generalhyperbolic equations from the Cauchy boundary data. In the first step, by truncating the Fourier seriesof the solution to this hyperbolic equation, we derive an approximation model, who solution directlyprovides the knowledge of the source. We then apply the quasi-reversibility method to solve this system.The convergence of the quasi-reversibility method is rigorously proved. Satisfactory numerical examplesillustrates the efficiency of our method.
Acknowledgments
Thuy Le and Loc Nguyen are supported by US Army Research Laboratory and US Army ResearchOffice grant W911NF-19-1-0044. The authors would like to thank Dr. Michael V. Klibanov for manyfruitful discussions.
References [1] R. A. Kruger, P. Liu, Y. R. Fang, and C. R. Appledorn. Photoacoustic ultrasound (PAUS)–reconstruction tomography.
Med. Phys. , 22:1605, 1995.[2] A. Oraevsky, S. Jacques, R. Esenaliev, and F. Tittel. Laser-based optoacoustic imaging in biologicaltissues.
Proc. SPIE , 2134A:122, 1994.[3] R. A. Kruger, D. R. Reinecke, and G. A. Kruger. Thermoacoustic computed tomography: technicalconsiderations.
Med. Phys. , 26:1832, 1999.[4] N. Do and L. Kunyansky. Theoretically exact photoacoustic reconstruction from spatially andtemporally reduced data.
Inverse Problems , 34(9):094004, 2018.[5] M. Haltmeier. Inversion of circular means and the wave equation on convex planar domains.
Comput.Math. Appl. , 65:1025–1036, 2013.[6] F. Natterer. Photo-acoustic inversion in convex domains.
Inverse Probl. Imaging , 6:315–320, 2012.[7] L. V. Nguyen. A family of inversion formulas in thermoacoustic tomography.
Inverse Probl. Imaging ,3:649–675, 2009.[8] V. Katsnelson and L. V. Nguyen. On the convergence of time reversal method for thermoacoustictomography in elastic media.
Applied Mathematics Letters , 77:79–86, 2018.[9] Y. Hristova. Time reversal in thermoacoustic tomography–an error estimate.
Inverse Problems ,25:055008, 2009.[10] Y. Hristova, P. Kuchment, and L. V. Nguyen. Reconstruction and time reversal in thermoacoustictomography in acoustically homogeneous and inhomogeneous media.
Inverse Problems , 24:055006,2008.[11] P. Stefanov and G. Uhlmann. Thermoacoustic tomography with variable sound speed.
InverseProblems , 25:075011, 2009. 2112] P. Stefanov and G. Uhlmann. Thermoacoustic tomography arising in brain imaging.
Inverse Prob-lems , 27:045004, 2011.[13] C. Clason and M. V. Klibanov. The quasi-reversibility method for thermoacoustic tomography in aheterogeneous medium.
SIAM J. Sci. Comput. , 30:1–23, 2007.[14] C. Huang, K. Wang, L. Nie, L. V. Wang, and M. A. Anastasio. Full-wave iterative image recon-struction in photoacoustic tomography with acoustically inhomogeneous media.
IEEE Trans. Med.Imaging , 32:1097–1110, 2013.[15] G. Paltauf, R. Nuster, M. Haltmeier, and P. Burgholzer. Experimental evaluation of reconstructionalgorithms for limited view photoacoustic tomography with line detectors.
Inverse Problems , 23:S81–S94, 2007.[16] G. Paltauf, J. A. Viator, S. A. Prahl, and S. L. Jacques. Iterative reconstruction algorithm foroptoacoustic imaging.
J. Opt. Soc. Am. , 112:1536–1544, 2002.[17] H. Ammari, E. Bretin, E. Jugnon, and V. Wahab. Photoacoustic imaging for attenuating acous-tic media. In H. Ammari, editor,
Mathematical Modeling in Biomedical Imaging II , pages 57–84.Springer, 2012.[18] H. Ammari, E. Bretin, J. Garnier, and V. Wahab. Time reversal in attenuating acoustic media.
Contemp. Math. , 548:151–163, 2011.[19] M. Haltmeier and L. V. Nguyen. Reconstruction algorithms for photoacoustic tomography in het-erogeneous damping media.
Journal of Mathematical Imaging and Vision , 61:1007–1021, 2019.[20] S. Acosta and B. Palacios. Thermoacoustic tomography for an integro-differential wave equationmodeling attenuation.
J. Differential Equations , 5:1984–2010, 2018.[21] P. Burgholzer, H. Gr¨un, M. Haltmeier, R. Nuster, and G. Paltauf. Compensation of acoustic at-tenuation for high-resolution photoa- coustic imaging with line detectors.
Proc. SPIE , 6437:643724,2007.[22] A. Homan. Multi-wave imaging in attenuating media.
Inverse Probl. Imaging , 7:1235–1250, 2013.[23] R. Kowar. On time reversal in photoacoustic tomography for tissue similar to water.
SIAM J.Imaging Sci. , 7:509–527, 2014.[24] R. Kowar and O. Scherzer. Photoacoustic imaging taking into account attenuation. In H. Ammari,editor,
Mathematics and Algorithms in Tomography II, Lecture Notes in Mathematics , pages 85–130.Springer, 2012.[25] A. I. Nachman, J. F. Smith III, and R.C. Waag. An equation for acoustic propagation in inhomo-geneous media with relaxation losses.
J. Acoust. Soc. Am. , 88:1584–1595, 1990.[26] B. Cox, S. Arridge, and P. Beard. Photoacoustic tomography with a limited-aperture planar sensorand a reverberant cavity.
Inverse Problems , 23:S95, 2007.[27] B. Cox and P. Beard. Photoacoustic tomography with a single detector in a reverberant cavity.
TheJournal of the Acoustical Society of America , 123:3371–3371, 2008.2228] L. Kunyansky, B. Holman, and B. Cox. Photoacoustic tomography in a rectangular reflecting cavity.
Inverse Problems , 29:125010, 2013.[29] L. V. Nguyen and L. Kunyansky. A dissipative time reversal technique for photo-acoustic tomographyin a cavity.
SIAM J. Imaging Sci. , 9:748–769, 2016.[30] R. Latt`es and J. L. Lions.
The Method of Quasireversibility: Applications to Partial DifferentialEquations . Elsevier, New York, 1969.[31] E. B´ecache, L. Bourgeois, L. Franceschini, and J. Dard´e. Application of mixed formulations of quasi-reversibility to solve ill-posed problems for heat and wave equations: The 1d case.
Inverse Problems& Imaging , 9(4):971–1002, 2015.[32] L. Bourgeois. A mixed formulation of quasi-reversibility to solve the Cauchy problem for Laplace’sequation.
Inverse Problems , 21:1087–1104, 2005.[33] L. Bourgeois. Convergence rates for the quasi-reversibility method to solve the Cauchy problem forLaplace’s equation.
Inverse Problems , 22:413–430, 2006.[34] L. Bourgeois and J. Dard´e. A duality-based method of quasi-reversibility to solve the Cauchy problemin the presence of noisy data.
Inverse Problems , 26:095016, 2010.[35] J. Dard´e. Iterated quasi-reversibility method applied to elliptic and parabolic data completionproblems.
Inverse Problems and Imaging , 10:379–407, 2016.[36] M. V. Klibanov , A. V. Kuzhuget, S. I. Kabanikhin, and D. Nechaev. A new version of the quasi-reversibility method for the thermoacoustic tomography and a coefficient inverse problem.
ApplicableAnalysis , 87:1227–1254, 2008.[37] M. V. Klibanov. Carleman estimates for global uniqueness, stability and numerical methods forcoefficient inverse problems.
J. Inverse and Ill-Posed Problems , 21:477–560, 2013.[38] M. V. Klibanov and F. Santosa. A computational quasi-reversibility method for Cauchy problemsfor Laplace’s equation.
SIAM J. Appl. Math. , 51:1653–1675, 1991.[39] M. V. Klibanov and J. Malinsky. Newton-Kantorovich method for 3-dimensional potential inversescattering problem and stability for the hyperbolic Cauchy problem with time dependent data.
Inverse Problems , 7:577–596, 1991.[40] M. V. Klibanov. Carleman estimates for the regularization of ill-posed Cauchy problems.
AppliedNumerical Mathematics , 94:46–74, 2015.[41] L. H. Nguyen, Q. Li, and M. V. Klibanov. A convergent numerical method for a multi-frequencyinverse source problem in inhomogenous media.
Inverse Problems and Imaging , 13:1067–1094, 2019.[42] Q. Li and L. H. Nguyen. Recovering the initial condition of parabolic equations from lateral Cauchydata via the quasi-reversibility method.
Inverse Problems in Science and Engineering , 28:580–598,2020.[43] M. V. Klibanov and L. H. Nguyen. PDE-based numerical method for a limited angle X-ray tomog-raphy.
Inverse Problems , 35:045009, 2019. 2344] V. A. Khoa, M. V. Klibanov, and L. H. Nguyen. Convexification for a 3D inverse scattering problemwith the moving point source.
SIAM J. Imaging Sci. , 13(2):871–904, 2020.[45] M. V. Klibanov, T. T. Le, and L. H. Nguyen. Convergent numerical method for a linearized traveltime tomography problem with incomplete data. to appear on SIAM Journal on Scientific Comput-ing, see also https://arxiv.org/abs/1911.04581 , 2020.[46] L. C. Evans.
Partial Differential Equations . Graduate Studies in Mathematics, Volume 19. Amer.Math. Soc., 2010.[47] M. V. Klibanov and D-L. Nguyen. Convergence of a series associated with the convexification methodfor coefficient inverse problems. preprint, arXiv:2004.05660 , 2020.[48] M. V. Klibanov. Convexification of restricted Dirichlet to Neumann map.
J. Inverse and Ill-PosedProblems , 25(5):669–685, 2017.[49] M. M. Lavrent’ev, V. G. Romanov, and S. P. Shishat · ski˘i. Ill-Posed Problems of MathematicalPhysics and Analysis . Translations of Mathematical Monographs. AMS, Providence: RI, 1986.[50] V. A. Khoa, G. W. Bidney, M. V. Klibanov, L. H. Nguyen, L. Nguyen, A. Sullivan,and V. N. Astratov. Convexification and experimental data for a 3D inverse scatter-ing problem with the moving point source.
Inverse Problems, accepted for publication, seehttps://iopscience.iop.org/article/10.1088/1361-6420/ab95aa/meta , 2020.[51] L. H. Nguyen. A new algorithm to determine the creation or depletion term of parabolic equationsfrom boundary measurements. preprint, arXiv:1906.01931 , 2019.[52] V. A. Khoa, G. W. Bidney, M. V. Klibanov, L. H. Nguyen, L. Nguyen, A. Sullivan, and V. N. Astra-tov. An inverse problem of a simultaneous reconstruction of the dielectric constant and conductivityfrom experimental backscattering data. to appear on Inverse Problems in Science and Engineering,see also arXiv:2006.00913 , 2020.[53] A. V. Smirnov, M. V. Klibanov, and L. H. Nguyen. On an inverse source problem for the fullradiative transfer equation with incomplete data.