Analysis of the error in constitutive equation approach for time-harmonic elasticity imaging

TL;DR

This paper develops a theoretical foundation for the MECE approach in time-harmonic elasticity imaging, explaining its robustness and properties, especially under incomplete boundary conditions, with implications for practical inverse problems.

Contribution

It proves the existence and uniqueness of solutions for the coupled system in MECE formulations, even with incomplete boundary data, and links MECE to classical least squares and ECE methods.

Findings

Unique and stable solutions at any frequency with sufficient data

Applicability to partial interior data in elastography

Convergence and differentiability of finite element discretizations

Abstract

We consider the identification of heterogeneous linear elastic moduli in the context of time-harmonic elastodynamics. This inverse problem is formulated as the minimization of the modified error in constitutive equation (MECE), an energy-based cost functional defined as an weighted additive combination $\mathcal{E}+\kappa\mathcal{D}$ of the error in constitutive equation (ECE) $\mathcal{E}$, expressed using an energy seminorm, and a quadratic error term $\mathcal{D}$ incorporating the kinematical measurements. MECE-based identification are known from existing computational evidence to enjoy attractive properties such as improved convexity, robustness to resonant frequencies, and tolerance to incompletely specified boundary conditions (BCs). The main goal of this work is to develop theoretical foundations, in a continuous setting, allowing to explain and justify some of the aforementioned beneficial properties, in particular addressing the general case where BCs may be underspecified. A specific feature of MECE formulations is that forward and adjoint solutions are governed by a fully coupled system, whose mathematical properties play a fundamental role in the qualitative and computational aspects of MECE minimization. We prove that this system has a unique and stable solution at any frequency, provided data is abundant enough (in a sense made precise therein) to at least compensate for any missing information on BCs. As a result, our formulation leads in such situations to a well-defined solution even though the relevant forward problem is not \emph{a priori} clearly defined. This result has practical implications such as applicability of MECE to partial interior data (with important practical applications including ultrasound elastography), convergence of finite element discretizations and differentiability of the reduced MECE functional. In addition, we establish that usual least squares and pure ECE formulations are limiting cases of MECE formulations for small and large values of $\kappa$, respectively. For the latter case, which corresponds to exact enforcement of kinematic data, we furthermore show that the reduced MECE Hessian is asymptotically positive for any parameter perturbation supported on the measurement region, thereby corroborating existing computational evidence on convexity improvement brought by MECE functionals. Finally, numerical studies that support and illustrate our theoretical findings, including a parameter reconstruction example using interior data, are presented.

Equations

$$\text{mbox}div\mathbf{\sigma }+\mathbf{b}=-\rho {\omega }^2\mathbf{u}\text{in }\Omega ,\mathbf{\sigma }\cdot \mathbf{n}=\mathbf{g}\text{on }{\Gamma }_{\text{N}},$$

$$\bm{u}=\mathbf{0}\quad\text{on }\Gamma_{\text{D}},\qquad\bm{\varepsilon}[\bm{u}]=\mbox{$\frac{1}{2}$}(\bm{\nabla}\bm{u}+\bm{\nabla}\bm{u}^{\text{\scriptsize T}})\quad\text{in }\Omega,$$

$$\mathbf{\sigma }=\text{boldmathC}:\mathbf{\epsilon }[\mathbf{u}]\text{in }\Omega ,$$

$$\big{(}\hskip 1.00006pt\bm{a},\bm{b}\hskip 1.00006pt\big{)}:=\int_{\Omega}\bm{a}\!:\!\overline{\bm{b}}\;\text{d}V=\int_{\Omega}a_{ij}\overline{b}_{ij}\;\text{d}V,$$

$$\big{(}\hskip 1.00006pt\bm{\sigma},\bm{\varepsilon}[\widetilde{\bm{w}}{}]\hskip 1.00006pt\big{)}-\omega^{2}\big{(}\hskip 1.00006pt\rho\bm{u},\widetilde{\bm{w}}{}\hskip 1.00006pt\big{)}-\big{(}\hskip 1.00006pt\bm{\sigma}\!\cdot\!\bm{n},\widetilde{\bm{w}}{}\hskip 1.00006pt\big{)}_{\Gamma\setminus\Gamma_{\text{N}}}=\big{\langle}\hskip 1.00006pt\bm{f},\widetilde{\bm{w}}{}\hskip 1.00006pt\big{\rangle}_{\mathcal{V}^{\prime},\mathcal{V}}\qquad\text{for all }\widetilde{\bm{w}}{}\in\mathcal{V}:=H^{1}(\Omega;\mathbb{R}^{d}),$$

$$\mathcal{U}$$

$$\mathcal{S}(\mathbf{u})$$

$$\mathcal{Q}$$

$$\mathcal{E}(\mathbf{u},\mathbf{\sigma },\text{boldmathC}):=\frac{1}{2}{\int }_{\Omega }(\mathbf{\sigma }-\text{boldmathC}:\mathbf{\epsilon }[\mathbf{u}]):{\text{boldmathC}}^{-1}:(\mathbf{\sigma }-\text{boldmathC}:\mathbf{\epsilon }[\mathbf{u}])\text{d}V$$

$${\Lambda }_{\kappa }(\mathbf{u},\mathbf{\sigma },\text{boldmathC}):=\mathcal{E}(\mathbf{u},\mathbf{\sigma },\text{boldmathC})+\kappa \mathcal{D}(\mathbf{u}-{\mathbf{u}}^{\text{m}},\mathbf{u}-{\mathbf{u}}^{\text{m}}),$$

$$(\mathbf{u},\mathbf{\sigma },\text{boldmathC}):=\underset{\mathbf{v}\in \mathcal{U},\mathbf{\tau }\in \mathcal{S}(\mathbf{v}),\text{boldmathA}\in \mathcal{Q}}{\mathrm{argmin}}{\Lambda }_{\kappa }(\mathbf{v},\mathbf{\tau },\text{boldmathA}).$$

$$\mathcal{L}(\bm{u},\bm{w},\bm{\sigma},\text{\boldmath$\mathcal{C}$}):=\Lambda_{\kappa}(\bm{u},\bm{\sigma},\text{\boldmath$\mathcal{C}$})-\Re\big{\{}\hskip 1.00006pt\big{(}\hskip 1.00006pt\bm{\sigma},\bm{\varepsilon}[\bm{w}]\hskip 1.00006pt\big{)}-\omega^{2}\big{(}\hskip 1.00006pt\rho\bm{u},\bm{w}\hskip 1.00006pt\big{)}-\big{(}\hskip 1.00006pt\bm{\sigma}\!\cdot\!\bm{n},\bm{w}\hskip 1.00006pt\big{)}_{\Gamma\setminus\Gamma_{\text{N}}}-\big{\langle}\hskip 1.00006pt\bm{f},\bm{w}\hskip 1.00006pt\big{\rangle}_{\mathcal{V}^{\prime},\mathcal{V}}\hskip 1.00006pt\big{\}},$$

$${\partial }_{\sigma }\mathcal{L}(\mathbf{u},\mathbf{w},\mathbf{\sigma },\text{boldmathC})[\tilde{\mathbf{\sigma }}]$$

$${\partial }_w\mathcal{L}(\mathbf{u},\mathbf{w},\mathbf{\sigma },\text{boldmathC})[\tilde{\mathbf{w}}]$$

$${\partial }_u\mathcal{L}(\mathbf{u},\mathbf{w},\mathbf{\sigma },\text{boldmathC})[\tilde{\mathbf{u}}]$$

$${\partial }_{\text{boldmathC}}\mathcal{L}(\mathbf{u},\mathbf{w},\mathbf{\sigma },\text{boldmathC})[\widehat{\text{boldmathC}}]$$

$$\Re\big{\{}\hskip 1.00006pt\big{(}\hskip 1.00006pt\tilde{\bm{\sigma}},\text{\boldmath$\mathcal{C}$}^{-1}\!:\!\bm{\sigma}-\bm{\varepsilon}[\bm{u}+\bm{w}]\hskip 1.00006pt\big{)}\hskip 1.00006pt\big{\}}+\Re\big{\{}\hskip 1.00006pt\big{(}\hskip 1.00006pt\tilde{\bm{\sigma}}\!\cdot\!\bm{n},\bm{w}\hskip 1.00006pt\big{)}_{\Gamma\setminus\Gamma_{\text{N}}}\hskip 1.00006pt\big{\}}=0\qquad\text{for all }\tilde{\bm{\sigma}}\in\mathcal{S}.$$

$$\mathbf{\sigma }=\text{boldmathC}:\mathbf{\epsilon }[\mathbf{u}+\mathbf{w}]\text{in }\Omega ,$$

$$\mathbf{w}=0\text{on }\Gamma \setminus {\Gamma }_{\text{N}},$$

$$\mathcal{W}=\big{\{}\hskip 1.00006pt\bm{v}\in\mathcal{V},\;\bm{v}=\mathbf{0}\text{ on }\Gamma_{\text{D}}\cup\Gamma_{c}\hskip 1.00006pt\big{\}}.$$

$$\big{(}\hskip 1.00006pt\bm{\varepsilon}[\bm{u}+\bm{w}],\text{\boldmath$\mathcal{C}$}\!:\!\bm{\varepsilon}[\widetilde{\bm{w}}{}]\hskip 1.00006pt\big{)}-\omega^{2}\big{(}\hskip 1.00006pt\rho\bm{u},\widetilde{\bm{w}}{}\hskip 1.00006pt\big{)}=\mathcal{F}(\widetilde{\bm{w}}{}),\quad\text{for all }\widetilde{\bm{w}}{}\in\mathcal{W}.$$

$$\partial_{u}\mathcal{L}(\bm{u},\bm{w},\bm{\sigma},\text{\boldmath$\mathcal{C}$})[\widetilde{\bm{u}}{}]=\Re\big{\{}\hskip 1.00006pt\big{(}\hskip 1.00006pt\bm{\varepsilon}[\bm{u}]-\text{\boldmath$\mathcal{C}$}^{-1}\!:\!\bm{\sigma}),\text{\boldmath$\mathcal{C}$}\!:\!\bm{\varepsilon}[\widetilde{\bm{u}}{}]\hskip 1.00006pt\big{)}+\kappa\mathcal{D}\big{(}\hskip 1.00006pt\bm{u}-\bm{u}^{\text{\scriptsize m}},\widetilde{\bm{u}}{}\hskip 1.00006pt\big{)}+\omega^{2}\big{(}\hskip 1.00006pt\rho\widetilde{\bm{u}}{},\bm{w}\hskip 1.00006pt\big{)}\hskip 1.00006pt\big{\}},$$

$$\big{(}\hskip 1.00006pt\bm{\varepsilon}[\bm{w}],\text{\boldmath$\mathcal{C}$}\!:\!\bm{\varepsilon}[\widetilde{\bm{u}}{}]\hskip 1.00006pt\big{)}-\omega^{2}\big{(}\hskip 1.00006pt\rho\widetilde{\bm{u}}{},\bm{w}\hskip 1.00006pt\big{)}-\kappa\mathcal{D}\big{(}\hskip 1.00006pt\bm{u}-\bm{u}^{\text{\scriptsize m}},\widetilde{\bm{u}}{}\hskip 1.00006pt\big{)}=0\quad\text{for all }\widetilde{\bm{u}}{}\in\mathcal{U},$$

$$\big{(}\hskip 1.00006pt\bm{\varepsilon}[\bm{u}]\otimes\bm{\varepsilon}[\bm{u}]-(\text{\boldmath$\mathcal{C}$}^{-1}\!:\!\bm{\sigma})\otimes(\text{\boldmath$\mathcal{C}$}^{-1}\!:\!\bm{\sigma})\,,\,\hat{\text{\boldmath$\mathcal{C}$}}\hskip 1.00006pt\big{)}=0\qquad\text{for all }\hat{\text{\boldmath$\mathcal{C}$}}\in\mathcal{Q}$$

$$\tilde{\mathcal{E}}(\text{\boldmath$\mathcal{C}$}):=\mathcal{E}(\bm{u},\bm{\sigma},\text{\boldmath$\mathcal{C}$})=\big{(}\hskip 1.00006pt\bm{\varepsilon}[\bm{w}],\text{\boldmath$\mathcal{C}$}\!:\!\bm{\varepsilon}[\bm{w}]\hskip 1.00006pt\big{)},\qquad\tilde{\mathcal{D}}(\text{\boldmath$\mathcal{C}$}):=\mathcal{D}(\bm{u}-\bm{u}^{\text{\scriptsize m}},\bm{u}-\bm{u}^{\text{\scriptsize m}}),$$

$$\underset{\text{boldmathC}\in \mathcal{Q}}{\min }\tilde{\Lambda }{}_{\kappa }(\text{boldmathC}),\tilde{\Lambda }{}_{\kappa }(\text{boldmathC}):={\Lambda }_{\kappa }(\mathbf{u},\mathbf{\sigma },\text{boldmathC})=\tilde{\mathcal{E}}(\text{boldmathC})+\kappa \tilde{\mathcal{D}}(\text{boldmathC}),$$

$$\tilde{\Lambda }{}_{\kappa }(\text{boldmathC})=\underset{\mathbf{u}\in \mathcal{U},\mathbf{\sigma }\in \mathcal{S}(\mathbf{u})}{\mathrm{argmin}}{\Lambda }_{\kappa }(\mathbf{u},\mathbf{\sigma },\text{boldmathC}).$$

$$\tilde{\Lambda }{}_{\kappa }({\text{boldmathC}}^{\star })={\Lambda }_{\kappa }({\mathbf{u}}^{\star },{\mathbf{w}}^{\star },{\text{boldmathC}}^{\star })\le {\Lambda }_{\kappa }({\mathbf{u}}_{\mathcal{C}},{\mathbf{w}}_{\mathcal{C}},\text{boldmathC})\le {\Lambda }_{\kappa }(\mathbf{u},\mathbf{w},\text{boldmathC})$$

$$\text{Find }(\bm{w},\bm{u})\in\mathcal{W}\times\mathcal{U},\quad\left\{\begin{aligned} \text{(a)}&&\mathcal{A}(\bm{w},\widetilde{\bm{w}}{},\text{\boldmath$\mathcal{C}$})+\mathcal{B}(\bm{u},\widetilde{\bm{w}}{},\text{\boldmath$\mathcal{C}$})&=\big{\langle}\hskip 1.00006pt\bm{f},\widetilde{\bm{w}}{}\hskip 1.00006pt\big{\rangle}_{\mathcal{W}^{\prime},\mathcal{W}}&&\text{for all }\widetilde{\bm{w}}{}\in\mathcal{W}\\ \text{(b)}&&\mathcal{B}(\widetilde{\bm{u}}{},\bm{w},\text{\boldmath$\mathcal{C}$})-\kappa\mathcal{D}(\bm{u},\widetilde{\bm{u}}{})&=-\kappa\big{\langle}\hskip 1.00006pt\bm{d},\widetilde{\bm{u}}{}\hskip 1.00006pt\big{\rangle}_{\mathcal{U}^{\prime},\mathcal{U}}&&\text{for all }\widetilde{\bm{u}}{}\in\mathcal{U}\end{aligned}\right.$$

$$(\mathbf{v},\tilde{\mathbf{v}}{)}_{\omega }=\mathcal{A}(\mathbf{v},\tilde{\mathbf{v}})+\rho {\omega }^2(\mathbf{v},\tilde{\mathbf{v}}),\Vert \mathbf{v}{\Vert }_{\omega }^2:=(\mathbf{v},\mathbf{v}{)}_{\omega }.$$

$$\mathcal{A}(\mathbf{w},\tilde{\mathbf{w}})\le a\Vert \mathbf{w}{\Vert }_{\omega }\Vert \tilde{\mathbf{w}}{\Vert }_{\omega },\mathcal{B}(\mathbf{u},\tilde{\mathbf{w}})\le b\Vert \mathbf{u}{\Vert }_{\omega }\Vert \tilde{\mathbf{w}}{\Vert }_{\omega },\mathcal{D}(\mathbf{u},\tilde{\mathbf{u}})\le d\Vert \mathbf{u}{\Vert }_{\omega }\Vert \tilde{\mathbf{u}}{\Vert }_{\omega }$$

$$B:\mathcal{U}\to\mathcal{W}{}^{\prime},\;\;\big{\langle}\hskip 1.00006ptB\bm{u},\widetilde{\bm{w}}{}\hskip 1.00006pt\big{\rangle}_{\mathcal{W}^{\prime},\mathcal{W}}=\mathcal{B}(\bm{u},\widetilde{\bm{w}}{})\;\;\text{for all }(\bm{u},\widetilde{\bm{w}}{})\in\mathcal{U}\times\mathcal{W},$$

$$\big{\langle}\hskip 1.00006pt\widetilde{\bm{u}}{},B^{\text{t}}\bm{w}\hskip 1.00006pt\big{\rangle}_{\mathcal{U},\mathcal{U}^{\prime}}=\big{\langle}\hskip 1.00006ptB\widetilde{\bm{u}}{},\bm{w}\hskip 1.00006pt\big{\rangle}_{\mathcal{W}^{\prime},\mathcal{W}}=\mathcal{B}(\widetilde{\bm{u}}{},\bm{w})\qquad\text{for all }(\widetilde{\bm{u}}{},\bm{w})\in\mathcal{U}\times\mathcal{W},$$

Full text

1. Introduction

Energy-based objective functionals are strong alternatives to conventional least square methods for various parameter identification problems. Such functionals, often called error in constitutive equation (ECE) functionals in the area of solid mechanics, were initially introduced in [22] for error estimation in the linear elastic FEM and in e.g. [20, 19] for electrical impedance tomography. ECE functionals have since been successfully used in various mechanical parameter identification problems under linear static [15], nonlinear quasistatic [24], time-harmonic [23, 21, 5, 4] and, more recently, transient conditions  [1, 14, 26, 8]. Mathematical and numerical issues are also discussed in e.g. [10, 17].

The main idea behind ECE approaches is to relax the constitutive equations connecting fluxes (e.g. stresses) and gradients of state variables (e.g strains). To this end, stresses (assumed to satisfy dynamic equilibrium) and strains (assumed to satisfy kinematical constraints) are treated as independent quantities. ECE functionals then measure residuals in constitutive equations evaluated on given stresses and strains, thereby assuming a physical meaning directly connected to constitutive parameter identification. In their original form, ECE functionals assume the measured displacements or strains to be strictly enforced as part of the kinematical admissibility constraints, but this is often undesirable as real data is usually polluted with noise. For this reason, ECE-based identification is nowadays rather formulated by means of so-called modified error in constitutive equation (MECE) functionals [23], where reliable and unreliable informations are treated differently. Equilibrium equations, initial conditions and known boundary conditions (BCs) are deemed reliable and enforced strictly (e.g. using Lagrange multipliers). By contrast, measured data, constitutive properties and (when applicable) imperfectly known BCs are deemed unreliable and incorporated as constitutive or observation residuals using an additive combination $\mathcal{E}+\kappa\mathcal{D}$ of ECE and least-squares components $\mathcal{E}$ and $\mathcal{D}$, with $\kappa$ a positive weight parameter.

Optimization problems that arise from MECE-based identification have been observed to display very attractive properties such as improved convexity, robustness to resonant frequencies, and tolerance to partially or completely unknown BCs. For instance, the MECE objective for transient elastodynamics was observed in [14] to be convex over a wide region of the parameter space. Convexity improvement relative to least-squares functionals was also reported in [9] in the context of linear elastostatics. Moreover, for time-harmonic conditions, the FEM-discretized coupled stationarity problem was shown in [4] to remain uniquely solvable at resonant prescribed frequencies. In addition, as shown in [12, 8], ECE approaches can naturally accommodate configurations with partially or completely unknown BCs, a feature also used in [6] for parameter identification using elastostatic interior data (i.e. unknown BCs). This makes such formulations perfectly suited to situations where interior data is abundant over a subset of the medium being probed. Important practical applications include elasticity imaging [16] where interior displacements are tracked inside soft tissue using ultrasound, but information about region boundaries is difficult to ascertain.

These advantages of MECE formulations over conventional least squares counterparts are backed by abundant numerical and experimental evidence, but to our best knowledge lack theoretical support in infinite-dimensional settings . Some of the existing analyses address discretized MECE formulations [4, 8], while [17] studied existence of solutions and convexity of a continuous ECE formulation using complete internal data. Accordingly, the main goal of this work is to develop theoretical foundations allowing to explain and justify some of the aforementioned beneficial properties of MECE formulations. We focus on elasticity imaging under time-harmonic conditions and adopt a Hilbert space setting. Moreover, our analysis addresses the general case where BCs may be underdetermined.

A specific feature of MECE formulations is that first-order optimality conditions lead to fully coupled forward and adjoint problems, rather than unidirectionally coupled problems arising in conventional least squares approaches. The mathematical properties of this coupled system play a fundamental role in the qualitative (e.g. existence of solutions, continuity and differentiability of the objective with respect to parameters) and numerical aspects of MECE-based imaging. In particular, our formulation leads to a well-defined solution in the underspecified BC case, for which a relevant forward problem is not a priori clearly defined. Hence, establishing the well-posedness of the coupled system is our first main goal. Treating this system as a perturbed mixed problem, we prove that it has a unique and stable solution at any frequency, subject to conditions that ensure that the available data more than compensates insufficient information on BCs. This result has important practical implications such as applicability of MECE to partial interior data, convergence of finite element discretizations and differentiability of the reduced MECE objective functional. Secondly, establish that least squares and pure ECE formulations are limiting forms for $\kappa\to 0$ and $\kappa\to\infty$, respectively, of MECE formulations. Thirdly, for the latter case (corresponding to exact enforcement of kinematic data), we show that the reduced MECE hessian becomes positive for any parameter perturbation supported on the measurement region. The latter result also has strong practical implications as convexity translates into robustness with respect to the initial guess in gradient-based optimization.

This work is organized as follows. In Section 2, we define the elasticity imaging problem and set notations, introduce the relevant MECE functional and state the relevant first-order stationarity conditions. Section 3 then addresses the well-posedness of the coupled stationarity problem, allowing for underdetermined BCs. Then, the limiting forms of the stationarity solution as $\kappa\to 0$ or $\kappa\to\infty$ are established in Section 4. The reduced MECE Hessian, and in particular its asymptotic convexity in the $\kappa\to\infty$ limit, is studied in Sec. 5. Then, numerical studies that support and illustrate our theoretical findings, including a parameter reconstruction example using interior data, are presented in Section 6. Finally, most of the proofs are deferred to Section 7.

2. Problem setting

Let a solid elastic body occupy a bounded and connected domain $\Omega\subset\mathbb{R}^{d}$ $(1\leq d\leq 3)$ with boundary $\Gamma$. The time-harmonic motion of this body is governed by (i) the balance equations

[TABLE]

where $\bm{u}$ is the displacement field, $\omega$ represents the specified angular frequency, $\rho$ denotes the known mass density, $\bm{b}$ is a given body force density, $\bm{\sigma}$ represents the stress tensor, $\bm{g}$ and $\Gamma_{\text{N}}\subseteq\Gamma$ are the given surface force density (traction) and its support, respectively, and $\bm{n}$ is the outward unit vector normal to $\Gamma$; (ii) the kinematic compatibility equations

[TABLE]

where $\bm{\varepsilon}[\bm{u}]$ denotes the linearized strain tensor associated with $\bm{u}$ and $\Gamma_{\text{D}}\subseteq\Gamma$ is the portion of the boundary where the displacement is prescribed; and (iii) the (linear elastic) constitutive relation

[TABLE]

where $\mathcal{C}$ is the fourth-order elasticity tensor field. For simplicity, the kinematical boundary condition specified in (2) is homogeneous; the case of a non-homogeneous boundary condition can be treated with minor modifications, as can loading arrangements other than those appearing in (1).

Here, the boundary subsets $\Gamma_{\text{N}}$ and $\Gamma_{\text{D}}$ are only required to not overlap (i.e., $\Gamma_{\text{N}}\cap\Gamma_{\text{D}}=\emptyset$) and may cover $\Gamma$ only partially (i.e., $\Gamma_{\text{N}}\cup\Gamma_{\text{D}}\subseteq\Gamma$). In other words, we may have $\Gamma_{c}\not=\emptyset$, where $\Gamma_{c}:=\Gamma\setminus(\Gamma_{\text{N}}\cup\Gamma_{\text{D}})$. In such event, equations (1)-(3) admit multiple solutions, whereas a unique solution exists (except for a countable set of eigenfrequencies $\omega$) when $\Gamma_{\text{N}}\cup\Gamma_{\text{D}}=\Gamma$. This non-standard boundary condition setting is adopted here to model experimental situations where full-field interior data (to be introduced thereafter) is available while boundary conditions are underspecified.

Measurements

In addition to the fundamental equations (1)-(3), we assume availability for the prescribed frequency $\omega$ of measured time-harmonic displacements $\bm{u}^{\text{\scriptsize m}}$ in $\Omega_{\text{\scriptsize m}}\subset\Omega$ (interior data), which is for example usual in elastography applications [16]. The forthcoming analyses can be adapted to accommodate measured displacements on a portion of $\Gamma\setminus\Gamma_{\text{D}}$.

Inverse problem

We address the inverse problem of reconstructing the elasticity tensor field $\mathcal{C}$ such that (i) the governing equations of motion (1)-(3) are satisfied, and (ii) $\mathcal{C}$ is consistent with the measurement $\bm{u}^{\text{\scriptsize m}}$, under prescribed time-harmonic conditions. Equality $\bm{u}=\bm{u}^{\text{\scriptsize m}}\text{ in }\Omega_{\text{\scriptsize m}}$ between the data $\bm{u}^{\text{\scriptsize m}}$ and its model counterpart $\bm{u}$, implicit in requirement (ii), should hold for error-free model and measurements but will be relaxed in the upcoming formulation to allow for expected uncertainties.

2.1. Weak formulation of motion

Let \big{(}\hskip 1.00006pt\bm{a},\bm{b}\hskip 1.00006pt\big{)} denote the $L^{2}(\Omega)$ scalar product of second-order tensor fields $\bm{a},\bm{b}\in L^{2}(\Omega;\mathbb{R}^{d\times d})$:

[TABLE]

where the overline means complex conjugation. Repeated indices imply summation wherever indicial notation is used. The $L^{2}(\Omega)$ scalar product of vector or scalar fields follows the same notational style with suitable adjustments; so does the $L^{2}(\Gamma)$ scalar product of fields defined on a surface $\Gamma$, denoted as \big{(}\hskip 1.00006pt\bm{a},\bm{b}\hskip 1.00006pt\big{)}_{\Gamma}. The weak formulation of the balance equations (1) then reads

[TABLE]

where the continuous linear functional $\bm{f}\in\mathcal{V}^{\prime}$ embodies all given excitations; for example \big{\langle}\hskip 1.00006pt\bm{f},\widetilde{\bm{w}}{}\hskip 1.00006pt\big{\rangle}_{\mathcal{V}^{\prime},\mathcal{V}}:=\big{(}\hskip 1.00006pt\bm{b},\bm{w}\hskip 1.00006pt\big{)}+\big{(}\hskip 1.00006pt\bm{g},\bm{w}\hskip 1.00006pt\big{)}_{\Gamma_{\text{N}}} for the given force densities $\bm{b},\bm{g}$ appearing in (1). In (5) and thereafter, $\langle\cdot,\cdot\rangle_{X^{\prime},X}$ denotes the duality pairing between a Hilbert space $X$ an its dual $X^{\prime}$ (i.e. \big{\langle}\hskip 1.00006pt\bm{f},\widetilde{\bm{w}}{}\hskip 1.00006pt\big{\rangle}_{X^{\prime},X} evaluates the linear functional $\bm{f}\in X^{\prime}$ at $\widetilde{\bm{w}}{}\in X$). In addition, let the spaces, $\mathcal{U}$, $\mathcal{S}(\bm{u})$ and $\mathcal{Q}$ of kinematically admissible displacements, dynamically admissible stresses and admissible elasticity tensor fields, respectively, be defined as

[TABLE]

(where $\mathbbm{Q}$ denotes the finite-dimensional vector space of fourth-order tensors $\mathcal{C}$ with major and minor symmetries, i.e. $\mathcal{C}_{ijkl}=\mathcal{C}_{klij}=\mathcal{C}_{jilk}$, and $c_{0}$ is some positive constant). Finally, the mass density field $\rho\in L^{\infty}(\Omega)$ must be bounded below by a positive constant.

2.2. MECE optimization problem

We follow the inversion approach initiated in [4, 27], whereby the foregoing inverse problem is formulated as an optimization problem in which the unknown elasticity tensor field $\mathcal{C}$ is estimated by minimizing an objective function that additively combines two error terms: (i) an error in constitutive equation (ECE) functional [22] defined by

[TABLE]

that measures (in units of energy) the discrepancy in the constitutive equation (3), and (ii) a quadratic error term $\mathcal{D}(\bm{u}-\bm{u}^{\text{\scriptsize m}},\bm{u}-\bm{u}^{\text{\scriptsize m}})$, where $\mathcal{D}$ is a positive bilinear form, that quantifies the mismatch between the predicted (or model) displacements and the measured ones. This objective function, hereafter called the modified ECE (MECE) functional, is defined by

[TABLE]

where $\kappa>0$ is a weight parameter. For a given triple (\bm{u},\bm{\sigma},\text{\boldmath\mathcal{C}})\in\mathcal{U}\times\mathcal{S}(\bm{u})\times\mathcal{Q} (with these sets as defined in (6a,6b,6c)), \Lambda_{\kappa}(\bm{u},\bm{\sigma},\text{\boldmath\mathcal{C}}) defines a quantitative measure of the consistency of these variables with (i) the constitutive equation, and (ii) the available measurements $\bm{u}^{\text{\scriptsize m}}$. Accordingly, the elasticity imaging inverse problem is formulated as the PDE-constrained optimization problem

[TABLE]

2.3. Stationarity conditions

We now collect the first-order stationarity conditions for the minimization problem (9). To this end, let the Lagrangian $\mathcal{L}:\mathcal{U}\times\mathcal{V}\times\mathcal{S}\times\mathcal{Q}\to\mathbb{R}$ be defined as

[TABLE]

where (i) $\bm{w}\in\mathcal{V}$ plays the role of the Lagrange multiplier, (ii) the constraint is the dynamic balance equation (5), and (iii) the term \big{(}\hskip 1.00006pt\bm{\sigma}\!\cdot\!\bm{n},\bm{w}\hskip 1.00006pt\big{)}_{\Gamma\setminus\Gamma_{\text{N}}} is crucial for the case in which $\Gamma_{\text{N}}\cup\Gamma_{\text{D}}\not=\Gamma$ (i.e., boundary conditions are not prescribed over the entire boundary). The first-order optimality conditions for the minimization problem  (9) are then given, in terms of first-order Gâteaux derivatives of $\mathcal{L}$, by

[TABLE]

where the test functions $\widetilde{\bm{u}}{}$ are constrained, consistently with the assumption $\bm{u}\in\mathcal{U}$, whereas $\widetilde{\bm{w}}{}$ are for now unconstrained. We begin by exploiting conditions (11a) and (11b). Condition (11a) yields

[TABLE]

while condition (11b) simply restates the balance constraint (5). Using first the subspace of $\mathcal{S}$ containing all $\tilde{\bm{\sigma}}$ with vanishing trace on $\Gamma$, equation (12) provides

[TABLE]

The second term in (12) then enforces the essential condition

[TABLE]

the stress-type unknown $\bm{\sigma}\!\cdot\!\bm{n}\mid_{\Gamma\setminus\Gamma_{\text{N}}}$ in the fourth term of (5) being the associated Lagrange multiplier, see e.g. [3]. We elect in this work to treat condition (14) via elimination, and hence seek $\bm{w}$ in the function space $\mathcal{W}$ defined by

[TABLE]

For the same reason, we now restrict equation (5) to test functions $\widetilde{\bm{w}}{}\in\mathcal{W}$, thereby also eliminating the Lagrange multiplier $\bm{\sigma}\bm{n}\mid_{\Gamma\setminus\Gamma_{\text{N}}}$, and obtain, after using (13):

[TABLE]

We note that $\mathcal{W}\subseteq\mathcal{U}$; more precisely, $\mathcal{W}\subset\mathcal{U}$ (with strict inclusion) when $\Gamma_{N}\cup\Gamma_{D}\not=\Gamma$ (insufficient boundary data), while $\mathcal{W}=\mathcal{U}$ when $\Gamma_{N}\cup\Gamma_{D}=\Gamma$ (sufficient boundary data).

Then, moving to the third condition (11c), we find

[TABLE]

so that, on substituting (13) and rearranging, (11c) yields

[TABLE]

Finally, condition (11d) yields the following equation, which is nonlinear in $\mathcal{C}$:

[TABLE]

Concluding, the first-order stationarity conditions for the minimization problem (9) consist of the coupled weak equations (16), (18), (19) governing (\bm{u},\bm{w},\text{\boldmath\mathcal{C}}), with $\bm{\sigma}$ then given explicitly by (13).

2.4. Reduced optimization problem

In a full-space approach, the stationarity system (13,16,18,19) is solved (iteratively) as a whole, as done e.g. in [13]. Alternatively, a reduced-space approach can be formulated from observing that equations (16,18) define for given $\mathcal{C}$ a linear (coupled) problem for $(\bm{u},\bm{w})$. Then, letting $(\bm{u},\bm{w})=(\bm{u}_{\mathcal{C}},\bm{w}_{\mathcal{C}})$ solve equations (16,18) for a given \text{\boldmath\mathcal{C}}\in\mathcal{Q}, reduced versions of the ECE and data-misfit components of $\Lambda_{\kappa}$ that depend only on $\mathcal{C}$ can be defined by

[TABLE]

where we have used \bm{\sigma}=\text{\boldmath\mathcal{C}}\!:\![\bm{u}+\bm{w}], see (13). Then, problem (9) can be recast in reduced minimization form as

[TABLE]

We observe that the functional \Lambda_{\kappa}(\bm{u},\bm{\sigma},\text{\boldmath\mathcal{C}}) is, for any fixed \text{\boldmath\mathcal{C}}\in\mathcal{Q}, differentiable and convex in $(\bm{u},\bm{\sigma})$ due to the requisite ellipticity of $\mathcal{C}$ (see (6c)). Problem (16,18) is hence equivalent to solving the partial minimization of \Lambda_{\kappa}(\bm{u},\bm{\sigma},\text{\boldmath\mathcal{C}}) with $\mathcal{C}$ given, and we have

[TABLE]

For that reason, we will henceforth refer to problem (16,18) as the stationarity system.

In this work, we adopt and study this reduced-space approach, which thanks to the characterization (22) can in fact be shown to be equivalent to the full-space approach in the following sense:

Lemma 2.1**.**

Problems (21) and (9) are equivalent: a solution to (21) also solves (9) and vice versa.

Proof 2.2**.**

First, let \text{\boldmath\mathcal{C}}^{\star} solve the reduced problem (21), and $(\bm{u}^{\star},\bm{w}^{\star})$ solve (16-18) for \text{\boldmath\mathcal{C}}=\text{\boldmath\mathcal{C}}^{\star}. Then:

[TABLE]

for any (\bm{u},\bm{w},\text{\boldmath\mathcal{C}})\in\mathcal{U}\times\mathcal{W}\times\mathcal{Q}, where the first inequality stems from (21) and the second from (22). Hence, (\bm{u}^{\star},\bm{w}^{\star},\text{\boldmath\mathcal{C}}^{\star}) solves (9).

Conversely, let (\bm{u}^{\sharp},\bm{w}^{\sharp},\text{\boldmath\mathcal{C}}^{\sharp}) solve problem (9). This triple then verifies (16-18) with \text{\boldmath\mathcal{C}}=\text{\boldmath\mathcal{C}}^{\sharp}, i.e. has the form (\bm{u}_{\mathcal{C}^{\sharp}},\bm{w}_{\mathcal{C}^{\sharp}},\text{\boldmath\mathcal{C}}^{\sharp}), so that \Lambda_{\kappa}(\bm{u}^{\sharp},\bm{w}^{\sharp},\text{\boldmath\mathcal{C}}^{\sharp})=\tilde{\Lambda}{}_{\kappa}(\text{\boldmath\mathcal{C}}^{\sharp}). Since in addition problem (21) consists in minimizing $\Lambda_{\kappa}$ over a subset of $\mathcal{U}\times\mathcal{W}\times\mathcal{Q}$, \text{\boldmath\mathcal{C}}^{\sharp} is a minimizer of $\tilde{\Lambda}{}_{\kappa}$. The proof is complete.

2.5. Coupled stationarity system: a key component of MECE-based imaging

The stationarity system (16,18) plays for several reasons a fundamental role in MECE-based imaging:

1. (a)

The definition (21) of a reduced-space approach needs well-posedness of the stationarity system; 2. (b)

The present MECE framework aims at treating situations with underspecified BCs, for which a relevant forward problem is not a priori clearly defined, in contrast with usual PDE-constrained inversion methods. Instead, the field $\bm{u}$ acts as a forward solution, provided problem (16,18) is well-posed. It is in particular important to determine conditions on the interior data ensuring that it compensates the insufficient BC information and make problem (16,18) well-posed. 3. (c)

In cases involving underspecified BCs, the relevant forward problem is not a priori clearly defined. Instead, the stationarity system (16,18) acts as a combination of the forward and adjoint problems. The latter are coupled in the present framework, while they are usually uncoupled in inverse problems solved using standard $L^{2}$ minimization [13, 25, 18] (see Sec. 5.5 for further elaboration). 4. (d)

Well-posedness of the stationarity system implies that the solution mapping \text{\boldmath\mathcal{C}}\mapsto(\bm{u},\bm{w}) is Fréchet differentiable (by virtue of the implicit function theorem, e.g. [11, Thm. 7.13-1], whose applicability can then readily be verified). Hence, \tilde{\Lambda}{}_{\kappa}(\text{\boldmath\mathcal{C}}) is also differentiable in that case, since \Lambda_{\kappa}(\bm{u},\bm{\sigma},\text{\boldmath\mathcal{C}}) is, upon expressing $\sigma$ with (13), quadratic in $(\bm{u},\bm{w})$ and affine in $\mathcal{C}$. 5. (e)

In turn, continuity of \tilde{\Lambda}{}_{\kappa}(\text{\boldmath\mathcal{C}}) is one of two key ingredients needed to establish existence of solutions to Problem  (21). The other key ingredient would be establishing the existence of minimizing sequences in $\mathcal{Q}$ that contain convergent subsequences; it is outside of the scope of this work. Interested readers are referred to [17] for an in-depth study of existence of solutions in problems involving functionals similar to those studied herein.

The above considerations show the importance of establishing the well-posednedess of the stationarity system (16,18); this is the goal of the next section.

3. Analysis of the stationarity problem

In preparation for the forthcoming analysis, we rewrite the stationarity problems  (18-16) in the form

[TABLE]

where $\mathcal{A}$ is the elastic stiffness bilinear form given for any $(\bm{v},\widetilde{\bm{v}}{})\in\mathcal{U}\times\mathcal{U}$ by \mathcal{A}(\bm{v},\widetilde{\bm{v}}{},\text{\boldmath\mathcal{C}}):=\big{(}\hskip 1.00006pt\bm{\varepsilon}[\bm{v}],\text{\boldmath\mathcal{C}}\!:\!\bm{\varepsilon}[\widetilde{\bm{v}}{}]\hskip 1.00006pt\big{)}, $\mathcal{B}:=\mathcal{A}-\omega^{2}(\rho\cdot,\cdot)$ is the dynamic stiffness bilinear form, the linear functional $\bm{d}\in\mathcal{U}^{\prime}$ defined by \big{\langle}\hskip 1.00006pt\bm{d},\widetilde{\bm{u}}{}\hskip 1.00006pt\big{\rangle}_{\mathcal{U}^{\prime},\mathcal{U}}=\mathcal{D}(\bm{u}^{\text{\scriptsize m}},\widetilde{\bm{u}}{}) incorporates the kinematic data, the linear functional $\bm{f}\in\mathcal{W}^{\prime}$ synthesizes all given applied loads, and the function spaces $\mathcal{U},\mathcal{W}$ are respectively defined by (6a) and (15). Moreover, we endow $\mathcal{U}$ and $\mathcal{W}$ with the inner product $(\cdot,\cdot)_{\omega}$ and norm $\|\cdot\|_{\omega}$ defined by

[TABLE]

The norm $\|\cdot\|_{\omega}$ is equivalent (for fixed $\omega>0$) to the standard $H^{1}$ norm (if $|\Gamma_{\text{D}}|>0$, the alternative definition $\|\widetilde{\bm{u}}{}\|^{2}=\mathcal{A}(\widetilde{\bm{u}}{},\widetilde{\bm{u}}{})$ is also suitable). The bilinear forms $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{D}$ are assumed to be continuous for this norm, i.e. there exist constants $a>0$, $b>0$ and $d>0$ such that

[TABLE]

for all $\bm{w},\widetilde{\bm{w}}{}\in\mathcal{W}$ and $\bm{u},\widetilde{\bm{u}}{}\in\mathcal{U}$. In what follows, we will denote by $A,B,D$ the linear operators implicitly defined by the bilinear forms $\mathcal{A},\mathcal{B},\mathcal{D}$, e.g.:

[TABLE]

Moreover, $B^{\text{t}}:\mathcal{W}\to\mathcal{U}^{\prime}$ will denote the transposed operator associated to $B$, defined by

[TABLE]

and $N(F)$ will denotes the null space of a linear operator $F$.

3.1. Well-posedness of the coupled problem

We regard problem (24) as a perturbed mixed problem [7, Sec. 4.3]. It can in fact be given the form

[TABLE]

where $\mathbb{U}$ is the Hilbert space $\mathbb{U}:=\mathcal{U}\times\mathcal{W}$ equipped with the inner product defined by $(\mathsf{U},\widetilde{\mathsf{U}})_{\mathbb{U}}=(\bm{u},\widetilde{\bm{u}}{})_{\omega}+(\bm{w},\widetilde{\bm{w}}{})_{\omega}$ (the associated norm being given by $\|\mathsf{U}\|^{2}_{\mathbb{U}}=(\mathsf{U},\mathsf{U})_{\mathbb{U}}$) and with the (continuous, symmetric) bilinear form $G:\mathbb{U}\times\mathbb{U}\to\mathbb{R}$ and the (continuous) linear form $\mathsf{F}\in\mathbb{U}^{\prime}$ defined by

[TABLE]

Remark 3.1**.**

The quadratic form $\mathsf{U}\mapsto G(\mathsf{U},\mathsf{U})$ with $\kappa\geq 0$ is not sign-definite: $G(\mathsf{U},\mathsf{U})=\mathcal{A}(\bm{w},\bm{w})\geq 0$ for $\mathsf{U}=(\bm{w},\mathbf{0})$, and $G(\mathsf{U},\mathsf{U})=-\mathcal{A}(\bm{w},\bm{w})-\kappa\mathcal{D}(\bm{u},\bm{u})\leq 0$ for any $(\bm{w},\bm{u})$ satisfying (24a) with $\bm{f}=\mathbf{0}$.

Let $\mathsf{G}:\mathbb{U}\to\mathbb{U}^{\prime}$ be the (bounded) linear operator associated to the bilinear form $G$. The above definition of $G$ implies $\mathsf{G}=\mathsf{G}^{\text{t}}$. The variational problem (24) can be shown to be well-posed by checking the applicability of the closed range theorem (see e.g. [11, Thm. 5.11-6]), which here has the form:

Lemma 3.2**.**

Let $\mathbb{U}$ be a Hilbert space and $\mathsf{G}:\mathbb{U}\to\mathbb{U}^{\prime}$ a bounded linear operator such that $\mathsf{G}=\mathsf{G}^{\text{t}}$. $\mathsf{G}$ is invertible with bounded inverse if and only if there is a constant $\eta>0$ such that $\|\mathsf{G}\mathsf{U}\|_{\mathbb{U}^{\prime}}\geq\eta\|\mathsf{U}\|_{\mathbb{U}}$ for any $\mathsf{U}\in\mathbb{U}$, in which case $\|\mathsf{G}^{-1}\|\leq\eta^{-1}$.

The null spaces $\mathcal{H}:=N(B)$ and $\mathcal{K}:=N(B^{\text{t}})$ of $B$ and $B^{\text{t}}$ play an important role for studying the well-posedness of problem (24), so we now characterize them. To this end, let $\mathcal{Z}$ denote the subspace of $\mathcal{U}$ such that any $\bm{u}\in\mathcal{Z}$ solves the homogeneous problem

[TABLE]

As is well-known, $\mathcal{Z}=\{\mathbf{0}\}$ unless $\omega$ belongs to the countable set of eigenvalues for problem (27), in which case $\mathcal{Z}$ is finite-dimensional.

An element $\bm{u}$ of $\mathcal{H}\subset\mathcal{U}$ is characterized by $\mathcal{B}(\bm{u},\widetilde{\bm{w}}{})=0$ for all $\widetilde{\bm{w}}{}\in\mathcal{W}$. On applying the first Green identity in $\mathcal{A}(\cdot,\cdot)$, the strong form of this variational problem is found to be

[TABLE]

in strong form. Then, two cases arise:

Case (i). If $\mathcal{W}=\mathcal{U}$ (i.e. $\Gamma_{\text{D}}\cup\Gamma_{\text{N}}=\Gamma$), problems (28a) and (28b) are identical and coincide with problem (27); therefore $\mathcal{H}=\mathcal{K}=\mathcal{Z}$.

Case (ii). $\mathcal{W}\subsetneq\mathcal{U}$ (i.e. $\Gamma_{c}\not=\emptyset$), the situation is completely different as the boundary conditions are undetermined in (28a) and overdetermined in (28b). In the latter case, homogeneous Dirichlet and Neumann data is simultaneously imposed on $\Gamma_{c}$, therefore problem (28b) has only the trivial solution by virtue of the unique continuation principle [2, Corollary], i.e. $\mathcal{K}=\{\mathbf{0}\}$. By contrast, $\mathcal{H}$ now includes forced responses for any excitation applied on $\Gamma_{c}$ and eigenfunctions when $\omega$ is an eigenvalue for any kind of homogeneous data on $\Gamma_{c}$ (see problem (28a)), and is thus infinite-dimensional. We do not attempt to characterize $\mathcal{H}$ more precisely, since in this case we will only use the fact that $\mathcal{K}=\{\mathbf{0}\}$.

The following property of the bilinear form $\mathcal{B}$, whose proof is given in Sec. 7.1, is crucial for establishing the well-posedness of the stationarity problem (24):

Lemma 3.3**.**

There exists $\beta>0$ such that the bilinear form $\mathcal{B}:\mathcal{U}\times\mathcal{W}\to\mathbb{R}$ introduced in problem (24) satisfies the inf-sup condition

[TABLE]

Moreover, let $0\leq\omega_{1}\leq\omega_{2}\leq\ldots$ be the eigenvalues associated with the homogeneous problem (27). Then, there exists $0<\xi\leq 1$ such that the inf-sup constant $\beta$ satisfies

[TABLE]

If either $\mathcal{U}=\mathcal{W}$ or $\mathcal{H}=\{\mathbf{0}\}$, the above inequality holds with $\xi=1$.

In addition, the following assumptions are made on the bilinear forms $\mathcal{D}$ and $\mathcal{A}$:

Assumption 1**.**

The measurement bilinear form $\mathcal{D}$ is coercive on $\mathcal{H}\times\mathcal{H}$: there exists $\delta>0$ such that $\mathcal{D}(\bm{v},\bm{v})\geq\delta\|\bm{v}\|^{2}_{\omega}$ for all $\bm{v}\in\mathcal{H}$.

Assumption 2**.**

The experimental conditions are such that $|\Gamma_{\text{D}}\cup\Gamma_{c}|>0$. The stiffness bilinear form $\mathcal{A}$ is therefore coercive on $\mathcal{W}\times\mathcal{W}$.

Remark 3.4**.**

Let $\mathcal{N}:=N(D)\subset\mathcal{U}$ be the null space of $D$. Assumption 1 implies that $\mathcal{H}\cap\mathcal{N}=\{\mathbf{0}\}$. If it is not verified, there exists $\bm{z}\not=\mathbf{0},\ \bm{z}\in\mathcal{H}\cap\mathcal{N}$. Then, $(\bm{w},\bm{u})=(\mathbf{0},\bm{z})$ solves (24) with $\bm{f}=\bm{d}=\mathbf{0}$, i.e. uniqueness fails for the original stationarity problem. Assumption 1 is therefore necessary for the well-posedness of problem (24), and will be seen to be also sufficient. It means that any nontrivial elastodynamic state satisfying the (possibly incomplete) homogeneous boundary conditions on $\Gamma_{\text{D}},\Gamma_{\text{N}}$ must register on the measurement apparatus.

For the case of incomplete BCs (for which $\mathcal{H}$ is infinite-dimensional), assumption 1 is a stringent requirement, as it cannot be met with a $L^{2}$ norm on measurement residuals (since $D$ would then be compact for the $\|\cdot\|_{\omega}$ norm). In this case, $\mathcal{D}$ may for example be defined so as to be equivalent to the $H^{1}(\Omega_{\text{\scriptsize m}})$ norm. By contrast, for the complete BC case (for which $\mathcal{H}$ is at most finite-dimensional), coercivity on $\mathcal{H}\times\mathcal{H}$ can be achieved with $\mathcal{D}$ defined in terms of a $L^{2}$ norm.

Assumption 2 excludes the case $\Gamma_{\text{N}}=\partial\Omega$, i.e. pure Neumann boundary conditions. The resulting coercivity of $\mathcal{A}$ over the whole $\mathcal{W}\times\mathcal{W}$ will make the forthcoming analyses simpler and is not detrimental in practice: under usual experimental conditions, the sample will not undergo known surface forces on the whole boundary $\partial\Omega$ while being completely “unsupported.”

We are now in a position to establish the well-posedness of the stationarity problem (24). To do so, still following [7], we use the decomposition $\mathcal{U}=\mathcal{H}\oplus\mathcal{H}^{\perp}$ and split any $\bm{u}\in\mathcal{U}$ and $\bm{d}\in\mathcal{U}^{\prime}$ according to

[TABLE]

to cater for $\mathcal{H}$ being potentially non-trivial, noting that $\mathcal{B}(\bm{u},\widetilde{\bm{w}}{})=\mathcal{B}(\bm{u}_{1},\widetilde{\bm{w}}{})$. We therefore have

[TABLE]

Using these definitions and notations, we obtain the following main result, whose proof is given in Section 7.2:

Theorem 3.5**.**

Let Assumptions 1 and 2 hold. Then, for every $\bm{f}\in\mathcal{W}^{\prime}$ and $\bm{d}\in\mathcal{U}^{\prime}$ and for any $\kappa>0$, problem (24) has a unique solution $(\bm{w},\bm{u})\in\mathcal{W}\times\mathcal{U}$, that moreover satisfies

[TABLE]

where the constant $C$ depends only on $\kappa$, the stability constants $\alpha,\beta,\delta$ and the continuity constants $a,\,d$. More precisely, with reference to the form (36a,b) of problem (24), the following estimates hold (with all dependences on $\kappa$ in the continuity constants made explicit):

[TABLE]

where the non-dimensional constants $q_{a},q_{d},r_{a},r_{d}$ are defined by $q_{a}:=a/\alpha$, $q_{d}:=d/\delta$, $r_{a}:=a/\beta$ and $r_{d}:=d/\beta$ and the constant $Q$ is given by 2Q:=\kappa\alpha^{-1}\big{(}\hskip 1.00006ptq_{d}r_{a}+\sqrt{4\alpha(\kappa\delta)^{-1}+q^{2}_{c}r^{2}_{a}}\hskip 1.00006pt\big{)}.

Remark 3.6**.**

Consider finite-dimensional subspaces $\mathcal{W}_{h}\subset\mathcal{W}$ and $\mathcal{U}_{h}\subset\mathcal{W}$ (e.g. from finite element discretizations), whose dimension depends on a discretization parameter $h$. If the conditions of Theorem 3.5 hold uniformly with $h$ (i.e. for all discretization levels), the discrete systems arising from (24) are well-posed. Moreover, if the approximability property holds (i.e. elements of $\mathcal{U},\mathcal{W}$ can be approximated arbitrarily closely by elements of $\mathcal{U}_{h},\mathcal{W}_{h}$ for some small enough $h$), then the sequence of solutions of the discrete systems converges (as $h\to 0$) to the solution of (24) in the given norm [7].

3.2. Supplementary assumption for identification feasibility

Assumptions 1 and 2 ensure the well-posedness of the stationarity problem (24). Additional requirements are however needed to ensure that the data $\bm{u}^{\text{\scriptsize m}}$ also provides useful information towards the original elastic imaging inverse problem. To see this, consider the case where $N(D)=:\mathcal{N}=\mathcal{H}^{\perp}$, for which the experimental data is just sufficient to satisfy Assumption 1. In that case, the stationarity problem (24) reads

[TABLE]

having written equation (24b) separately for $\widetilde{\bm{u}}{}=\widetilde{\bm{u}}{}_{1}\in\mathcal{H}^{\perp}$ and $\widetilde{\bm{u}}{}=\widetilde{\bm{u}}{}_{0}\in\mathcal{H}$. Equation (c) determines $\bm{u}_{0}$ solely from the measurements (so $\bm{u}_{0}$ depends neither on $\bm{f}$ nor on the assumed elastic properties $\mathcal{C}$), while equations (a), (b) show that $\bm{w},\bm{u}_{1}$ do not depend on the measurements. The case $\mathcal{N}=\mathcal{H}^{\perp}$ is therefore a limiting situation where the measurement $\bm{u}^{\text{\scriptsize m}}$ carries no information on $\mathcal{C}$. We therefore introduce the following additional assumption on the measurement configuration:

Assumption 3**.**

$\mathcal{N}$* is a proper subspace of $\mathcal{H}^{\perp}$, i.e. has a nontrivial orthogonal complement in $\mathcal{H}^{\perp}$.*

4. Stationarity solution asymptotics

To gain insight into the effect of the adjustable weight parameter $\kappa$ on the minimization problem (9), we now seek the limiting forms of the solution of the stationarity problem (24) in the limiting situations $\kappa\to 0$ and $\kappa\to+\infty$. Consequences of the results of this section, in particular regarding sign properties of the Hessian of the MECE reduced functional, are discussed in Sections 5.4 and 5.5.

We start by recasting problem (24) using the splitting (30) and operator notation, to obtain

[TABLE]

The above block form (36a,b) of the stationarity problem (24) is such that its first row is equation (24a) whereas the remaining two rows are equation (24b) with $\widetilde{\bm{u}}{}=\widetilde{\bm{u}}{}_{0}\in\mathcal{H}$ and $\widetilde{\bm{u}}{}=\widetilde{\bm{u}}{}_{1}\in\mathcal{H}^{\perp}$, in that order. The operators $D_{00},D_{10},D_{01},D_{11}$ are the $\mathcal{H}\to\mathcal{H}^{\prime}$, $\mathcal{H}^{\perp}\to\mathcal{H}^{\prime}$, $\mathcal{H}\to(\mathcal{H}^{\perp})^{\prime}$ and $\mathcal{H}^{\perp}\to(\mathcal{H}^{\perp})^{\prime}$ restrictions of $D:\mathcal{U}\to\mathcal{U}^{\prime}$, respectively (these restrictions being defined by $D_{ij}=P_{j}DE_{i}$ in terms of the extension operators $E_{0}:\mathcal{H}\to\mathcal{U}$, $E_{1}:\mathcal{H}^{\perp}\to\mathcal{U}$ and the orthogonal projectors $P_{0}:\mathcal{U}\to\mathcal{H}$, $P_{1}:\mathcal{U}\to\mathcal{H}^{\perp}$). The zero blocks in (36b) account for $\mathcal{H}$ being the null space of the operator $B$.

4.1. Small-$\kappa$ expansion

We begin by deriving the leading expansion of the stationarity solution $(\bm{w},\bm{u})=(\bm{w}^{\kappa},\bm{u}^{\kappa})$ about $\kappa=0$, assuming a priori the expansion to have the form

[TABLE]

(i.e. to follow the format $\mathsf{U}_{\kappa}=\mathsf{U}^{(0)}+\kappa\mathsf{U}^{(1)}+\ldots$), inserting the ansatz (37) into problem (36a), (36b) and writing the resulting $O(1),\,O(\kappa)\ldots$ equations. The $O(1)$ equations are readily obtained as:

[TABLE]

Proposition 1**.**

The stationarity solution admits the small-$\kappa$ expansion $\mathsf{U}=\mathsf{U}^{(0)}+O(\kappa)$, in the sense of the $\|\cdot\|_{\mathbb{U}}$ norm. The components $\bm{w},\bm{u}_{0},\bm{u}_{1}$ of $\mathsf{U}^{(0)}$ solve the well-posed problems (38a) and (38b). Moreover, we have $\bm{w}^{(0)}=\mathbf{0}$ in the “usual” case where $\mathcal{K}=\{\mathbf{0}\}$.

Proof 4.1**.**

Define the truncation error $\Delta\mathsf{U}=\{\bm{w}-\bm{w}^{(0)},\,\bm{u}_{0}-\bm{u}_{0}^{(0)},\,\bm{u}_{1}-\bm{u}_{1}^{(0)}\}=\mathsf{U}^{\kappa}-\mathsf{U}^{(0)}$. On setting $\mathsf{U}^{\kappa}=\mathsf{U}^{(0)}+\Delta\mathsf{U}$ in (36a) and using equations (38a)-(38c), the governing system for $\Delta\mathsf{U}$ is

[TABLE]

and therefore has the form (36a,b) with $\bm{f}=\bm{d}_{0}=\mathbf{0}$ and $\bm{d}_{1}=-B^{\text{t}}\bm{w}^{(1)}$. Theorem 3.5 hence applies to problem (39), yielding the bounds

[TABLE]

Consequently there exists $\kappa_{0}>0$ and a constant $C(\kappa_{0})$ such that we have $\|\Delta\mathsf{U}\|_{\mathbb{U}}\leq C(\kappa_{0})\kappa$ for all $\kappa<\kappa_{0}$. Finally, if $\mathcal{K}=\{\mathbf{0}\}$, the second equation of (38a) implies $\bm{w}^{(0)}=\mathbf{0}$. The proof is complete.

4.2. Large-$\kappa$ expansion

We now seek an expansion of $\mathsf{U}_{\kappa}$ about $\kappa=+\infty$, assumed to have the form $\mathsf{U}_{\kappa}=\mathsf{U}^{(0)}+\kappa^{-1}\mathsf{U}^{(1)}+\ldots$. Using this ansatz into the block form (36a,b) of problem (24), the arising leading-order $O(\kappa)$ equation is

[TABLE]

The chosen ansatz requires that this equation be well-posed. Assumption 1 on $\mathcal{D}$, which postulates only the coercivity of $D_{00}$, is insufficient in this respect and must be replaced by a stronger requirement:

Assumption 4**.**

$\mathcal{D}$* is coercive on $\mathcal{N}^{\perp}\times\mathcal{N}^{\perp}$ (with $\mathcal{N}=N(D)$): there exists $\delta>0$ such that $\mathcal{D}(\bm{v},\bm{v})\geq\delta\|\bm{v}\|^{2}_{\omega}$ for all $\bm{v}\in\mathcal{N}^{\perp}$.*

This also makes the splitting (30) unsuitable for studying the large-$\kappa$ limiting case. We instead split $\bm{u}$ according to

[TABLE]

which implies $\bm{d}_{1}=\mathbf{0}$. The operator form (36a) of the coupled problem now uses the definitions

[TABLE]

based on the splitting (41), instead of (36a). Its solution satisfies (as shown in Sec. 7.3) the estimates

[TABLE]

with $s_{a}:=b/\alpha$ and $s_{d}:=b/\delta$. Problem (36a), (42) therefore has a unique solution $\mathsf{U}$ that depends continuously on the data $\bm{f},\bm{d}_{0}$ (as already known from Theorem 3.5). Moreover, and importantly, all continuity constants in estimates (43) are bounded in the limit $\kappa\to\infty$. Choosing $\kappa_{0}>0$, we obtain the existence of a constant $C(\kappa_{0})>0$, independent on $\kappa$, such that $\|\mathsf{U}\|_{\mathbb{U}}\leq C(\kappa_{0})$ for all $\kappa\geq\kappa_{0}$.

We now proceed by assuming the expansion of the stationarity solution $(\bm{w},\bm{u})=(\bm{w}^{\kappa},\bm{u}^{\kappa})$ about $\kappa=+\infty$ to have, in terms of the new splitting (41), the form

[TABLE]

Inserting this ansatz into (42) results in $O(\kappa),\,O(1)\ldots$ equations. The $O(\kappa)$ equation is

[TABLE]

it is satisfied by $\bm{u}^{(0)}_{0}=\bm{u}^{\text{\scriptsize m}}$ since the definition of $\bm{d}$ in (24) and decomposition (41) imply that $\bm{d}_{0}=D\bm{u}^{\text{\scriptsize m}}$. Then, the $O(1)$ equations yield the two uncoupled well-posed systems

[TABLE]

Proposition 2**.**

The stationarity solution admits the large-$\kappa$ expansion $\mathsf{U}=\mathsf{U}^{(0)}+O(\kappa^{-1})$, in the sense of the $\|\cdot\|_{\mathbb{U}}$ norm, with the components $\bm{w}^{(0)},\bm{u}^{(0)}_{0},\bm{u}^{(0)}_{1}$ of $\mathsf{U}^{(0)}$ solving the well-posed problems (46a,b).

Proof 4.2**.**

We examine the truncation error $\Delta\mathsf{U}$, defined as for Prop. 1. The governing system for $\Delta\mathsf{U}$ still has the format (39), its right-hand side $\mathsf{Y}$ being now given, from using (45) and (46a,b), by

[TABLE]

The problem (36a), (47) for $\Delta\mathsf{U}$ has the form (36a), (42), with $\bm{f}=\mathbf{0}$ and $\bm{d}_{0}=\kappa^{-1}B^{\text{t}}_{0}\bm{w}^{(0)}$. Its solution therefore satisfies estimates (43), which for this particular right-hand side give

[TABLE]

Consequently, there exists $C(\kappa_{0})>0$ such that $\|\Delta\mathsf{U}\|_{\mathbb{U}}\leq\kappa^{-1}C(\kappa_{0})$ for all $\kappa\geq\kappa_{0}$.

5. Derivatives of the reduced MECE functional

In this section, we derive general formulas for the gradient and Hessian of \tilde{\Lambda}{}_{\kappa}(\text{\boldmath\mathcal{C}}) at any \text{\boldmath\mathcal{C}}\in\mathcal{Q}. Then, taking advantage of the results of Sec. 4, we show that the Hessian of \tilde{\Lambda}{}_{\kappa}(\text{\boldmath\mathcal{C}}) is asymptotically positive in the large-$\kappa$ case but sign-indefinite in the small-$\kappa$ case connected to least-squares mininimization.

5.1. First-order derivative of $\tilde{\Lambda}{}_{\kappa}$

Recall, from remark (d) of Section 2.5, that the reduced objective $\tilde{\Lambda}{}_{\kappa}$ is continuously differentiable with respect to $\mathcal{C}$. The first-order derivative is a priori given by

[TABLE]

where $D_{\mathcal{C}}$ denotes a total derivative w.r.t. $\mathcal{C}$ and the prime $(\cdot)^{\prime}$ symbol is used as a shorthand notation for a derivative w.r.t. $\mathcal{C}$ in a given direction \hat{\text{\boldmath\mathcal{C}}} (e.g. \tilde{\Lambda}{}_{\kappa}^{\prime}(\text{\boldmath\mathcal{C}})=\tilde{\Lambda}{}_{\kappa}^{\prime}(\text{\boldmath\mathcal{C}})[\hat{\text{\boldmath\mathcal{C}}}]). The last two equalities follow from definition (21) of $\tilde{\Lambda}{}_{\kappa}$ implying verification of the stationarity equations (11a-c), and in particular of the balance constraint, for any $\mathcal{C}$. An explicit expression of \tilde{\Lambda}{}_{\kappa}^{\prime}(\text{\boldmath\mathcal{C}}) is then obtained as

[TABLE]

In particular, any \text{\boldmath\mathcal{C}}^{\star} solving the reduced minimization problem (21) must satisfy the first-order optimality condition

[TABLE]

5.2. First-order derivative of the stationarity solution

The second-order derivative of $\tilde{\Lambda}{}_{\kappa}$ will involve (from applying the total derivative operator $D_{\mathcal{C}}$ to (50)) the first-order derivative $(\bm{u}^{\prime},\bm{w}^{\prime})$ of the stationarity solution $(\bm{u},\bm{w})$ (see remark (d) of Section 2.5). Upon differentiating w.r.t. $\mathcal{C}$ the stationarity problem (24), the stationarity solution derivative $(\bm{u}^{\prime},\bm{w}^{\prime})\in\mathcal{U}\times\mathcal{W}$ solves the variational problem

[TABLE]

which is well posed since its governing operator is identical to that of the stationarity problem (24).

5.3. Second-order derivative of $\tilde{\Lambda}{}_{\kappa}$

For convenience, we focus in the sequel on the $\mathcal{Q}\to\mathbb{R}$ quadratic form defined by \hat{\text{\boldmath\mathcal{C}}}\mapsto\tilde{\Lambda}{}_{\kappa}^{\prime\prime}(\text{\boldmath\mathcal{C}})[\hat{\text{\boldmath\mathcal{C}}},\hat{\text{\boldmath\mathcal{C}}}], denoted \tilde{\Lambda}{}_{\kappa}^{\prime\prime}(\text{\boldmath\mathcal{C}}) for short, the corresponding $\mathcal{Q}\times\mathcal{Q}\to\mathbb{R}$ bilinear mapping being then given by the usual polarization identity 4\tilde{\Lambda}{}_{\kappa}^{\prime\prime}(\text{\boldmath\mathcal{C}})[\hat{\text{\boldmath\mathcal{C}}}_{1},\hat{\text{\boldmath\mathcal{C}}}_{2}]=\tilde{\Lambda}{}_{\kappa}^{\prime\prime}(\text{\boldmath\mathcal{C}})[\hat{\text{\boldmath\mathcal{C}}}_{1}+\hat{\text{\boldmath\mathcal{C}}}_{2},\hat{\text{\boldmath\mathcal{C}}}_{1}+\hat{\text{\boldmath\mathcal{C}}}_{2}]-\tilde{\Lambda}{}_{\kappa}^{\prime\prime}(\text{\boldmath\mathcal{C}})[\hat{\text{\boldmath\mathcal{C}}}_{1}-\hat{\text{\boldmath\mathcal{C}}}_{2},\hat{\text{\boldmath\mathcal{C}}}_{1}-\hat{\text{\boldmath\mathcal{C}}}_{2}]. The second-order derivative \tilde{\Lambda}{}_{\kappa}^{\prime\prime}(\text{\boldmath\mathcal{C}}) can be derived by differentiating (50) (e.g. [11, Thm. 7.8-2]) and rearranging terms, to obtain

[TABLE]

This expression of \tilde{\Lambda}{}_{\kappa}^{\prime\prime}(\text{\boldmath\mathcal{C}}) can be recast in a convenient alternative form as follows. Setting $\widetilde{\bm{w}}{}=\bm{w}^{\prime}$ in (51a) and $\widetilde{\bm{u}}{}=\bm{u}^{\prime}$ in (51b), we find the identity

[TABLE]

so that (52) takes the symmetric form

[TABLE]

Moreover, \tilde{\Lambda}{}_{\kappa}^{\prime\prime}(\text{\boldmath\mathcal{C}}) can alternatively be given a sign-revealing form. To this aim, invoking again the splitting (41), the derivative problem (51) can be written, in operator form, as

[TABLE]

where e.g. $A^{\prime}$ is the operator associated with \mathcal{A}(\cdot,\cdot,\hat{\text{\boldmath\mathcal{C}}}). Using this in (54) produces the following result (whose proof is given in Sec. 7.4), which gives \tilde{\Lambda}{}_{\kappa}^{\prime\prime}(\text{\boldmath\mathcal{C}}) as an algebraic sum of positive quadratic forms:

Proposition 3**.**

Let $(\bm{w},\bm{u}_{0},\bm{u}_{1})$ solve the stationarity problem (36a), (42) and $(\bm{w}^{\prime},\bm{u}^{\prime}_{0},\bm{u}^{\prime}_{1})$ solve the derivative problem (55). The second-order derivative of the reduced MECE functional is given by

[TABLE]

with the (symmetric, positive, coercive) operator $Z$ defined by $Z:=A+\kappa^{-1}B_{0}D^{-1}B^{\text{t}}_{0}$ in terms of the operators appearing in problem (55).

5.4. Reduced MECE functional: large-$\kappa$ limiting case

Regarding the leading-order term $(\bm{u}^{(0)},\bm{w}^{(0)})$ of the stationarity solution expansion (44), we observe that equation (46a) for $(\bm{u}_{1}^{(0)},\bm{w}^{(0)})$ coincides with the stationarity problem for the minimization of the pure ECE functional (7) with the kinematic measurement enforced strictly, the latter being achieved by the solution $\bm{u}^{(0)}_{0}$ of the remaining equation (45). Moreover, using expansion (44) in \tilde{\mathcal{E}}(\text{\boldmath\mathcal{C}}),\tilde{\mathcal{D}}(\text{\boldmath\mathcal{C}}) defined by (21), we find

[TABLE]

Consistently with the fact that the data is enforced exactly in the $\kappa\to\infty$ limit, we therefore observe that the value of the reduced objective \tilde{\Lambda}{}_{\kappa}(\text{\boldmath\mathcal{C}}) is for large $\kappa$ dominated by its ECE component.

Moreover, exploiting the large-$\kappa$ expansion of the stationarity solution in the reduced objective \tilde{\Lambda}{}_{\kappa}^{\prime\prime}(\text{\boldmath\mathcal{C}}) as given by Proposition 3 for situations where elastic moduli are kept fixed outside the measurement region, \tilde{\Lambda}{}_{\kappa}^{\prime\prime}(\text{\boldmath\mathcal{C}}) is found to be asymptotically convex in the large-$\kappa$ limit:

Theorem 5.1**.**

Assume that \text{supp}(\hat{\text{\boldmath\mathcal{C}}})\subset\Omega^{\text{\scriptsize m}} (i.e. that the perturbation \hat{\text{\boldmath\mathcal{C}}} vanishes outside the measurement region), so that $B_{1}^{\prime}=0$. Then, the second-order derivative \tilde{\Lambda}{}_{\kappa}^{\prime\prime}(\text{\boldmath\mathcal{C}}) admits the large-$\kappa$ expansion

[TABLE]

with its leading term \tilde{\Lambda}{}_{\kappa}^{\prime\prime}{}^{(0)}(\text{\boldmath\mathcal{C}}) given by

[TABLE]

In the above expression, the operator $P_{0}:\mathcal{W}^{\prime}\to\mathcal{W}^{\prime}$ is defined by $P_{0}:=B_{1}(B_{1}^{\text{t}}A^{-1}B_{1}^{-1})^{-1}B_{1}^{\text{t}}A^{-1}$ (so $P_{0}P_{0}=P_{0}$, i.e. $P_{0}$ is a projector, and $A^{-1}P_{0}=P_{0}^{\text{t}}A^{-1}$), while $\bm{w}^{(0)}$ is the (leading) zeroth-order term of the large-$\kappa$ expansion of $\bm{w}$, see Prop. 2. The above leading-order coefficient \tilde{\Lambda}{}_{\kappa}^{\prime\prime}{}^{(0)}(\text{\boldmath\mathcal{C}}) is therefore positive. Moreover, the subsequent coefficient \tilde{\Lambda}{}_{\kappa}^{\prime\prime}{}^{(1)}(\text{\boldmath\mathcal{C}}) is sign-indefinite.

Proof 5.2**.**

See Section 7.5.

5.5. Reduced MECE functional: small-$\kappa$ limiting case

We restrict this discussion to the “usual” case where $\mathcal{K}=\{\mathbf{0}\}$, for which we have $\bm{w}^{(0)}=\mathbf{0}$ (see Prop. 1). Using this and expansion (37) of the stationarity solution in (21) and (50), the reduced MECE functional \tilde{\Lambda}{}_{\kappa}(\text{\boldmath\mathcal{C}}) and its derivative \tilde{\Lambda}{}_{\kappa}^{\prime}(\text{\boldmath\mathcal{C}}) have the $O(\kappa)$ expansions

[TABLE]

Consistently with the constitutive equation is enforced exactly in the $\kappa\to 0$ limit, we see that the value (60a) of the reduced objective \tilde{\Lambda}{}_{\kappa}(\text{\boldmath\mathcal{C}}) is for small $\kappa$ dominated by its data-misfit component $\mathcal{D}$.

Moreover, equations (38a,c) give $B\bm{u}^{(0)}_{1}=\bm{f}$ and $B^{\text{t}}\bm{w}^{(1)}=-\bm{d}_{1}+D_{01}\bm{u}^{(0)}_{0}+D_{11}\bm{u}^{(0)}_{1}$. If in addition $\mathcal{H}=\{\mathbf{0}\}$ (i.e. boundary conditions are well-posed and $\omega$ is not an eigenvalue for problems (27) and (28a)), we have that (i) $\bm{u}_{0}^{(0)}=\mathbf{0}$, (ii) $\bm{u}_{1}^{(0)}$ is the forward solution, (iii) $\bm{w}^{(1)}$ is the adjoint solution for the objective function \tilde{\mathcal{D}}(\text{\boldmath\mathcal{C}}), (iv) the leading contribution to \kappa^{-1}\tilde{\Lambda}{}_{\kappa}(\text{\boldmath\mathcal{C}}) as $\kappa\to 0$ is the reduced quadratic misfit functional \tilde{\mathcal{D}}(\text{\boldmath\mathcal{C}}):=\mathcal{D}(\bm{u}^{(0)}-\bm{u}^{\text{\scriptsize m}},\bm{u}^{(0)}-\bm{u}^{\text{\scriptsize m}}) commonly used for solving PDE-constrained inverse problems, and (v) the leading term in equation (60b) coincides with the known expression for \tilde{\mathcal{D}}^{\prime}(\text{\boldmath\mathcal{C}}). Consequently, the MECE-based inversion in reduced form becomes the minimization of the non-regularized least-squares misfit $\mathcal{D}$ in the limit $\kappa\to 0$. Similar remarks apply when $\mathcal{H}\not=\{\mathbf{0}\}$, the (now nontrivial) field $\bm{u}^{(0)}_{0}=D_{00}^{-1}(\bm{d}_{0}-D_{10}\bm{u}^{(1)})$ being such that the leading-order field $\bm{u}^{(0)}_{0}+\bm{u}^{(0)}_{1}$ minimizes $\mathcal{D}(\bm{u}^{(0)}-\bm{u}^{\text{\scriptsize m}},\bm{u}^{(0)}-\bm{u}^{\text{\scriptsize m}})$ for given $\mathcal{C}$.

Then, the stationarity solution expansion (37) yields expansions $\bm{w}^{\prime}=\kappa\bm{w}^{\prime\,(1)}+o(\kappa)$ and $\bm{u}^{\prime}=\bm{u}^{\prime\,(0)}+\kappa\bm{u}^{\prime\,(1)}+o(\kappa)$ for the solution derivatives (since $\bm{w}^{(0)}=\mathbf{0}$, which implies $\bm{w}^{\prime\,(0)}=\mathbf{0}$). Consequently, (54) provides

[TABLE]

whose leading term is a priori sign-indefinite (since it can be recast as an algebraic sum of positive quadratic expressions), by contrast with the corresponding result of Theorem 5.1 for the large-$\kappa$ case.

6. Numerical results

In this section, we present numerical studies that support our theoretical findings. We first show (Sec. 6.1) that the inf-sup constant remains strictly positive over a wide range of frequencies and displays convergent behavior upon mesh discretization. We then show (Sec. 6.2) that finite element discretizations are convergent as long as the two key necessary conditions are met. The next example (Sec. 6.3) illustrates how the reduced objective becomes convex as $\kappa$ increases. Finally, we show in Sec. 6.4 a parameter reconstruction example using interior data, which demonstrates the capability of the method to accommodate underspecified boundary conditions.

6.1. Stability of coupled system: 1D example

In this example, we study the behavior of the inf-sup constant $\beta$ defined in (29) for the operator $\mathcal{B}$ corresponding to longitudinal vibrations at frequency $f=\omega/2\pi$ of a one-dimensional bar fixed at one end, and whose length, mass density, and Young’s modulus are taken as one. The boundary condition at the other end is unspecified. The relevant spaces are \mathcal{W}:=\big{\{}\hskip 1.00006ptw\in H^{1}(\Omega),\;w(0)=w(1)=0\hskip 1.00006pt\big{\}} and \mathcal{U}:=\big{\{}\hskip 1.00006ptu\in H^{1}(\Omega),\;u(0)=0\hskip 1.00006pt\big{\}}, and we have $\mathcal{W}\subsetneq\mathcal{U}$. The bar is discretized using linear finite elements.

The inf-sup constant $\beta$ was computed for frequencies $f$ in the $[1,20]$ Hz range, in increments of $0.1$ Hz, the finite element mesh was refined until there was negligible change in $\beta$ over the latter frequency range. The numerical evaluation of $\beta$ uses a discretized version of (65). The requisite approximation $\widehat{\bm{u}}_{h}$ of the Riesz representative $\widehat{\bm{u}}$ of $B^{\text{t}}\bm{w}\in\mathcal{U}^{\prime}$ is obtained by setting $\widehat{\bm{u}}_{h}=\mathbf{N}_{h}\bm{d}$ and $\bm{w}_{h}=\mathbf{N}_{h}\bm{v}$, where $\mathbf{N}_{h}$ is a matrix of finite element shape functions for a given discretization with characteristic mesh size $h$. We then take $\mathbf{S}_{h}$ to be the positive-definite matrix such as \bm{d}^{T}\mathbf{S}_{h}\bm{d}:=\big{(}\hskip 1.00006pt\widehat{\bm{u}}_{h},\widehat{\bm{u}}_{h}\hskip 1.00006pt\big{)}_{\omega} and $\mathbf{B}_{h}$ the matrix associated with the bilinear form $\mathcal{B}$ after discretization with finite elements. Then, the nodal values $\bm{d}$ associated with $\widehat{\bm{u}}_{h}$ can be computed as $\bm{d}=\mathbf{S}_{h}^{-1}\mathbf{B}_{h}\bm{v}$. The inf-sup constant is obtained from

[TABLE]

(the first inequality following from the fact that $\mathcal{K}_{h}^{\perp}\subset\mathcal{W}_{h}$), i.e. by finding the smallest eigenvalue of the generalized eigenvalue problem $\mathbf{B}_{h}^{T}\mathbf{S}_{h}^{-T}\mathbf{B}_{h}\bm{v}=\lambda\mathbf{S}_{h}\bm{v}$. If that eigenvalue is positive, so is $\beta_{h}$.

Figure 1(a) shows that indeed $\beta_{h}$ remains positive for all frequencies investigated in this example. Thus, if the measured data satisfies Assumption 1, the coupled system corresponding to this 1D bar is well-posed in the considered frequency range, in spite of the boundary conditions being underspecified, as expected per Theorem 3.5. Figure 1(b) then shows how $\beta_{h}$ changes with element size $h$ for a frequency of $20$ Hz, and clearly demonstrates the expected convergence of $\beta$ as the mesh is refined.

6.2. Well-posedness and convergence of coupled system

In this example, we demonstrate the implications of Assumption 1 (i.e. $\mathcal{H}\cap\mathcal{N}=\{\mathbf{0}\}$), which is required for problem (24) to be well-posed, by showing its incidence on the convergence behavior of finite element approximations. We consider again a 1D bar, this time with a Dirichlet condition at both ends, so that we now have \mathcal{W}=\mathcal{U}=\big{\{}\hskip 1.00006ptu\in H^{1}(\Omega),\;u(0)=u(1)=0\hskip 1.00006pt\big{\}}. We use an assumed solution $\bar{u}=\sin(\omega x)$ and $\bar{w}=\bar{u}$, which corresponds to an excitation $f=-\omega^{2}\sin(\omega x)$ in problem (24). The observation operator $D$ (as introduced in (36b)) is of the form $D\bm{u}=\sum_{k=1}^{M}u(x_{k})$, i.e. involves pointwise measurements at $M$ locations to be specified (pointwise values of elements of $\mathcal{U}$ being well-defined in this 1D setting).

We study in this example whether finite element discretizations of problem (24) converge and if so at what rate, given that the exact assumed solution $\bar{u}$ is in $\mathcal{H}$. We define (piecewise-linear, conforming) finite element approximations $u_{h}$ and $w_{h}$, with $f_{h}$ being an interpolant of $f$, and consider two versions $D_{1}$ and $D_{2}$ of $D$, such that the $M$ measurement locations $x_{k}$ are taken as $M_{1}:=\{\frac{\pi}{\omega},\frac{2\pi}{\omega},\frac{3\pi}{\kappa},...,1\}$ for $D=D_{1}$ while $D_{2}$ uses randomly generated measurement abscissae. Consequently, we have $\bar{u}\in\mathcal{N}$ for $D=D_{1}$ but $\bar{u}\notin\mathcal{N}$ for $D=D_{2}$. The case $D=D_{1}$ thus makes $\mathcal{N}\cap\mathcal{H}$ non-trivial, in conflict with Assumption 1, whereas $\mathcal{N}\cap\mathcal{H}=\{0\}$ if $D_{2}$ is used. We thus expect the coupled system (24) to be ill-posed in the former case, but well-posed in the latter case.

Figure 2 shows the relative error $e_{ls}$ between the manufactured solution $(\bar{u},\bar{w})$ of the coupled system (24) and its finite element approximation $(u_{h},w_{h})$, defined as

[TABLE]

as a function of the element size $h$ and for either measurement set. The expected convergence of $(u_{h},w_{h})$ to $(\bar{u},\bar{w})$ is indeed observed, with the expected $O(h^{2})$ rate, if $D=D_{2}$. By contrast, using $D=D_{1}$ leads to rapid divergence of $e_{ls}$ as the mesh is refined.

6.3. Convexity of the reduced MECE functional

We now show an example that demonstrates the convexification of the reduced objective (21) as $\kappa$ increases, predicted by part 2 of Theorem 5.1. We consider again a 1D bar problem, this time of the form

[TABLE]

with the inhomogeneous Young modulus taken as $E(x)=1_{]0,0.5[}(x)\,E_{1}+1_{]0.5,1[}(x)\,E_{2}$ and the excitation frequency set to $f=3$ Hz. The bar is discretized using $100$ linear finite elements. Displacements are assumed to be measured at all nodes. The reduced objective $\tilde{\Lambda}{}_{\kappa}(E_{1},E_{2})$ was computed for $(E_{1},E_{2})\in[0.5,2.5]^{2}$ and several values of $\kappa$, by solving the coupled system (24) for each combination $(E_{1},E_{2})$ and then computing (21).

Figures 3(a)-3(d) illustrate how the reduced objective $(E_{1},E_{2})\mapsto\tilde{\Lambda}{}_{\kappa}(E_{1},E_{2})$ changes with $\kappa$, and in particular show clearly that this function is progressively smoother as $\kappa$ increases and becomes convex for large $\kappa$ (Figure 3(d)). Furthermore, for very small values of $\kappa$ (Figure 3(a)) the objective resembles that of a least squares functional, as predicted by the asymptotic analysis of Section 5.5.

These results have strong practical implications. The least-squares objective function (small-$\kappa$ limiting case) has many local minima, making the inversion results strongly dependent on the initial guess. On the other hand, the MECE objective for intermediate and large $\kappa$ is much smoother and convex, which translates into robustness with respect to the initial guess.

It is important to mention that the value of $\kappa$ used for a given problem has to be set according to the level of noise in the data. Hence, although our analysis indicates that choosing $\kappa$ as large as possible is beneficial for convexity, doing so without regard to noise level would likely produce poor reconstructions. This results from the fact that, as $\kappa$ increases, $\bm{u}$ becomes closer to the measured data $\bm{u}^{\text{\scriptsize m}}$ and may over-fit the noise. See [27] for some discussion on how to choose $\kappa$ according to the Morozov discrepancy principle and a heuristic approach termed error balance.

6.4. 2D identification example

In this last example, we demonstrate that we can successfully reconstruct material properties with MECE when boundary conditions are completely unknown but full field measurements are available in part of the domain. We consider a square elastic domain with an unknown circular inclusion, under 2D plane-strain conditions; see Figure  4(a). The synthetic experiment consists in loading the body with a uniform time-harmonic pressure of $1$ kPa applied on the top side at a frequency of $10$ Hz, the bottom side being kept fixed, and with the bulk and shear moduli set to $B=8\text{\,kPa},\ G=1.5\text{\,kPa}$ (background) and $B=20\text{\,kPa},\ G=4\text{\,kPa}$ (inclusion). The mass density is uniform, with $\rho=1000\,\text{kg/m}^{3}$. Both components of the displacement are assumed to be measured over a dense grid of points in the subdomain $\Omega^{\text{\scriptsize m}}$ (delineated with a dashed line in Figure 4(a)). Corresponding synthetic measurements were generated with a fine finite element mesh and interpolated onto a coarser, regular, reconstruction mesh made of $69\times 69$ square 4-noded elements. Each element of that mesh supports two material unknowns $B,G$. The boundary conditions on all four sides of $\Omega^{\text{\scriptsize m}}$ are taken as unknown.

The reconstruction results shown in Figure 4(b) demonstrate that material properties can be imaged accurately even without the knowledge of boundary conditions; in particular, the recovered values of $B$ and $G$ are close on average to their target values. Interestingly, the reconstructed inclusion displays clear and sharp edges despite the fact that no additional regularization is used for the fields $B,G$.

7. Proofs

7.1. Proof of Lemma 3.3

Let $\widehat{\bm{u}}\in\mathcal{U}$ solve the variational problem \big{(}\hskip 1.00006pt\widetilde{\bm{u}}{},\widehat{\bm{u}}\hskip 1.00006pt\big{)}_{\omega}=\mathcal{B}(\widetilde{\bm{u}}{},\bm{w}) for all $\widetilde{\bm{u}}{}\in\mathcal{U}$ (i.e. $\widehat{\bm{u}}$ is the $\bm{w}$-dependent Riesz representant of the linear functional $B^{\text{t}}\bm{w}\in\mathcal{U}^{\prime}$). We then have

[TABLE]

Let $(\bm{\psi}_{n})_{n\geq 1}$ and $(\omega_{n})_{n\geq 1}$ denote the countable sets of eigenfunctions and eigenvalues associated with problem (27), i.e. such that each $\bm{\psi}_{n}\in\mathcal{U}$ verifies $\mathcal{A}(\bm{\psi}_{n},\widetilde{\bm{u}}{})-\omega_{n}^{2}(\rho\bm{\psi}_{n},\widetilde{\bm{u}}{})=0$ for all $\widetilde{\bm{u}}{}\in\mathcal{U}$ (we assume for definiteness the normalization \big{(}\hskip 1.00006pt\rho\bm{\psi}_{n},\bm{\psi}_{n}\hskip 1.00006pt\big{)}=1 and the ordering $0\leq\omega_{1}\leq\omega_{2}\leq\ldots$); then $(\bm{\psi}_{n})_{n\geq 1}$ is a Hilbert basis of $\mathcal{U}$. Let $I:=\{n\in\mathbb{N},\;\omega=\omega_{n}\}$, so that \mathcal{Z}=\text{span}\big{(}\hskip 1.00006pt\bm{\psi}_{n,n\in I}\hskip 1.00006pt\big{)}, noting that $\mathcal{Z}$ is finite-dimensional. With this convention, $I=\emptyset$, i.e. $\mathcal{Z}=\{\mathbf{0}\}$, if $\omega$ is not an eigenvalue for problem (27). Expanding the fields $\widehat{\bm{u}}\in\mathcal{U}$ and $\bm{w}\in\mathcal{W}$ on the Hilbert basis $(\bm{\psi}_{n})_{n\geq 1}$ (noting for the latter that $\mathcal{W}\subset\mathcal{U}$) as

[TABLE]

the weak problem linking $\widehat{\bm{u}}$ to $\bm{w}$ then implies (using $\widetilde{\bm{u}}{}=\psi_{n}$ as test functions)

[TABLE]

We therefore obtain

[TABLE]

To establish the existence of an inf-sup constant $\beta$ such that (29) holds, we separately examine two cases: (i) $\mathcal{W}=\mathcal{U}$ and (ii) $\mathcal{W}\subsetneq\mathcal{U}$ (recall that $\mathcal{W}\subset\mathcal{U}$).

Case (i): $\mathcal{W}=\mathcal{U}$.

This case corresponds to $\Gamma=\Gamma_{\text{D}}\cup\Gamma_{\text{N}}$, and we have $\mathcal{H}=\mathcal{K}=\mathcal{Z}$. Let $P$ be the orthogonal projection on $\mathcal{Z}$, so that

[TABLE]

Using (66), we have

[TABLE]

Moreover, $\|(I-P)\bm{w}\|^{2}_{\omega}=\|\bm{w}\|^{2}_{\omega}$ since $\bm{w}\in\mathcal{K}{}^{\perp}=\mathcal{Z}^{\perp}$. Therefore

[TABLE]

Since $\beta$ does not depend on $\bm{w}$, the inf-sup condition (29) holds true with the above value of $\beta$.

Case (ii): $\mathcal{W}\subsetneq\mathcal{U}$.

In this case, $\mathcal{K}=\{\mathbf{0}\}$, so that the infimum in (29) is taken over $\bm{w}\in\mathcal{W}$. If $I=\emptyset$, then $\mathcal{Z}=\{\mathbf{0}\}$ and we again have $\|(I-P)\bm{w}\|^{2}_{\omega}=\|\bm{w}\|^{2}_{\omega}$; consequently, the argument of case (i) still applies, equations (67) and (68) remain valid (with infimums taken over all integers), and the inf-sup condition (29) again holds with $\beta$ as given by (68).

On the other hand, a distinct approach is needed for the case $I\not=\emptyset$ because $\|(I-P)\bm{w}\|^{2}_{\omega}=\|\bm{w}\|^{2}_{\omega}$ no longer holds for any $\bm{w}\in\mathcal{W}$. Instead, we will show below that

[TABLE]

implying

[TABLE]

The lemma results from (70), since that value of $\beta$ does not depend on $\bm{w}$, provided (69) holds.

Proof of (69).

We first note that $\mathcal{W}\cap\mathcal{Z}=\{\mathbf{0}\}$: $\bm{w}\in\mathcal{W}$ requires $\bm{w}=\mathbf{0}$ on $\Gamma_{c}$, whereas $\bm{w}\in\mathcal{Z}$ implies $\bm{L}(\omega)\bm{w}=\mathbf{0}$ and $\bm{t}[\bm{w}]=\mathbf{0}$ on $\Gamma_{c}$, and the unique continuation principle implies that only $\bm{w}=\mathbf{0}$ can fulfill all requirements.

The proof then proceeds by contradiction. Assume that (69) is false. In that case, we have

[TABLE]

Choose a sequence $(\xi_{n})>0$ such that $\xi_{n}\to 0$, $n\to\infty$, and for each $\xi_{n}$ choose $\bm{w}_{n}$ such that $\|\bm{w}_{n}-P\bm{w}_{n}\|_{\omega}<\xi_{n}\|\bm{w}_{n}\|_{\omega}$. $P$ being linear, $\|\bm{w}_{n}\|=1$ may be assumed for each $n$ without detriment. Then, $P$ being a projection, $\|P\bm{w}_{n}\|\leq 1$: each $P\bm{w}_{n}$ belongs to the unit ball $Z=\{\bm{z}\in\mathcal{Z},\;\|\bm{z}\|\leq 1\}$ of $\mathcal{Z}$. As $\text{dim}(\mathcal{Z})<\infty$, (i) $Z$ is compact, so the sequence $(P\bm{w}_{n})$ contains a convergent subsequence (still denoted $(P\bm{w}_{n})$), whose limit is denoted $\bm{z}$, and (ii) $\mathcal{Z}$ is closed, so $\bm{z}\in\mathcal{Z}$. On the other hand, we have

[TABLE]

implying that $\bm{w}_{n}\to\bm{z}$ as $n\to\infty$. Since $\bm{w}_{n}\in\mathcal{W}$ and the Dirichlet trace on $\Gamma_{u}\cup\Gamma_{c}$ is continuous, $\bm{z}\in\mathcal{W}$; moreover we have that $\|\bm{z}\|_{\omega}=1$ (as the limit of a sequence of elements of unit norm), and hence $\bm{w}\not=\mathbf{0}$. Summarizing, we simultaneously have $\bm{w}\in\mathcal{W}\cap\mathcal{Z}$ and $\bm{w}\not=\mathbf{0}$, which leads to a contradiction because $\mathcal{W}\cap\mathcal{Z}=\{\mathbf{0}\}$. Concluding, (69) is true.

7.2. Proof of Theorem 3.5

The proof methods follows that of [7, Thm. 4.3.1]. We first observe that setting $\widetilde{\bm{u}}{}=-\bm{u}$ and $\widetilde{\bm{w}}{}=\bm{w}$ in (24) and adding the resulting equalities yields \mathcal{A}(\bm{w},\bm{w})+\kappa\mathcal{D}(\bm{u},\bm{u})=\kappa\big{\langle}\hskip 1.00006pt\bm{d},\bm{u}\hskip 1.00006pt\big{\rangle}_{\mathcal{U}^{\prime},\mathcal{U}}+\big{\langle}\hskip 1.00006pt\bm{f},\bm{w}\hskip 1.00006pt\big{\rangle}_{\mathcal{W}^{\prime},\mathcal{W}} which, by virtue of the assumed coercivity of $\mathcal{A}$ on $\mathcal{W}\times\mathcal{W}$, implies the inequality

[TABLE]

Then, equation (24b) with $\widetilde{\bm{u}}{}=\bm{u}_{0}\in\mathcal{H}$ gives

[TABLE]

which implies the inequality \delta\|\bm{u}_{0}\|^{2}\leq\|\bm{d}_{0}\|\big{(}\hskip 1.00006pt\|\bm{u}_{0}\|+d\,\|\bm{u}_{1}\|\hskip 1.00006pt\big{)}, where the left-hand side results from $\mathcal{D}$ being coercive on $\mathcal{H}\times\mathcal{H}$ by Assumption 1, and the right-hand side from assumed continuity of $\mathcal{D}$ and $\bm{d}_{0}$. This results in

[TABLE]

Finally, using equation (24a) in operator form (i.e. $A\bm{w}+B\bm{u}_{1}=\bm{f}$, since $B\bm{u}_{0}=\mathbf{0}$), we have

[TABLE]

which, using the inf-sup condition (29) (which implies that $\|B\bm{u}_{1}\|\geq\beta\|\bm{u}_{1}\|$ for any $\bm{u}_{1}\in\mathcal{H}^{\perp}$), yields the inequality

[TABLE]

We now exploit inequalities (71), (72) and (74), considering separately each of the three possible cases where only one of $\bm{f},\bm{d}_{0},\bm{d}_{1}$ is nonzero. Considering first the case $\bm{f}\not=\mathbf{0}$, inequalities (71), (72) and (74) readily provide the estimates

[TABLE]

with $q_{a},q_{d}$ as defined in the statement of the proposition.

Consider next the case $\bm{d}_{1}\not=\mathbf{0}$. Inequalities (71), (72) and (74) then become

[TABLE]

(with $r_{a}$ as defined in the statement of the proposition), from which we obtain the estimates

[TABLE]

Finally, for the case $\bm{d}_{0}\not=\mathbf{0}$, inequalities (71), (72) and (74) provide

[TABLE]

Concatenating the above three bounds gives the inequality

[TABLE]

Being of the form $\|\bm{w}\|^{2}-A\|\bm{w}\|-B\leq 0$ with $A,B>0$, it holds for all $\|\bm{w}\|\leq W$ with $W$ the positive root of $w^{2}-Aw-B=0$, i.e.:

[TABLE]

The remaining sought estimates are then

[TABLE]

The stationarity problem (24) being linear, the estimates for $\bm{u}_{0},\,\bm{u}_{1},\,\bm{w}$ for general right-hand sides $\bm{d}=\bm{d}_{0}+\bm{d}_{1}$ and $\bm{f}$ follow from the triangle inequality.

Finally, since $\mathcal{U}=\mathcal{H}\oplus\mathcal{H}^{\perp}$, we have $\|\mathsf{F}\|_{\mathbb{U}}^{2}=\|\kappa\bm{d}_{0}\|^{2}+\|\kappa\bm{d}_{1}\|^{2}+\|\bm{f}\|^{2}$, implying $\|\kappa\bm{d}_{0}\|\leq\|\mathsf{F}\|_{\mathbb{U}}$, $\|\kappa\bm{d}_{1}\|\leq\|\mathsf{F}\|_{\mathbb{U}}$ and $\|\bm{f}\|\leq\|\mathsf{F}\|_{\mathbb{U}}$. The obtained estimates of $\|\bm{u}_{0}\|$, $\|\bm{u}_{1}\|$ and $\|\bm{w}\|$ therefore collectively imply that there exists a constant $C>0$ such that $\|\mathsf{W}\|_{\mathbb{U}}<C\|\mathsf{G}\mathsf{W}\|_{\mathbb{U}}$ for any $\mathsf{W}\in\mathbb{U}$. Well-posedness of problem (24) follows by lemma 3.2 (with $\eta=C^{-1}$).

7.3. Proof of estimates (43)

Let $\bm{z}_{0}\in\mathcal{N}^{\perp}$ solve the problem

[TABLE]

which is well-posed ($\mathcal{D}$ being coercive on $\mathcal{N}^{\perp}\times\mathcal{N}^{\perp}$ by Assumption 4); moreover, $\bm{z}_{0}$ obeys the estimate $\|\bm{z}_{0}\|\leq\delta^{-1}\|\bm{d}_{0}\|$. Introducing the new unknown $\bm{y}:=\bm{u}-\bm{z}_{0}$, the coupled problem (24) becomes

[TABLE]

with $\bm{h}:=\bm{f}-B_{0}\bm{z}_{0}$ (since $\bm{z}_{0}\in\mathcal{N}^{\perp}$). Setting $\widetilde{\bm{u}}{}=-\bm{y}$ and $\widetilde{\bm{w}}{}=\bm{w}$ in the above system and adding the resulting equalities yields \mathcal{A}(\bm{w},\bm{w})+\kappa\mathcal{D}(\bm{y},\bm{y})=\big{\langle}\hskip 1.00006pt\bm{h},\bm{w}\hskip 1.00006pt\big{\rangle}_{\mathcal{W}^{\prime},\mathcal{W}}, which (by coercivity of $\mathcal{A}$ on $\mathcal{W}\times\mathcal{W}$ and applying the triangle inequality to the right-hand side) implies the inequality

[TABLE]

Then, equation (82b) with $\widetilde{\bm{u}}{}=\bm{u}_{0}\in\mathcal{N}^{\perp}$ gives (in operator form) $\kappa D\bm{y}_{0}=B^{\text{t}}_{0}\bm{w}$, and hence, using (83), implies

[TABLE]

Finally, using equation (82a) in operator form (i.e. $A\bm{w}+B_{0}\bm{y}_{0}+B_{1}\bm{u}_{1}=\bm{h}$), we have

[TABLE]

which, using the inf-sup condition (29) (which provides $\|B_{1}\bm{u}_{1}\|\geq\beta\|\bm{u}_{1}\|$ for any $\bm{u}_{1}\in\mathcal{N}$ since $\mathcal{N}\subset\mathcal{H}^{\perp}$), yields the inequality

[TABLE]

Estimates (43) finally stem from recalling that $\bm{u}_{0}=\bm{y}_{0}+\bm{z}_{0}$.

7.4. Proof of Proposition 3

First, recasting (54) using operator notation, we have

[TABLE]

We then observe that solving the derivative problem (55) gives

[TABLE]

with the operator $Z:\mathcal{W}\to\mathcal{W}^{\prime}$ as given in the proposition statement by $Z:=A+\kappa^{-1}B_{0}D^{-1}B^{\text{t}}_{0}$ and having set $\bm{H}:=-A^{\prime}\bm{w}-B^{\prime}_{0}\bm{u}_{0}-B^{\prime}_{1}\bm{u}_{1}$, $\bm{F}:=\bm{H}-\kappa^{-1}B_{0}D^{-1}B_{0}^{\prime}{}^{\text{t}}\bm{w}$, $\bm{G}:=-B_{1}^{\text{t}}{}^{\prime}\bm{w}$. Using (88) and noting that $\kappa^{-1}B_{0}D^{-1}B^{\text{t}}_{0}=Z-A$, we obtain the identities

[TABLE]

with the symmetric, positive operator $\Delta:\mathcal{W}\to\mathcal{W}^{\prime}$ given by $\Delta:=B_{0}^{\prime}D^{-1}B_{0}^{\prime}{}^{\text{t}}$. Upon substitution into (87), the above formulas provide

[TABLE]

Since (88) implies $\bm{F}=Z\bm{w}^{\prime}+B_{1}\bm{u}^{\prime}_{1}$, the above expression of $\tilde{\Lambda}{}_{\kappa}^{\prime\prime}$ yields the claimed formula.

7.5. Proof of Theorem 5.1

To begin, the expression (88) of $\bm{u}_{1}^{\prime}$ gives

[TABLE]

with the operator $P:\mathcal{W}^{\prime}\to\mathcal{W}^{\prime}$ defined by $P:=B_{1}(B_{1}^{\text{t}}Z^{-1}B_{1})^{-1}B_{1}^{\text{t}}Z^{-1}$. It is easy to see that $P$ verifies $PP=P$ (i.e. $P$ is a projection) and $Z^{-1}P=P^{\text{t}}Z^{-1}$.

For the case where $B_{1}^{\prime}=0$, we have $\bm{G}=\mathbf{0}$, and therefore $B_{1}\bm{u}^{\prime}_{1}=P\bm{F}$. The expression (93) of \tilde{\Lambda}{}_{\kappa}^{\prime\prime}(\text{\boldmath\mathcal{C}}) then becomes (with the above definition of $P$)

[TABLE]

Next, applying the Sherman-Morrison-Woodbury formula to $Z=A+\kappa^{-1}B_{0}D^{-1}B^{\text{t}}_{0}$ gives

[TABLE]

wherein $A_{11}:=B_{1}^{\text{t}}A^{-1}B_{1}$, $A_{01}:=B_{0}^{\text{t}}A^{-1}B_{1}$, $A_{10}:=B_{1}^{\text{t}}A^{-1}B_{0}$. Using these identities and performing straightforward algebra, we find

[TABLE]

with the (projection) operator $P_{0}$ as defined in the theorem statement. We also note that $R^{-1}=\kappa^{-1}D^{-1}+o(\kappa^{-1})$, so that

[TABLE]

We now substitute the above expansion into (95) and set $\bm{F}=\bm{F}^{(0)}+\kappa^{-1}\bm{F}^{(1)}+o(\kappa^{-1})$ (noting that $\bm{F}^{(0)}=-A^{\prime}\bm{w}^{(0)}-B^{\prime}_{0}\bm{u}^{\text{\scriptsize m}}$) and $\bm{w}=\bm{w}^{(0)}+\kappa^{-1}\bm{w}^{(1)}+o(\kappa^{-1})$, to obtain

[TABLE]

with

[TABLE]

Concluding, the above value of $\tilde{\Lambda}{}_{\kappa}^{\prime\prime}{}^{(0)}$ is that claimed in the theorem statement, and is clearly positive, whereas $\tilde{\Lambda}{}_{\kappa}^{\prime\prime}{}^{(1)}$ is obtained as an algebraic sum of positive quadratic expressions (so is a priori sign-indefinite). The proof of Theorem 5.1 is complete.

8. Conclusions

In this work we studied some of the most salient mathematical properties of the Modified Error in Constitutive Equations (MECE) approach for inverse problems in the context of frequency-domain elastodynamics. In particular, we proved (under conditions on the available interior data) well-posedness of the coupled system that arises as part of the first order optimality conditions. The coupled problem remains well posed even when boundary conditions are partially or completely unknown and at resonant frequencies. The latter findings have strong practical implications in inverse problems where interior data is abundant and boundary conditions are difficult to ascertain. We have exploited this benefit in our recent work in elastography [16].

We also showed that the reduced MECE functional becomes convex in the limit where the weight given to the data misfit component goes to infinity. Moreover, convexification of the reduced objective occurs continuously as the parameter increases, as demonstrated in the numerical examples. This characteristic of the MECE functional also has strong practical implications as solutions to the inverse problem are less sensitive to initial guesses. This fact has also been exploited in our work on elasticity and viscoelasticity imaging [12, 16]. Future work includes exploring the possibility of devising a unified MECE formulation that can be applied to a wide range of physics such as nonlinear elasticity, electromagnetism, plasticity, fluid dynamics, etc.