Theoretical stability in coefficient inverse problems for general hyperbolic equations with numerical reconstruction

TL;DR

This paper establishes theoretical local stability results for determining spatial coefficients in hyperbolic equations and proposes a numerical iterative method using Tikhonov regularization and Poisson equations, demonstrated through 1D examples.

Contribution

It provides new theoretical stability estimates for coefficient inverse problems in hyperbolic equations and introduces a practical iterative reconstruction method with numerical validation.

Findings

Local Hölder stability under geometric conditions

Effective iterative reconstruction method based on Tikhonov regularization

Numerical examples demonstrating the method's performance

Abstract

In this article, we investigate the determination of the spatial component in the time-dependent second order coefficient of a hyperbolic equation from both theoretical and numerical aspects. By the Carleman estimates for general hyperbolic operators and an auxiliary Carleman estimate, we establish local H\"older stability with both partial boundary and interior measurements under certain geometrical conditions. For numerical reconstruction, we minimize a Tikhonov functional which penalizes the gradient of the unknown function. Based on the resulting variational equation, we design an iteration method which is updated by solving a Poisson equation at each step. One-dimensional prototype examples illustrate the numerical performance of the proposed iteration.

Tables

Table 1. Table 1: The numerical performance of Algorithm 4.4 for various choices of true solutions in Example 5.1 .

Case	$p_{\mathrm{true}}\left(x\right)$	$N$	$\mathrm{err}$	Elapsed time
(a)	$\frac{1}{2}\sin \pi x+1$	$21$	$0.86\%$	$0.51s$
(b)	$x\left(x-1\right)\left(x-\frac{3}{2}\right)+1$	$35$	$0.61\%$	$0.97s$
(c)	$\frac{1}{2}\min (x,1-x)+\frac{1}{4}\sin \pi x+1$	$9$	$1.72\%$	$0.25s$

Table 2. Table 2: Parameters and corresponding numerical performance in Example 5.2 under various combinations of the relative noise levels δ 0 subscript 𝛿 0 \delta_{0} and the observation subdomains ω 𝜔 \omega .

$\omega$	${\delta }_0$	$K$	$\alpha$	$N$	$\mathrm{err}$	Elapsed time
$\Omega \setminus [0.2,0.8]$	$1\%$	$4\times {10}^{-5}$	${10}^{-7}$	$17$	$0.66\%$	$0.43s$
$\Omega \setminus [0.1,0.9]$	$1\%$	$2\times {10}^{-5}$	${10}^{-7}$	$21$	$0.86\%$	$0.51s$
$\Omega \setminus [0.05,0.95]$	$1\%$	${10}^{-5}$	${10}^{-7}$	$39$	$0.94\%$	$0.80s$
$\Omega \setminus [0.1,0.9]$	$0\%$	$2\times {10}^{-5}$	${10}^{-9}$	$9$	$0.54\%$	$0.31s$
$\Omega \setminus [0.1,0.9]$	$2\%$	$2\times {10}^{-5}$	$2\times {10}^{-7}$	$33$	$1.39\%$	$0.96s$
$\Omega \setminus [0.1,0.9]$	$4\%$	$2\times {10}^{-5}$	$4\times {10}^{-7}$	$128$	$3.06\%$	$2.98s$
$\Omega \setminus [0.1,0.9]$	$8\%$	$2\times {10}^{-5}$	$8\times {10}^{-7}$	$121$	$18.70\%$	$2.23s$

Equations

$${\mathcal{H}}_{p}w:={\partial }_t^{2}w-\mathrm{div}(pa\nabla w)-b\cdot \nabla w-cw\text{mbox}inQ:=\Omega \times (-T,T),$$

$$\mathcal{B}w:=pa\nabla w\cdot \nu \text{mbox}on\partial \Omega \times (-T,T).$$

$$\{\begin{array}{ll}{\mathcal{H}}_{p}u=F& \text{mbox}inQ,\\ u=u_0,{\partial }_{t}u=u_1& \text{mbox}in\Omega \times \{0\}.\end{array}$$

$$a\in (W^{3,\infty }(Q){)}^{n\times n},\exists {\kappa }_0>0\text{mbox}suchthata\xi \cdot \xi \ge {\kappa }_0|\xi {|}^2\text{mbox}in\overline{Q},\forall \xi \in {\mathbb{R}}^n,$$

$$b\in (W^{2,\infty }(-T,T;L^{\infty }(\Omega )){)}^n,c\in W^{2,\infty }(-T,T;L^{\infty }(\Omega )),p\in W^{2,\infty }(\Omega ),$$

$$F\in H^2(-T,T;L^2(\Omega )),u_0\in W^{3,\infty }(\Omega ),u_1\in H^2(\Omega ).$$

$$u(p)\in \bigcap _{k=0}^2{H}^{4-k}(-T,T;H^k(\Omega )),$$

$$\sum _{k=0}^2\Vert u(p){\Vert }_{H^{3-k}(-T,T;H^k(\Omega ))}\le M_0$$

$${\partial }_{t}a={\partial }_{t}b={\partial }_{t}c={\partial }_{t}F=u_1=0\text{mbox}in\Omega \times \{0\}.$$

$$\psi (x,t):=d(x)-\beta t^2,0<\beta <1.$$

$$\Omega (\delta ):=\{x\in \Omega ;d(x)>\delta \},Q(\delta ):=\{(x,t)\in Q;\psi (x,t)>\delta \}.$$

$$\phi (x,t):={\mathrm{e}}^{\lambda \psi (x,t)}.$$

$${\int }_{Q(\delta )}s\left(|{\partial }_{t}w{|}^2+|\nabla w{|}^2+s^2|w{|}^2\right){\mathrm{e}}^{2s\phi }\mathrm{d}x\mathrm{d}t\le C_0{\int }_{Q(\delta )}|{\mathcal{H}}_{p}w{|}^2{\mathrm{e}}^{2s\phi }\mathrm{d}x\mathrm{d}t$$

$${\int }_{Q(\delta )}\sum _{i,j=1}^n|{\partial }_i{\partial }_{j}w{|}^2{\mathrm{e}}^{2s\phi }\mathrm{d}x\mathrm{d}t$$

$${\int }_{Q(\delta )}|{\partial }_t^{2}w{|}^2{\mathrm{e}}^{2s\phi }\mathrm{d}x\mathrm{d}t$$

$${\partial }_{t}w\in H^2(Q),\mathrm{supp}w\subset \overline{Q(\delta )},{\partial }_t^{k}w(\cdot ,\pm T)=0,k=0,1,2.$$

$$p(x)a_{ij}(x,t)=\{\begin{array}{ll}{a}_0(x,t),& i=j,\\ 0,& i\ne j,\end{array}$$

$$d(x)=x_1^2+\cdots +x_{n-1}^2+(x_n-1{)}^2,$$

$$a_0>\sqrt{\beta },\beta \left(2-\frac{t{\partial }_t{a}_0}{a_0}+\frac{2|t\nabla a_0|}{\sqrt{a_0}}\right)<2a_0-\nabla a_0\cdot \nabla d+{\partial }_n{a}_0\text{mbox}in\overline{Q},$$

$$\overline{\Omega ({\delta }_0)}\subset \Omega \cup \Gamma ,\partial \omega \supset \Gamma ,T^2>\Vert d{\Vert }_{C(\overline{\Omega })},$$

$$\mathcal{U}:=\{p\in W^{2,\infty }(\Omega );\Vert p{\Vert }_{H^1(\Omega )}\le M_1,p\ge {\kappa }_1\text{mbox}in\overline{\Omega },p=h_0,{\partial }_{\nu }p=h_1\text{mbox}on\Gamma \}.$$

$$|a(\cdot ,0)\nabla u_0\cdot \nabla d|\ge {\kappa }_2\text{mbox}a.e.in\Omega .$$

$$\Vert p-q{\Vert }_{H^1(\Omega (\delta ))}\le CD+C{M}^{1-\theta }{D}^{\theta },$$

$$D:=\{\begin{array}{rlrl}& \sum _{k=0}^2\Vert u(p)-u(q){\Vert }_{H^{4-k}(-T,T;H^{k-1/2}(\Gamma ))}& & \\ & +\Vert \mathcal{B}u(p)-\mathcal{B}u(q){\Vert }_{H^2(-T,T;H^{1/2}(\Gamma ))},& & \text{mbox}Type(\mathrm{I}),\\ & \sum _{k=0}^2\Vert u(p)-u(q){\Vert }_{H^{4-k}(-T,T;H^k(\omega ))},& & \text{mbox}Type(\mathrm{II}).\end{array}$$

$$\mathrm{div}((p-q)a\nabla u(q))$$

$$\mathcal{B}u(p)-\mathcal{B}u(q)=h_{0}a\nabla (u(p)-u(q))\cdot \nu =\mathcal{B}(u(p)-u(q))\text{mbox}on\Gamma \times (-T,T)$$

$$f\in H^2(\Omega ),f=\nabla f=0\text{mbox}on\Gamma ,\Vert f{\Vert }_{H^1(\Omega )}\le 2M_1\le 2M.$$

$$\{\begin{array}{ll}\mathcal{H}v=\mathcal{L}f:=\mathrm{div}(fa\nabla u(q))& \text{mbox}inQ=\Omega \times (-T,T),\\ v={\partial }_{t}v=0& \text{mbox}in\Omega \times \{0\},\end{array}$$

$$v\in \bigcap _{k=0}^2{H}^{4-k}(-T,T;H^k(\Omega )),\sum _{k=0}^2\Vert v{\Vert }_{H^{3-k}(-T,T;H^k(\Omega ))}\le 2M_0\le 2M.$$

$$\{\begin{array}{ll}\tilde{v}=v,\mathcal{B}\tilde{v}=\mathcal{B}v\text{mbox}on\Gamma \times (-T,T),& \text{mbox}Type(I),\\ \tilde{v}=v\text{mbox}in\omega \times (-T,T),& \text{mbox}Type(II),\end{array}$$

$$\sum _{k=0}^2\Vert \tilde{v}{\Vert }_{H^{3-k}(-T,T;H^k(\Omega ))}\le CM,\sum _{k=0}^2\Vert \tilde{v}{\Vert }_{H^{4-k}(-T,T;H^k(\Omega ))}\le CD.$$

$$y_{\ell }:={\partial }_t^{\ell }(v-\tilde{v}),\ell =0,1,2.$$

$$\mathcal{H}{y}_0=\mathcal{L}f-\mathcal{H}\tilde{v},\mathcal{H}{y}_1={\mathcal{A}}^{\prime }{y}_0+{\mathcal{L}}^{\prime }f-{\partial }_t(\mathcal{H}\tilde{v}),$$

Full text

[TABLE]

**AMS Subject Classifications ** 35L20, 35R30, 35A23, 65M32, 65F22

††footnotetext:

Manuscript last updated: .

†

School of Mathematical Sciences, Fudan University, No.​ 220 Handan Road, Shanghai 200433, China.

‡

Graduate School of Mathematical Sciences, The University of Tokyo, 3-8-1 Komaba, Meguro-ku, Tokyo 153-8914, Japan.

∗

Corresponding author. E-mail: ykliu@ms.u-tokyo.ac.jp

@hangfrom@seccntformat

sectionIntroduction

@xsect

2.3ex plus.2ex

Let $T>0$ and $\Omega\subset\mathbb{R}^{n}$ ($n=1,2,3$) be an open bounded domain whose boundary $\partial\Omega$ is of $C^{3}$ class. Let $\nu=(\nu_{1},\ldots,\nu_{n})$ be the outward unit normal vector to $\partial\Omega$, and denote the normal derivative on $\partial\Omega$ by $\partial_{\nu}w:=\nabla w\cdot\nu$. Introduce the hyperbolic operator

[TABLE]

where $a=(a_{ij}(x,t))_{1\leq i,j\leq n}$ is a symmetric matrix, $b=(b_{i}(x,t))_{1\leq i\leq n}$ is a vector, and $p=p(x),c=c(x,t)$ are scalar functions. For the non-degeneracy, we assume that $p$ is strictly positive in $\overline{\Omega}\,$, and $a$ is strictly positive-definite in $\overline{Q}\,$. For later use, we denote the normal derivative associated with the second order coefficient $p\,a$ as

[TABLE]

In this paper, we consider the initial value problem for a hyperbolic equation

[TABLE]

From the physical point of view, (1.3) models the general acoustic wave in a highly anisotropic medium depending on both space and time. In a recent paper [22], it turns out that (1.3) also describes the one-dimensional time cone model for phase transformation. Note that for the well-posedness of the forward problem, we should attach (1.3) with a boundary condition. However, for the inverse problem proposed later, it suffices to access e.g.​ only the partial boundary value, which is regarded as a part of observation data. Therefore, for the moment we do not include any boundary condition in (1.3). To emphasize the dependency, throughout this paper we will write any solution satisfying (1.3) with the coefficient $p$ as $u(p)$. Detailed assumptions on the coefficients $p,a,b,c,F,u_{0},u_{1}$ involved in (1.3) will be given later in Section Theoretical Stability in Coefficient Inverse Problems for General Hyperbolic Equations with Numerical Reconstruction.

This paper is mainly concerned with the following coefficient inverse problem on the determination of the spatial component $p$ in the principal part of the hyperbolic operator (1.1).

Type (I)*  the partial boundary observation of $u(p)$ and $\mathcal{B}u(p)$ on $\Gamma\times(-T,T)$, or*

Type (II)*  the partial interior observation of $u(p)$ in $\omega\times(-T,T)$.*

In the formulation of Problem 1.1, the second order coefficient $p(x)a(x,t)$ takes the form of incompletely separated variables, where the unknown spatial component $p(x)$ contributes an important part in the wave propagation speed. In view of the acoustic equation, Problem 1.1 stands for the identification of the bulk modulus, which is of practical significance. Hence, we will not only investigate the theoretical aspect of Problem 1.1 due to our interest in mathematics, but also consider the reconstruction method to solve $p(x)$ numerically.

In retrospect, researches on coefficient inverse problems for hyperbolic equations started soon after the pioneering work of Bukhgeim and Klibanov [5] which discovered the potential of Carleman estimates. We refer e.g.​ to [17, 18, 11] for some early results mainly on the uniqueness. Around 2000s, [24, 8, 9] established the global Lipschitz stability for determining the zeroth order coefficient $c(x)$ in $(\partial_{t}^{2}-\triangle-c)u=0$ by the same types of data in Problem 1.1. For Problem 1.1 with $a=I_{n\times n}$ and $b=c\equiv 0$ in (1.1), Imanuvilov and Yamamoto [10] employed an $H^{-1}$ Carleman estimate to obtain the global Hölder stability by partial boundary observation. Later, this result was improved to Lipschitz in Bellassoued and Yamamoto [4] with time-independent coefficients, i.e., $a=(a_{ij}(x))$ in (1.1). For other references on this direction, see also [3, 12]. It reveals that all of the above literature more or less imposes some geometrical conditions because of the finite wave propagation speed. Meanwhile, the above results mostly rely on a linearization approach, which reduces the problem to a corresponding inverse source problem. In addition, we emphasize the difference between e.g.​ $\mathrm{div}(p\nabla u)$ and $p\,\triangle u$, where the former is more physical and the latter is technically easier (see [2]). However, there seems no publication treating Problem 1.1 in the case of time-dependent principal parts due to the essential difficulty. Recently, Jiang, Liu and Yamamoto [15] established the second order Carleman estimate for (1.1) and proved the local Hölder stability for a related inverse source problem. Motivated by [15], we first attempt to generalize the result in [10] with a more general hyperbolic operator $\mathcal{H}_{p}$ with time-dependent coefficients, which is the first focus of this article.

Simultaneously, we have witnessed the recent applications of the iterative thresholding algorithm to inverse problems for partial differential equations. For the abstract formulation and convergence analysis of the algorithm, we refer to [6, 7, 23]. Attracted by its efficiency and robustness in many image processing problems, [13] first utilized the iterative thresholding algorithm to solve inverse problems for elliptic and parabolic equations. In [20, 15, 16], similar iteration methods were implemented to treat inverse source problems for hyperbolic-type equations with different types of observation data. Following the same line, we also attempt to develop the same class of iteration method to solve Problem 1.1 numerically, which is the second focus of this article. However, we should realize the underlying ill-posedness as well as the nonlinearity of Problem 1.1, which differs considerably from inverse source problems. For the numerical reconstruction of a time-dependent principal coefficient, we refer to [21].

Based on the newly established Carleman estimate in [15], we first prove the local Hölder stability of Problem 1.1 for both types of observation data (see Theorem 2.3). Due to the lack of an $H^{-1}$ Carleman estimate for $\mathcal{H}_{p}$ as that in [10], we have to argue in an alternative way to evaluate the $H^{1}$-norm of the difference in unknown functions, which results from the divergence form $\mathrm{div}(p\,a\nabla u)$ in (1.3). For the numerical reconstruction, we reformulate Problem 1.1 as a minimization problem with the Tikhonov regularization penalizing the $L^{2}$-norm of $\nabla p$. Deriving the variational equation of the minimizer, we arrive at a novel iteration method which needs to solve a Poisson equation at each step.

The rest of this article is organized as follows. Preparing the necessary ingredients including the key Carleman estimates, in Section Theoretical Stability in Coefficient Inverse Problems for General Hyperbolic Equations with Numerical Reconstruction we state the main result on the theoretical stability of Problem 1.1. Then we give the proof of the main result in Section Theoretical Stability in Coefficient Inverse Problems for General Hyperbolic Equations with Numerical Reconstruction. Next, Section Theoretical Stability in Coefficient Inverse Problems for General Hyperbolic Equations with Numerical Reconstruction is devoted to the derivation of an iteration method for Problem 1.1, followed by Section Theoretical Stability in Coefficient Inverse Problems for General Hyperbolic Equations with Numerical Reconstruction illustrating several one-dimensional numerical examples. Finally, we provide some concluding remarks in Section Theoretical Stability in Coefficient Inverse Problems for General Hyperbolic Equations with Numerical Reconstruction.

@hangfrom@seccntformat

sectionPreliminaries and Main Results

@xsect

2.3ex plus.2ex

In this section, we start from the general settings and assumptions concerning the governing equation (1.3), and prepare the key Carleman estimates for the hyperbolic operator (1.1). Then we state the main result of this paper, which gives the stability estimate for Problem 1.1.

Throughout this paper, we write $\partial_{i}=\frac{\partial}{\partial x_{i}}$ and $\partial_{i}\partial_{j}=\frac{\partial^{2}}{\partial x_{i}\partial x_{j}}$ ($1\leq i,j\leq n$) for the partial derivatives in space. We recall the definition $Q:=\Omega\times(-T,T)$ and the notations $W^{k,\infty}(-T,T;W^{\ell,\infty}(\Omega))$, $H^{k}(\Omega)$, $H^{k-1/2}(\partial\Omega)$, etc.​ ($k,\ell=0,1,\ldots$) for the usual Sobolev spaces (see Adams [1]). For the various coefficients appearing in (1.3), we basically make the following assumptions:

[TABLE]

With some suitably given boundary condition and the compatibility condition, it is well known that the initial-boundary value problem governed by (1.3) admits a unique solution $u(p)$ which depends continuously on the involved coefficients (see e.g.​ [19, 15]). In order to prove the theoretical stability for Problem 1.1, we have to assume

[TABLE]

which satisfies the a priori estimate

[TABLE]

with a constant $M_{0}>0$. As was mentioned in [15], the smoothness of coefficients assumed in (2.1) may not guarantee the regularity in (2.2). Instead, we only understand (2.1)–(2.3) as the minimum necessary assumptions for the stability.

If the initial value problem (1.3) is only formulated in $\Omega\times(0,T)$, we should additionally impose

[TABLE]

Then we can perform even reflections of $a,b,c,F$ and the solution $u(p)$ with respect to $t$, so that they preserve the same regularity in $Q$ as assumed in (2.1) and (2.2). To circumvent the above technical assumption, we simply consider (1.3) in $Q=\Omega\times(-T,T)$ without loss of generality.

Next, we recall the Carleman estimates for the hyperbolic operator $\mathcal{H}_{p}$ defined in (1.1), which play an essential role in both the statement and the proof of the main result. Similarly to [15], we pick $d\in C^{2}(\overline{\Omega})$ such that $d>0$ in $\overline{\Omega}$ , and set

[TABLE]

By $d(x)$ and $\psi(x,t)$, we introduce the level sets with a parameter $\delta\geq 0$ as

[TABLE]

With a sufficiently large parameter $\lambda>0$, we define the weight function as

[TABLE]

We collect the Carleman estimates concerning (1.1) in the following lemma.

[TABLE]

holds for all $s\geq s_{0}$ and $w\in H^{2}(Q)$ satisfying $\mathrm{supp}\,w\subset\overline{Q(\delta)}\,$, where $Q(\delta)\ (\delta\geq 0)$ is defined in (2.5). Then for any $\delta\geq 0$, there exist constants $s_{1}>0$ and $C_{1}>0$ such that

[TABLE]

holds for all $s\geq s_{1}$ and $w\in H^{2}(Q)$ satisfying

[TABLE]

[TABLE]

under certain geometrical assumptions and the choice

[TABLE]

a sufficient condition for (2.7) can be (see Isakov [12, Theorem 3.4.3])

[TABLE]

where $\beta\in(0,1)$ is the parameter in (2.4). On the other hand, the second order Carleman estimates (2.8)–(2.9) was established in [15, Lemma 3.1] on the basis of (2.7).

Now we turn our attention to Problem 1.1. For the subboundary $\Gamma\subset\partial\Omega$, the subdomain $\omega\subset\Omega$ and $T>0$ where the observation is taken, we make the following assumption. Given a constant $\delta_{0}\geq 0$ and a function $d\in C^{2}(\overline{\Omega})$ such that $d>0$ in $\overline{\Omega}\,$, we assume

[TABLE]

where $\Omega(\delta_{0})$ is the level set defined in (2.5). Although the above assumption looks technical, it generalizes the similar assumption in treating the same inverse problems for simpler wave equations (see [9, 20]), where the function $d$ was simply taken as $d(x)=|x-x_{0}|^{2}$ with $x_{0}\not\in\overline{\Omega}\,$.

Finally, for the unknown function $p(x)$ to be determined, we define an admissible set $\mathcal{U}$ in the following way. For given constants $M_{1}>0,\kappa_{1}>0$ and given functions $h_{0}\in W^{2,\infty}(\Gamma)$ and $h_{1}\in W^{1,\infty}(\Gamma)$ where $\Gamma$ satisfies (2.10), we restrict $p$ in

[TABLE]

Collecting all necessary ingredients, now we are in a position to state the main theoretical result, which answers the stability issue of Problem 1.1.

[TABLE]

Then for any $\delta>\delta_{0}$, there exist constants $C>0$ and $\theta\in(0,1)$ depending on $a,b,c,d,u_{0},\delta,T,\Gamma$ or $\omega$ such that

[TABLE]

where $\Omega(\delta)$ is defined in (2.5), $M:=\max\{M_{0},M_{1}\}$ with the a priori bounds $M_{0},M_{1}$ in (2.3) and (2.11) respectively, and

[TABLE]

(1)  At a first glance, the statement of Theorem 2.3 resembles that of [15, Theorem 2.3] to a large extent. Indeed, as the majority of existing literature did, here we also linearize the system to consider the governing equation of the difference $u(p)-u(q)$, which partially reduces Problem 1.1 to an inverse source problem. Even though, Theorem 2.3 is nontrivial because one should deal with the source term

[TABLE]

as that in [10]. However, in our case there seems no $H^{-1}$ Carleman estimate and thus we cannot evaluate (2.15) simply by $\|p-q\|_{L^{2}(\Omega)}$. As a result, we should treat (2.15) in the $L^{2}$ sense, so that we have to argue more to dominate the $H^{1}$-norm of $p-q$. Hence, although (2.13) gives an $H^{1}$ estimate which looks stronger than that in [10], the fact is that we fail to provide an $L^{2}$ one instead.

(2)  In Theorem 2.3, the unknown coefficient $p$ is restricted in the admissible set $\mathcal{U}$. According to the definition (2.11), it means that the partial Cauchy data of $p$ is known on the subboundary $\Gamma$, which seems reasonable because the observation is taken there. Moreover, the boundary operator $\mathcal{B}$ defined in (1.2) becomes linear on $\Gamma\times(-T,T)$ within the admissible set $\mathcal{U}$, which allows us to write

[TABLE]

in the observation data (2.14).

(3)  We explain the geometry of the observation data in the assumption (2.10). Since $d$ is strictly positive in $\overline{\Omega}\,$, by definition we know $\Omega(\delta)=\Omega$ for $0\leq\delta<\min_{x\in\overline{\Omega}}d(x)$. Therefore, if we fix $\delta_{0}=0$, then (2.10) requires $\Gamma=\partial\Omega$ and $\partial\omega\supset\partial\Omega$, that is, the observation should be taken on the whole boundary. Nevertheless, in this case (2.13) becomes a global Hölder estimate with sufficiently small $\delta>0$. On the other hand, if one attempt to shrink the observation region, then the estimate (2.13) tends to be local. On the opposite extreme, since $\Omega(\delta)=\emptyset$ for any $\delta\geq\|d\|_{C(\overline{\Omega})}$, our result becomes trivial with large $\delta$. As a result, the choice of $d$ in the weight function plays a delicate role in the balance between the cost of measurements and the stability. The optimal choice of $d$ seems to be an interesting topic, but we will not discuss it in this paper.

@hangfrom@seccntformat

sectionProof of Theorem 2.3

@xsect

2.3ex plus.2ex

In this section, we prove Theorem 2.3 under the same assumptions therein. The proof basically follows the same line as that in [15], which relies heavily on the Carleman estimates in Lemma 2.1. However, many details are different from each other, and in our problem we should especially argue more for the $H^{1}$ estimate in (2.13).

Henceforth, by $C>0$ we denote generic constants independent of the parameter $s>0$ in Carleman estimates and the observation data $D>0$, which may change line by line.

For clearness, we divide the proof into five steps.

Step 1  We start from the general preparation of notations and auxiliary functions. For simplicity, we abbreviate the hyperbolic operator $\mathcal{H}_{p}$ as $\mathcal{H}$ throughout this section. Fixing any $p,q\in\mathcal{U}$, we set $f:=p-q$. Then it follows from the definition (2.11) that $f\in H^{2}(\Omega)$ and $f=\partial_{\nu}f=0$ on $\Gamma$. Together with the a priori bound, we conclude

[TABLE]

Correspondingly, we set $v:=u(p)-u(q)$. According to (1.3), (2.2) and (2.3), $v$ satisfies

[TABLE]

We note that $\mathcal{L}$ is a first order differential operator with respect to $f$.

By the Sobolev extension theorem, there exists $\widetilde{v}\in\bigcap_{k=0}^{2}H^{4-k}(-T,T;H^{k}(\Omega))$ such that

[TABLE]

For later use, we further introduce

[TABLE]

Then it is readily seen that $y_{\ell}\in\bigcap_{k=0}^{2}H^{4-\ell-k}(-T,T;H^{k}(\Omega))$ ($\ell=0,1,2$), and $y_{0},y_{1},y_{2}$ satisfy

[TABLE]

where we define

[TABLE]

Moreover, it follows from (3.4), (3.3) and (3.5) that for $\ell=0,1,2$,

[TABLE]

Note that $y_{\ell}$ ($\ell=0,1,2$) may not vanish outside the observation region, which prevents us from applying the Carleman estimates in Lemma 2.1. To this end, it is necessary to introduce a cutoff function as that in [15]. Recall the constant $\delta_{0}\geq 0$ and the function $d\in C^{2}(\overline{\Omega})$ given in advance. For any $\delta>\delta_{0}$, we fix a constant $\delta_{1}\in(\delta_{0},\delta)$, e.g., $\delta_{1}=\frac{\delta_{0}+\delta}{2}$. With the level set $Q(\delta)$ defined in (2.5), we define $\mu\in C^{\infty}(\overline{Q})$ such that

[TABLE]

Especially, we see that the condition (3.10) is possible by the definition (2.4) of $\psi$. Meanwhile, owing to the assumption $T^{2}>\|d\|_{C(\overline{\Omega})}$ in (2.10), we can choose $\beta\in(0,1)$ such that $\beta T^{2}>\|d\|_{C(\overline{\Omega})}$, indicating $\psi(\,\cdot\,,\pm T)<0$ in $\overline{\Omega}\,$. Together with the assumption $\overline{\Omega(\delta_{0})}\subset\Omega\cup\Gamma$ in (2.10), we conclude

[TABLE]

Now we further set

[TABLE]

Employing (3.7), (3.11) and the assumption $\partial\omega\supset\Gamma$ in (2.10), we obtain for $\ell=0,1,2$ that

[TABLE]

On the other hand, it is obvious that $z_{\ell}$ share the same regularity as that of $y_{\ell}$, i.e., $z_{\ell}\in H^{2}(Q)$ ($\ell=0,1,2$) and especially $\partial_{t}z_{0},\partial_{t}z_{1}\in H^{2}(Q)$. By (3.6) and direct calculations, we deduce

[TABLE]

where

[TABLE]

Similarly to [15], it turns out that $[\mathcal{H},\mu]$, $[\mathcal{A}^{\prime},\mu]$ and $[\mathcal{A}^{\prime\prime},\mu]$ are all first order differential operators which only involve derivatives of $\mu$. Then the definition (3.9) of $\mu$ implies

[TABLE]

Step 2  Now that $z_{\ell}$ ($\ell=0,1,2$) satisfy (3.12), we are well prepared to apply the Carleman estimates in Lemma 2.1 to $z_{\ell}$. We will utilize estimates (3.1), (3.5) and (3.8) repeatedly throughout the proof, and take advantage of the properties of the cutoff function $\mu$ as well as the weight function $\varphi$.

Applying the first order Carleman estimate (2.7) in Lemma 2.1 to (3.15), we have

[TABLE]

Our aim in this step is to give an estimate for $\|F_{2}\,\mathrm{e}^{s\varphi}\|_{L^{2}(Q(\delta_{0}))}$. To this end, we apply (2.7)–(2.8) in Lemma 2.1 to (3.13)–(3.14) to obtain for $\ell=0,1$ that

[TABLE]

To proceed, we should estimate $\|(\partial_{t}^{k}F_{\ell})\mathrm{e}^{s\varphi}\|_{L^{2}(Q(\delta_{0}))}$ for $k,\ell=0,1$. First, we combine the properties (3.9)–(3.10) with (3.16), (3.5) and (3.8) to dominate

[TABLE]

where we used $\varphi=\mathrm{e}^{\lambda\psi}\leq\mathrm{e}^{\lambda\delta_{1}}=:\eta_{1}$ in $\overline{Q(\delta_{0})}\setminus Q(\delta_{1})$ by the definition (2.5). Similarly, by

[TABLE]

we further estimate

[TABLE]

where we turned to the a priori estimate (3.1) to treat the term $(\partial_{t}\mu)\mathcal{L}f$. Hence, substituting (3.19)–(3.20) into (3.18) with $\ell=0$ yields

[TABLE]

Next, we deal with $F_{1}$ defined in (3.14). Noting the fact that $\mathcal{A}^{\prime}$ is a second order differential operator in space, we apply (3.21), (3.5) and (3.8) to derive

[TABLE]

for all $s\gg 1$. To bound $\|(\partial_{t}F_{1})\mathrm{e}^{s\varphi}\|_{L^{2}(Q(\delta_{0}))}$, we calculate

[TABLE]

Similarly to the treatment for $\partial_{t}F_{0}$, we employ (3.21) again to estimate

[TABLE]

Substituting (3.22)–(3.23) into (3.18) with $\ell=1$, we deduce

[TABLE]

Choosing $s>0$ sufficiently large, we can absorb the first term on the right-hand side into the left-hand side and conclude

[TABLE]

Using (3.5), (3.8) and (3.16) again, we apply (3.21) and (3.24) to the definition (3.15) of $F_{2}$ and arrive at the estimate

[TABLE]

for all $s\gg 1$.

Step 3  We further introduce

[TABLE]

Then simple calculations yield

[TABLE]

We attempt to give upper and lower estimates for

[TABLE]

First, since $\mathcal{P}z_{2}=F_{2}+b\cdot\nabla z_{2}+c\,z_{2}$, the application of (3.17) immediately gives

[TABLE]

for all $s\gg 1$. Noticing $\mathrm{supp}\,z=\mathrm{supp}\,z_{2}\subset\overline{Q(\delta_{0})}$ and using (3.17) again, we employ (3.26)–(3.27) to estimate $I_{0}$ from above as

[TABLE]

On the other hand, it follows from (3.12) that

[TABLE]

which allows us to perform integration by parts to estimate $I_{0}$ from below as

[TABLE]

where we applied the lower bounds in (2.1) and (2.11) to obtain (3.29). Utilizing (3.26) and (3.17) again, we further estimate

[TABLE]

Combining the above inequality with (3.29), (3.28) and (3.25), we obtain for all $s\gg 1$ that

[TABLE]

Next, we shall relate the above estimate of $\nabla\widetilde{z}(\,\cdot\,,0)$ with that of $f$. In fact, by the definition of $\widetilde{z}$, we take $t=0$ in (3.2) and find

[TABLE]

where

[TABLE]

Hence, by $\mathcal{L}_{0}(\mu_{0}f)=\mu_{0}\,\mathcal{L}_{0}f+f\,a(\,\cdot\,,0)\nabla u_{0}\cdot\nabla\mu_{0}$, we obtain

[TABLE]

Recalling the definition (2.5) of $\Omega(\delta)$, we see $\mathrm{supp}(\nabla\mu_{0})\subset\overline{\Omega(\delta_{0})}\setminus\Omega(\delta_{1})$ and $\varphi_{0}\leq\eta_{1}$ in $\overline{\Omega(\delta_{0})}\setminus\Omega(\delta_{1})$. By the same argument as before, we apply (3.30) to bound

[TABLE]

where we used the Sobolev embedding $C[-T,T]\subset H^{1}(-T,T)$ and (3.5) to dominate

[TABLE]

Step 4  In order to relate the above estimate (3.33) with the $H^{1}$-norm of $f$, we need the following Carleman estimate for the first order differential operator $\mathcal{L}_{0}$ in (3.32).

[TABLE]

holds for all $s\geq s_{2}$ and $g\in H_{0}^{2}(\Omega)$.

[TABLE]

In fact, by setting $\widetilde{g}:=\mathrm{e}^{s\varphi_{0}}g$, we calculate

[TABLE]

By the assumption (2.12) and $d>0$ in $\overline{\Omega}\,$, we obtain the lower bound

[TABLE]

Then we can estimate

[TABLE]

where

[TABLE]

By $g\in H_{0}^{1}(\Omega)$, we perform integration by parts to treat $I_{1}$ as

[TABLE]

According to the assumption (2.1), both $\|a(\,\cdot\,,0)\|_{W^{1,\infty}(\Omega)}$ and $\|u_{0}\|_{W^{2,\infty}(\Omega)}$ are bounded, which indicates $|I_{1}|\leq C\|\widetilde{g}\|_{L^{2}(\Omega)}^{2}$ and thus

[TABLE]

Therefore, there exists a constant $s_{3}>0$ such that the right-hand side is strictly positive for all $s\geq s_{3}$, which implies (3.34).

Next, since $g\in H_{0}^{2}(\Omega)$ gives $\nabla g\in(H_{0}^{1}(\Omega))^{n}$, we apply (3.34) to $\nabla g$ to obtain

[TABLE]

To further estimate the right-hand side of (3.35), we calculate

[TABLE]

Since (2.1) also gives the boundedness of $\|a(\,\cdot\,,0)\|_{W^{2,\infty}(\Omega)}$ and $\|u_{0}\|_{W^{3,\infty}(\Omega)}$, we estimate

[TABLE]

Substituting the above inequality into (3.35), we obtain

[TABLE]

which, together with (3.34), yields

[TABLE]

Then there exists a constant $s_{2}>s_{3}$ such that

[TABLE]

Finally, it follows from $\nabla(e^{s\varphi_{0}}\mathcal{L}_{0}g)=\mathrm{e}^{s\varphi_{0}}\{\nabla(\mathcal{L}_{0}g)+s(\mathcal{L}_{0}g)\nabla\varphi_{0}\}$ that

[TABLE]

Applying the above estimate to (3.36), we obtain

[TABLE]

Since $\mathcal{L}_{0}$ is a first order differential operator, it reveals that $\mathrm{e}^{s\varphi_{0}}\mathcal{L}_{0}g\in H_{0}^{1}(\Omega)$, which allows us to apply the Poincaré inequality to conclude

[TABLE]

This completes the proof of Lemma 3.1. ∎

Step 5  We complete the proof of Theorem 2.3 in this step. By (2.10) and the definition of $\mu_{0}$, we see $\mu_{0}=\nabla\mu_{0}=0$ on $\partial\Omega(\delta_{0})\supset\partial\Omega\setminus\Gamma$. Together with (3.1), we obtain $\mu_{0}f\in H_{0}^{2}(\Omega)$, which allows us to take advantage of (3.33) and Lemma 3.1 with $g=\mu_{0}f$ to derive

[TABLE]

Substituting

[TABLE]

into the above inequality, we arrive at

[TABLE]

Owing to the choice of the weight function, we can employ Lebesgue’s dominated convergence theorem to absorb the first term in the right-hand side of the above estimate in the sense that

[TABLE]

Consequently, we have

[TABLE]

Recalling the facts that $\delta>\delta_{1}>\delta_{0}$ and $\mu_{0}\equiv 1$, $\varphi_{0}\geq\mathrm{e}^{\lambda\delta}=:\eta$ in $\Omega(\delta)$, we further estimate the left-hand side of the above inequality from below to deduce

[TABLE]

Since $\eta>\eta_{1}$ and $s^{2}\leq\mathrm{e}^{(\eta-\eta_{1})s}$ for all $s\gg 1$, we conclude

[TABLE]

Finally, discussing the two cases $D\geq M$ and $D<M$ as that in [15], we eventually obtain (2.13) and finish the proof.

@hangfrom@seccntformat

sectionIteration Method for Numerical Reconstruction

@xsect

2.3ex plus.2ex

As the theoretical stability is guaranteed by Theorem 2.3, in this section we study Problem 1.1 from the numerical viewpoint and aim at the derivation of an iteration method. Basically, the derivation is parallel to its counterpart in [15], where the corresponding inverse source problem was investigated. However, for Problem 1.1 we shall pay special attention to the nonlinearity and the ill-posedness in the recovery of the second order coefficient.

Instead of the general governing equation (1.3) whose coefficients are all assumed to be time-dependent, throughout this section we consider the initial-boundary value problem for a hyperbolic equation with the homogeneous Neumann boundary condition

[TABLE]

Here the boundary condition is simplified from $p\nabla u\cdot\nu=0$ on $\partial\Omega\times(0,T)$ because $p$ is assumed as strictly positive on $\overline{\Omega}$ for the non-degeneracy. We limit the numerical treatment to the time-independent case (4.1) not only for its simplicity, but also due to the belief that the ill-posedness are essentially the same.

For later use, we recall the classical theory on the well-posedness of (4.1).

[TABLE]

and the $k$th order compatibility condition is satisfied on $\partial\Omega\times\{0\}$. Then there exists a unique solution $u\in C([0,T];H^{k}(\Omega))\cap C^{1}([0,T];H^{k-1}(\Omega))$ to (4.1). Moreover, there exists a constant $C>0$ depending on $p,\Omega,T$ such that

[TABLE]

Henceforth, we basically assume

[TABLE]

and the second order compatibility condition is satisfied on $\partial\Omega\times\{0\}$. Then according to Lemma 4.1 with $k=2$, problem (4.1) admits a unique solution $u\in C([0,T];H^{2}(\Omega))\cap C^{1}([0,T];H^{1}(\Omega))$. As before, we still denote the unique solution to (4.1) as $u(p)$.

To deal with Problem 1.1 from the numerical aspect, we restrict ourselves to the following situation. Regarding the observation region, we consider the partial interior observation of $u(p)$ in $\omega\times(0,T)$ with a subdomain $\omega\subset\Omega$ satisfying $\partial\omega\supset\partial\Omega$. In other words, we require $\omega$ to cover the whole boundary $\partial\Omega$, which is special in Type (II) of Problem 1.1. In fact, although there seems no difference between boundary and interior measurements in the theoretical stability, the latter is definitely more informative and suitable for the numerical implementation. On the other hand, it follows from Remark 2.4 that $\partial\omega\supset\partial\Omega$ is a sufficient condition for the global stability of Problem 1.1, which is desirable for determining $p$ in the whole domain $\Omega$.

In accordance with the above setting, we restrict the unknown function $p$ in

[TABLE]

with given constants $M_{1}>0,\kappa_{1}>0$ and a given function $h_{0}\in W^{2,\infty}(\partial\Omega)$. Compared with the admissible set $\mathcal{U}$ defined in (2.11) for the theoretical stability, here we remove the restriction of $\partial_{\nu}p$ on $\partial\Omega$. Nevertheless, we still require that $p$ is known on the whole boundary due to the key assumption (2.10) for the stability. We refer to [10] for the same type of admissible sets as $\mathcal{U}_{1}$.

In practice, we are given the noisy observation data $u^{\delta}\in L^{2}(\omega\times(0,T))$ such that

[TABLE]

where $p_{\mathrm{true}}\in\mathcal{U}_{1}$ and $\delta>0$ stand for the true solution and the noise level respectively. Now we are well prepared to recast Problem 1.1 into a minimization problem with the Tikhonov regularization

[TABLE]

where $\alpha>0$ denotes the regularization parameter. Unlike the formulation in [20, 15], here we penalize the $L^{2}$-norm of $\nabla p$ because one can expect certain smoothness of $p$ as the second order coefficient. Meanwhile, there is no need to penalize the $H^{1}$-norm of $p$ due to the boundary condition $p=h_{0}$ on $\partial\Omega$.

As usual, we shall compute the Fréchet derivative of $J(p)$ in order to characterize its possible minimizer $p_{*}$. For arbitrarily fixed $p\in\mathcal{U}_{1}$, we may choose any $\widetilde{p}\in W^{2,\infty}(\Omega)$ such that

[TABLE]

holds for all sufficiently small $\varepsilon>0$. By (4.4), we directly calculate

[TABLE]

In order to pass $\varepsilon\downarrow 0$ in (4.6), we need the following technical lemma.

[TABLE]

where $w_{0}$ satisfies

[TABLE]

[TABLE]

Since $p+\varepsilon\,\widetilde{p}$ lies in the $\varepsilon$-neighborhood of $p$ by (4.5), it follows from Lemma 4.1 with $k=2$ that there exists a constant $C_{2}>0$ such that

[TABLE]

holds uniformly for all sufficiently small $\varepsilon\geq 0$. This indicates

[TABLE]

uniformly for all sufficiently small $\varepsilon\geq 0$, where we have $\|\widetilde{p}\,\|_{W^{1,\infty}(\Omega)}\leq 1$ by (4.5). Applying Lemma 4.1 with $k=1$ to (4.8), we obtain

[TABLE]

In the same manner, we further set $w_{\varepsilon}:=\varepsilon^{-1}v_{\varepsilon}$ and manipulate (4.7)–(4.8) to find

[TABLE]

Then we employ (4.9) and Lemma 4.1 with $k=0$ to conclude

[TABLE]

which finishes the proof. ∎

Now that the convergence is guaranteed by the above lemma, we can pass $\varepsilon\downarrow 0$ in (4.6) to deduce

[TABLE]

where $\chi_{\omega}$ denotes the characteristic function of $\omega$, and we utilized $\widetilde{p}=0$ on $\partial\Omega$ to obtain (4.10) by integration by parts.

In order to derive the explicit form of $J^{\prime}(p)$, we should further transform the first term on the right-hand side of (4.10). To this end, we follow the same line as that in [20, 15] to introduce the backward problem

[TABLE]

To clarify the dependency, we also denote the solution to (4.11) as $z(p)$. Since $\chi_{\omega}\left(u(p)-u^{\delta}\right)\in L^{2}(\Omega\times(0,T))$, Lemma 4.1 gives $z(p)\in H^{1}(\Omega\times(0,T))$. On the other hand, it can be inferred from the proof of Lemma 4.2 that the solution of (4.7) satisfies $w_{0}\in H^{1}(\Omega\times(0,T))$ and $w_{0}|_{t=0}=0$. Hence, in view of the weak solution of hyperbolic equations, we can regard $w_{0}$ and $z(p)$ as mutual test functions of each other, so that we can further treat

[TABLE]

Substituting the above identity into (4.10), we arrive at

[TABLE]

Since $\widetilde{p}$ was taken arbitrarily which satisfies (4.5), this suggests a characterization of the minimizer to the problem (4.4).

[TABLE]

where $u(p_{*})$ and $z(p_{*})$ solve the forward system (4.1) and the backward one (4.11) with the coefficient $p_{*}$, respectively.

On the basis of (4.12), we design the following iteration scheme

[TABLE]

where $K>0$ is a tuning parameter. In other words, given the result $p_{m}$ of the previous step, we have to solve the forward system (4.1), the backward system (4.11) and the boundary value problem (4.13) for a Poisson equation subsequently to obtain $p_{m+1}$. In comparison with the inverse source problems treated in [20, 15, 16, 14], we see that both solutions to forward and backward problems appear in (4.13) due to the nonlinearity of Problem 1.1. More importantly, here we should update $p_{m}$ indirectly by solving an extra Poisson equation since we penalize $\nabla p$ instead of $p$ itself. Such an additional procedure, however, does not affect the efficiency because the computational cost for solving (4.13) is rather minor compared with that for solving two time evolution equations. On the other hand, in view of the variational principle, it is readily seen that the solution $p_{m+1}$ of (4.13) coincides with the minimizer of the minimization problem

[TABLE]

for all $p\in H^{1}(\Omega)$ satisfying $p=h_{0}$ on $\partial\Omega$.

Concerning the convergence issue, we notice the relation between the iteration (4.13) and the minimization problem of a surrogate functional

[TABLE]

Indeed, let us fix $q$ and consider the minimization of $J^{s}(p,q)$ with respect to $p$ which is sufficiently close to $q$, e.g., $\|p-q\|_{W^{1,\infty}(\Omega)}\ll 1$. Separating the terms involving $p$ from others in (4.15), we treat $J^{s}(p,q)$ as

[TABLE]

where $v:=u(p)-u(q)$ satisfies

[TABLE]

Utilizing the backward problem (4.11), we take $z(q)$ and $v$ as mutual test functions to deduce

[TABLE]

where we used $p-q=0$ on $\partial\Omega$ and applied Lemma 4.1 with $k=1$ to $v$ to estimate

[TABLE]

within the admissible set $\mathcal{U}_{1}$. Substituting (4.17) into (4.16), we collect the constant component of $J^{s}(p,q)$ as $C_{3}$ and conclude

[TABLE]

Comparing the above expression with (4.14), we figure out that the iterative update (4.13) is almost equivalent to solving a series of minimization problem

[TABLE]

provided that $\|p-p_{m}\|_{W^{1,\infty}(\Omega)}$ is sufficiently small. On the other hand, it is well known that the convergence of (4.18) is guaranteed by the positivity of the surrogate functional $J^{s}(p,q)$ for all $p,q\in\mathcal{U}_{1}$ (see [23]). By definition, this is achieved by taking sufficiently large $K>0$ such that

[TABLE]

Consequently, it reveals that the convergence of (4.18) almost indicates the convergence of our proposed iteration (4.13). Unfortunately, due to the nonlinearity of the problem, we cannot remove the second order term $O(\|p-q\|_{W^{1,\infty}(\Omega)}^{2})$ in (4.17), which prevents us from proving the convergence rigorously.

We close this section by summarizing the main algorithm for the numerical reconstruction.

Compute $p_{m+1}$ according to the iterative update (4.13). 2. 2.

If $\|p_{m+1}-p_{m}\|_{L^{2}(\Omega)}/\|p_{m}\|_{L^{2}(\Omega)}\leq\epsilon$, then stop the iteration. Otherwise, update $m\leftarrow m+1$ and return to Step 1.

@hangfrom@seccntformat

sectionNumerical Examples

@xsect

2.3ex plus.2ex

In this section, we apply the iteration method proposed in the previous section to the numerical treatment for Problem 1.1, and evaluate its numerical performance. More precisely, we shall implement Algorithm 4.4 to reconstruct the principal coefficient $p$ in the hyperbolic equation (4.1). As the first attempt, we restrict ourselves to one spatial dimension, and simply set $\Omega=(0,1)$ and $T=1$. We divide $\overline{\Omega}\times[0,T]=[0,1]^{2}$ into $100\times 100$ equidistant meshes, and employ some unconditionally finite difference methods to solve the $3$ equations involved in Algorithm 4.4, namely, (4.1), (4.11) and (4.13).

We specify various coefficients and parameters to be used in the numerical tests as follows. For the source term $F$ and the initial value $u_{0}$ of (4.1), we fix

[TABLE]

For the boundary condition of $p$ required in (4.3), we simply set

[TABLE]

Given the true solution $p_{\mathrm{true}}$ and thus the noiseless data $u(p_{\mathrm{true}})$, we generate the noisy data $u^{\delta}$ by adding uniform random noises in such a way that

[TABLE]

where $\mathrm{rand}(-1,1)$ denotes the random number uniformly distributed in $[-1,1]$. For the noise level $\delta>0$, we choose it as a certain portion of the amplitude of the noiseless data, that is,

[TABLE]

For the tuning parameter $K$, it should be chosen sufficiently large to guarantee the convergence (see (4.19)). Roughly speaking, it depends on the operator norm of the forward operator which maps $p$ to $u(p)|_{\omega\times(0,T)}$, which is impossible to compute in practice. Hence, we have to postulate that $K$ is proportional to the size $|\omega|$ of the observation subdomain. Analogously, we also make empirical choices of the the regularization parameter $\alpha$ and the stopping criteria $\epsilon$ in Algorithm 4.4 in such a way that

[TABLE]

In all examples, we fix the initial guess as $p_{0}\equiv 1$. As usual, we evaluate the numerical performance of Algorithm 4.4 by the number $N$ of iterations, the relative $L^{2}$ error

[TABLE]

the elapsed time and sometimes the illustrative figures, and we recognize $p_{N}$ as the result of the numerical reconstruction.

1. (a)

A smooth and symmetric true solution $p_{\mathrm{true}}(x)=\frac{1}{2}\sin\pi x+1$. 2. (b)

An asymmetric true solution $p_{\mathrm{true}}(x)=x(x-1)(x-\frac{3}{2})+1$. 3. (c)

A non-smooth true solution $p_{\mathrm{true}}(x)=\frac{1}{2}\min(x,1-x)+\frac{1}{4}\sin\pi x+1$.

Various aspects of the numerical performance are listed in Table 1. The comparisons of true solutions with their reconstructed ones are illustrated in Figure 1.

[TABLE]

and evaluate the performance of Algorithm 4.4 with different combinations of noise levels and observation subdomains. In detail, we first fix the relative noise level as $\delta_{0}=1\%$ as that in Example 5.1, and change the observation subdomain $\omega$ as

[TABLE]

with decreasing sizes. Next, we fix $\omega=\Omega\setminus[0.1,0.9]$ and increase the relative noises $\delta_{0}$ as $0\%$, $1\%$, $2\%$, $4\%$ and $8\%$. In accordance with the above combinations of $\delta_{0}$ and $\omega$, we also change the parameters $M$ and $\alpha$ according to (5.1). The choices of parameters in the tests and the resulting numerical performance are listed in Table 2.

The above examples demonstrate the accuracy and robustness of Algorithm 4.4 as its previous applications e.g.​ in [20, 15]. Especially, in the reconstruction of the principal coefficient, it is obvious that the problem suffers from stronger ill-posedness and nonlinearity compared with the corresponding inverse source problem. Even though, the proposed method still provides satisfactory results with rather small observation subdomain and moderately large noise in observation data.

In Example 5.1, our method proves its feasibility for various choices of true solutions, even including a non-smooth one. This suggests the possibility of relaxing the assumption $p_{\mathrm{true}}\in W^{2,\infty}(\Omega)$ in the derivation of Algorithm 4.4. Unfortunately, it is shown in Figure 1(c) that our method fails to capture the local non-smoothness far away from $\omega$ because of its $L^{2}$-based formulation.

On the other hand, Example 5.2 illustrates the influence of the size of $\omega$ and the noise level upon the numerical performance. The results agree well with our common sense, namely, a smaller observation subdomain $\omega$ results in a larger relative error with more iteration steps until convergence. Still, we need the coverage $\partial\omega\supset\partial\Omega$ for the numerical stability. Meanwhile, both error and iteration steps also increase with larger noise level as expected. However, we see in Table 2 that an $8\%$ relative noise causes considerably large error in the reconstruction, possibly due to the strong nonlinearity of the problem.

@hangfrom@seccntformat

sectionConcluding Remarks

@xsect

2.3ex plus.2ex

The propose of this article is to investigate Problem 1.1, namely, the determination of the spatial component $p(x)$ in the second order coefficient of a hyperbolic equation, from both theoretical and numerical aspects. Theoretically, we are mainly motivated by the existing literature represented by [10] and formulate the problem within the general hyperbolic operator $\mathcal{H}_{p}$ with a time-dependent principal part. On the same direction of [15], we take advantage of the key Carleman estimates for $\mathcal{H}_{p}$ to establish a local Hölder stability result for Problem 1.1. The proof starts from the routine linearization, but unlike [10] we should turn to another Carleman estimate to dominate the $H^{1}$-norm of $p-q$ by that of $\mathrm{div}((p-q)a\nabla u(q))$. The reason traces back to our choice of including $p$ in the divergence in (1.1) for a concrete physical meaning. Instead, if we base the discussion on a nearly non-divergence form

[TABLE]

then the source term after linearization becomes $(p-q)\mathrm{div}(a\nabla u(q))$, and the $L^{2}$ estimate of $p-q$ reduces to an immediate corollary of [15, Theorem 2.3]. The same comment applies to the determination of any spatial components of lower order coefficients in $\mathcal{H}_{p}$, by which we can expect the identical stability result. In these cases, it suffices to replace (2.12) by some analogous non-vanishing assumptions, and we omit the details here.

In the numerical aspect, we adopt the orthodox Tikhonov regularization to interpret Problem 1.1 as a minimization problem. For the highest order coefficient $p$, we penalize the $L^{2}$-norm of $\nabla p$ with its information given on the whole boundary. Calculating the Fréchet derivative, we derive the variational equation for a minimizer of the Tikhonov functional, which involves a backward problem and the Laplacian of $p$. This suggests a novel iterative update (4.13), where one should solve a Poisson equation at each step. Moreover, by the variational principle we find a link between (4.13) and the minimization of a corresponding surrogate functional. Unfortunately, the convergence of the latter does not imply that of the former, because their equivalence is not rigorous due to the nonlinearity of Problem 1.1.

We conclude this paper with some possible future topics related to Problem 1.1. As was mentioned in Remark 2.4, the local stability in Theorem 2.3 relies heavily on the choice of the weight function $\varphi$ in Carleman estimates. We shall consider the possibility of a clever choice of $\varphi$ which optimizes the stability and reduces the observation cost. Meanwhile, another interesting issue is the simultaneous determination of several coefficients, e.g., finding $p,q$ in

[TABLE]

For such kind of problems, possibly one should take several measurements. As a similar but far more difficult case, one can study the same problem for linear anisotropic Lamé systems with time-dependent principal parts. Numerically, the idea of solving an auxiliary equation seems fresh in iteration methods to the best of our knowledge. We are interested in applying it to other inverse problems and analyze its properties, especially convergence.

Acknowledgement  The authors appreciate the valuable discussions with Jin Cheng and Shuai Lu (Fudan University). This work is supported by A3 Foresight Program “Modeling and Computation of Applied Inverse Problems”, Japan Society for the Promotion of Science (JSPS) and National Natural Science Foundation of China. The second and the third authors are partially supported by Grant-in-Aid for Scientific Research (S) 15H05740, JSPS. The second author is supported by JSPS Postdoctoral Fellowship for Overseas Researchers and Grant-in-Aid for JSPS Fellows 16F16319.