Lipschitz stability estimate and reconstruction of Lam\'e parameters in linear elasticity

TL;DR

This paper establishes a Lipschitz stability estimate for recovering isotropic elastic parameters from boundary measurements in linear elasticity and proposes a numerical iterative reconstruction method with numerical examples.

Contribution

It provides the first Lipschitz stability estimate for Lamé parameters under certain regularity and monotonicity assumptions, along with a practical reconstruction algorithm.

Findings

Lipschitz stability estimate proven for Lamé parameters

Numerical reconstruction algorithm developed and tested

Numerical examples demonstrating effectiveness

Abstract

In this paper, we consider the inverse problem of recovering an isotropic elastic tensor from the Neumann-to-Dirichlet map. To this end, we prove a Lipschitz stability estimate for Lam\'e parameters with certain regularity assumptions. In addition, we assume that the Lam\'e parameters belong to a known finite subspace with a priori known bounds and that they fulfill a monotonicity property. The proof relies on a monotonicity result combined with the techniques of localized potentials. To numerically solve the inverse problem, we propose a Kohn-Vogelius-type cost functional over a class of admissible parameters subject to two boundary value problems. The reformulation of the minimization problem via the Neumann-to-Dirichlet operator allows us to obtain the optimality conditions by using the Fr\'echet differentiability of this operator and its inverse. The reconstruction is then performed by means of an iterative algorithm based on a quasi-Newton method. Finally, we give and discuss several numerical examples.

Figures

Tables

Table 1. Table 1: Simulation results for Example 1: Reconstruction of constant Lamé parameters.

	$({\lambda }_i,{\mu }_i)$	$({\lambda }_c,{\mu }_c)$	$({\lambda }_e,{\mu }_e)$	$\frac{\left|{\lambda }_c-{\lambda }_e\right|}{\left|{\lambda }_e\right|}$	$\frac{\left|{\mu }_c-{\mu }_e\right|}{\left|{\mu }_e\right|}$
$ϵ=0.0,\rho =0.0$	(1,1)	(2.9999, 6.9999)	(3,7)	3.33e-07	1.42e-07
$ϵ=0.03,\rho =0.00001$	(1,1)	(2.6564, 6.8989)	(3,7)	0.1145	0.0144
$ϵ=0.05,\rho =0.00001$	(1,1)	( 2.6527,6.8869)	(3,7)	0.1157	0.0161

Equations

$$\partial \Omega =\overline{{\Gamma }_{\text{N}}}\cup \overline{{\Gamma }_{\text{D}}},{\Gamma }_{\text{N}}\cap {\Gamma }_{\text{D}}=\varnothing .$$

$$\{\begin{array}{rlrl}-\mathrm{div}(\mathbb{C}\hat{\nabla }u)& =0& & \text{ in }\Omega ,\\ (\mathbb{C}\hat{\nabla }u)\nu & =g& & \text{ on }{\Gamma }_{\text{N}},\\ u& =0& & \text{ on }{\Gamma }_{\text{D}},\end{array}$$

$$\hat{\nabla }u=\frac{1}{2}\left(\nabla u+(\nabla u{)}^T\right),\mathbb{C}\hat{\nabla }u={\left(\sum _{k,l=1}^d{\mathbb{C}}_{ijkl}\frac{\partial u_k}{\partial x_l}\right)}_{1\le i,j\le d},$$

$$\nabla u(x)=\left(\frac{\partial }{\partial x_1}u(x),\cdots ,\frac{\partial }{\partial x_d}u(x)\right),\mathrm{div}F(x)=\left(\begin{array}{c}{\sum }_{i=1}^d\frac{\partial }{\partial x_i}{F}_{i1}(x)\\ ⋮\\ {\sum }_{i=1}^d\frac{\partial }{\partial x_i}{F}_{id}(x)\end{array}\right)$$

$${\mathbb{C}}_{ijkl}:=\lambda {\delta }_{ij}{\delta }_{kl}+\mu \left({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}\right),$$

$$\displaystyle\mu(x)\geq\delta_{0},\quad\lambda(x)+2\mu(x)\geq\delta_{0}\quad\mathrm{a.\,e.}\,\,x\in{\color[rgb]{0,0,0}\Omega},$$

$$\displaystyle\|\mu\|_{C^{0,1}({\color[rgb]{0,0,0}\Omega})}+\|\lambda\|_{L^{\infty}(\Omega)}\leq M_{0},$$

$$-\mathrm{div}(\mathbb{C}\hat{\nabla }u)$$

$$u$$

$$-\mathrm{div}(\mathbb{C}\hat{\nabla }u)$$

$$(\mathbb{C}\hat{\nabla }u)\nu$$

$$u$$

$$\mathcal{A}=\left\{\mathbb{C}=\mathbb{C}(\lambda ,\mu ):(\lambda ,\mu )\in L^{\infty }(\Omega )\times C^{0,1}(\Omega ),{\alpha }_1\le \lambda \le {\alpha }_2,{\beta }_1\le \mu \le {\beta }_2\right\}.$$

$${\int }_{\Omega }\mathbb{C}\hat{\nabla }u:\hat{\nabla }vdx={\int }_{{\Gamma }_{\text{N}}}g\cdot vds\text{ for all }v\in \mathcal{V},$$

$$\mathcal{V}:=\left\{v\in H^1(\Omega {)}^d:v_{{|}_{{\Gamma }_{\text{D}}}}=0\right\}.$$

$$\Lambda (\mathbb{C}):L^2({\Gamma }_{\text{N}}{)}^d\to L^2({\Gamma }_{\text{N}}{)}^d:g\mapsto u_{\mathbb{C}}^g{|}_{{\Gamma }_{\text{N}}},$$

$$⟨g,\Lambda (\mathbb{C})h⟩={\int }_{\Omega }\mathbb{C}\hat{\nabla }{u}_{\mathbb{C}}^g:\hat{\nabla }{u}_{\mathbb{C}}^{h}dx\text{ for all }g,h\in L^2({\Gamma }_{\text{N}}{)}^d,$$

$$\text{ Find }\mathbb{C}\text{ or }(\lambda ,\mu )\text{ knowing}⟨g,\Lambda (\mathbb{C})h⟩.$$

$${\int }_{\Omega }({\mathbb{C}}_1-{\mathbb{C}}_2)\hat{\nabla }{u}_2:\hat{\nabla }{u}_{2}dx\ge ⟨g,\Lambda ({\mathbb{C}}_2)g⟩-⟨g,\Lambda ({\mathbb{C}}_1)g⟩\ge {\int }_{\Omega }({\mathbb{C}}_1-{\mathbb{C}}_2)\hat{\nabla }{u}_1:\hat{\nabla }{u}_{1}dx.$$

$${\int }_{\Omega }{\mathbb{C}}_1\hat{\nabla }{u}_1:\hat{\nabla }{u}_{2}dx=⟨g,\Lambda ({\mathbb{C}}_2)g⟩={\int }_{\Omega }{\mathbb{C}}_2\hat{\nabla }{u}_2:\hat{\nabla }{u}_{2}dx.$$

$${\int }_{\Omega }{\mathbb{C}}_1\hat{\nabla }(u_1-u_2):\hat{\nabla }(u_1-u_2)dx$$

$$={\int }_{\Omega }{\mathbb{C}}_1\hat{\nabla }{u}_1:\hat{\nabla }{u}_{1}dx+{\int }_{\Omega }{\mathbb{C}}_1\hat{\nabla }{u}_2:\hat{\nabla }{u}_{2}dx-2{\int }_{\Omega }{\mathbb{C}}_1\hat{\nabla }{u}_1:\hat{\nabla }{u}_{2}dx$$

$$=⟨g,\Lambda ({\mathbb{C}}_1)g⟩-⟨g,\Lambda ({\mathbb{C}}_2)g⟩+{\int }_{\Omega }({\mathbb{C}}_1-{\mathbb{C}}_2)\hat{\nabla }{u}_2:\hat{\nabla }{u}_{2}dx.$$

$${\lambda }_1\le {\lambda }_2\text{ and }{\mu }_1\le {\mu }_2\text{ implies }\Lambda ({\mathbb{C}}_1)\ge \Lambda ({\mathbb{C}}_2).$$

$$\underset{n\to \infty }{\lim }{\int }_{D_1}{\left(\mathrm{div}{u}^{g_n}\right)}^{2}dx=\infty ,$$

$$\underset{n\to \infty }{\lim }{\int }_{D_2}{\left(\mathrm{div}{u}^{g_n}\right)}^{2}dx=0,$$

$$\underset{n\to \infty }{\lim }{\int }_{D_1}\hat{\nabla }{u}^{g_n}:\hat{\nabla }{u}^{g_n}dx=\infty ,$$

$$\underset{n\to \infty }{\lim }{\int }_{D_2}\hat{\nabla }{u}^{g_n}:\hat{\nabla }{u}^{g_n}dx=0.$$

$$\displaystyle A_{j}:{\color[rgb]{0,0,0}L^{2}(D_{j})}\rightarrow L^{2}(\Gamma_{\textup{N}})^{d},\quad F\mapsto v|_{\Gamma_{\textup{N}}},$$

$${\int }_{\Omega }\mathbb{C}\hat{\nabla }v:\hat{\nabla }wdx={\int }_{D_j}F\mathrm{div}wdx\text{ for all }w\in \mathcal{V},$$

$$\displaystyle B_{j}:{\color[rgb]{0,0,0}L^{2}(D_{j})^{d\times d}}\rightarrow L^{2}(\Gamma_{\textup{N}})^{d},\quad G\mapsto v|_{\Gamma_{\textup{N}}},$$

$${\int }_{\Omega }\mathbb{C}\hat{\nabla }v:\hat{\nabla }wdx={\int }_{D_j}G:\hat{\nabla }wdx\text{ for all }w\in \mathcal{V}.$$

$$\displaystyle A^{\prime}_{j}:L^{2}(\Gamma_{\textup{N}})^{d}\rightarrow{\color[rgb]{0,0,0}L^{2}(D_{j})},\quad j=1,2,$$

$$\displaystyle B^{\prime}_{j}:L^{2}(\Gamma_{\textup{N}})^{d}\rightarrow{\color[rgb]{0,0,0}L^{2}(D_{j})^{d\times d}},\quad j=1,2,$$

$$\displaystyle{\color[rgb]{0,0,0}\int_{\Omega}FA^{\prime}_{j}g\,dx=\int_{\Gamma_{\textup{N}}}g\cdot A_{j}F\,ds=\int_{\Omega}\mathbb{C}\hat{\nabla}v:\hat{\nabla}u\,dx=\int_{D_{j}}F\mathrm{div}(u)\,dx.}$$

Full text

1 Introduction

In this paper, we consider the inverse problem of recovering the elastic tensor $\mathbb{C}$ of a linear isotropic elastic body from the Neumann-to-Dirichlet operator $\Lambda(\mathbb{C})$. The main motivations of this problem are non-destructive testing for elastic bodies in order to detect and reconstruct material inclusions (as presented in [1] and the references therein), geophysical (see, e.g., [2]), and medical applications (as considered in [3]), in particular localization of potential tumors via a medical imaging modality called elastography. Elastography is concerned with the reconstruction of the elastic properties in biological tissues and the present article aims at giving access to these features. For a topical review of inverse problems in elasticity and their applications the reader is, e.g., referred to [4].

The paper is split into two parts and gives a twofold perspective on determining Lamé parameters in linear elasticity. Part one is on proving a Lipschitz stability estimate based on the Neumann-to-Dirichlet map for the corresponding tensors $\mathbb{C}_{1}$ and $\mathbb{C}_{2}$, depending on the Lamé parameters $(\lambda_{1},\mu_{1})$ and $(\lambda_{2},\mu_{2})$, respectively. The estimate is given under the following a priori assumptions.

• (i)

The Lamé parameters $(\lambda_{j},\mu_{j})$ for $j=1,2$ belong to a known finite dimensional subspace, with $\lambda$ being piecewise continuous and $\mu$ being Lipschitz continuous, and with fixed uniform lower and upper bounds.

• (ii)

The two pairs of parameters satisfy a monotonicity assumption that $(\lambda_{1},\mu_{1})$ is, as a pair, either a lower or an upper bound to $(\lambda_{2},\mu_{2})$.

Part two deals with the reconstruction of the parameters $(\lambda,\mu)$ based on minimizing a Kohn-Vogelius type functional. Let us stress that this numerical part does not build on the theoretical results presented in Section 3 but rather approaches the problem from a heuristic numerical side to demonstrate that useful numerical reconstructions are indeed possible. It remains a challenging open task how to unite the theoretical and numerical approaches in order to find rigorously justified reconstruction methods that work well in practically relevant settings.

From the theoretical point of view, the inverse problem of recovering $\mathbb{C}$ (or the Lamé moduli $\lambda,\mu$; cf. (2)) has been studied by several authors. In the two dimensional case Ikehata [6] proves that the deflection $h$ between $(\lambda+h,\mu+h)$ and $(\lambda,\mu)$ can be uniquely determined by the first-order approximation of the Dirichlet-to-Neumann operator. Akamatsu, Nakamura and Steinberg [7], give an inversion formula for the normal derivatives at the boundary of the Lamé coefficients $\lambda,\mu\in C^{\infty}$ from the Dirichlet-to-Neumann map. At the same time they present stability estimates for the boundary values of $\lambda,\mu$. Nakamura and Uhlmann [8] established that the Lamé coefficients are uniquely determined from the Dirichlet-to-Neumann operator, assuming that they are sufficiently close to a pair of positive constants. Imanuvilov and Yamamoto [9] proved that the Lamé coefficient $\lambda$ can be recovered from partial Cauchy data if the coefficient $\mu$ is some positive constant. Boundary determination of Lamé coefficients can be found in the work by Lin and Nakamura [11]. In addition, we want to mention a global uniqueness results in 2D by Imanuvilov and Yamamoto [12].

For the three dimensional case, Nakamura and Uhlmann [13, 14] and Eskin and Ralston [15] proved uniqueness results for both Lamé coefficients when $\mu$ is assumed to be close to a positive constant. The proofs in the above papers rely on the construction of complex geometric optics solutions. For a partial data version, uniqueness for recovering piecewise constant Lamé parameters and a quantitative Lipschitz stability result was proved in [16, 17], and some boundary determination results were shown in [18, 13, 11]. For fully anisotropic $\mathbb{C}$, uniqueness was proved by C$\hat{\mathrm{a}}$rstea, Honda and Nakamura [19] for a piecewise homogeneous medium. Isakov, Wang and Yamamoto [20], proved Hölder and Lipschitz stability estimates of determining all coefficients of a dynamical Lamé system with residual stress, including the density Lamé parameters, and the residual stress, by three pairs of observations from the whole boundary or from a part of it.

In this paper, we prove a Lipschitz stability result when the Lamé coefficient $\lambda$ is piecewise continuous, $\mu$ is Lipschitz and a monotonicity assumption holds. In more detail, we assume that the Lamé parameters belong to a known finite subspace with a priori known lower and upper bounds. Our approach relies on the monotonicity of the Neumann-to-Dirichlet operator with respect to the elastic tensor and the techniques of localized potentials [26, 29, 30, 31, 32, 34, 35]. For the numerical solution, we reformulate the inverse problem into a minimization problem using a Kohn-Vogelius type cost functional, and use a quasi-Newton method which employs the analytic gradient of the cost function and the approximation of the inverse Hessian is updated by a BFGS (Broyden, Fletcher, Goldfarb, Shanno) scheme [39].

Let us give some more remarks on the relation of this work to previous results. Stability for inverse coefficient problems are derived in general from technically challenging approaches involving quantitative unique continuation estimates and the study of special solutions [40, 41, 13, 14]. Our approach on proving a Lipschitz stability result is relatively simple and easy to extend to other settings, and has already led to new results on uniqueness and Lipschitz stability in electrical impedance tomography (EIT) with finitely many electrodes [42] as well as for the inverse Robin transmission problem [43, 44] and on the stability in machine learning reconstruction algorithms [45] under a definiteness assumption.

The paper is organized as follows. In Section 2, we introduce the forward as well as the inverse problem and the Neumann-to-Dirichlet operator. In Section 3, we show a monotonicity result between the Lamé parameters and the Neumann-to-Dirichlet operator and deduce the existence of localized potentials. Then, we prove the Lipschitz stability estimate for a finite dimensional subset with bounded Lamé parameters and a monotonicity assumption. In Section 4, we present a numerical approach for solving the inverse problem and in Section 5 we show some numerical results.

2 Problem formulation

Let $\Omega\subset\mathbb{R}^{d}$ $(d\geq 2)$, be a bounded and connected open set, occupied by an isotropic material with linear stress-strain relation. The boundary $\partial\Omega$, is assumed to be $C^{1,1}$ and consists of two non-empty disjoint open parts, the fixed ”Dirichlet-boundary” $\Gamma_{\textup{D}}$ and the ”Neumann-boundary” $\Gamma_{\textup{N}}$ with

[TABLE]

The choice of mixed boundary conditions is based on the physical treatment of the elasticity problem. The Neumann-to-Dirichlet operator with fixed Dirichlet part is an idealized model for fixing an elastic object in place on one part of the boundary, applying different pressure patterns to the remaining part and measuring the resulting displacements.

We denote a given surface load by $g\in L^{2}(\Gamma_{\textup{N}})^{d}$. Then the displacement vector $u:\overline{\Omega}\rightarrow\mathbb{R}^{d}$ satisfies the following boundary value problem

[TABLE]

where $\nu$ is the outer unit normal vector to $\partial\Omega$, which is similar to the boundary value problem considered, e.g., in [46]. The linearized strain tensor $\hat{\nabla}u$ and the stress tensor $\mathbb{C}\hat{\nabla}u$ are given by

[TABLE]

where

[TABLE]

for $F:\Omega\mapsto\mathbb{R}^{d\times d}$.

The isotropic elastic tensor is defined as

[TABLE]

where $\lambda,\mu$ are the Lamé coefficients.

Next, we state a unique continuation principle from [47]:

Theorem 2.1** (Weak Unique Continuation Principle).**

Let $\Omega$ be a connected open set in $\mathbb{R}^{d}$ for $d\geq 2$. Let $\mu(x)\in C^{0,1}(\Omega)$ and $\lambda(x)\in L^{\infty}(\Omega)$ satisfy

[TABLE]

with positive constants $\delta_{0}$, $M_{0}$, where we define \|f\|_{C^{0,1}(\Omega)}=\|f\|_{L^{\infty}({\color[rgb]{0,0,0}\Omega})}+\|\nabla f\|_{L^{\infty}(\Omega)}. If $u\in H^{1}(\Omega)^{d}$ satisfies

[TABLE]

then $u=0$ in $\Omega$.

Proof.

The reader is referred to Theorem 1.2 from [47]. ∎

Corollary 2.2** (Unique Continuation Principle for Local Cauchy Data).**

Let the same assumptions as in Theorem 2.1 hold and let $\Omega$ have a $C^{1,1}$ boundary and let $\Gamma$ be a nonempty open subset of $\partial\Omega$. If $u\in H^{1}(\Omega)^{d}$ satisfies

[TABLE]

*then $u=0$ in $\Omega$. *

Proof.

A solution with vanishing Cauchy data in $\Omega$ can be extended by zero to a solution in an extended open set $\tilde{\Omega}$, where $\Omega\subsetneq\tilde{\Omega}$. Hence, the weak UCP (Theorem 2.1) can be applied to show that $u=0$ in $\Omega$. ∎

Next, for given constants \alpha_{1},\alpha_{2},\beta_{1},{\color[rgb]{0,0,0}\beta_{2}} satisfying $0<\alpha_{1}\leq\alpha_{2}$, $0<\beta_{1}\leq\beta_{2}$, we define the set of admissible elastic tensors by

[TABLE]

Hence, the Lamé parameters of every $\mathbb{C}(\lambda,\mu)\in\mathcal{A}$ satisfy the conditions of the unique continuation principles Theorem 2.1 and Corollary 2.2.

In what follows, we denote $A:B=\sum_{i,j=1}^{d}a_{ij}b_{ij}$, for matrices $A=(a_{ij})$ and $B=(b_{ij})$.

The weak formulation of problem (1) is given by

[TABLE]

where

[TABLE]

It is easy to see that for each $\mathbb{C}\in\mathcal{A}$, problem (3) has a unique solution $u\in\mathcal{V}$, which follows by the Lax-Milgram theorem and is shown, e.g., in [49].

We introduce the Neumann-to-Dirichlet operator $\Lambda(\mathbb{C})$:

[TABLE]

where here and in the following $u_{\mathbb{C}}^{g}\in\mathcal{V}$ always denotes the unique solution of (1) with surface load $g\in L^{2}(\Gamma_{\textup{N}})^{d}$ and elastic tensor $\mathbb{C}\in\mathcal{A}$.

It is well known (and an easy consequence from the variational formulation (3) and the compactness of the trace operator) that $\Lambda(\mathbb{C})$ is a self-adjoint compact linear operator that fulfills

[TABLE]

where $\langle\cdot,\cdot\rangle$ denotes the $L^{2}(\Gamma_{\textup{N}})$-inner product.

The inverse problem we consider here is the following:

[TABLE]

3 Monotonicity, localized potentials and Lipschitz stability

In this section, we show a monotonicity estimate between the elastic tensor and the Neumann-to-Dirichlet operator and the existence of localized potentials. Then we deduce a Lipschitz stability estimate for a finite dimensional subset with bounded Lamé parameters and a definiteness assumption.

Lemma 3.1** (Monotonicity estimate).**

Let $\mathbb{C}_{1},\mathbb{C}_{2}\in\mathcal{A}$ be elastic tensors, and $g\in L^{2}(\Gamma_{\textup{N}})^{d}$ be an applied boundary load. The corresponding solutions of (1) are denoted by $u_{1}:=u^{g}_{\mathbb{C}_{1}},\ u_{2}:=u^{g}_{\mathbb{C}_{2}}\in\mathcal{V}$. Then

[TABLE]

Proof.

Since $\Lambda(\mathbb{C}_{2})g=u_{2}|_{\Gamma_{\textup{N}}}$ we can use the variational formulation (3) for $\mathbb{C}_{1}$ and $\mathbb{C}_{2}$ with $v:=u_{2}$ and obtain

[TABLE]

Thus

[TABLE]

Since the left-hand side is nonnegative, the first asserted inequality follows.

Interchanging $\mathbb{C}_{1}$ and $\mathbb{C}_{2}$, the second inequality follows. ∎

Based on the previous lemma and the definition of $\mathbb{C}$, we are led to the following monotonicity property.

Corollary 3.2** (Monotonicity).**

For $\mathbb{C}_{1}:=\mathbb{C}(\lambda_{1},\mu_{1}),\mathbb{C}_{2}:=\mathbb{C}(\lambda_{2},\mu_{2})\in\mathcal{A}$

[TABLE]

Theorem 3.3** (Localized potentials).**

Let $\mathbb{C}\in\mathcal{A}$ and $D_{1},D_{2}$ be two open sets with $\overline{D}_{1},\overline{D}_{2}\Subset\Omega$, $\overline{D}_{1}\cap\overline{D}_{2}=\emptyset$ and let $\Omega\setminus(\overline{D}_{1}\cup\overline{D}_{2})$ be connected. Then there exists a sequence $(g_{n})_{n\in\mathbb{N}}\subset L^{2}(\Gamma_{\textup{N}})^{d}$, such that the corresponding solutions $(u^{g_{n}})_{n\in\mathbb{N}}$ of (1) fulfill

[TABLE]

Proof.

This proof is based on the UCP for local Cauchy data (Corollary 2.2). First, we define the virtual measurement operators

• (a)

$A_{j}$ ($j=1,2$) by

[TABLE]

where $v\in\mathcal{V}$ solves

[TABLE]

• (b)

$B_{j}$ ($j=1,2$) by

[TABLE]

where $v\in\mathcal{V}$ solves

[TABLE]

First, we show that the dual operators

[TABLE]

are given by $A^{\prime}_{j}g=\operatorname{div}(u)|_{D_{j}}$ and $B^{\prime}_{j}g=\hat{\nabla}u|_{D_{j}}$, where $u$ solves problem (1).

• To (a):

Let $F\in L^{2}(\Omega)$, $g\in L^{2}(\Gamma_{\textup{N}})^{d}$, $u$, $v\in\mathcal{V}$ solve (1) and (11), respectively. Then,

[TABLE]

• To (b):

Let $G\in L^{2}(\Omega)^{d\times d}$, $g\in L^{2}(\Gamma_{\textup{N}})^{d}$, $u$, $v\in\mathcal{V}$ solve (1) and (12), respectively. Then,

[TABLE]

Next, we will prove that

[TABLE]

Let $\varphi\in\mathcal{R}(A_{1})\cap\mathcal{R}(B_{2})$. Then there exist $v_{1},v_{2}\in\mathcal{V}$ such that $v_{1}|_{\Gamma_{\textup{N}}}=v_{2}|_{\Gamma_{\textup{N}}}=\varphi$, and

[TABLE]

for all $w\in\mathcal{V}$ with $\text{supp}(w)\subset\overline{\Omega}\setminus\overline{D}_{j}$, $j=1,2$. Hence,

[TABLE]

and $(\mathbb{C}\hat{\nabla}v_{1})\nu|_{\Gamma_{\textup{N}}}=(\mathbb{C}\hat{\nabla}v_{2})\nu|_{\Gamma_{\textup{N}}}=0$. The unique continuation principle for Cauchy data (Corollary 2.2) yields that $v_{1}=v_{2}$ in $\Omega\setminus(\overline{D}_{1}\cup\overline{D}_{2})$. Hence $v:=v_{1}\chi_{D_{2}}+v_{2}\chi_{\Omega\setminus\overline{D}_{2}}\in\mathcal{V}$ and satisfies

[TABLE]

It follows that $v=0$ and thus $\varphi=v|_{\Gamma_{\textup{N}}}=0$, and consequently $\mathcal{R}(A_{1})\cap\mathcal{R}(B_{2})=\left\{0\right\}$.

Next, we take a closer look at the operator $A_{1}:L^{2}(D_{1})\to L^{2}(\Gamma_{\textup{N}})^{d}$, $F\mapsto v|_{\Gamma_{\textup{N}}}$. Let $O$ be an open ball with $\overline{O}\subseteq D_{1}$ and let $v\in\mathcal{V}$ solve (11) for $F=\chi_{O}$. We will show that $v|_{\Gamma_{\textup{N}}}=A_{1}\chi_{O}\neq 0$. We argue by contradiction and assume $v|_{\Gamma_{\textup{N}}}=0$. From (11) we obtain that

[TABLE]

for all $w\in\mathcal{V}$ with $\mathrm{supp}(w)\subseteq\overline{\Omega}\setminus\overline{O}$. Thus, $v$ fulfills

[TABLE]

Now, we apply Corollary 2.2 on $\Omega\setminus\overline{O}$ which results in $v=0$ on $\Omega\setminus\overline{O}$, so that $v|_{O}\in H_{0}^{1}(O)^{d}$. Also, for all $w\in H_{0}^{1}(O)^{d}$, we obtain from (11) with $F=\chi_{O}$

[TABLE]

where $\tilde{w}\in H_{0}^{1}(\Omega)^{d}$ is the zero extension of $w\in H_{0}^{1}(O)^{d}$. Since (14) is uniquely solvable by the Lax-Milgram-Theorem it follows that $v|_{O}=0$, so that $v=0$ in all of $\Omega$ which contradicts (11) with $F=\chi_{O}$. Hence, $A_{1}\chi_{O}\neq 0$ which shows that $\mathcal{R}(A_{1})\neq\{0\}$.

This proves (13) and this implies $\mathcal{R}(A_{1})\not\subseteq\mathcal{R}(B_{2})$. Using [[34], Corollary 2.6] it follows that there exists a sequence $(g_{n})_{n\in\mathbb{N}}\subset L^{2}(\Gamma_{\textup{N}})^{d}$:

[TABLE]

and

[TABLE]

i.e. (7) and (10) hold. Since

[TABLE]

this also implicates (8) and (9).

∎

Next, we go over to the background of the Lipschitz stability and introduce the definition of piecewise continuous functions.

Definition 3.4**.**

A function $f\in L^{\infty}(\Omega)$ is called piecewise continuous, if there exists a finite decomposition of $\Omega$ into non-empty open subsets $\Omega_{i}\subseteq\Omega$, $i=1,...,n$, so that $\Omega\setminus\bigcup_{i=1}^{n}\Omega_{i}$ is a Lebesgue null set, $\Omega_{i}\cap\Omega_{j}=\emptyset$ ($i\neq j$), and $f|_{\Omega_{i}}$ is continuous for all $i=1,...,n$.

Inverse elliptic coefficient problems are known to be ill-posed and stability results can only be obtained under a-priori assumptions, cf. the works cited in the introduction. For our problem, we will prove a stability result under the assumption that the coefficients belong to an a-priori known finite-dimensional subspace (e.g., stemming from the parameter parametrization or a desired finite resolution), that upper and lower bounds are a-priori known, and that a definiteness condition holds. More precisely, let $\mathcal{F}$ be a finite dimensional subspace of $\check{C}(\Omega)\times C^{0,1}(\Omega)$, where $\check{C}(\Omega)$ is the space of piecewise continuous functions. We consider four constants $0<a\leq b$ and $0<c\leq d$ which are the lower and upper bounds of the Lamé parameter and define the sets

[TABLE]

[TABLE]

In the following main result of this paper, the domain $\Omega$, the finite-dimensional subspace $\mathcal{F}$ and the bounds $0<a\leq b$ and $0<c\leq d$ are fixed, and the constant in the Lipschitz stability result will depend on them.

Theorem 3.5** (Lipschitz stability).**

There exists a positive constant $C>0$ such that for all $\mathbb{C}_{1}:=\mathbb{C}(\lambda_{1},\mu_{1}),\mathbb{C}_{2}:=\mathbb{C}(\lambda_{2},\mu_{2})\in\mathcal{E}$ with either

[TABLE]

we have

[TABLE]

*Here $\|.\|_{*}$ is the natural norm of $\|.\|_{\mathcal{L}(L^{2}(\Gamma_{\textup{N}}))}$. *

Proof.

It suffices to prove the theorem for assumption (b) since the other case follows from interchanging $(\lambda_{1},\mu_{1})$ and $(\lambda_{2},\mu_{2})$. For the sake of brevity, we write $\|\cdot\|$ for $\|\cdot\|_{L^{2}(\Gamma_{\textup{N}})^{d}}$. We start with the reformulation of the right-hand side of estimate (15). Since $\Lambda(\mathbb{C}_{1})$ and $\Lambda(\mathbb{C}_{2})$ are self-adjoint, and assumption (b) implies $\Lambda(\mathbb{C}_{2})\geq\Lambda(\mathbb{C}_{1})$ by Corollary 3.2, we have that

[TABLE]

Next, we apply the second inequality in the monotonicity relation (5) in Lemma 3.1 and thus obtain for all $g\in L^{2}(\Gamma_{\textup{N}})^{d}$

[TABLE]

where $u_{\mathbb{C}_{1}}^{g}=u_{(\lambda_{1},\mu_{1})}^{g}\in\mathcal{V}$ denotes the solution of (1) with Neumann data $g$ and elastic tensor $\mathbb{C}_{1}:=\mathbb{C}(\lambda_{1},\mu_{1})$.

The estimate (3) contains the linear differences $\lambda_{2}-\lambda_{1}$ and $\mu_{2}-\mu_{1}$, but it also contains the solution $u_{(\lambda_{1},\mu_{1})}^{g}$ that depends non-linearly on the coefficients. Following the ideas of [43, 42], we separate these dependencies by introducing

[TABLE]

Thus, we obtain for $\mathbb{C}_{1}\neq\mathbb{C}_{2}$,

[TABLE]

By assumption (b), and the definition of $d_{\Omega}(\mathbb{C}_{1},\mathbb{C}_{2})$, it follows that the second argument of $\Psi$ in (17) belongs to the compact set

[TABLE]

Hence, (17) yields that

[TABLE]

The assertion of Theorem 3.5 follows if we can show that the right-hand side of (18) is strictly positive. Since $\Psi$ is continuous, we can conclude that the function

[TABLE]

is semi-lower continuous, so that it attains its minimum on the compact set

$\mathcal{C}\times\mathcal{F}_{[a,b]\times[c,d]}$. Hence, to prove Theorem 3.5, it suffices to show that

[TABLE]

In order to prove that (19) holds true, let $\left((\zeta_{1},\zeta_{2}),(\kappa,\tau)\right)\in\mathcal{C}\times\mathcal{F}_{[a,b]\times[c,d]}$. Then there exist an open subset $\emptyset\neq D_{1}\subset\Omega$ and a constant $0<\delta<1$, such that either

[TABLE]

We use the localized potentials sequence in Theorem 3.3 to obtain an open subset $D_{2}\subset\Omega$ with $\overline{D_{1}}\cap\overline{D_{2}}=\emptyset$, and a boundary load $\tilde{g}\in L^{2}(\Gamma_{\textup{N}})^{d}$ with

[TABLE]

In case (i), this leads to

[TABLE]

and in case (ii), we have

[TABLE]

Hence, in both cases,

[TABLE]

so that Theorem 3.5 is proven. ∎

Remark 1*.*

As a consequence of Theorem 3.5, we end up with the following uniqueness result for all $\mathbb{C}_{1}:=\mathbb{C}(\lambda_{1},\mu_{1}),\mathbb{C}_{2}:=\mathbb{C}(\lambda_{2},\mu_{2})$ from the finite dimensional space $\mathcal{E}$, provided that either, $\lambda_{1}-\lambda_{2}\leq 0$ and $\mu_{1}-\mu_{2}\leq 0$, or $\lambda_{1}-\lambda_{2}\geq 0$ and $\mu_{1}-\mu_{2}\geq 0$:

[TABLE]

Remark 2*.*

All of the results in this section stay valid (with the same proofs) when the Neumann-to-Dirichlet operator $\Lambda(\mathbb{C})$ is extended to the spaces $H^{-\frac{1}{2}}(\Gamma_{\textup{N}})\to H^{\frac{1}{2}}(\Gamma_{\textup{N}})$, where $H^{\frac{1}{2}}(\Gamma_{\textup{N}})=\{u|_{\Gamma_{\textup{N}}}:u\in\mathcal{V}\}$ and $H^{-\frac{1}{2}}(\Gamma_{\textup{N}})$ is its dual. On these spaces, it is easily shown that $\Lambda(\mathbb{C})$ is bijective, and its inverse is the Dirichlet-to-Neumann operator $\Lambda_{\textup{D}}:f\mapsto u_{\textup{D}}|_{\Gamma_{\textup{N}}}$, where $u_{\textup{D}}$ solves

[TABLE]

4 Numerical approach to solve the inverse problem

In this section, we consider Lamé parameters $(\lambda,\mu)\in\tilde{\mathcal{F}}$, where $\tilde{\mathcal{F}}$ is a finite dimensional subset of $L^{\infty}(\Omega)\times L^{\infty}(\Omega)$-functions with positive minima.

We first take a look at the inverse problem

[TABLE]

where $f_{k}\in L^{2}(\Gamma_{\textup{N}})^{d}$ is a measurement of the displacement corresponding to the input surface load $g_{k}$, and $K\in\mathbb{N}$ is the number of measurements.

To solve the inverse problem (22) numerically, we first consider a single measurement $(f,g)$ and define a minimization problem of Kohn-Vogelius type:

[TABLE]

Here $u_{\textup{N}}$ and $u_{\textup{D}}$ solve the following problems:

[TABLE]

[TABLE]

Theorem 4.1**.**

The functional $J:L^{\infty}_{+}(\Omega)\times L^{\infty}_{+}(\Omega)\to\mathbb{R}$, defined in (23) is Fréchet differentiable, and its Fréchet derivative at $(\lambda,\mu)\in L^{\infty}_{+}(\Omega)\times L^{\infty}_{+}(\Omega)$ in the direction $(\tilde{\lambda},\tilde{\mu})\in L^{\infty}_{+}(\Omega)\times L^{\infty}_{+}(\Omega)$ is given by

[TABLE]

Proof.

From the definition of the functional $J$, and applying Green’s formula once, we have

[TABLE]

Since $\Lambda(\mathbb{C}(\lambda,\mu))$ is Fréchet differentiable with

[TABLE]

(c.f., e.g. [10] for a proof that uses only the variational formulation) and

[TABLE]

implies

[TABLE]

(see Remark 2), the assertion follows.

∎

In the next section, we treat the case of several measurements $(f_{1},g_{1}),\ldots,(f_{K},g_{K})$ by adding the respective functionals of the form (23), together with a regularization term.

5 Implementation details and numerical examples

As stated in the introduction, our aim is to obtain first results for the solution of the inverse problem of elasticity. Thus, we examine the intuitive 2 dimensional setting as a first descriptive and meaningful example. In doing so, we consider the following setup for our numerical example: The domain $\Omega$ under consideration is the two dimensional unit disk centered at the origin. We use a Delaunay triangular mesh and a standard finite element method with piecewise finite elements to numerically compute the states for our problem. The exact data $f$ are computed synthetically by solving the direct problem (1). In the real-world, the data $f$ are experimentally acquired and thus can be corrupted by noise such as arising from quantization errors.

In our numerical examples the simulated noise data are generated using the following formula:

[TABLE]

where $\delta$ is a uniform distributed random variable and $\varepsilon$ indicates the level of noise. For our examples, the random variable $\delta$ is realized using the Matlab function rand$()$. We use the BFGS algorithm [5], from the optimization toolbox of Matlab, to minimize the cost function defined in (23). This quasi-Newton method is well adapted to such problem.

5.1 Numerical examples

For the following numerical examples, we use four measurements

$f_{1},\ldots,f_{4}$ corresponding to the following surface loads:

[TABLE]

We reconstruct the Lamé parameters by minimizing the regularized cost functional

[TABLE]

where $u^{g_{k}}_{\textup{N}}$ and $u^{f_{k}}_{\textup{D}}$ are the solutions to problem (24) and (25) respectively with respect to the boundary load $g_{k}$ and the corresponding measurement data $f_{k}$. Here, $\rho$ is a heuristically chosen regularization parameter. The derivative of (27) is easily obtained from Theorem 4.1.

5.1.1 Example 1

In this example, the Lamé parameters to be recovered are assumed to be constant on $\Omega$ and we consider the minimization problem of recovering two scalar parameters. Let $(\lambda_{i},\mu_{i}){\in\mathbb{R}^{2}}$ denote the initialization, $(\lambda_{e},\mu_{e}){\in\mathbb{R}^{2}}$ the exact parameters to be recovered and $(\lambda_{c},\mu_{c}){\in\mathbb{R}^{2}}$ the computed parameters.

Table 1 summarizes the computational results of the algorithm. Figures 1-3, show the decrease of the cost function $J$ and the $L^{\infty}$-norm of $J^{\prime}$ in the course of the optimization process. The numerical solution represents a good approximation and it is stable with respect to small amounts of noise.

5.1.2 Example 2

In this example, the exact Lamé parameters (Figure 4) to be reconstructed are given by

[TABLE]

We reconstruct $\mu$ and $\lambda$ by minimizing the functional (27) in the space of piecewise constant functions on the FEM mesh. The choice $(\lambda_{0},\mu_{0})=(0.3,0.5)$ was made as an initial guess.

The resulting reconstructions and the absolute errors are depicted in Figures 5-6 and Figures 7-8. The history of cost function $J$ and the $L^{\infty}$-norm of $J^{\prime}$ in the course of the optimization process are depicted in Figures 9-10. The parameter $\lambda_{e}$ is well reconstructed, while the parameter $\mu_{e}$ is less well reconstructed.

5.1.3 Example 3

In this example, the exact Lamé parameters (Figure 11) to be recovered are given by

[TABLE]

and

[TABLE]

As in Example 2, we reconstruct $\mu$ and $\lambda$ by minimizing the functional (27) in the space of piecewise constant functions on the FEM mesh, and use the initial guess $(\lambda_{0},\mu_{0})=(0.3,0.5)$. The resulting reconstructions and the absolute errors are depicted in Figures 12-13 and Figures 14-15. The history of cost function $J$ and the $L^{\infty}$-norm of $J^{\prime}$ in the course of the optimization process are depicted in Figures 16-17. The parameter $\mu_{e}$ is well reconstructed, while the parameter $\lambda_{e}$ is less well reconstructed. However, the location of the peaks are obtained.

Remark 3*.*

The accuracy of the numerical results presented in this paper depend on the initialization of the proposed algorithm and the choice of the regularization parameter. Meta-heuristic algorithms can be adapted to select a good initialization and avoid being trapped in a local minima. The optimal choice of the regularization parameter is based on some selection methods (such as the discrepancy principle, the generalized cross validation, the L-curve criterion) which is beyond the scope of this paper.

6 Conclusion and Outlook

In this paper, we dealt with the identification of Lamé parameters in linear elasticity. We introduced the inverse problem and the corresponding Neumann-to-Dirichlet operator. Based on this, we analyzed the connection between the Lamé parameters and the Neumann-to-Dirichlet operator which led to a monotonicity result. In order to prove a Lipschitz stability estimate for Lamé parameters which belong to a known finite subspace with a priori known bounds as well as certain regularity and monotonicity properties, we applied the monotonicity result combined with the localized potentials. The numerical solution of the inverse problem itself, was obtained via the minimization of a Kohn-Vogelius-type cost functional. In more detail, the reconstruction was performed via an iterative algorithm based on a quasi-Newton method. Finally, we presented our numerical examples and discussed them. We want to remark that the monotonicity properties of the Neumann-to-Drichlet operator as well as the results of the localized potentials build the basis for the monotonicty methods for linear elasticity (see [53]) .