733
Views
14
CrossRef citations to date
0
Altmetric
Research Article

Lipschitz stability estimate and reconstruction of Lamé parameters in linear elasticity

ORCID Icon, ORCID Icon, &
Pages 396-417 | Received 05 Jun 2019, Accepted 01 Jul 2020, Published online: 27 Jul 2020

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é parameters with certain regularity assumptions. In addition, we assume that the Lamé parameters belong to a known finite subspace with a priori known bounds and that they fulfil 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échet 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.

2010 Mathematics Subject Classifications:

1. Introduction

In this paper, we consider the inverse problem of recovering the elastic tensor C of a linear isotropic elastic body from the Neumann-to-Dirichlet operator Λ(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 [Citation1] and the references therein), geophysical (see, e.g. [Citation2]), and medical applications (as considered in [Citation3]), in particular localization of potential tumours 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 [Citation4].

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 C1 and C2, depending on the Lamé parameters (λ1,μ1) and (λ2,μ2), respectively. The estimate is given under the following a priori assumptions.

  1. The Lamé parameters (λj,μj) for j = 1, 2 belong to a known finite-dimensional subspace, with λ being piecewise continuous and μ being Lipschitz continuous, and with fixed uniform lower and upper bounds.

  2. The two pairs of parameters satisfy a monotonicity assumption that (λ1,μ1) is, as a pair, either a lower or an upper bound to (λ2,μ2).

Part two deals with the reconstruction of the parameters (λ,μ) 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 C (or the Lamé moduli λ,μ; cf. (Equation2)) has been studied by several authors. In the two dimensional case Ikehata [Citation5] proves that the deflection h between (λ+h,μ+h) and (λ,μ) can be uniquely determined by the first-order approximation of the Dirichlet-to-Neumann operator. Akamatsu et al. [Citation6] give an inversion formula for the normal derivatives at the boundary of the Lamé coefficients λ,μC from the Dirichlet-to-Neumann map. At the same time they present stability estimates for the boundary values of λ,μ. Nakamura and Uhlmann [Citation7] 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 [Citation8] proved that the Lamé coefficient λ can be recovered from partial Cauchy data if the coefficient μ is some positive constant. Boundary determination of Lamé coefficients can be found in the work by Lin and Nakamura [Citation9]. In addition, we want to mention a global uniqueness results in 2D by Imanuvilov and Yamamoto [Citation10].

For the three dimensional case, Nakamura and Uhlmann [Citation11,Citation12] and Eskin and Ralston [Citation13] proved uniqueness results for both Lamé coefficients when μ 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 [Citation14,Citation15], and some boundary determination results were shown in [Citation9,Citation11,Citation16]. For fully anisotropic C, uniqueness was proved Caˆrstea et al. [Citation17] for a piecewise homogeneous medium. Isakov et al. [Citation18] 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. Further on, we want to mention [Citation19], where Hölder stability for non-linear inverse problems is considered.

In this paper, we prove a Lipschitz stability result when the Lamé coefficient λ is piecewise continuous, μ 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 [Citation20–26]. 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 [Citation27].

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 [Citation11,Citation12,Citation28,Citation29]. 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 [Citation30] as well as for the inverse Robin transmission problem [Citation31,Citation32] and on the stability in machine learning reconstruction algorithms [Citation33] under a definiteness assumption. In addition, we want to mention the works [Citation34–38], [Citation39–41] which are concerned with monotonicity-based methods for EIT as well as [Citation42] for monotonicity in inverse medium scattering and [Citation43] for monotonicity methods for the p-Laplace equation.

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 ΩRd (d2), be a bounded and connected open set, occupied by an isotropic material with linear stress-strain relation. The boundary Ω, is assumed to be C1,1 and consists of two non-empty disjoint open parts, the fixed ‘Dirichlet-boundary’ ΓD and the ‘Neumann-boundary’ ΓN with Ω=ΓN¯ΓD¯,ΓNΓD=. 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 gL2(ΓN)d. Then the displacement vector u:Ω¯Rd satisfies the following boundary value problem (1) div(Cˆu)=0in Ω,(Cˆu)ν=gon ΓN,u=0on ΓD,(1) where ν is the outer unit normal vector to Ω, which is similar to the boundary value problem considered, e.g. in [Citation44]. The linearized strain tensor ˆu and the stress tensor Cˆu are given by ˆu=12u+(u)T,Cˆu=k,l=1dCijklukxl1i,jd, where u(x)=x1u(x),,xdu(x),divF(x)=i=1dxiFi1(x)i=1dxiFid(x) for F:ΩRd×d. The isotropic elastic tensor is defined as (2) Cijkl:=λδijδkl+μδikδjl+δilδjk,(2) where λ,μ are the Lamé coefficients.

Next, we state a unique continuation principle from [Citation45] (see also [Citation46] for a recent related result):

Theorem 2.1

Weak Unique Continuation Principle

Let Ω be a connected open set in Rd for d2. Let μ(x)C0,1(Ω) and λ(x)L(Ω) satisfy μ(x)δ0,λ(x)+2μ(x)δ0a.e. xΩ,μC0,1(Ω)+λL(Ω)M0, with positive constants δ0, M0, where we define fC0,1(Ω)=fL(Ω)+fL(Ω). If uH1(Ω)d satisfies div(Cˆu)=0in Ω,u=0 in some ball BΩ, then u = 0 in Ω.

Proof.

The reader is referred to Theorem 1.2 from [Citation45].

Corollary 2.2

Unique Continuation Principle for Local Cauchy Data

Let the same assumptions as in Theorem 2.1 hold and let Ω have a C1,1 boundary and let Γ be a non-empty open subset of Ω. If uH1(Ω)d satisfies div(Cˆu)=0in Ω,(Cˆu)ν=0on Γ,u=0on Γ, then u = 0 in Ω.

Proof.

A solution with vanishing Cauchy data in Ω can be extended by zero to a solution in an extended open set Ω~, where ΩΩ~. Hence, the weak UCP (Theorem 2.1) can be applied to show that u = 0 in Ω.

Next, for given constants α1,α2,β1,β2 satisfying 0<α1α2, 0<β1β2, we define the set of admissible elastic tensors by A=C=C(λ,μ): (λ,μ)L(Ω)×C0,1(Ω),α1λα2,β1μβ2. Hence, the Lamé parameters of every C(λ,μ)A satisfy the conditions of the unique continuation principles Theorem 2.1 and Corollary 2.2.

In what follows, we denote A:B=i,j=1daijbij, for matrices A=(aij) and B=(bij).

The weak formulation of problem (Equation1) is given by (3) ΩCˆu:ˆvdx=ΓNgvdsfor all vV,(3) where V:=vH1(Ω)d:v|ΓD=0. It is easy to see that for each CA, problem (Equation3) has a unique solution uV, which follows by the Lax-Milgram theorem and is shown, e.g. in [Citation47].

We introduce the Neumann-to-Dirichlet operator Λ(C): Λ(C):L2(ΓN)dL2(ΓN)d:guCg|ΓN, where here and in the following uCgV always denotes the unique solution of (Equation1) with surface load gL2(ΓN)d and elastic tensor CA.

It is well known (and an easy consequence from the variational formulation (Equation3) and the compactness of the trace operator) that Λ(C) is a self-adjoint compact linear operator that fulfils g,Λ(C)h=ΩCˆuCg:ˆuChdxfor all g,hL2(ΓN)d, where , denotes the L2(ΓN)-inner product.

The inverse problem we consider here is the following: (4) Find C or (λ,μ) knowing g,Λ(C)h.(4)

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 C1,C2A be elastic tensors, and gL2(ΓN)d be an applied boundary load. The corresponding solutions of (Equation1) are denoted by u1:=uC1g, u2:=uC2gV. Then (5) Ω(C1C2)ˆu2:ˆu2dxg,Λ(C2)gg,Λ(C1)gΩ(C1C2)ˆu1:ˆu1dx.(5)

Proof.

Since Λ(C2)g=u2|ΓN we can use the variational formulation (Equation3) for C1 and C2 with v:=u2 and obtain ΩC1ˆu1:ˆu2dx=g,Λ(C2)g=ΩC2ˆu2:ˆu2dx. Thus ΩC1ˆ(u1u2):ˆ(u1u2)dx=ΩC1ˆu1:ˆu1dx+ΩC1ˆu2:ˆu2dx2ΩC1ˆu1:ˆu2dx=g,Λ(C1)gg,Λ(C2)g+Ω(C1C2)ˆu2:ˆu2dx. Since the left-hand side is nonnegative, the first asserted inequality follows.

Interchanging C1 and C2, the second inequality follows.

Based on the previous lemma and the definition of C, we are led to the following monotonicity property.

Corollary 3.2

Monotonicity

For C1:=C(λ1,μ1),C2:=C(λ2,μ2)A (6) λ1λ2 and μ1μ2 implies Λ(C1)Λ(C2).(6)

Theorem 3.3

Localized potentials

Let CA and D1,D2 be two open sets with D¯1,D¯2Ω, D¯1D¯2= and let Ω(D¯1D¯2) be connected. Then there exists a sequence (gn)nNL2(ΓN)d, such that the corresponding solutions (ugn)nN of (Equation1) fulfil (7) limnD1divugn2dx=,(7) (8) limnD2divugn2dx=0,(8) (9) limnD1ˆugn:ˆugndx=,(9) (10) limnD2ˆugn:ˆugndx=0.(10)

Proof.

This proof is based on the UCP for local Cauchy data (Corollary 2.2).

First, we define the virtual measurement operators

  1. Aj (j = 1, 2) by Aj:L2(Dj)L2(ΓN)d,Fv|ΓN,where vV solves (11) ΩCˆv:ˆwdx=DjFdivwdxfor all wV,(11)

  2. Bj (j = 1, 2) by Bj:L2(Dj)d×dL2(ΓN)d,Gv|ΓN,where vV solves (12) ΩCˆv:ˆwdx=DjG:ˆwdxfor all wV.(12)

First, we show that the dual operators Aj:L2(ΓN)dL2(Dj),j=1,2,Bj:L2(ΓN)dL2(Dj)d×d,j=1,2, are given by Ajg=div(u)|Dj and Bjg=ˆu|Dj, where u solves problem (Equation1).

To (a):

Let FL2(Ω), gL2(ΓN)d, u, vV solve (Equation1) and (Equation11), respectively. Then, ΩFAjgdx=ΓNgAjFds=ΩCˆv:ˆudx=DjFdiv(u)dx.

To (b):

Let GL2(Ω)d×d, gL2(ΓN)d, u, vV solve (Equation1) and (Equation12), respectively. Then, ΩG:Bjgdx=ΓNgBjGds=ΩCˆv:ˆudx=DjG:ˆudx.

Next, we will prove that (13) R(A1)R(B2)={0}andR(A1){0}.(13)

Let ϕR(A1)R(B2). Then there exist v1,v2V such that v1|ΓN=v2|ΓN=ϕ, and ΩCˆvj:ˆwdx=0 for all wV with supp(w)Ω¯D¯j, j = 1, 2.

Hence, div(Cˆv1)=0in ΩD¯1,div(Cˆv2)=0in ΩD¯2, and (Cˆv1)ν|ΓN=(Cˆv2)ν|ΓN=0. The unique continuation principle for Cauchy data (Corollary 2.2) yields that v1=v2 in Ω(D¯1D¯2). Hence v:=v1χD2+v2χΩD¯2V and satisfies div(Cˆv)=0in Ω,(Cˆv)ν=0on ΓN. It follows that v = 0 and thus ϕ=v|ΓN=0, and consequently R(A1)R(B2)=0.

Next, we take a closer look at the operator A1:L2(D1)L2(ΓN)d, Fv|ΓN. Let O be an open ball with O¯D1 and let vV solve (Equation11) for F=χO. We will show that v|ΓN=A1χO0. We argue by contradiction and assume v|ΓN=0. From (Equation11) we obtain that ΩOCˆv:ˆwdx=0 for all wV with supp(w)Ω¯O¯. Thus, v fulfils div(Cˆv)=0 in ΩO¯ and (Cˆv)ν=0 on ΓN. Now, we apply Corollary 2.2 on ΩO¯ which results in v = 0 on ΩO¯, so that v|OH01(O)d. Also, for all wH01(O)d, we obtain from (Equation11) with F=χO (14) OCˆv:ˆwdx=ΩCˆv:ˆw~dx=Odiv(w)dx=Owνds=0,(14) where w~H01(Ω)d is the zero extension of wH01(O)d. Since (Equation14) is uniquely solvable by the Lax-Milgram-Theorem it follows that v|O=0, so that v = 0 in all of Ω which contradicts (Equation11) with F=χO. Hence, A1χO0 which shows that R(A1){0}.

This proves (Equation13) and this implies R(A1)R(B2). Using [Citation25, Corollary 2.6] it follows that there exists a sequence (gn)nNL2(ΓN)d: limnA1gnL2(Ω)2=limnD1divugn2dx= and limnB2gnL2(Ω)d×d2=limnD2ˆugn:ˆugndx=0, i.e. (Equation7) and (Equation10) hold. Since trˆugn=divugn this also implicates (Equation8) and (Equation9).

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

Definition 3.4

A function fL(Ω) is called piecewise continuous, if there exists a finite decomposition of Ω into non-empty open subsets ΩiΩ, i=1,,n, so that Ωi=1nΩi is a Lebesgue null set, ΩiΩj= (ij), and f|Ω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 F be a finite-dimensional subspace of Cˇ(Ω)×C0,1(Ω), where Cˇ(Ω) is the space of piecewise continuous functions. We consider four constants 0<ab and 0<cd which are the lower and upper bounds of the Lamé parameter and define the sets F[a,b]×[c,d]=(λ,μ)F: aλ(x)b, cμ(x)d for all xΩ,E=C(λ,μ): (λ,μ)F[a,b]×[c,d]. In the following main result of this paper, the domain Ω, the finite-dimensional subspace F and the bounds 0<ab and 0<cd 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 C1:=C(λ1,μ1),C2:=C(λ2,μ2)E with either (a)λ1λ2 and μ1μ2or(b)λ1λ2 and μ1μ2, we have (15) dΩ(C1,C2):=maxλ1λ2L(Ω),μ1μ2L(Ω)CΛ(C1)Λ(C2).(15) Here . is the natural norm of .L(L2(ΓN)).

Proof.

It suffices to prove the theorem for assumption (b) since the other case follows from interchanging (λ1,μ1) and (λ2,μ2). For the sake of brevity, we write for L2(ΓN)d. We start with the reformulation of the right-hand side of estimate (Equation15). Since Λ(C1) and Λ(C2) are self-adjoint, and assumption (b) implies Λ(C2)Λ(C1) by Corollary 3.2, we have that Λ(C2)Λ(C1)=supg=1g,Λ(C2)Λ(C1)g=supg=1g,Λ(C2)Λ(C1)g. Next, we apply the second inequality in the monotonicity relation (Equation5) in Lemma 3.1 and thus obtain for all gL2(ΓN)d (16) g,Λ(C2)Λ(C1)gΩ(C1C2)ˆuC1g:ˆuC1gdx,=Ω(λ1λ2)divu(λ1,μ1)g2dx+2Ω(μ1μ2)ˆu(λ1,μ1)g:ˆu(λ1,μ1)gdx(16) where uC1g=u(λ1,μ1)gV denotes the solution of (Equation1) with Neumann data g and elastic tensor C1:=C(λ1,μ1).

The estimate (Equation16) contains the linear differences λ2λ1 and μ2μ1, but it also contains the solution u(λ1,μ1)g that depends non-linearly on the coefficients. Following the ideas of [Citation30,Citation31], we separate these dependencies by introducing Ψ: L2(ΓN)d×F×F[a,b]×[c,d]RΨg,(ζ1,ζ2),(κ,τ):=Ωζ1divu(κ,τ)g2dx+2Ωζ2ˆu(κ,τ)g:ˆu(κ,τ)gdx. Thus, we obtain for C1C2, (17) Λ(C2)Λ(C1)dΩ(C1,C2)supg=1Ψg,λ1λ2dΩ(C1,C2),μ1μ2dΩ(C1,C2),(λ1,μ1).(17) By assumption (b), and the definition of dΩ(C1,C2), it follows that the second argument of Ψ in (Equation17) belongs to the compact set C:=(ζ1,ζ2)F:ζ1,ζ20 and maxζ1L(Ω),ζ2L(Ω)=1. Hence, (Equation17) yields that (18) Λ(C2)Λ(C1)dΩ(C1,C2)inf(ζ1,ζ2)C,(κ,τ)F[a,b]×[c,d]supg=1Ψg,(ζ1,ζ2),(κ,τ).(18) The assertion of Theorem 3.5 follows if we can show that the right-hand side of (Equation18) is strictly positive. Since Ψ is continuous, we can conclude that the function (ζ1,ζ2),(κ,τ)supg=1Ψg,(ζ1,ζ2),(κ,τ) is semi-lower continuous, so that it attains its minimum on the compact set C×F[a,b]×[c,d]. Hence, to prove Theorem 3.5, it suffices to show that (19) supg=1Ψg,(ζ1,ζ2),(κ,τ)>0for all (ζ1,ζ2)C,(κ,τ)F[a,b]×[c,d].(19) In order to prove that (Equation19) holds true, let ((ζ1,ζ2),(κ,τ))C×F[a,b]×[c,d]. Then there exist an open subset D1Ω and a constant 0<δ<1, such that either (i) ζ1|D1δ, and ζ20,or(ii) ζ2|D1δ, and ζ10.

We use the localized potentials sequence in Theorem 3.3 to obtain an open subset D2Ω with D1¯D2¯=, and a boundary load g~L2(ΓN)d with (20) D1divu(κ,τ)g~2dx1δ and D1ˆu(κ,τ)g~:ˆu(κ,τ)g~dx12δ.(20) In case (i), this leads to Ψg~,(ζ1,ζ2),(κ,τ)Ωζ1divu(κ,τ)g~2dx+2Ωζ2ˆu(κ,τ)g~:ˆu(κ,τ)g~dxD1ζ1divu(κ,τ)g~2dxδD1divu(κ,τ)g~2dx1 and in case (ii), we have Ψg~,(ζ1,ζ2),(κ,τ)Ωζ1divu(κ,τ)g~2dx+2Ωζ2ˆu(κ,τ)g~:ˆu(κ,τ)g~dx2D1ζ2ˆu(κ,τ)g~:ˆu(κ,τ)g~dx2δD1ˆu(κ,τ)g~:ˆu(κ,τ)g~dx1. Hence, in both cases, supg=1Ψ(g,(ζ1,ζ2),(κ,τ))Ψg~g~,(ζ1,ζ2),(κ,τ)=1g~2Ψ(g~,(ζ1,ζ2),(κ,τ))>0, so that Theorem 3.5 is proven.

Remark 3.6

As a consequence of Theorem 3.5, we end up with the following uniqueness result for all C1:=C(λ1,μ1),C2:=C(λ2,μ2) from the finite-dimensional space E, provided that either, λ1λ20 and μ1μ20, or λ1λ20 and μ1μ20: Λ(C1)=Λ(C2)if and only ifC1=C2.

Remark 3.7

All of the results in this section stay valid (with the same proofs) when the Neumann-to-Dirichlet operator Λ(C) is extended to the spaces H1/2(ΓN)H1/2(ΓN), where H1/2(ΓN)={u|ΓN:uV} and H1/2(ΓN) is its dual. On these spaces, it is easily shown that Λ(C) is bijective, and its inverse is the Dirichlet-to-Neumann operator ΛD:fuD|ΓN, where uD solves (21) div(C(λ,μ)ˆuD)=0in Ω,uD=fon ΓN,uD=0on ΓD.(21)

4. Numerical approach to solve the inverse problem

In this section, we consider Lamé parameters (λ,μ)F~, where F~ is a finite-dimensional subset of L(Ω)×L(Ω)-functions with positive minima.

We first take a look at the inverse problem (22) Find C(λ,μ) knowing measurements fk=Λ(C(λ,μ))gk,k=1,K,(22) where fkL2(ΓN)d is a measurement of the displacement corresponding to the input surface load gk, and KN is the number of measurements.

To solve the inverse problem (Equation22) numerically, we proceed similarly to [Citation31], see also [Citation48–50]. We first consider a single measurement (f,g) and define a minimization problem of Kohn-Vogelius type: (23) min(λ,μ)F~J(λ,μ)=ΩC(λ,μ)ˆ(uNuD):ˆ(uNuD)dx.(23) Here uN and uD solve the following problems: (24) div(C(λ,μ)ˆuN)=0in Ω,(C(λ,μ)ˆuN)ν=gon ΓN,uN=0on ΓD,(24) (25) div(C(λ,μ)ˆuD)=0in Ω,uD=fon ΓN,uD=0on ΓD.(25)

Theorem 4.1

The functional J:L+(Ω)×L+(Ω)R, defined in (Equation23) is Fréchet differentiable, and its Fréchet derivative at (λ,μ)L+(Ω)×L+(Ω) in the direction (λ~,μ~)L+(Ω)×L+(Ω) is given by (26) J(λ,μ)(λ~,μ~)=ΩC(λ~,μ~)ˆuD:ˆuDdxΩC(λ~,μ~)ˆuN:ˆuNdx.(26)

Proof.

From the definition of the functional J, and applying Green's formula once, we have J(λ,μ)=ΩC(λ,μ)ˆuN:ˆuNdx+ΩC(λ,μ)ˆuD:ˆuDdx2ΓNgfds=g,Λ(C(λ,μ))g+ΛD(C(λ,μ))f,f2ΓNgfds.

Since Λ(C(λ,μ)) is Fréchet differentiable with g,Λ(C(λ,μ))(λ~,μ~)g=ΩC(λ~,μ~)ˆuN:ˆuNdx (c.f., e.g. [Citation51] for a proof that uses only the variational formulation) and ΛD(C(λ,μ))=Λ(C(λ,μ))1 implies ΛD(C(λ,μ))(λ~,μ~)f,f=Λ(C(λ,μ))1f,Λ(C(λ,μ))(λ~,μ~)Λ(C(λ,μ))1f=ΩC(λ~,μ~)ˆuD:ˆuDdx (see Remark 3.7), the assertion follows.

In the next section, we treat the case of several measurements (f1,g1),,(fK,gK) by adding the respective functionals of the form (Equation23), 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 two dimensional setting as a first descriptive and meaningful example. In doing so, we consider the following setup for our numerical example: The domain Ω under consideration is the two dimensional unit disk centred 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 (Equation1). 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: f~(x1,x2)=f(x1,x2)1+εδon ΓN, where δ is an uniform distributed random variable and ϵ indicates the level of noise. For our examples, the random variable δ is realized using the Matlab function rand(). We use the BFGS algorithm [Citation52], from the optimization toolbox of Matlab, to minimize the cost function defined in (Equation23). This quasi-Newton method is well adapted to such problem.

5.1. Numerical examples

For the following numerical examples, we use four measurements f1,,f4 corresponding to the following surface loads: g1=(0.1,0.1),g2=(0.1,0.2), g3=(0.2,0.1), and g4=(0.3,0.5) on ΓN. We reconstruct the Lamé parameters by minimizing the regularized cost functional (27) J(λ,μ):=k=14ΩC(λ,μ)ˆ(uNgkuDfk):ˆ(uNgkuDfk)dx+ρ2Ωλ2+μ2dx,(27) where uNgk and uDfk are the solutions to problem (Equation24) and (Equation25) respectively with respect to the boundary load gk and the corresponding measurement data fk. Here, ρ is a heuristically chosen regularization parameter. The derivative of (Equation27) 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 Ω and we consider the minimization problem of recovering two scalar parameters. Let (λi,μi)R2 denote the initialization, (λe,μe)R2 the exact parameters to be recovered and (λc,μc)R2 the computed parameters. Table  summarizes the computational results of the algorithm. Figures  show the decrease of the cost function J and the L-norm of J 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.

Figure 1. Simulation results for Example 1: History of the cost function J and the L-norm of J in the case of ε=0.0 and ρ=0.0.

Figure 1. Simulation results for Example 1: History of the cost function J and the L∞-norm of J′ in the case of ε=0.0 and ρ=0.0.

Figure 2. Simulation results for Example 1: History of the cost function J and the L-norm of J in the case of ε=0.03 and ρ=0.00001.

Figure 2. Simulation results for Example 1: History of the cost function J and the L∞-norm of J′ in the case of ε=0.03 and ρ=0.00001.

Figure 3. Simulation results for Example 1: History of the cost function J and the L-norm of J in the case of ε=0.05 and ρ=0.00001.

Figure 3. Simulation results for Example 1: History of the cost function J and the L∞-norm of J′ in the case of ε=0.05 and ρ=0.00001.

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

5.1.2. Example 2

In this example, the exact Lamé parameters (Figure ) to be reconstructed are given by μe(x1,x2)=x12+x22 and λe(x1,x2)=1. We reconstruct μ and λ by minimizing the functional (Equation27) in the space of piecewise constant functions on the FEM mesh. The choice (λ0,μ0)=(0.3,0.5) was made as an initial guess.

Figure 4. Simulation results for Example 2: The exact Lamé parameters.

Figure 4. Simulation results for Example 2: The exact Lamé parameters.

The resulting reconstructions and the absolute errors are depicted in Figures . The history of cost function J and the L-norm of J in the course of the optimization process are depicted in Figures  and . The parameter λe is well reconstructed, while the parameter μe is less well reconstructed.

Figure 5. Simulation results for Example 2: The computed Lamé parameters (ϵ=0.0 and ρ=0.0).

Figure 5. Simulation results for Example 2: The computed Lamé parameters (ϵ=0.0 and ρ=0.0).

Figure 6. Simulation results for Example 2: The computed Lamé parameters (ϵ=0.03 and ρ=0.0001).

Figure 6. Simulation results for Example 2: The computed Lamé parameters (ϵ=0.03 and ρ=0.0001).

Figure 7. Simulation results for Example 2: Error for λ and μ (ϵ=0.0 and ρ=0.0).

Figure 7. Simulation results for Example 2: Error for λ and μ (ϵ=0.0 and ρ=0.0).

Figure 8. Simulation results for Example 2: Error for λ and μ (ϵ=0.03 and ρ=0.0001).

Figure 8. Simulation results for Example 2: Error for λ and μ (ϵ=0.03 and ρ=0.0001).

Figure 9. Simulation results for Example 2: History of the cost function J and the L-norm of J in the ϵ=0.0 and ρ=0.0.

Figure 9. Simulation results for Example 2: History of the cost function J and the L∞-norm of J′ in the ϵ=0.0 and ρ=0.0.

Figure 10. Simulation results for Example 2: History of the cost function J and the L-norm of J in the case of ϵ=0.03 and ρ=0.0001.

Figure 10. Simulation results for Example 2: History of the cost function J and the L∞-norm of J′ in the case of ϵ=0.03 and ρ=0.0001.

5.1.3. Example 3

In this example, the exact Lamé parameters (Figure ) to be recovered are given by μe=x12+x22,

and λe=exp(5((x11/2)2+(x21/2)2))+exp(5((x1+1/2)2+(x2+1/2)2)). As in Example 2, we reconstruct μ and λ by minimizing the functional (Equation27) in the space of piecewise constant functions on the FEM mesh, and use the initial guess (λ0,μ0)=(0.3,0.5). The resulting reconstructions and the absolute errors are depicted in Figures . The history of cost function J and the L-norm of J in the course of the optimization process are depicted in Figures  and . The parameter μe is well reconstructed, while the parameter λe is less well reconstructed. However, the location of the peaks are obtained.

Figure 11. Simulation results for Example 3: The exact Lamé parameters.

Figure 11. Simulation results for Example 3: The exact Lamé parameters.

Figure 12. Simulation results for Example 3: The computed Lamé parameters (ϵ=0.0 and ρ=0.0).

Figure 12. Simulation results for Example 3: The computed Lamé parameters (ϵ=0.0 and ρ=0.0).

Figure 13. Simulation results for Example 3: The computed Lamé parameters (ϵ=0.03 and ρ=0.0001).

Figure 13. Simulation results for Example 3: The computed Lamé parameters (ϵ=0.03 and ρ=0.0001).

Figure 14. Simulation results for Example 3: Error for λ and μ (ϵ=0.0 and ρ=0.0).

Figure 14. Simulation results for Example 3: Error for λ and μ (ϵ=0.0 and ρ=0.0).

Figure 15. Simulation results for Example 3: Error for λ and μ (ϵ=0.03 and ρ=0.0001).

Figure 15. Simulation results for Example 3: Error for λ and μ (ϵ=0.03 and ρ=0.0001).

Figure 16. Simulation results for Example 3: History of the cost function J and the L-norm of J (ϵ=0.0 and ρ=0.0).

Figure 16. Simulation results for Example 3: History of the cost function J and the L∞-norm of J′ (ϵ=0.0 and ρ=0.0).

Figure 17. Simulation results for Example 3: History of the cost function J and the L-norm of J (ϵ=0.03 and ρ=0.0001).

Figure 17. Simulation results for Example 3: History of the cost function J and the L∞-norm of J′ (ϵ=0.03 and ρ=0.0001).

Remark 5.1

The accuracy of the numerical results presented in this paper depends 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 analysed 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 [Citation53]).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Correction Statement

This article has been republished with minor changes. These changes do not impact the academic content of the article.

References

  • Ammari H, Bretin E, Garnier J. Mathematical methods in elasticity imaging. Princeton (NJ): Princeton University Press; 2015.
  • Weglein AB, Araújo FV, Carvalho PM, et al. Inverse scattering series and seismic exploration. Inverse Probl. 2003;19:R27–R83.
  • Barbone E, Gokhale NH. Elastic modulus imaging: on the uniqueness and nonuniqueness of the elastography inverse problem in two dimensions. Inverse Probl. 2004;20:203–96.
  • Bonnet M, Constantinescu A. Topical review: inverse problems in elasticity. Inverse Probl. 2015;21(2):R1–R50.
  • Ikehata M. Inversion formulas for the linearized problem for an inverse boundary value problem in elastic prospection. SIAM J Appl Math. 1990;50(6):1635–1644.
  • Akamatsu M, Nakamura G, Steinberg S. Identification of Lamé coefficients from boundary observations. Inverse Probl. 1991;7(3):335.
  • Nakamura G, Uhlmann G. Identification of Lamé parameters by boundary measurements. Am J Math. 1993;115:1161–1187.
  • Imanuvilov OY, Yamamoto M. On reconstruction of Lamé coefficients from partial Cauchy data. J Inverse Ill-Posed Probl. 2011;19(6):881–891.
  • Lin YH, Nakamura G. Boundary determination of the Lamé moduli for the isotropic elasticity system. Inverse Probl. 2017;33(12):125004.
  • Imanuvilov OY, Yamamoto M. Global uniqueness in inverse boundary value problems for the Navier-Stokes equations and Lamé system in two dimensions. Inverse Probl. 2015;31(3):46.
  • Nakamura G, Uhlmann G. Inverse problems at the boundary for an elastic medium. SIAM J Math Anal. 1995;26(2):263–279.
  • Nakamura G, Uhlmann G. Global uniqueness for an inverse boundary value problem arising in elasticity. Inventiones Math. 2003;152(1):205–207.
  • Eskin G, Ralston J. On the inverse boundary value problem for linear isotropic elasticity. Inverse Probl. 2002;18(3):907.
  • Beretta E, Francini E, Morassi A, et al. Lipschitz continuous dependence of piecewise constant Lamé coefficients from boundary data: the case of non-flat interfaces. Inverse Probl. 2014;30(12):125005.
  • Beretta E, Francini E, Vessella S. Uniqueness and Lipschitz stability for the identification of Lamé parameters from boundary measurements. Inverse Probl Imaging. 2014;8(3):611–644.
  • Nakamura G, Tanuma K, Uhlmann G. Layer stripping for a transversely isotropic elastic medium. SIAM J Appl Math. 1999;59(5):1879–1891.
  • Cârstea CI, Honda N, Nakamura G. Uniqueness in the inverse boundary value problem for piecewise homogeneous anisotropic elasticity. SIAM J Math Anal. 2018;50(3):3291–3302.
  • Isakov V, Wang JN, Yamamoto M. An inverse problem for a dynamical Lamé system with residual stress. SIAM J Math Anal. 2007;39(4):1328–1343.
  • de Hoop MV, Qiu L, Scherzer O. Local analysis of inverse problems: Hölder stability and iterative reconstruction. Inverse Probl. 2012;28(4):045001.
  • Harrach B, Pohjola V, Salo M. Monotonicity and local uniqueness for the Helmholtz equation. Anal. PDE 2019; 12 (7), 1741–1771.
  • Harrach B, Lin YH, Liu H. On localizing and concentrating electromagnetic fields. SIAM J Appl Math. 2018;78(5):2558–2574.
  • Harrach B. Simultaneous determination of the diffusion and absorption coefficient from boundary data. Inverse Probl Imaging. 2012;6(4):663–679.
  • Harrach B, Lin YH. Monotonicity-based inversion of the fractional Schrödinger equation I. positive potentials. SIAM J Math Anal. 2019; 51 (4), 3092–3111.
  • Harrach B, Lin YH. Monotonicity-based inversion of the fractional Schrödinger equation II. general potentials and stability. SIAM J Math Anal. 2020; 52 (1), 402–436.
  • Gebauer B. Localized potentials in electrical impedance tomography. Inverse Probl Imaging. 2008;2(2):251–269.
  • Harrach B. On uniqueness in diffuse optical tomography. Inverse Probl. 2009;25(5):055010.
  • Murea C. The BFGS algorithm for a nonlinear least squares problem arising from blood flow in arteries. Comput Math Appl. 2005;49(2-3):171–186.
  • Bellassoued M, Imanuvilov O, Yamamoto M. Inverse problem of determining the density and two Lamé coefficients by boundary data. SIAM J Math Anal. 2008;40(1):238–265.
  • Bellassoued M, Yamamoto M. Lipschitz stability in determining density and two Lamé coefficients.  Mathematical Analysis and Applications J Math Anal Appl. 2007;329(2):1240–1259.
  • Harrach B. Uniqueness and Lipschitz stability in electrical impedance tomography with finitely many electrodes. Inverse Probl. 2019;35(2):024005.
  • Harrach B, Meftahi H. Global uniqueness and Lipschitz-stability for the inverse Robin transmission problem. SIAM J Appl Math. 2019;79(2):525–550.
  • Harrach B. Uniqueness, stability and global convergence for a discrete inverse elliptic Robin transmission problem. arXiv preprint arXiv:1907.02759; 2020.
  • Seo JK, Kim KC, Jargal A. A learning-based method for solving ill-posed nonlinear inverse problems: a simulation study of lung EIT. SIAM J Imaging Sci. 2019; 12 (3), 1275–1295.
  • Harrach B, Minh MN. Enhancing residual-based techniques with shape reconstruction features in electrical impedance tomography. Inverse Probl. 2016;32(12):125002.
  • Harrach B, Ullrich M. Resolution guarantees in electrical impedance tomography.  Transactions on Medical Imaging IEEE Trans Med Imaging. 2015;34(7):1513–1521.
  • Harrach B, Lee E, Ullrich M. Combining frequency-difference and ultrasound modulated electrical impedance tomography. Inverse Probl. 2015;31(9):095003.
  • Harrach B, Minh MN. Monotonicity-based regularization for phantom experiment data in electrical impedance tomography. In: Hofmann B, Leitão A, Zubelli J (eds): New Trends in Parameter Identification for Mathematical Models, pp 107–120. Trends Math., Birkhäuser, Cham, 2018.
  • Barth A, Harrach B, Hyvönen N, et al. Detecting stochastic inclusions in electrical impedance tomography. Inverse Probl. 2017;33(11):115012.
  • Harrach B, Seo JK. Exact shape-reconstruction by one-step linearization in electrical impedance tomography. SIAM J Math Anal. 2010;42(4):1505–1518.
  • Harrach B, Ullrich M. Monotonicity-based shape reconstruction in electrical impedance tomography. SIAM J Math Anal. 2013;45(6):3382–3403.
  • Harrach B, Ullrich M. Local uniqueness for an inverse boundary value problem with partial data. Proc Am Math Soc. 2017;145(3):1087–1095.
  • Griesmaier R, Harrach B. Monotonicity in inverse medium scattering on unbounded domains. SIAM J Appl Math. 2018;78(5):2533–2557.
  • Brander T, Harrach B, Kar M, et al. Monotonicity and enclosure methods for the p-Laplace equation. SIAM J Appl Math. 2018;78(2):742–758. doi: 10.1137/17M1128599.
  • Ciarlet PG. The finite element method for elliptic problems. Amsterdam: North Holland; 1978.
  • Lin CL, Nakamura G, Uhlmann G, et al. Quantitative strong unique continuation for the Lamé system with less regular coefficients. Methods Appl Anal. 2011;18(1):85–92.
  • Gosh A, Gosh T. Unique continuation for a non bi-Laplacian fourth order elliptic operator. arXiv preprint arXiv:1908.05882; 2019.
  • Ciarlet PG. Mathematical elasticity, volume I: Three-dimensional elasticity. Series “Studies in Mathematics and its Applications”. Amsterdam: North Holland; 1988.
  • Delfour MC, Zolésio JP. Shapes and geometries. 2nd ed. Philadelphia (PA): Society for Industrial and Applied Mathematics (SIAM); 2011. (Advances in Design and Control; Vol. 22). Metrics, analysis, differential calculus, and optimization. doi: 10.1137/1.9780898719826.
  • Ekeland I, Temam R. Analyse convexe et problèmes variationnels. Dunod: Collection Études Mathématiques; 1974.
  • Correa R, Seeger A. Directional derivative of a minimax function. Nonlinear Anal. 1985;9(1):13–22. http://dx.doi.org/10.1016/0362-546X(85)90049-5
  • Lechleiter A, Rieder A. Newton regularizations for impedance tomography: convergence by local injectivity. Inverse Probl. 2008;24(6):065009. 18.
  • Dai Y-H. Convergence properties of the BFGS algoritm. SIAM J. Optim.. 2002;13(3):693–701.
  • Eberle S, Harrach B. Shape reconstruction in linear elasticity: standard and linearized monotonicity method. arXiv preprint arXiv:2003.02598; 2020.

Reprints and Corporate Permissions

Please note: Selecting permissions does not provide access to the full text of the article, please see our help page How do I view content?

To request a reprint or corporate permissions for this article, please click on the relevant link below:

Academic Permissions

Please note: Selecting permissions does not provide access to the full text of the article, please see our help page How do I view content?

Obtain permissions instantly via Rightslink by clicking on the button below:

If you are unable to obtain permissions via Rightslink, please complete and submit this Permissions form. For more information, please visit our Permissions help page.