Reconstruction of a Time-dependent Potential from Wave Measurements

TL;DR

This paper addresses the challenging problem of reconstructing a time-dependent potential in a wave equation from wave measurements, providing theoretical guarantees and a practical numerical algorithm.

Contribution

It offers a new existence and uniqueness result for the inverse problem, computes the Fréchet derivative, and develops a regularized Newton-like method for reconstruction.

Findings

Successful numerical reconstructions demonstrating feasibility

High reconstruction quality in simulated examples

Efficient algorithm with favorable computational performance

Abstract

We add a time-dependent potential to the inhomogeneous wave equation and consider the task of reconstructing this potential from measurements of the wave field. This dynamic inverse problem becomes more involved compared to static parameters, as, e.g. the dimensions of the parameter space do considerably increase. We give a specifically tailored existence and uniqueness result for the wave equation and compute the Fr\'echet derivative of the solution operator, for which also show the tangential cone condition. These results motivate the numerical reconstruction of the potential via successive linearization and regularized Newton-like methods. We present several numerical examples showing feasibility, reconstruction quality, and time efficiency of the resulting algorithm.

Figures

Tables

Table 1. Table 1: Errors of the reconstruction from u ε superscript 𝑢 𝜀 u^{\varepsilon} for ε = 10 − 2 𝜀 superscript 10 2 \varepsilon=10^{-2} .

	$c_{\text{hat}}$		$c_{\text{plateau}}$
	$L^2([0,T],L^2)$	$H^2([0,T],L^2)$	$L^2([0,T],L^2)$	$H^2([0,T],L^2)$
$n=1$	$26.28\%$	$57.93\%$	$29.47\%$	$33.94\%$
$n=2$	$44.68\%$	$71.90\%$	$37.06\%$	$41.15\%$
$n=3$	$58.85\%$	$81.66\%$	$40.71\%$	$43.50\%$

Table 2. Table 2: Numerically observed orders of convergence when reconstructing from u ε superscript 𝑢 𝜀 u^{\varepsilon} .

	$c_{\text{hat}}$		$c_{\text{plateau}}$
	$L^2([0,T],L^2)$	$H^2([0,T],L^2)$	$L^2([0,T],L^2)$	$H^2([0,T],L^2)$
$n=1$	$0.25$	$0.00$	$0.29$	$0.30$
$n=2$	$0.31$	$0.07$	$0.25$	$0.22$
$n=3$	$0.32$	$0.08$	$0.33$	$0.29$

Table 3. Table 3: Errors of the reconstruction from g ε superscript 𝑔 𝜀 g^{\varepsilon} for ε = 10 − 2 𝜀 superscript 10 2 \varepsilon=10^{-2} .

	$c_{\text{hat}}$		$c_{\text{plateau}}$
	$L^2([0,T],L^2)$	$H^2([0,T],L^2)$	$L^2([0,T],L^2)$	$H^2([0,T],L^2)$
$n=1$	$26.55\%$	$58.74\%$	$27.28\%$	$43.67\%$
$n=2$	$46.11\%$	$73.02\%$	$34.48\%$	$39.24\%$
$n=3$	$66.17\%$	$84.72\%$	$49.30\%$	$51.49\%$

Equations

$$u^{\prime \prime }-\Delta u+cu=f\text{in }[0,T]\times \Omega .$$

$$\{\begin{array}{ll}{u}^{\prime \prime }(t,x)-\Delta u(t,x)+c(t,x)u(t,x)=f(t,x)& (t,x)\in (0,T)\times \Omega ,\\ u(0,x)=u^{\prime }(0,x)=0& x\in \Omega ,\\ u(t,x)=0& (t,x)\in (0,T)\times \partial \Omega .\end{array}$$

$$⟨u^{\prime \prime }(t),\phi ⟩+{\int }_{\Omega }\left[\nabla u(t)\cdot \nabla \phi +c(t)u(t)\phi \right]\mathrm{d}x=\left(f(t),\phi \right)$$

$$a:[0,T]\times H_0^1\times H_0^1\to \mathbb{R},a(t,u,v):={\int }_{\Omega }\left[\nabla u\cdot \nabla v+c(t)uv\right]\mathrm{d}x,$$

$$⟨u^{\prime \prime }(t),\phi ⟩+a(t,u(t),\phi )=\left(f(t),\phi \right)\text{for all }\phi \in H_0^1\text{ and a.e. t∈[0,T]}.$$

$$⟨u_m^{\prime \prime }(t),\phi ⟩+a(t,u_m(t),\phi )=\left(f(t),\phi \right)\text{for all }\phi \in \mathrm{lin}\{{\phi }_j|j=1,\dots ,m\}$$

$${\Vert c-c_b\Vert }_{L^{\infty }([0,T],L^2)}<\delta (\Omega ):={\left(2(1+C_{\mathrm{P}}^2)C_{H_0^1\hookrightarrow L^4}^2\right)}^{-1},$$

$$\underset{t\in [0,T]}{\mathrm{esssup}}{\Vert u_m(t)\Vert }_{H_0^1}+$$

$$\le \left(1+{\Vert c\Vert }_{L^2([0,T],L^2)}\right)e^{C_1\left(1+{\Vert c^{\prime }\Vert }_{L^{\infty }([0,T],L^2)}\right)}{\Vert f\Vert }_{L^2([0,T],L^2)}$$

$$\underset{t\in [0,T]}{\mathrm{esssup}}{\Vert u_m(t)\Vert }_{H_0^1}+\underset{t\in [0,T]}{\mathrm{esssup}}{\Vert u_m^{\prime }(t)\Vert }_{L^2}+{\Vert u_m^{\prime \prime }\Vert }_{L^2([0,T],H^{-1})}\le e^{C_2\left(1+{\Vert c\Vert }_{H^2([0,T],L^2)}\right)}{\Vert f\Vert }_{L^2([0,T],L^2)},$$

$$⟨u_m^{\prime \prime }(t),u_m^{\prime }(t)⟩+a(t,u_m(t),u_m^{\prime }(t))=\left(f(t),u_m^{\prime }(t)\right)\text{for a.e. t∈[0,T]}.$$

$$⟨u_m^{\prime \prime }(t),u_m^{\prime }(t)⟩=\sum _{i=1}^m{\alpha }_i^{\prime \prime }(t){\alpha }_i^{\prime }(t)=\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\sum _{i=1}^m{\alpha }_i^{\prime }(t{)}^2=\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}{\Vert u_m^{\prime }(t)\Vert }_{L^2}^2,$$

$$a^{\prime }(t,v,w):=\frac{\mathrm{d}}{\mathrm{d}t}a(t,v,w)={\int }_{\Omega }{c}^{\prime }(t)vw\mathrm{d}x\text{for }v,w\in H_0^1,$$

$$\frac{\mathrm{d}}{\mathrm{d}t}a(t,u_m(t),u_m(t))$$

$$=2a(t,u_m(t),u_m^{\prime }(t))+a^{\prime }(t,u_m(t),u_m(t)).$$

$$2\left(f(t),u_m^{\prime }(t)\right)+a^{\prime }(t,u_m(t),u_m(t))=\frac{\mathrm{d}}{\mathrm{d}t}a(t,u_m(t),u_m(t))+\frac{\mathrm{d}}{\mathrm{d}t}{\Vert u_m^{\prime }(t)\Vert }_{L^2}^2.$$

$$\int_{0}^{t}\big{[}2\,{\left({f(s)},{u_{m}^{\prime}(s)}\right)}+a^{\prime}(s,u_{m}(s),u_{m}(s))\big{]}{\,\mathrm{d}s}=a(t,u_{m}(t),u_{m}(t))+{{{{\left\|u_{m}^{\prime}(t)\right\|}}}}^{2}_{L^{2}}.$$

$$\displaystyle a(t,u_{m}(t),u_{m}(t))\ +\$$

$$\ge (1+C_{\mathrm{P}}^2{)}^{-1}{\Vert u_m(t)\Vert }_{H_0^1}^2+{\int }_{\Omega }c(t)|u_m(t){|}^2\mathrm{d}x+{\Vert u_m^{\prime }(t)\Vert }_{L^2}^2.$$

$${\int }_{\Omega }c(t)u_m(t{)}^2\mathrm{d}x$$

$$\ge c_0{\Vert u_m(t)\Vert }_{L^2}^2-{\left(2(1+C_{\mathrm{P}}^2)C_{H_0^1\hookrightarrow L^4}^2\right)}^{-1}{\Vert u_m(t)\Vert }_{L^4}^2$$

$$\ge c_0{\Vert u_m(t)\Vert }_{L^2}^2-\frac{1}{2(1+C_{\mathrm{P}}^2)}{\Vert u_m(t)\Vert }_{H_0^1}^2.$$

$$a(t,u_m(t),u_m(t))+{\Vert u_m^{\prime }(t)\Vert }_{L^2}^2\ge \frac{1}{2(1+C_{\mathrm{P}}^2)}{\Vert u_m(t)\Vert }_{H_0^1}^2+c_0{\Vert u_m(t)\Vert }_{L^2}^2+{\Vert u_m^{\prime }(t)\Vert }_{L^2}^2.$$

$$\frac{1}{2(1+C_{\mathrm{P}}^2)}{\Vert u_m(t)\Vert }_{H_0^1}^2+{\Vert u_m^{\prime }(t)\Vert }_{L^2}^2+c_0{\Vert u_m(t)\Vert }_{L^2}^2$$

$$\le {\Vert f\Vert }_{L^2([0,T],L^2)}^2+{\int }_0^t\left[{\Vert u_m^{\prime }(s)\Vert }_{L^2}^2+{\Vert c^{\prime }\Vert }_{L^{\infty }([0,T],L^2)}{C}_{H_0^1\hookrightarrow L^4}^2{\Vert u_m(s)\Vert }_{H_0^1}^2\right]\mathrm{d}s.$$

$$c_0{\Vert u_m(t)\Vert }_{L^2}^2$$

$$\ge \min (0,c_0){\int }_{\Omega }t{\int }_0^t|u_m^{\prime }(s){|}^2\mathrm{d}s\mathrm{d}x\ge \min (0,T{c}_0){\int }_0^t{\Vert u_m^{\prime }(s)\Vert }_{L^2}^2\mathrm{d}s,$$

$$\frac{1}{2(1+C_{\mathrm{P}}^2)}\left({\Vert u_m(t)\Vert }_{H_0^1}^2+{\Vert u_m^{\prime }(t)\Vert }_{L^2}^2\right)\le$$

$$\times \left({\Vert f\Vert }_{L^2([0,T],L^2)}^2+{\int }_0^t{\Vert u_m(s)\Vert }_{H_0^1}^2+{\Vert u_m^{\prime }(s)\Vert }_{L^2}^2\mathrm{d}s\right).$$

$${\Vert u_m(t)\Vert }_{H_0^1}^2+{\Vert u_m^{\prime }(t)\Vert }_{L^2}^2\le {\Vert f\Vert }_{L^2([0,T],L^2)}^2\exp \left(C(T,c_0,\Omega )\left(1+{\Vert c^{\prime }\Vert }_{L^{\infty }([0,T],L^2)}\right)\right)$$

$$⟨u_m^{\prime \prime }(t),v⟩$$

$$\le {\Vert f(t)\Vert }_{L^2}{\Vert v\Vert }_{L^2}+{\Vert \nabla u_m(t)\Vert }_{L^2}{\Vert \nabla v\Vert }_{L^2}+{\Vert c(t)\Vert }_{L^2}{\Vert u_m(t)\Vert }_{L^4}{\Vert v\Vert }_{L^4}$$

$$\le {\Vert f(t)\Vert }_{L^2}+\left(1+C_{H_0^1\hookrightarrow L^4}^2{\Vert c(t)\Vert }_{L^2}\right){\Vert u_m(t)\Vert }_{H_0^1}.$$

$${\Vert u_m^{\prime \prime }\Vert }_{L^2([0,T],H^{-1})}^2={\int }_0^T\underset{\Vert v{\Vert }_{H_0^1}=1}{\sup }|⟨u_m^{\prime \prime }(t),v⟩{|}^2\mathrm{d}t$$

$$\le {\int }_0^T{\left({\Vert f(t)\Vert }_{L^2}+\left(1+C_{H_0^1\hookrightarrow L^4}^2{\Vert c(t)\Vert }_{L^2}\right){\Vert u_m(t)\Vert }_{H_0^1}\right)}^2\mathrm{d}t$$

$$\le {\left(1+{\Vert c\Vert }_{L^2([0,T],L^2)}\right)}^2\exp \left(C(T,c_0,\Omega )\left(1+{\Vert c^{\prime }\Vert }_{L^{\infty }([0,T],L^2)}\right)\right){\Vert f\Vert }_{L^2([0,T],L^2)}^2.$$

Full text

1 Introduction

We consider the inhomogeneous wave equation in a bounded time-space domain $[0,T]\times\Omega$ with a time- and space-dependent potential $c$ and a source $f$,

[TABLE]

For this setting, we tackle the dynamic inverse problem to reconstruct $c$ from measurements of $u$ at specific measurement points for many time steps. This inverse problem provides a simplified model for the non-destructive testing via time-dependent waves in dynamic environments of, e.g. complex carbon-fiber-reinforced polymers under loadings; it is additionally crucial for the detection of non-linear terms in the wave equations from merely an approximate linear model.

Our aim is to show that this dynamic inverse problem for the time-varying quantity $c$ can be mathematically rigorously formulated, analyzed, and stably solved by successive linearization in reasonable computation time. Thus, we first construct suitable function spaces for the coefficients and the solutions to solve the latter partial differential equation with homogeneous initial and boundary conditions in $n=1$, $2$, or $3$ dimensions. Then we show that the parameter-to-solution map is Fréchet differentiable on a suitable domain of definition.

The inverse problem to determine $c$ from full measurements of $u$, as well as its linearization, both turn out to be ill-posed. As we can however show that for our setting the tangential cone condition of Scherzer [Sch95] is satisfied, we consider successive linearization as the starting point for an inversion algorithm.

To be able to cope with the more important case of reduced point measurements of the wave field, too, we compute the necessary operator adjoints and finally detail several numerical experiments computed by the so-called REGINN algorithm of Rieder [Rie05] when applied to the inverse problem. Roughly speaking, these experiments show that the observable region in space is determined by the excitations and the sensor positions up to errors due to the noise level.

There are not so many papers in the literature tackling inverse problems for space- and time-dependent parameters of a wave equation. On the theoretical side, there are a couple of papers proving uniqueness result for various types of coefficients and data, see, e.g. [Ste89, RS91], together with several more recent works that particularly indicate the rising interest in the topic, see [Kia16, Ben15, Sal13, Esk07]. The main tool of many of these papers are geometric optics solutions. Concerning numerical algorithms, there does not seem to exist a similar variety of results, apart from the detection of time-dependent (point) sources for the wave equation, see, e.g. [EH01]. We would like to further note reconstruction results for non-linear elastic materials in [BSS15] indicating future potential fields of application for the algorithms from this paper.

Our solution theory for (1.1) follows the weak solution theory of Lions und Magenes [LM72] as the latter can also be used for more complicated problems than considered in this paper. This weak solution theory shows existence and uniqueness of solution to (1.1) for all $c\in H^{2}([0,T],L^{2}(\Omega))$ that are bounded from below by some (possibly negative) $c_{0}\in\mathbb{R}$ and all square-integrable $f$.

A similar framework provided in Evans’ book [Eva10, Chapter 7.2] requires $c\in C^{1}([0,T]\times\overline{\Omega})$. If one aims to embed the latter space in some Hilbert- or at least in some reflexive Banach space, one typically ends up in a high-order Sobolev space. Thus, resulting reconstructed parameters will then typically possess some extra spatial smoothness, which is somewhat inconvenient from the point of view of applications.

The above-mentioned set of suitable parameters for our solution theory is not yet ready to show, e.g. Fréchet differentiability of the solution operator $c\mapsto u$, but by well-known analytic tools we prove that a suitable open domain for such a derivative exists.

Of course, more general settings than (1.1) should include more general variable coefficients, which is however out of this paper’s scope. Similarly, regularization methods that can be rigorously formulated in Banach spaces to the inversion problem under consideration is a natural continuation of this work that will be considered in a future work.

The remainder of this paper is structured as follows: Sections 2 and 3 treat weak solutions and Fréchet derivatives for the wave equation. Ill-posedness of the resulting inverse problems is shown in Section 4. Sections 5 and 6 consider the discretization of all operators and their adjoints that are involved in the setting. Finally, Section 7 details the inversion algorithm and Section 8 presents numerical examples.

Notation: If there is no danger of confusion, we write $L^{p}$ instead of $L^{p}(\Omega)$, and analogously $H^{k}$ instead of $H^{k}(\Omega)$ for Sobolev spaces; corresponding spaces of functions with zero traces get an additional index [math]. Further, ${\langle{\cdot},{\cdot}\rangle}$ is the duality product between $H^{-1}\coloneqq(H^{1}_{0}){\vphantom{1pt}}^{\prime}$ and $H^{1}_{0}$ and ${\left({\cdot},{\cdot}\right)}$ is the scalar product of $L^{2}$. We identify $L^{2}$ with its dual space such that we typically work in the Gelfand triple $H^{1}_{0}\subset L^{2}\subset H^{-1}$. We do not distinguish between scalar and vector-valued functions and, e.g. also write ${\left({\nabla u},{\nabla v}\right)}$ for the ${L^{2}(\Omega)}^{n}$-scalar product of $\nabla u$ and $\nabla v$.

2 Existence of solution to the wave equation for time-dependent parameters

In this section, we show existence theory for the wave equation with time-dependent parameter $c=c(t,x)$ with $t\in[0,T]$ for some $T>0$ and $x$ in a bounded Lipschitz domain $\Omega\subset\mathbb{R}^{n}$ for $n\in\{1,2,3\}$. The system excitation is modeled by a source $f\in L^{2}([0,T]\times\Omega)$, such that the initial boundary value problem for the wave field $u$ reads

[TABLE]

We initially require $c\in L^{2}([0,T],L^{2})$, but we finally will require more regularity of this parameter. Multiplication of the wave equation in (2.1) by a space-dependent test function $\varphi\in H^{1}_{0}$, integration over $\Omega$, and partial integration yields the following definition of a weak solution to (2.1).

Definition 2.1**.**

A function $u\in H^{2}([0,T],H^{-1})\cap H^{1}([0,T],L^{2})\cap L^{2}([0,T],H^{1}_{0})$ is a weak solution to (2.1) if

[TABLE]

holds for almost every (a.e.) $t\in[0,T]$ and all $\varphi\in H^{1}_{0}$, and if $u$ satisfies the initial conditions $u(0)=0$ in $L^{2}$ and $u^{\prime}(0)=0$ in $H^{-1}$.

Note that the integral in (2.2) involving $c$ is well defined by the smoothness of both $u(t)$ and $\varphi(t)\in H^{1}_{0}$: For $n\leq 3$ the embedding $H^{1}_{0}\hookrightarrow L^{4}$ is continuous, such that the integrand is at least in $L^{1}$. Further, the initial conditions for $u(0)$ and $u^{\prime}(0)$ are well defined because $u\in H^{1}([0,T],L^{2})\hookrightarrow C([0,T],L^{2})$ and $u^{\prime}\in H^{1}([0,T],H^{-1})\hookrightarrow C([0,T],H^{-1})$, such that $u(0)$ naturally belongs to $L^{2}$ and $u^{\prime}(0)\in H^{-1}$.

If we set

[TABLE]

then the weak formulation from the last Definition 2.1 is equivalent to

[TABLE]

To construct such a weak solution, we proceed by Galerkin approximation in finite-dimensional subspaces of $H^{1}_{0}$. To this end, choose some orthogonal basis $(\varphi_{j})_{j\in\mathbb{N}}$ of $H^{1}_{0}$ that is at the same time an orthonormal basis of $L^{2}$ (e.g. via the eigenfunctions of the Laplacian). Working with the Gelfand triple $H^{1}_{0}\subset L^{2}\subset H^{-1}$, the equalities ${\langle{\varphi_{i}},{\varphi_{j}}\rangle}=\delta_{ij}$ imply that the $\varphi_{j}$ are also dual and normalized to each other for the duality product between $H^{-1}$ and $H^{1}_{0}$.

Plugging the finite-dimensional ansatz $u_{m}(t)\coloneqq\sum_{j=1}^{m}\alpha_{j}(t)\varphi_{j}\in\operatorname{lin}\{\varphi_{j}\,|\,j=1,\dots,m\}$ for some $m\in\mathbb{N}$ into (2.3), we note that $u_{m}$ needs to solve

[TABLE]

with zero initial conditions. As is well-known for parameters that are constant in time, this yields $m$ ordinary differential equations for the coefficients $\alpha_{j}$ for $j=1,\dots,m$ with right-hand sides $t\mapsto{\left({f(t)},{\varphi_{i}}\right)}$ in $L^{2}([0,T])$ that are all uniquely solvable in $H^{2}([0,T])$.

Lemma 2.2**.**

For $c\in L^{2}([0,T],L^{2})$ and $m\in\mathbb{N}$ there is $u_{m}\in H^{2}([0,T],H^{1}_{0}(\Omega))$ with zero initial conditions that solves (2.4) for a.e. $t\in[0,T]$.

We next compute an explicit constant that bounds the norms of all $u_{m}$. As $u_{m}(t)$ belongs to the dense subset $\operatorname{lin}\{\varphi_{j}\,|\,j\in\mathbb{N}\}$ of $H^{1}_{0}(\Omega)$, this will finally show that $u_{m}$ converges to a weak solution of the wave equation.

Lemma 2.3** (Energy estimates).**

For $c\in H^{1}([0,T],L^{2})$ with $c^{\prime}\in L^{\infty}([0,T],L^{2})$ we assume that there is $c_{b}\in L^{\infty}([0,T],L^{2})$ with $c_{b}\geq c_{0}\in\mathbb{R}$ such that

[TABLE]

where $C_{H^{1}_{0}\hookrightarrow L^{4}}$ is the operator norm of the embedding $H^{1}_{0}(\Omega)\hookrightarrow L^{4}(\Omega)$ and $C_{\mathrm{P}}$ denotes the Poincaré constant of $H^{1}_{0}(\Omega)$. Then $u_{m}$ from (2.4) satisfies for all $m\in\mathbb{N}$ that

[TABLE]

with $C_{1}>0$ depending only on $T$, $\Omega$ and $c_{0}$. If, additionally, $c\in H^{2}([0,T],L^{2})$, then

[TABLE]

where $C_{1,2}>0$ merely differ by a constant depending on $\Omega$ and $T$.

The energy estimate depends on the lower bound $c_{0}\in\mathbb{R}$ of the comparison parameter $c_{b}$. Even if $c_{0}$ is completely arbitrary, we fix this constant from now on to avoid technicalities.

Proof.

(1) As $u_{m}^{\prime}(t)=\sum_{j=1}^{m}\alpha_{j}^{\prime}(t)\varphi_{j}$ belongs to $\operatorname{lin}\{\varphi_{j}\,|\,j\in\mathbb{N}\}\subset H^{1}_{0}$, we plug $u_{m}^{\prime}$ as test function into the weak formulation of $u_{m}$,

[TABLE]

Due to the representation of $u_{m}^{\prime}(t)$, we directly get that

[TABLE]

and, analogously, ${\left({\nabla u_{m}^{\prime}(t)},{\nabla u_{m}(t)}\right)}=\tfrac{{\,\mathrm{d}}}{{\,\mathrm{d}t}}{{{{\left\|\nabla u_{m}\right\|}}}}^{2}_{L^{2}(\Omega)^{n}}/2$. Setting

[TABLE]

one computes, again by finiteness of the sum representation of $u_{m}^{\prime}$, that

[TABLE]

Together with (2.6), the latter equality shows that for a.e. $t\in[0,T]$ there holds

[TABLE]

Thus, the fundamental theorem of analysis and the zero initial conditions for $u_{m}$ imply that

[TABLE]

Bounding the right-hand side from below by Poincaré’s estimate ${{{{\left\|u_{m}\right\|}}}}_{L^{2}}\leq C_{\mathrm{P}}{{{{\left\|\nabla u_{m}\right\|}}}}_{L^{2}}$ we hence arrive at

[TABLE]

(2) As $c$ may take negative values, it is now crucial to bound the second integral on the right via the auxiliary function $c_{b}(t)$ that is by assumption bounded from below by $c_{0}\in\mathbb{R}$,

[TABLE]

As the constant in front of ${{{{\left\|u_{m}(t)\right\|}}}}^{2}_{H^{1}_{0}}$ is positive, we conclude that

[TABLE]

In combination with (2.8), we have hence shown that

[TABLE]

(3) As we aim to apply Gronwall’s inequality (see, e.g. [Eva10, Appendix B.2]), we need to control the term $c_{0}{{{{\left\|u_{m}(t)\right\|}}}}_{L^{2}}^{2}$. This is simple if $c_{0}$ is non-negative, as that term can then be dropped. More generally,

[TABLE]

such that

[TABLE]

Now, Gronwall’s inequality implies that

[TABLE]

where $C(T,c_{0},\Omega)$ is a placeholder for a constant depending only on $T,c_{0}$ and $\Omega$.

(4) To obtain $H^{-1}$-bounds for $u_{m}^{\prime\prime}(t)$, let us finally choose any $v\in H^{1}_{0}(\Omega)$ with ${{{{\left\|v\right\|}}}}_{H^{1}_{0}}=1$ and note that for a.e. $t\in[0,T]$,

[TABLE]

Representing the $H^{-1}$-norm as a dual norm shows that

[TABLE]

Together with (2.9) this shows the lemma’s first claimed bound. For $c\in H^{2}([0,T],L^{2})$ the estimate simplifies as $c^{\prime}\in H^{1}([0,T],L^{2})$, a space which is continuously embedded in $L^{\infty}([0,T],L^{2})$. The arising norms of $c$ und $c^{\prime}$ can hence be bounded by ${{{{\left\|c\right\|}}}}_{H^{2}([0,T],L^{2})}$. ∎

As is well-known, the estimate of the last lemma cannot be shown transferred to the (weak) limit $u$ of the bounded sequence $u_{m}$, because the regularity of $u(t)$ is too low to test its variational formulation by $u^{\prime}(t)$. The energy estimates, however, do allow to prove existence of at least one solution to (2.3). Uniqueness of this solution is then achieved by a standard regularity trick.

Theorem 2.4**.**

Under the assumptions of Lemma 2.3 there is a unique weak solution $u$ to (2.1) that satisfies the energy estimates of Lemma 2.3.

Proof.

Concerning existence we note that the sequence $(u_{m})_{m\in\mathbb{N}}$ is bounded in $L^{2}([0,T],H^{1}_{0})$, $(u_{m}^{\prime})_{m\in\mathbb{N}}$ is bounded in $L^{2}([0,T],L^{2})$ and $(u_{m}^{\prime\prime})_{m\in\mathbb{N}}$ is bounded in $L^{2}([0,T],H^{-1})$. Due to reflexivity of these spaces we obtain a subsequence, that we denote for notational simplicity also by $(u_{m})_{m\in\mathbb{N}}$, such that

[TABLE]

It is easy to show that $w=u^{\prime}$ as well as $z=u^{\prime\prime}$. The argument showing that $u$ is in fact a weak solution to (2.1) and that it is unique can be found in [LM72]. Finally, we take a look at the energy estimates for $u$. Because of the weak convergences in (2.10) we only obtain $L^{2}$-estimates in time at first, but it also follows that $(u_{m})_{m\in\mathbb{N}}$ is bounded in $L^{\infty}([0,T],H^{1}_{0})$ and $(u_{m}^{\prime})_{m\in\mathbb{N}}$ is bounded in $L^{\infty}([0,T],L^{2})$. Up to extraction of a further subsequence (that we do again not denote explicitly), this shows that

[TABLE]

The energy estimates of Lemma 2.3 hence transfer to $u$ which in particular belongs to the space $H^{2}([0,T],H^{-1})\cap W^{1,\infty}([0,T],L^{2})\cap L^{\infty}([0,T],H^{1}_{0})$. ∎

The last result states that for parameters that are first smooth enough and second close enough to a function bounded from below by some constant, there is a unique solution to the wave equation that satisfies an energy estimate. Before stating this as a corollary, recall the number $\delta(\Omega)$ from (2.5) and the arbitrary, but fixed, constant $c_{0}\in\mathbb{R}$ from Lemma 2.3.

Corollary 2.5**.**

The wave equation (2.1) possesses for every $f\in L^{2}([0,T],L^{2})$ and every

[TABLE]

a unique weak solution. This weak solution satisfies the energy estimate

[TABLE]

with $C$ depending only on $T$, $\Omega$ and $c_{0}$.

3 Fréchet differentiability with respect to the parameter

The solution theory from the last section allows to define a solution operator $S:\,c\,\mapsto\,u$ mapping the time-dependent parameter $c$ to the wave $u$. This operator is obviously non-linear since, e.g. $c=0$ is not mapped to the trivial solution. For inversion, we will hence exploit that $S$ can be locally linearized by its Fréchet derivative.

Differentiability of $S$ follows from Lipschitz continuity that we derive, by and large, via the energy estimate from Corollary 2.5. The precise setting is fixed in the following formal definition of $S$. For simplicity, we fix the source term $f\in L^{2}([0,T],L^{2})$ for a moment and recall a last time the fixed constants $\delta(\Omega)$ from (2.5) and $c_{0}\in\mathbb{R}$ from Lemma 2.3.

All constants $C$ we use in the sequel may change value from line to line but depend merely on $\Omega$, $c_{0}$, $T$, and an additional Lebesgue index $p$ that is fixed in the following definition.

Definition 3.1**.**

We consider ${S}\,\colon\,{\mathcal{D}(S)\subset X}\to{Y}$ mapping the parameter $c$ to the solution $u$ of (2.3) as an operator with domain $\mathcal{D}(S)$ embedded in $X\coloneqq H^{2}([0,T],L^{2})\cap L^{2}([0,T],L^{p})$ for a fixed $p=p(n)$ with $p(1)\geq 2$, $p(2)>2$, and $p(3)>3$,

[TABLE]

Due to Corollary 2.5, we further set the solution space $Y$ to

[TABLE]

Both $X$ and $Y$ are equipped with their natural norms.

Note that the last definition requires parameters $c$ to belong to $L^{2}([0,T],L^{p})$ for some index $p$. The reason behind is that the product $h(t)w(t)$ of an $L^{p}$- and an $H^{1}_{0}$-function on a bounded domain belongs to $L^{2}$, which becomes fundamental for Lemma 3.2. This is due to the continuous embedding $H^{1}_{0}\hookrightarrow L^{q}$ for $1\leq q\leq\infty$ and $n=1$, $1\leq q<\infty$ and $n=2$, as well as $1\leq q<6$ and $n=3$. Precisely, for $h\in X$, $w\in Y$ and a.e. $t\in[0,T]$,

[TABLE]

which follows from the Hölder inequality, that is, $2p/(p-2)<\infty$ for $p>2$ and $2p/(p-2)<6$ for $p>3$. Squaring and integrating (3.2) finally results in

[TABLE]

Note further that the norm on $X$ is stronger than the $L^{\infty}([0,T],L^{2})$-norm, such that $\mathcal{D}(S)$ is always open in $X$, which is crucial for Fréchet differentiability.

When treating inverse problems, we of course work with noisy data in $L^{2}([0,T],L^{2})$; nevertheless, the following results profit from the somewhat more involved image space $Y$ from (3.1).

Lemma 3.2**.**

The forward map $S$ is locally Lipschitz continuous.

Proof.

To $c_{1}$ and $c_{2}\in\mathcal{D}(S)$ we assign $u_{1}\coloneqq Sc_{1}$ and $u_{2}\coloneqq Sc_{2}$. The difference $w\coloneqq u_{1}-u_{2}$ then satisfies

[TABLE]

for a.e. $t\in[0,T]$ and all $\varphi\in H^{1}_{0}$, subject to homogeneous initial values. Corollary 2.5 shows that

[TABLE]

where we exploited that $(c_{1}-c_{2})u_{2}\in L^{2}([0,T],L^{2})$ can be estimated in norm as in (3.2) via the generalized Hölder’s inequality, our choice of $p$, and the Poincaré estimate by the $X$-norm of $c_{1}-c_{2}$ times the norm of $u_{2}$ in $L^{\infty}([0,T],H^{1}_{0})$. All arising constants can be uniformly bounded on every bounded set in $X$, which shows the claim. ∎

Formally computing the derivative of the weak formulation (2.3) of $u\in X$ with respect to $c\in\mathcal{D}(S)$ in direction $h\in X$ shows that $u_{h}\coloneqq(S^{\prime}c)[h]$ needs to solve the variational formulation

[TABLE]

with zero initial conditions, i.e. $u_{h}\in X$ is the weak solution to

[TABLE]

with zero initial and boundary conditions. The next theorem makes this formal argument rigorous.

Theorem 3.3**.**

The forward operator $S$ is Fréchet differentiable at $c\in\mathcal{D}(S)$: The derivative equals $S^{\prime}c\in\mathcal{L}(X,Y)$ and satisfies

[TABLE]

For every $h\in X$ the function $u_{h}=(S^{\prime}c)[h]\in Y$ is the unique weak solution of the weak formulation

[TABLE]

for all $\varphi\in H^{1}_{0}$, a.e. $t\in[0,T]$, and subject to homogeneous initial conditions $u_{h}(0)=u_{h}^{\prime}(0)=0$.

Proof.

As $\mathcal{D}(S)$ is open in $X$, there is for every $c\in\mathcal{D}(S)$ an open $X$-ball centered in $c$ such that all $h$ in this ball satisfy that $c+h\in\mathcal{D}(S)$. The solution $u_{h}$ from (3.4) is well-defined because the right-hand side belongs to $L^{2}([0,T],L^{2})$; $u_{h}$ further satisfies the energy estimate from Corollary (2.5) with $f$ replaced by $-hu$. The difference $w\coloneqq S(c+h)-Sc-(S^{\prime}c)[h]$ solves

[TABLE]

for all $\varphi\in H^{1}_{0}$ and a.e. $t\in[0,T]$, with homogeneous initial conditions. Thus, Corollary 2.5 and (3.2) imply that

[TABLE]

Lipschitz continuity of $S$ now implies that

[TABLE]

as $h\to 0$ in $X$. Clearly, $S^{\prime}c$ is a linear operator, such that the energy estimate from Corollary 2.5 shows that $S^{\prime}c$ belongs to $\mathcal{L}(X,Y)$, which finishes the proof. ∎

The last proof shows that the definition of the parameter space $X$ in Definition 3.1 as a subset of $L^{2}([0,T],L^{p})$ is crucial to be able to bound the right-hand side $-hu$ in the variational formulation for the derivative $u_{h}$.

4 Ill-posedness of the inverse problem

As we intend to reconstruct $c$ from wave measurements for several sources $f$ we still need to extend the forward operator $S$. To this end, we first discuss suitable measurement operators $\psi$ to construct vector-valued solution and measurement operators. For simplicity, we restrict ourselves to a particular linear and continuous measurement operator that we use later on for our numerical experiments. Second, we rigorously define the inverse problem and prove its ill-posedness in several settings.

In applications, measurements are typically taken by sensors at fixed spatial positions $x_{i}^{s}$ and many instances of time $t_{i}^{s}$, which yields measured values of the wave at points $(t_{1}^{s},x_{1}^{s}),\dots,(t_{l}^{s},x_{l}^{s})\in[0,T]\times\Omega$. Unfortunately, our solution space $Y$ does not imply that point measurements of $u\in Y$ at $(t_{i}^{s},x_{i}^{s})$ depend continuously on $u\in Y$. However, sensors anyway provide mean values of the wave in time and space over small regions around the introduced measurement points $(t_{i}^{s},x_{i}^{s})$. Thus, we model measurements as a convolution in time and space of the wave field against an integral kernel $k:\,\mathbb{R}\times\mathbb{R}^{n}\to\mathbb{R}$ in $L^{2}(\mathbb{R}^{n+1})$,

[TABLE]

and evaluate the convolved version $u_{k}$ of $u$ at the measurement points $(t_{i}^{s},x_{i}^{s})$. Precisely, for numbers $r_{t}>0$ and $r_{x}>0$ we define a particular kernel $k$ via the normalized auxiliary function $g:\,\mathbb{R}\to\mathbb{R}$ that equals $\sqrt{3}\,(1-s)$ if $s<1$ and zero otherwise as follows,

[TABLE]

This kernel clearly belongs to $L^{2}(\mathbb{R}^{n+1})\cap C(\mathbb{R}^{n+1})$ such that the convolution $u_{k}$ depends in the maximum norm continuously on $u\in L^{2}([0,T],L^{2})$. The associated evaluation operator

[TABLE]

hence roughly speaking models sensor measurements close to $u(t_{i}^{s},x_{i}^{s})$ if $r_{t}$ and $r_{x}$ are small. Due to the Cauchy-Schwarz inequality, ${{{{\left\|\psi\right\|}}}}_{\mathcal{L}(L^{2}([0,T],L^{2}),\mathbb{R}^{l})}\leq\sqrt{l}\,{{{{\left\|k\right\|}}}}_{L^{2}(\mathbb{R}^{n+1})}=\sqrt{l}$. Later on, we will require the adjoint $\psi^{*}$ of $\psi$ and hence note already here that $\psi^{*}:\,\mathbb{R}^{l}\to L^{2}([0,T],L^{2})$ can be characterized by

[TABLE]

Finally, we allow for several sources $f_{1},\dots,f_{d}\in L^{2}([0,T],L^{2})$ for some $d\in\mathbb{N}$ as excitation for the wave equation (2.3). For any parameter $c$ that satisfies the assumptions of Corollary 2.5, the source $f_{i}$ defines a wave $u_{i}\in Y$, such that for every $c\in X$ we can associate a forward operator $S_{i}$ to $f_{i}$ as in Definition 3.1. The vector-valued solution operator $\mathbf{S}$ simply collects all $S_{i}$ in a vector (and we do not distinguish column- and row vectors here).

Definition 4.1**.**

For the spaces $X$ and $\mathcal{D}(\mathbf{S})\coloneqq\mathcal{D}(S)$ as in Definition 3.1 we set the vector-valued solution operator $\mathbf{S}$ as

[TABLE]

Further, the evaluation operator $\Psi:\,L^{2}([0,T],L^{2})^{d}\to(\mathbb{R}^{l})^{d}$ with $\Psi(u_{1},\dots,u_{d})=(\psi\,u_{1},\dots,\psi\,u_{d})$ is continuous, and the non-linear measurement operator mapping parameters to measurements for $d\in\mathbb{N}$ sources $f_{1},\dots,f_{d}$ is

[TABLE]

Figure 1 illustrates connections between the vector-valued operators $\Phi$, $\Psi$, and $\mathbf{S}$. Of course, $\mathbf{S}$ and $\Phi$ are component-wise Fréchet differentiable.

The inverse problem we consider from now on is to determine the parameter $c\in X$ either from the solution $\mathbf{S}c\in L^{2}([0,T],L^{2})^{d}$ or from sensor measurements modeled by $\Phi c\in(\mathbb{R}^{l})^{d}$. We hence aim to determine the parameter of a differential equation from its impact on a wave and of course expect ill-posedness of this task. This holds particularly if the data $\Phi c$ stems from sensors placed on a surface inside $\Omega$, since the unknown parameter then depends on more variables than the data. Since $\Phi$ is by construction an operator with finite-dimensional range (and hence in particular features a closed range), the question of ill-posedness in the sense of the subsequent definition is anyway irrelevant for this measurement operator. (Note the comment in the end of this section on different kinds of measurement models that partly touches this point.) We prove in the rest of this section that the solution operator $S$ or its derivative $S^{\prime}$ yield locally ill-posed inversion problems. This directly implies ill-posedness of operator equations involving $\mathbf{S}$ or its linearization.

For a general operator ${F}\,\colon\,{\mathcal{D}(F)\subset V}\to{W}$ between Banach spaces $V$ and $W$ we recall from [Sch+12, Definition 3.15] that the equation $Fx=y$ is locally ill-posed in $x^{+}\in\mathcal{D}(F)$ if there exist for all $r>0$ sequences $(x_{n})\subset B_{r}(x^{+})\cap\mathcal{D}(F)$ such that $\|Fx_{n}-Fx^{+}\|_{W}\to 0$ but $\|x_{n}-x^{+}\|_{V}\not\to 0$ as $n\to\infty$. (See [HS98, Definition 1.1] for the corresponding definition in Hilbert spaces.) As the notion of locality is meaningless for linear problems, a linear operator equation is either everywhere locally ill-posed or else everywhere locally well-posed. For reflexive spaces $V$ and $W$ and $D(F)=V$, a linear operator equation is ill-posed if and only if the linear operator $F$ possesses a non-closed range or fails to be injective, see [Sch+12, Proposition 3.9].

This last point gets essential when one linearizes $Sc=u$ via the Fréchet derivative $S^{\prime}$ in $c^{+}\in\mathcal{D}(S)$ and tackles $(S^{\prime}c^{+})[h]=g-Sc^{+}$ as an equation for $h\in X$. For the following results, recall that $Y=H^{2}([0,T],H^{-1})\cap W^{1,\infty}([0,T],L^{2})\cap L^{\infty}([0,T],H^{1}_{0})$ is the natural image space for $S$ that results from the energy estimates.

Lemma 4.2**.**

If $f\neq 0$, then ${S^{\prime}c}\,\colon\,{X}\to{L^{2}([0,T],L^{2})}$ is for all $c\in\mathcal{D}(S)$ a compact operator with infinite-dimensional range. In particular, $\mathrm{Rg}(S^{\prime}c)$ is not closed in $L^{2}([0,T],L^{2})$ and the linearized operator equation $(S^{\prime}c^{+})[h]=g-Sc^{+}$ is locally ill-posed in every $h\in X$.

Proof.

We already know that $S^{\prime}c$ is bounded and linear from $X$ into $Y$. The embedding $Y\hookrightarrow L^{2}([0,T],H^{1}_{0})\cap H^{1}([0,T],L^{2})$ is continuous and from Simon [Sim86, p. 85] we know that the compact embedding of $H^{1}_{0}$ in $L^{2}$ implies that

[TABLE]

is compact as well. Thus, ${S^{\prime}c}\,\colon\,{X}\to{L^{2}([0,T],L^{2})}$ is compact and linear. If this operator possesses a finite-dimensional range, then the set of right-hand sides $-h\cdot Sc$ in the variational formulation (3.4) of $(S^{\prime}c)[h]$ must also belong to a finite-dimensional space by unique solvability of the wave propagation problem (2.3). This forces $Sc$ and hence also $f$ to vanish, what we excluded in the lemma, and hence proves by contradiction that $\mathrm{Rg}(S^{\prime}c)$ cannot have finite dimension. As an infinite-dimensional range of a compact linear operator cannot be closed, we have shown the lemma’s claim. ∎

The next lemma prepares a subsequent example on the ill-posedness of the linearized operator equation at $c=0$.

Lemma 4.3**.**

If $c\in\mathcal{D}(S)$ satisfies that $Sc\neq 0$ a.e. in $[0,T]\times\Omega$, then $S^{\prime}c$ is injective.

Proof.

If $(S^{\prime}c)[h]$ vanishes for some non-zero $h\in X$, then the right-hand side $-h\cdot Sc$ of the formulation (3.4) for $(S^{\prime}c)[h]$ vanishes, which contradicts the assumption that $Sc\neq 0$ a.e.. ∎

Example 4.4**.**

Set $n=1$, $\Omega=(0,\pi)$, $f(t,x)=2\sin(x)\sin(t)$, $c=0$, and fix $T>0$ arbitrarily. Then the range $\mathrm{Rg}(S^{\prime}c)$ is not closed in $Y$, i.e. the linearized equation at $c=0$ is everywhere locally ill-posed between $X$ and $Y$ (and a fortiori also between $X$ and every Banach space containing $Y$).

Proof.

Separation of variables shows that the solution to $u^{\prime\prime}(t,x)-\partial_{x}^{2}u(t,x)=2\sin(x)\sin(t)$ is given by $u(t,x)=\sin(x)\left(\sin(t)-t\cos(t)\right)$. Obviously $u\neq 0$ almost everywhere, which means that $S^{\prime}c$ is injective. We define $(h_{k})_{k\in\mathbb{N}}\subset X=H^{2}([0,T],L^{2})\cap L^{2}([0,T],L^{p})$ as $h_{k}(t,x)=k^{1/2}$ if $x<1/k$ and [math] else. This sequence converges point-wise to zero, but not in the $X$-norm as ${{{{\left\|h_{k}\right\|}}}}_{X}={{{{\left\|h_{k}\right\|}}}}_{L^{2}([0,T],L^{2})}=T^{1/2}.$ The energy estimates imply that $u_{h}^{k}\coloneqq(S^{\prime}c)[h_{k}]$ is bounded by

[TABLE]

Due to injectivity of $S^{\prime}c$ we conclude that $S^{\prime}$ cannot have a bounded generalized inverse that is continuous, which proves the claim. ∎

The last example exploits the zeros of the solution $u(t,\cdot)$ at the boundary of $(0,\pi)$; analogous examples can be constructed independent of dimension, the right-hand side $f$, or the chosen boundary conditions, as long as $u$ is at least continuous and possesses zeros in $[0,T]\times\overline{\Omega}$.

In a Hilbert space setting there are many results known that connect the ill-posedness of a non-linear equation with the ill-posedness of its linearization, see, e.g. Section 2 in [HS98]. One example is the so-called tangential cone condition, see (4.4) below or [Sch95], which furthermore straightforwardly extends to Banach spaces.

Theorem 4.5**.**

Assume that ${F}\,\colon\,{\mathcal{D}(F)\subset V}\to{W}$ is Fréchet differentiable between Banach spaces $V$ and $W$. If for some $x^{+}\in\mathcal{D}(F)$ there are $r>0$ and $0\leq\omega<1$ such that

[TABLE]

holds for all $v,w\in B(x^{+},r)\cap\mathcal{D}(F)$, then

[TABLE]

In this case, the non-linear problem $F(x)=y$ is locally ill-posed in $x^{+}$ if and only if the linearized problem at $x^{+}$ is locally ill-posed everywhere; that is, either $F^{\prime}(x^{+})$ is not injective or $\mathrm{Rg}(F^{\prime}(x^{+}))$ is not closed in $W$.

Proof.

It is well-known that condition (4.3) transforms into (4.4) by the reverse triangle inequality. If the latter condition holds, the linearized residual in the enumerator behaves as the non-linear residual in the denominator, such that the non-linear and the linearized operator equation can only be jointly locally ill-posed. ∎

Recall from Definition 3.1 that we have fixed a dimension-dependent Lebesgue index $p\in[2,\infty)$ when setting up the parameter space $X$. The next lemma exploits the conjugate index $q\in(2,\infty]$ defined by $1/p+1/q=1/2$.

Theorem 4.6**.**

Define $q\in(2,\infty]$ as above. Then both settings ${S}\,\colon\,{\mathcal{D}(S)\subset X}\to{L^{\infty}([0,T],H^{1}_{0})}$ and ${S}\,\colon\,{\mathcal{D}(S)\cap L^{4}([0,T],L^{p})}\to{L^{4}([0,T],L^{q})}$ for the solution operator $S$ allow to prove the non-linearity condition (4.3). As both $L^{4}([0,T],L^{p})$ and $L^{4}([0,T],L^{q})$ are reflexive, the conclusion of Theorem 4.5 holds for the second case, i.e. local ill-posedness of the non-linear equation $Sc=u$ at some $c^{+}\in\mathcal{D}(S)\cap L^{4}([0,T],L^{p})$ is equivalent to the local ill-posedness of the corresponding linearized equation $(S^{\prime}c^{+})[h]=u-Sc^{+}$ between $L^{4}([0,T],L^{p})$ and $L^{4}([0,T],L^{q})$.

The following proof clearly shows that there are many more interesting settings for the pre-image and image space of $S$ than announced in the theorem.

Proof.

We have already estimated the linearization error in (3.5): For $c_{1,2}\in\mathcal{D}(S)$,

[TABLE]

We estimate the last term on the right as in (3.2),

[TABLE]

such that

[TABLE]

As the $X$-norm is stronger than the $L^{\infty}([0,T],L^{q})$-norm,

[TABLE]

holds whenever $c_{1}-c_{2}$ is small enough in the $X$-norm. Thus, the claimed bound (4.3) holds if we choose $X$ as pre-image space and $L^{\infty}([0,T],H^{1}_{0})$ as image space for $S$, and $r>0$ small enough. (The $Y$-norm can always be bounded from below by the $L^{\infty}([0,T],H^{1}_{0})$-norm.)

To prove the analogous result for the pre-image space $X\cap L^{4}([0,T],L^{p})$, one uses Hölder’s inequality with twice the index $4$ instead of (4.5),

[TABLE]

For the rest of the proof it is then sufficient to reduce $X$ to $X\cap L^{4}([0,T],L^{p})$ and to set the image space to $L^{4}([0,T],L^{q})$; notably, these choices merely yield reflexive Banach spaces. ∎

Apart from condition (4.4), one can also show Lipschitz continuity of $c\mapsto S^{\prime}c$ in the operator norm between $X$ and $L^{2}([0,T],L^{2})$, that is, ${{{{\left\|S^{\prime}c-S^{\prime}c^{+}\right\|}}}}\leq L\,{{{{\left\|c-c^{+}\right\|}}}}_{X}$ holds for all $c\in B(c^{+},r)\cap\mathcal{D}(S)$. If one embeds the parameter space $X$ into a Hilbert space, this allows to show that local ill-posedness of $Sc=g$ at $c^{+}\in\mathcal{D}(S)$ implies local ill-posedness of the linearized equation at $c^{+}$, see [HS98, Section 2].

Let us finally mention a classical result on the ill-posedness of equations involving $S$ composed with a measurement operator. To this end, consider an operator $M$ such that the product $M\circ S$ is compact, continuous, and weakly sequentially closed from, roughly speaking, all parameters in $H^{2}([0,T],L^{2})$ that possess some lower real bound, into some separable Hilbert space $Z$. (Precisely, the domain of definition is the set defined in (2.11) and included in $H^{2}([0,T],L^{2})$.) If $M\circ S$ additionally possesses an infinite-dimensional range, then the operator equation $(M\circ S)c=g$ is locally ill-posed at any parameter $c$ in this domain of definition, see, e.g. [EKN89, Proposition A3]. This general result is independent of notions of derivatives and hence merely requires parameters $c$ in the set from (2.11).

5 Discretization of the wave equation

In this section we discuss the discretization of the wave equation (2.1) that we use to compute our numerical examples in Section (8). Recall that our existence theory treats the weak formulation of the wave equation,

[TABLE]

together with zero initial conditions $u(0)=u^{\prime}(0)=0$. We discretize the latter problem by Rothe’s method, i.e. we start by discretization in time, and consider the first-order system gained from $v\coloneqq u^{\prime}$ as additional unknown,

[TABLE]

Using a fixed step size $\Delta t>0$ we obtain time steps

[TABLE]

write $u^{i}\coloneqq u(t_{i})$ and analogously $v^{i}$ and $f^{i}$ for $i=0,\dots,N-1$, and set

[TABLE]

We approximate all time derivatives in (5.1) by a $\theta$-scheme, i.e. a weighted average of forward- and backward difference quotients in $u^{i}$ and $v^{i}$. Further elmininating the dependence of the first equation on $v^{i}$ shows that $(u^{i},v^{i})$ solves

[TABLE]

for $i=1,\dots,N$, with initial values $u^{0}=v^{0}=0$. Given $(u^{i-1},v^{i-1})$, the first equation (5.2a) can be used to compute $u^{i}$ by solving one elliptic problem and then plug the result into (5.2b) to compute $v^{i}$ via a second elliptic problem. We actually choose $\theta=1/2$ to obtain the Crank-Nicolson scheme, which converges in each time step of second order as $\Delta t\to 0$. The error of the last step and consequently the total error are hence of the first order in $\Delta t\to 0$. In addition, the scheme is unconditionally stable and does not exhibit energy loss, see Larsson and Thomée [LT03].

To transform the semi-discrete system (5.2) into a fully discrete one we rely on the finite element method. For technical simplicity we assume that $\Omega$ is a polygon and consider shape-regular and quasi-uniform triangulations $\mathcal{S}$ of $\Omega$ that satisfy $\overline{\Omega}=\bigcup_{S\in\mathcal{S}}\overline{S}$ and $\operatorname{diam}S\leq h$ for all $S\in\mathcal{S}$. For all simplexes $S\in\mathcal{S}$ we denote all affine mappings on $S$ as $\mathbb{P}_{1}(S)$ and introduce the finite-dimensional variational approximation spaces

[TABLE]

These finite-dimensional spaces define discrete approximations $u^{i}_{h}\in V_{h}$ and $v^{i}_{h}\in V_{h}$ to $u^{i}$ and $v^{i}$ by restricting the test function $\varphi$ in (5.2) to $V_{h}$, too. Standard error estimates for, e.g. the $L^{2}$-error between $u^{i}_{h}$ and $u^{i}$ indicate this error to be of second order in the diameter of the largest simplex of $\mathcal{S}$, see, e.g. Brenner and Scott [BS02].

If we denote the nodal basis of $V_{h}$ by $\{\varphi_{1},\dots,\varphi_{K}\}$, then both $u^{i}_{h}$ and $v^{i}_{h}$ are represented by $K$ coefficients,

[TABLE]

Linearity of (5.2) shows that the latter system is equivalent to the linear system of size $K\times K$ for \vec{u}^{i}=\big{(}{\vec{u}}^{i}_{k}\big{)}_{k} one gets by inserting (5.3) into (5.2). Let us define the mass matrix $M$ and stiffness matrix $A^{i}$ through

[TABLE]

and abbreviate expressions involving $f$ as a vector $F^{i}=({\left({\theta f^{i}+(1-\theta)f^{i-1}},{\varphi_{j}}\right)})_{j=1,\dots,K}$ in $\mathbb{R}^{K}$. This establishes the fully discrete system

[TABLE]

which is best solved using an iterative method like GMRES. If $v^{i}$ is not needed for further computations then the solution of (5.5b) for $v^{i}$ may be omitted if $Mv^{i}$ is stored instead. Since we regard $S$ as a map into $L^{2}([0,T],L^{2})$ this is the case for us.

For the triangulation of $\Omega$, the bookkeeping of the basis functions, and the assembly of (5.5) we use the finite element toolbox ALBERTA [SS06].

Although $c$ is not the solution of a PDE we nevertheless discretize it as an element of $V_{h}$ at each of the time steps $t_{i}$, in the very same way as $u$. This approach has the advantage of not requiring additional data structures for searched-for parameters.

6 Computation of adjoints of derivatives

The discretization scheme (5.5) allows to numerically approximate $S$ and its derivative $S^{\prime}$. The regularization scheme for $c$ that we present in Section 4 however also requires an approximation of the (complex-valued) transpose operator $(S^{\prime}c)^{*}$. For simplicity, we will rather require the adjoint operator later on since we artificially change into a Hilbert space framework in the next section. The importance of knowing such an operator is however already clear from linear regularization theory via filter functions.

From now on we consider $S^{\prime}$ to be a linear operator from $X$ into $L^{2}([0,T],L^{2})$ and compute its transpose operator mapping $L^{2}([0,T],L^{2})$ into $X^{\prime}$. Note that $S^{\prime}$ from (3.4) can be decomposed as

[TABLE]

where, first, ${M_{c}}\,\colon\,{X}\to{L^{2}([0,T],L^{2})}$, $h\mapsto-uh$, multiplies $h$ by $-u=-Sc$. Second, $L_{c}$ is a (weak) solution operator for the wave equation $w^{\prime\prime}-\Delta w+cw=g$ with variable right-hand side $g$,

[TABLE]

with zero initial and Dirichlet boundary conditions for $w$. By the above decomposition of $S^{\prime}c$ into two bounded linear operators we next compute the transpose $(S^{\prime})^{\ast}$. The resulting numerical schemes will actually carry over to $\Phi^{\prime}c$ and $\mathbf{S}^{\prime}c$, such that we beforehand note the following corollary of (6.1).

Corollary 6.1**.**

Assume that $c\in\mathcal{D}(\mathbf{S})$ and recall from (4.2) that $\Phi=\Psi\circ\mathbf{S}$.

1. (1)

$\mathbf{S}$* and $\Phi$ are Fréchet-differentiable in $c$ and $\Phi^{\prime}c=\Psi\,\circ\,\mathbf{S}^{\prime}c$. For the solution operator $\mathbf{L}_{c}$ mapping $g\in L^{2}([0,T],L^{2})^{d}$ to $(L_{c}\,g_{1},\dots,L_{c}\,g_{d})\in L^{2}([0,T],L^{2})^{d}$ and $\mathbf{M}_{c}\coloneqq\left(M_{c,1},\dots,M_{c,d}\right)\in\mathcal{L}(X,L^{2}([0,T],L^{2})^{d})$, with $M_{c,i}h\coloneqq-h\,S_{i}c$ for $h\in X$, there holds*

[TABLE] 2. (2)

Further, $(\mathbf{S}^{\prime}c)^{*}=\mathbf{M}_{c}^{*}\circ\mathbf{L}_{c}^{*}\in\mathcal{L}(L^{2}([0,T],L^{2})^{d},X)$ and for $h\in L^{2}([0,T],L^{2})^{d}$ there holds

[TABLE]

Let us now determine a numerically computable representation of the adjoint $L_{c}^{\ast}$ between $L^{2}([0,T],L^{2})$ that is as usual characterized for $z$ and $f\in L^{2}([0,T],L^{2})$ with $w\coloneqq L_{c}f$ by

[TABLE]

As $L_{c}^{*}z(t)$ belongs to $H^{1}_{0}$ we can replace the right-hand side by the weak formulation for $w$,

[TABLE]

Two partial integrations in time show by the initial conditions for $w$ that

[TABLE]

Hence, (6.2) is fulfilled if $h=L_{c}^{*}z\in Y$ is the weak solution to

[TABLE]

with zero Dirichlet boundary values and zero end conditions $h(T)=h^{\prime}(T)=0$. Due to the theory in Section 2 the latter differential equation is uniquely solvable and $z\mapsto h$ defines a bounded linear operator on $L^{2}([0,T],L^{2})$ that is numerically evaluated in the same way as $S$.

We now turn to the transpose $M_{c}^{\ast}:\,L^{2}([0,T],L^{2})\to X^{\prime}$ of the multiplication operator $M_{c}$. (The dual space $X^{\prime}$ is computed for a weighted inner product of $H^{2}([0,T];L^{2})$, see below.) Writing $M_{c}^{*}z=w_{z}\in X^{\prime}$ for $z\in L^{2}([0,T],L^{2})$, the function $w_{z}\in X$ has to satisfy the equality

[TABLE]

Loosely speaking, $w_{z}$ can hence be interpreted as a smoothed version of $-uz$. As $M_{c}^{*}$ is the only operator in our reconstruction scheme mapping into $X$, it actually controls smoothing of the searched-for parameter. It is practical to steer this smoothing by weights $\alpha,\beta>0$ that define the following duality product, extending the analogous weighted inner product of $L^{2}([0,T],L^{2})$,

[TABLE]

for all $g\in X$, the dual Lebesgue index $q\in[1,2]$ such that $1/p+1/{q}=1$, and $f\in X^{\prime}=H^{2}([0,T],L^{2})\cap L^{2}([0,T],L^{q})$. Thus, we obtain from (6.4) that

[TABLE]

which is a weak formulation of a fourth-order differential equation in time. Discretization of the last problem for $w_{z}$ via finite elements seems most natural but is indeed tedious as conforming finite element spaces need to be $H^{2}$-smooth in time, which is typically not pre-coded in open finite element packages.

Following the latter idea by via finite-dimensional subspaces of $X^{\prime}$ leads into Banach-space valued regularization schemes that we do merely for simplicity not consider in this paper. Instead, we formally use a simpler finite difference scheme in each individual spatial degree of freedom that arises by first discretizing the latter problem in space: As in the last section, we represent test functions $h\in X$ as $h(t,x)=\sum_{k=1}^{K}\vec{h}(t)_{k}\varphi_{k}(x)$ such that $h(t)\in V_{h}$ for every $t\in[0,T]$ with some $\vec{h}=(\vec{h}_{k})_{k=1,\dots,K}\in H^{2}([0,T])^{K}$. In the same way we define $\vec{u},\vec{z}\in L^{2}([0,T])^{K}$, $\vec{w}_{z}\in H^{2}([0,T])^{K}$, and recall the mass matrix $M$ from (5.4). We denote by \vec{u}\text{\raisebox{-1.72218pt}{\scalebox{1.75}{\cdot}}}\vec{z} the component-wise multiplication of $\vec{u}$ and $\vec{z}$, and deduce by (formal) partial integration that

[TABLE]

The above equation is fulfilled if $v_{k}\coloneqq(M\vec{w}_{z})_{k}$ solves for $k=1,\dots,K$ the one-dimensional ordinary differential equations

[TABLE]

In our numerical examples, we solve these systems at the time points $t_{i}=i\,\Delta t$, $i=0,\dots,N-1$, that we already fixed when solving for $u\in Y$ or representing $c\in X$. We continue to use the notation $v_{k}^{i}\coloneqq v_{k}(t_{i})$ and replace the appearing time derivatives by the standard central difference quotients up to order four, see, e.g. [For88]. The resulting fully discrete equations at the time points $t_{i}$ then read

[TABLE]

The latter equation requires \enquoteimaginary nodes at $t_{-2},t_{-1},t_{N}$ and $t_{N+1}$ that enforce the boundary conditions at $i=0$ (for $t=0$) and $i=N$ (for $t=T$) in the second line of (6.6), i.e.

[TABLE]

To sum up, our scheme for the numerical evaluation of $M_{c}^{*}z$ works as follows:

Compute the matrix-vector products M(\vec{u}^{i}\text{\raisebox{-1.72218pt}{\scalebox{1.75}{\cdot}}}\vec{z}^{i}) for $i=0,\dots,N-1$. 2. 2.

Solve the $N+4$-dimensional linear system (6.7)-(6.11) for every degree of freedom $k=1,\dots,K$ in space and as a result obtain the vector $(v_{k}^{i})_{i=0,\dots,N-1}$. 3. 3.

Determine the coefficients of $\vec{w}_{z}(t_{i})$ of $(M_{c}^{*}z)(t_{i})$ with respect to the basis of $V_{h}$ as the solution of $M\vec{w}_{z}(t_{i})=(v^{i}_{k})_{k=1,\dots,K}$ for every time step $t_{i}$.

In our numerical examples in Section 8 the number $K$ of degrees of freedom in space will be much higher than the number $N$ of time steps. (Typical orders of magnitudes for $n=3$ are $N\approx 200$ and $K\approx 5000$.) Thus, the execution of step 1 in the scheme above as well as the solution of $N$ linear systems with the sparse mass matrix in step 3 can be done fast compared to, e.g. the numerical solution of a forward problem. Since the matrix of the $(N+4)\times(N+4)$ dimensional system (6.7)-(6.11) does not change throughout the reconstruction it can be factored once and then used for the fast solution of the $K$ equations in step 2. We use the library SuperLU [Li05] for this, which exploits sparsity of the linear system.

Note that the computation of $M_{c}^{*}$ would greatly simplify if $c$ did not have to be so smooth in time. If, e.g. the solution operator $S$ turned out to be well-defined in some open subset of $H^{1}([0,T],L^{2})$, then the simpler equation v_{k}-\alpha\,v_{k}^{\prime\prime}=-(\,M\left(\vec{u}\text{\raisebox{-1.72218pt}{\scalebox{1.75}{\cdot}}}\vec{z}\right))_{k} would arise in (6.6) for $\alpha>0$, subject to homogeneous end conditions. This problem could be easily approximated with finite elements. If the entire setting even required no smoothness of $c$ at all, then the variant of the multiplication operator $M_{c}$ operating on $L^{2}([0,T],L^{2})$ would even become self-adjoint.

7 Regularization using inexact Newton iterations

In the preceding sections we showed well-definedness, continuity and differentiability of the solution operator $S$ and the measurement operator $\Phi$. We assume now that

[TABLE]

to take a look at a particular regularization scheme that stably approximate $c^{+}$ from data. More precisely, we propose inversion by the REGINN (\enquoteREGularization based on INexact Newton iteration) algorithm, which was stated and analyzed by Rieder [Rie99]. We give a brief reminder how REGINN works, by considering merely the first inverse problem in (7.1a).

In the entire section we actually neglect that the pre-image space $X$ is a Banach- instead of a Hilbert space. In our numerical experiments, we instead use $H^{2}([0,T],L^{2})\supset X$ as Hilbert space, as more involved schemes in Banach spaces are out of the scope of this paper.

Of course, we do not assume that data $g^{+}$ can be measured exactly but instead suppose to know some noisy version $g^{\varepsilon}$ with relative noise level $\varepsilon>0$, i.e. ${{{{\left\|g^{+}-g^{\varepsilon}\right\|}}}}\leq\varepsilon{{{{\left\|g^{+}\right\|}}}}\approx\varepsilon{{{{\left\|g^{\varepsilon}\right\|}}}}$. As is customary, we assume to know $\varepsilon>0$ a-priori. We already mentioned that REGINN relies on successive linearization of (7.1a) starting with some initial guess $c_{0}\in\mathcal{D}(S)$ to generate a sequence $(c_{k})_{k\in\mathbb{N}_{0}}\subset\mathcal{D}(S)$ of approximations of $c^{+}$. Writing $c^{+}=c_{k}+s_{k}^{+}$ for each $k\in\mathbb{N}_{0}$, the best update $s_{k}^{+}$ solves

[TABLE]

Because of the linearization error $E(c^{+},c_{k})$ and the exact data $g^{+}$ we only know a perturbed right-hand side $b_{k}^{\delta}=g^{\varepsilon}-\Phi c_{k}$. Its noise level $\delta$ is also unknown, as we only have $\|b_{k}^{\delta}-b_{k}^{+}\|\leq\varepsilon{{{{\left\|g^{\varepsilon}\right\|}}}}+\mathcal{O}(\|c^{+}-c_{k}\|^{2})$. REGINN applies a regularization method for linear inverse problems to this problem and stops it when the relative linear residuum is smaller than a tolerance times the non-linear residuum. In our case the former is done via the method of conjugate gradients (CG), which creates an inner iteration that computes a sequence of approximations $(s_{k,i})_{i\in\mathbb{N}}$ of $s_{k}^{+}$. The stopping is done by choosing tolerances $\mu_{k}\in(0,1)$ and picking $s_{k}\coloneqq s_{k,i_{k}}$ with

[TABLE]

Afterwards we can set $c_{k+1}\coloneqq c_{k}+s_{k}$ and continue the iteration, which we stop using the discrepancy principle by a fixed parameter $\tau>1$,

[TABLE]

The combination of REGINN with CG as inner regularization method was also analyzed by Rieder [Rie05]. Convergence is only guaranteed if the $\mu_{k}$ stay in the interval $[a,b]\subset(0,1)$ where $a$ and $b$ depend on unknown constants, like $\eta$ in the non-linearity condition. Due to the shrinking linearization error we want to be able to reduce $\mu_{k}$ during the outer iteration. On the other hand this reduction should not increase the computing time (number of CG-steps) of the next outer step too much. Rieder proposes the following strategy in [Rie99]: Start with $\mu_{1}=\mu_{2}=\mu_{\text{start}}\in(0,1)$ and for $k\geq 3$ define

[TABLE]

The tolerance $\mu_{k}$ is then set to

[TABLE]

where $\mu_{\text{max}}\in(\mu_{\text{start}},1)$. This achieves a $\gamma$-linear reduction of $\mu_{k}$ with $\gamma\in(0,1)$ if the number of inner steps is decreasing. We use $\mu_{\text{start}}=0.7$, $\gamma=0.9$ and $\mu_{\text{max}}=0.99$. Algorithm 1 lists a pseudo-code for the whole reconstruction procedure.

8 Numerical examples

We want to show in this last section that REGINN with the CG-iteration is indeed able to provide an estimate of a time- and space-dependent parameter $c$ in acceptable time, especially for $n=3$. To this end, we set $T\coloneqq 2$, $\mathds{1}\coloneqq(1)_{i=1,\dots,n}$, and reconstruct two different parameters from artificial data measured in $\Omega\coloneqq(0,1)^{n}$ with $n\in\{1,2,3\}$. The first parameter is hat-shaped and moves in time from $1/4\cdot\mathds{1}$ to $3/4\cdot\mathds{1}$,

[TABLE]

This parameter is smooth in time and space, at least in theory. Since spatial smoothness is actually not required by the parameter space $X$, we also test a parameter with discontinuities,

[TABLE]

In all calculations we use the finite element interpolation of $c_{\text{plateau}}$ in $V_{h}$, which is of course continuous but has a sharp edge at $\partial B(\tfrac{1}{2}\mathds{1},\tfrac{1}{4})$.

Motivated by possible applications we set $d\coloneqq 2^{n}$ at positions into the domain $\Omega$; each of the elements of their position vectors $x^{a}_{k}$ for $k\in\{1,\dots,d\}$ either equals to $1/3$ or 2/3. We further consider that each actuators excites a wave in $\Omega$ that we model by $d$ right-hand sides $f_{1},\dots,f_{d}$. Precisely, for frequency $\omega=8\pi$ and actuator radius $r^{a}=0.1$ we define ${f_{k}}\,\colon\,{\mathbb{R}\times\mathbb{R}^{3}}\to{\mathbb{R}}$ by

[TABLE]

and $f_{k}(t,x)=0$ else.

For the discretization we define $\Delta t\coloneqq 10^{-2}$ and employ a spatial grid consisting of $6$, $5$ or $4$ global refinements of the trivial triangulation of $(0,1)^{n}$ for $n=1,2$ or $3$, respectively. The finite element interpolations of both parameters evaluated at $t=1.7$ are shown in Figure 2. (Here and in subsequent figures we restrict ourselves to $n=2$.)

In (6.5) we introduced numbers $\alpha$ and $\beta$ in front of the first and second order terms of the $H^{2}([0,T],L^{2})$ scalar product in order to control the smoothness of the reconstruction. Our primary goal is to minimize the $L^{2}([0,T],L^{2})$-error to the exact parameter. For this we chose $\alpha$, $\beta$ in such a way that numerical approximations of $\alpha\|c^{\prime}\|^{2}$ and $\beta\|c^{\prime\prime}\|^{2}$ are one order of magnitude smaller than $\|c\|^{2}$. Tests with both parameters led us to define $\alpha\coloneqq 2\cdot 10^{-2}$ and $\beta\coloneqq 2\cdot 10^{-3}$.

We start our numerical experiments by checking whether the reconstruction $c_{k^{*}(\varepsilon)}$ converges to the exact parameter in the $X$- or the $L^{2}([0,T],L^{2})$-norm when $\varepsilon$ tends to [math]. We do so by applying REGINN to artificial data $u^{\varepsilon}\in L^{2}([0,T],L^{2})^{d}$ with relative noise level $\varepsilon>0$, i.e. ${\|u^{\varepsilon}-u\|}_{L^{2}([0,T],L^{2})^{d}}=\varepsilon{{{{\left\|u\right\|}}}}$ for $u=\mathbf{S}c$. The additive noise consists of a scaled vector of uniformly distributed pseudo-random numbers in $[-1,1]$. The stopping index $k^{*}(\varepsilon)$ is determined by the discrepancy principle with $\tau=2$, see (7.3).

For $\varepsilon=10^{-2}$ both reconstructions are satisfactory, as can be seen in Figure 3. Although the time dependence can already be deduced from these reconstructions, the $L^{2}$-errors are relatively high and amount to $45\%$ for $c_{\text{hat}}$ and $37\%$ when estimating $c_{\text{plateau}}$. The corresponding values in one and three spatial dimensions as well as the $H^{2}$-errors are listed in Table 1. In all dimensions the $H^{2}$-error for the moving hat is much higher than the $L^{2}$-error; for the other parameter both norms yield similar values.

Figure 4 shows the dependence of these errors on $\varepsilon$ in a logarithmic scale. While both errors clearly converge for $c_{\text{plateau}}$, this is at least questionable for the $H^{2}$-error of the moving hat. This leads to the hypothesis that $c_{\text{hat}}$ is not sufficiently smooth in time, which also seems to be the case in one and three space dimensions, see Table 2. The $L^{2}$-error is approximately of order $\mathcal{O}(\varepsilon^{0.3})$, at least for $\varepsilon\in[2.5\cdot 10^{-4},5\cdot 10^{-2}]$. For very small $\varepsilon$ we expect a saturation of the error due to the fixed discretization. In the case of $c_{\text{hat}}$ the behavior of the $L^{2}$-error in Figure 4 already hints at this effect for $\varepsilon\leq 10^{-3}$.

Now we turn to reconstructing $c$ from incomplete noisy measurements $g^{\varepsilon}\approx\Phi c\in(\mathbb{R}^{l})^{d}$, considering two measurement setups. The first one consists of $5^{n}$ sensors which are grid-like distributed in $\Omega$, as shown in Figure 5a. Each sensor generates measurements for $20$ equidistant times in $(0,T)$. This defines the space-time-positions $(t_{i}^{s},x_{i}^{s})_{i=1,\dots,l}$ of $l=20\cdot 5^{n}$ measurement points. We furthermore set $r_{x}\coloneqq 0.05$ and $r_{t}\coloneqq 0.02$.

In a real application it might be impossible (or inaffordable) to fill the whole domain with sensors. For $n=2$ we simulate this by placing the $25$ sensors on the left and lower edges of the domain. To avoid overlap we reduce the sensor radius in this case from $0.05$ to $0.035$.

The reconstructions for $n=2$ from $2000$ values ($20\cdot 5^{2}$ values for each of the $2^{2}$ right-hand sides) in the grid-like setting with $1\%$ artificial noise are shown in Figure 6. They look very similar to the reconstructions in Figure 3, where the whole wave (about $1.6\cdot 10^{6}$ degrees of freedom) was available. The errors are listed in Table 3 and most of them are only slightly higher than the corresponding entry of Table 1. For $1\%$ noise this measurement setup seems to be sufficient to obtain roughly the same reconstruction quality as from $u^{\varepsilon}$.

When repositioning the sensors as shown in Figure 5b the quality of the reconstruction decreases, as can be seen in Figure 7. The reconstruction of the moving hat $c_{\text{hat}}$ only achieves an $L^{2}([0,T],L^{2})$ error of $66.29\%$, which is significantly higher than the error of $46.11\%$ when using the $5\times 5$ grid. For $c_{\text{plateau}}$ the results are more encouraging, $41.79\%$ compared to $34.48\%$. The reconstruction quality suffers in particular for $t>1.5$, when the reconstruction vanishes in a neighborhood of the corner $(1,1)$ farthest from the sensors. This is of course a consequence of the finite speed of propagation. For smaller $\varepsilon=10^{-4}$, $c_{\text{plateau}}$ is better approximated ($26.9\%$ error), but the error for $c_{\text{hat}}$ remains at $53.6\%$.

We wish to remark that the discretization was chosen in such a way that computations for $n=3$ and small $\varepsilon$ can be done in affordable time. Typical computing times range for $\varepsilon=10^{-2}$ ranged from $5$ seconds ($n=1$) to $25$ minutes ($n=3$) when measured on an Intel i7–2600 CPU and required between $50$ MiB and $750$ MiB of memory.