- Email: [email protected]
An equation error approach for the elasticity imaging inverse problem for predicting tumor location
Computers and Mathematics with Applications 67 (2014) 122–135
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
An equation error approach for the elasticity imaging inverse problem for predicting tumor location E. Crossen a , M.S. Gockenbach b , B. Jadamba a , A.A. Khan a,∗ , B. Winkler a a
Center for Applied and Computational Mathematics, School of Mathematical Sciences, Rochester Institute of Technology, 85 Lomb Memorial Drive, Rochester, NY, 14623, USA b Department of Mathematical Sciences, Michigan Technological University, Houghton, MI, USA
article
info
Article history: Received 1 April 2013 Received in revised form 20 August 2013 Accepted 21 October 2013 Keywords: Inverse problems Parameter identification Tumor identification Equation error approach Regularization Total variation
abstract The primary objective of this work is to study the elasticity imaging inverse problem of identifying cancerous tumors in the human body. This nonlinear inverse problem not only represents an important and interesting application, it also brings forth noteworthy mathematical challenges since the underlying model is a system of elasticity equations involving incompressibility. Due to the locking effect, classical finite element methods are not effective for incompressible elasticity equations. Therefore, special treatment is necessary for both the direct and inverse problems. To study the inverse problem in an optimization framework, we propose an extension of the equation error approach. We focus on two cases, namely when the material parameter is sufficiently smooth and when it is may be discontinuous. For the latter case, we extend the total variation regularization method to the elasticity imaging inverse problem. We give the existence results for the proposed equation error approach and give the convergence analysis for the discretized problem. We give sufficient details on the discrete formulas as well as on the implementation issues. Numerical examples for smooth and discontinuous coefficients are given. © 2013 Elsevier Ltd. All rights reserved.
1. Introduction In living soft tissue, differences in molecular makeup as well as in microscopic and macroscopic structure result in significant differences in tissue stiffness. Moreover, changes in the tissue stiffness generally correlate with pathological changes in the tissue’s state. Many cancers appear as hard nodules within the surrounding softer tissue. Such diseases are often detectable by the standard medical practice of palpation, but palpation remains a subjective technique and is mostly limited to the detection of large, stiff tissue abnormalities that lie near the skin’s surface. In contrast, a tool like ultrasound can be used to diagnose tumors deeper within the body, but it can also fail to find lesions that lack certain acoustical properties. The potential for using varying elastic properties to differentiate between healthy and diseased tissue in a more quantitative manner is clear and has been recognized by many authors. Recently, several researchers have proposed elastic imaging methods that operate on just such principals. The elasticity imaging inverse problem is an innovative approach that makes use of the varying elastic features of healthy and diseased tissue to detect lesions. The basic principle is to exert a small, external, and quasistatic compression force to the tissue and to measure its axial displacement field or, more indirectly, the tissue’s overall motion. From this measurement, a tumor can be identified by solving the inverse problem of determining the tissue’s underlying elasticity. Although several
∗
Corresponding author. E-mail addresses: [email protected] (E. Crossen), [email protected] (M.S. Gockenbach), [email protected] (B. Jadamba), [email protected] (A.A. Khan), [email protected] (B. Winkler). 0898-1221/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.camwa.2013.10.006
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
123
researchers have proposed elastic imaging methods, Raghavan and Yagle [1] were among the first authors to realize that this study can be best done in a framework of an inverse problem of reconstructing elasticity from measured strains and the equations of equilibrium. In their contribution, the equations of equilibrium were solved using finite difference schemes. Kallel and Bertrand [2] used the Newton–Raphson method to vary a finite-element-based model of the elasticity equations to fit, in an output least-squares sense, a set of axial tissue displacement fields estimated using a correlation technique applied to ultrasound signals. The ill-conditioning of the Hessian matrix in this method was eliminated employing the Tikhonov regularization technique. In the same vein, Doyley et al. [3], used a modified Newton–Raphson based iterative scheme for solving the inverse problem of computing the spatial distribution of the Young’s modulus. Computer simulations were conducted to determine the relative merits of reconstructing tissue elasticity using the knowledge of known displacement and stress boundary conditions. Some other important developments can be found in [4–9]. The most popular optimization approach for parameter identification in PDE models is output least-squares (OLS). In the context of the following elliptic problem
− ∇ · ( a∇ u) = f
in Ω ,
(1)
augmented with suitable boundary conditions, the OLS approach consists of minimizing the functional JOLS (a) = ∥u(a) − z ∥2 ,
(2)
where z is the data (the measurement of u), ∥ · ∥ is a suitable norm and u(a) solves the variational problem associated with (1). A drawback of OLS-based approaches is that each iteration of the optimization algorithm requires a solution of the weak form of (1), making these methods relatively expensive. The equation error approach is a general approach in which the unknown parameter is chosen to minimize the residual error. In the context of (1), the equation error approach chooses a to minimize JEE (a) =
1 2
∥∇ · (a∇ z ) + f ∥2H −1 (Ω ) .
(3)
The functional JEE is quadratic in a, so that (after discretization) minimizing JEE reduces to solving a positive (semi)definite linear system. This method has been analyzed by Acar [10], Kärkkäinen [11], Gockenbach and Khan [12], Gockenbach, Jadamba and Khan [12] and recently by Al-Jamal and Gockenbach [13] for (1), and by Gockenbach, Jadamba and Khan [14] for general elliptic inverse problems. See [15–22] for further details. At this juncture, we would like to discuss the so-called adjoint-weighted equation (AWE) approach that has been employed recently by Albocher, Oberai, Barbone and Harari [23] to solve the inverse problem of identifying the shear modulus in incompressible plane-stress elasticity. The AWE approach, introduced by Barbone, Oberai and Harari [24] to study the heat conductivity problem, is a variational framework that is weighted by the adjoint operator and is designed to obtain the direct solution of the inverse problem. The AWE approach has the advantage that it produces a well-posed problem under conditions which are typically milder than the ones required for the well-posedness of the strong form. It is known that under a suitable constraint, the AWE approach is identical to a least-squares formulation of the direct problem. In the above, (3) corresponds to the least-squares formulation for (1), where the dual norm is being used. As compared to the AWE approach, in the equation error approach, a variational formulation of the underlying direct problem is used to define the residual error which is then to be minimized. The uniqueness in the equation error approach is obtained via a regularization whereas an additional constraint leads to the uniqueness in the AWE approach. The AWE for the general linear elasticity problem is given in [25]. Recently, the so-called error in constitutive equation approach in studied by Banerjee et al. [26] to solve an inverse problem in frequency-domain elastodynamics. Pierron et al. [27] presented an extension of the virtual fields method (VFM) to elasto-plastic material identification. The VFM has a similar flavor and builds upon minimizing the error in the constitutive equations. We remark that the equation error approach has two distinct advantages over the OLS approach. Firstly, as mentioned above, it leads to a convex optimization problem and hence it only possesses global minimizers. Secondly, the equation approach is computationally inexpensive as there is no underlying variational problem to be solved. To be more specific, the optimization problem resulting from the equation error approach is an unconstrained optimization problem whereas the optimization problem resulting from the OLS approach is a constrained optimization problem where the PDE is the constraint. (Of course, in both formulations, we must take into account some simple constraints on the parameter.) On the other hand, a deficiency of the equation error approach is that it relies on differentiating the data and hence it is quite sensitive to noise. Furthermore, the equation error approach requires that the data be available on the whole domain. 1.1. Elasticity imaging inverse problem. Associated difficulties and our approach The elasticity imaging inverse problem builds upon a combination of mathematical, computational, and experimental techniques based on linear elasticity models. The underlying mathematical model for this nonlinear inverse problem is the following system of partial differential equations that describes the response of an isotropic elastic object to certain body forces and traction applied to its boundary:
−∇ · σ = f
in Ω ,
(4a)
124
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
σ = 2µϵ(u) + λdiv u I , u = f0 on Γ1 , σ n = f1 on Γ2 .
(4b) (4c) (4d)
Here the domain Ω is a subset of R or R and ∂ Ω = Γ1 ∪ Γ2 is its boundary. In (4), the vector-valued function u = u(x) represents the displacement of the elastic object, f is the applied body force, n is the unit outward normal, and ϵ(u) = 21 (∇ u + ∇ uT ) is the linearized strain tensor. The resulting stress tensor σ and the stress–strain law (4b) are obtained under the assumption that the elastic object is isotropic and the displacement is small enough so that a linear relationship holds. The coefficients µ and λ are the Lamé parameters that quantify the elastic properties of the elastic material. In this setting, the direct problem for (4) is to find the displacement u when functions f0 , f1 , the variable coefficients µ and λ, and the body force f are known. In the elasticity imaging inverse problem, the focus, however, is on the inverse problem of identifying the parameter µ when a certain measurement z of the displacement u is available. As mentioned above, using ultrasound, one can measure displacements in human tissue. Since cancerous tumors differ markedly in their elastic properties from the surrounding healthy tissue, these tumors are then located by solving the inverse problem of identifying the variable parameters that describe the elastic properties of the tissue. We remark that in many engineering applications the corresponding inverse problem is to find both µ and λ which typically vary in a small range (see [28,29] and the cited references therein). The primary objective of this work is to extend the equation error approach to the elasticity imaging inverse problem. The main technical challenge in the study of elasticity imaging inverse problem stems from the fact that human tissue is nearly incompressible, that is, the elasticity modulus λ is significantly large and the standard finite element methods become inefficient for the direct as well as the inverse problem. To have a better understanding of the difficulties associated with a large λ, we introduce some notation. The dot product of two tensors Υ1 and Υ2 will be denoted by Υ1 · Υ2 . Given a sufficiently smooth domain Ω ⊂ R2 , the L2 -norm of a tensor-valued function Υ = Υ (x) is given by 2
∥Υ ∥2L2 = ∥Υ ∥2L2 (Ω ) =
Υ ·Υ =
Ω
3
Ω
2 2 2 2 Υ11 + Υ12 + Υ21 + Υ22 .
On the other hand, for a vector-valued function u(x) = (u1 (x), u2 (x))T , the L2 -norm is given by
∥u∥2L2 = ∥u∥2L2 (Ω ) =
u21 + u22 ,
Ω
1
whereas the H -norm by
∥u∥2H 1 = ∥u∥2H 1 (Ω ) = ∥u∥2L2 + ∥∇ u∥2L2 . In the following discussion, for the sake of simplicity, in (4) we set f0 = 0. For this case, the space of test functions, denoted by V¯ , is given by V¯ = {¯v ∈ H 1 (Ω ) × H 1 (Ω ) : v¯ = 0 on Γ1 }.
(5)
By using the Green’s identity and the boundary conditions (4c) and (4d), we obtain the following weak form of the elasticity system (4): find u¯ ∈ V¯ such that
Ω
2µϵ(¯u) · ϵ(¯v ) +
Ω
λ(div u¯ )(div v¯ ) =
Ω
f v¯ +
Γ2
v¯ f1 ,
for every v¯ ∈ V¯ .
(6)
We define Ψ : V¯ × V¯ → R by
Ψ (¯u, v¯ ) =
Ω
2µϵ(¯u) · ϵ(¯v ) +
Ω
λ(div u¯ )(div v¯ ).
Provided that µ and µ + λ are both away from zero, it can be shown that there are two positive constants ψ1 > 0 and ψ2 > 0 with ψ1 ≤ µ and ψ2 ≥ λ + µ such that ψ1 ∥¯v ∥2¯ ≤ Ψ (¯v , v¯ ) ≤ ψ2 ∥¯v ∥2¯ , for every v¯ ∈ V¯ . V
V
Since λ ≫ µ, the ratio ψ3 = ψ2 /ψ1 is very large. However, since the constant ψ3 determines the error estimates defined by Cea’s lemma, it is natural to suspect that the error estimates could be much larger than the actual approximation error. In fact, this phenomenon is often observed and is one of the reasons that is attributed to the so-called locking effect (see Braess [30]). In the literature, several approaches have been proposed to overcome the locking effect. This includes the least-squares finite element methods (see [31]), the discontinuous Galerkin methods (see [32]), and the mixed finite element methods (see [33]), among others. One of the most popular of these methods is the mixed finite elements, which, in the present context, consists of introducing a pressure term p ∈ Q = L2 (Ω ) by p = λdiv u¯ .
(7)
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
125
The weak formulation of (7) reads:
Ω
(div u¯ )q −
1 Ω
λ
pq = 0,
for every q ∈ Q .
(8)
By using relation (7), the weak form (6) can be expressed as follows: find u¯ ∈ V¯ such that
Ω
2µϵ(¯u) · ϵ(¯v ) +
Ω
p(div v¯ ) =
Ω
f v¯ +
Γ2
v¯ f1 ,
for every v¯ ∈ V¯ ,
(9)
where the pressure p is also an unknown. Therefore, the problem of finding u¯ ∈ V¯ satisfying (6) has now been converted into the problem of finding (¯u, p) ∈ V¯ × Q satisfying (8) and (9). It turns out that (8) and (9) can conveniently be studied in the framework of saddle point problems, as is done in this work. We conclude this section, by briefly discussing an aspect of combining (8) and (9). Let us now define a bilinear form : V × V → R by Ψ
(u, v) = Ψ
Ω
2µϵ(¯u) · ϵ(¯v ) +
Ω
(div v¯ )p −
Ω
(div u¯ )q +
1 Ω
λ
pq,
where u = (¯u, p) and v = (¯v , q). By imposing suitable conditions on the elasticity parameters µ, it can be shown that there exists a constant ψ4 > 0 such that
(u, v) ≤ ψ4 ∥u∥V ∥v∥V , Ψ
for every u, v ∈ V .
Furthermore, assuming uniform positivity of µ and using the Korn’s inequality, we ensure that for every v = (¯v , q) ∈ V , we have
(v, v) = Ψ
Ω
2µϵ(¯v ) · ϵ(¯v ) +
1 Ω
λ
qq ≥ ψ5 (∥¯v ∥2V¯ + ∥q∥2Q )
(10)
is coercive and continuous. Therefore, the Lax–Milgram where ψ5 > 0 is a positive constant. That is, the bilinear map Ψ lemma ensures the existence of a unique solution u = (¯u, p) ∈ V of the saddle point problem (8) and (9). However it should (·, ·), although sufficient to give the unique solvability, is be noted that, unfortunately, the coercivity of the bilinear form Ψ influenced by the value Ψ5 ≈ λ1 , which for a large λ, is very small and serves as an indicator of the numerical instability. In conclusion, although (8) and (9) can be combined into a single variational problem, such a strategy does not alone alleviate the locking effect. In this work, we propose a new equation error approach for the elasticity imaging inverse problem of identifying the elasticity parameter µ in the saddle point problem (8) and (9). The proposed framework leads to a convex optimization problem with no PDEs constraints. Besides giving theoretical details, we also present encouraging numerical results to show the efficiency of the proposed approach. The contents of this paper are organized into five sections. Section 2 introduces the equation error formulation for the identification of the smooth parameter. The existence theorem and convergence results for the discrete problem are given. Discrete formulas are given in Section 3 and numerical results are given in Section 4. The paper concludes with some remarks and open questions. 2. Equation error approach for the elasticity imaging inverse problem We consider the following saddle point problem: given a suitable µ, find (¯u, p) ∈ V¯ × Q such that
Ω
Ω
2µϵ(¯u) · ϵ(¯v ) +
(div u¯ )q −
1 Ω
λ
Ω
p(div v¯ ) =
pq = 0,
Ω
f v¯ +
Γ2
v¯ f1 ,
for every v¯ ∈ V¯ ,
for every q ∈ Q ,
(11a)
(11b)
where Q = L2 (Ω ) and the space V¯ is given by (5). Given all the data, the direct problem in the context of the above saddle point problem is to find (¯u, p). However, our interest is in the inverse problem of finding a parameter µ that makes (11) nearly true for a measurement (¯z , zˆ ) of (¯u, p). For u = (¯u, p) ∈ V and v = (¯v , q) ∈ V , we define the following trilinear form: T (µ, u, v) =
Ω
2µϵ(¯u) · ϵ(¯v ) +
Ω
p(div v¯ ) −
Ω
(div u¯ )q +
1 Ω
λ
pq.
(12)
126
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
By (10), the above trilinear form is coercive. However, since λ1 is very small, for stable numerical approximation, it is necessary to obtain estimates independent of λ. Such estimates are achieved by using the fact that Q and V satisfy the following Babushka–Brezzi condition (see [33]): there exists a constant β > 0 such that
β ≤ inf sup q∈Q
(div v¯ ) q . ∥q∥ ∥¯v ∥
Ω
v¯ ∈V¯
(13)
The above condition is an abstract condition of the angle between the spaces V¯ and Q . Among other things, it provides a guidance for the choice of the finite element spaces used to discretize the saddle point problem. For every (µ, u) ∈ L∞ × V and m(v) := Ω f v¯ + Γ v¯ f1 , we have T (µ, u, ·) − m(·) ∈ V ∗ . We define e(µ, u) to be the 2 preimage of this element under the Riesz map. That is,
⟨e(µ, u), v⟩ =
Ω
2µϵ(¯u) · ϵ(¯v ) +
Ω
p(div v¯ ) −
Ω
(div u¯ )q +
1 Ω
λ
pq − Ω
f v¯ −
Γ2
v¯ f1 .
(14)
We have e(µ, u) = (e1 (µ, u), e2 (µ, u)) ∈ V¯ × Q , where e1 (µ, u) and e2 (µ, u) are defined by
⟨e1 (µ, u), v¯ ⟩ =
⟨e2 (µ, u), q⟩ =
Ω
Ω
2µϵ(¯u) · ϵ(¯v ) +
(div u¯ )q −
1 Ω
λ
Ω
pq,
p(div v¯ ) −
Ω
f v¯ −
Γ2
v¯ f1 ,
for every v¯ ∈ V¯ ,
for every q ∈ Q .
(15)
(16)
2.1. Equation error approach for the identification of smooth coefficient µ To define the equation error approach, we choose B = H 2 (Ω ) to be the parameter space and by A we denote the set of admissible coefficients. We assume that the set A is closed and convex in B. For a given z = (¯z , zˆ ) ∈ V and µ ∈ A, we define the following regularized equation error functional JE (µ) =
1 2
∥e(µ, z )∥2V +
κ 2
∥µ∥2H 2 (Ω ) ,
where κ > 0 is the regularization parameter and ∥µ∥2H 2 (Ω ) is the regularizing term. We now pose the following minimization problem: find µ ∈ A by solving: min JE (µ).
(17)
µ∈A
Notice that if µ = µ∗ and u = u∗ = (¯u∗ , p∗ ) satisfy (11), then e(µ∗ , u∗ ) = 0. The idea of the equation error method is that if z is a measurement of u∗ , then we estimate µ∗ by choosing µ to make the ‘‘residual’’ e(µ, z ) as small as possible in V ∗. We begin with the following existence result. Theorem 2.1. The convex minimization problem (17) is uniquely solvable. Proof. The convexity of the functional JE (·) is a direct consequence of the linearity of e(µ, z ) in µ. Due to the fact that JE (µ) ≥ 0 for every µ ∈ A, we can find a minimizing sequence {µn } ⊂ A such that lim JE (µn ) = inf J (µ).
n→∞
µ∈A
The inequality
κ 2
∥µn ∥2H 2 (Ω ) ≤
1 2
∥e(µn , z )∥2V +
κ 2
∥µn ∥2H 2 (Ω ) ,
confirms that the sequence {µn } is bounded in ∥ · ∥H 2 (Ω ) . By the compact embedding of H 2 (Ω ) into L∞ (Ω ), there exists a subsequence that converges weakly in H 2 (Ω ) and strongly in L∞ (Ω ). By using the same notation for the subsequences as well, we obtain that µn ⇀ µ ˜ ∈ A in H 2 Ω and µn → µ ˜ in L∞ (Ω ). By using the definition of e(·, ·), we have the following two equations
1 ⟨e(µn , z ), v⟩V = 2µn ϵ(¯z ) · ϵ(¯v ) + zˆ (div v¯ ) − (div z¯ )q + zˆ q − f v¯ − v¯ f1 , ∀ v ∈ V , Ω Ω Ω Ω λ Ω Γ2 1 ⟨e(µ, ˜ z ), v⟩V = 2µϵ(¯ ˜ z ) · ϵ(¯v ) + zˆ (div v¯ ) − (div z¯ )q + zˆ q − f v¯ − v¯ f1 , ∀ v ∈ V . Ω Ω Ω Ω λ Ω Γ2
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
127
By subtracting the above two equations, we obtain
⟨e(µn , z ) − e(µ, ˜ z ), v⟩ =
2(µn − µ)ϵ(¯ ˜ z ) · ϵ(¯v ),
Ω
for every v ∈ V .
We set v = e(µn , z ) − e(µ, z ) and perform a simple manipulation to conclude that e(µn , z ) → e(µ, ˜ z ) in V . Finally, using the lower-semicontinuity of the norm ∥ · ∥H 2 (Ω ) , we get JE (µ) ˜ =
1 2
∥e(µ, ˜ z )∥2V +
κ 2
∥µ∥ ˜ 2H 2 (Ω )
κ ∥e(µn , z )∥2V + lim inf ∥µn ∥2H 2 (Ω ) n →∞ 2 2 1 κ 2 = lim inf ∥e(µn , z )∥V + ∥µn ∥2H 2 (Ω ) ≤ lim
1
n→∞
n→∞
2
2
= inf JE (µ), µ∈A
which confirms that µ ˜ ∈ A is a solution of (17) and the proof is complete.
To develop a computational framework for the direct as well as for the inverse problem, discretization is necessary. In this work, we employ the finite element discretization for the direct and the inverse problems. To continue with an abstract treatment of the saddle point problem, we assume that we are given a parameter h converging to 0 and a family V¯ h of finite-dimensional subspaces of V¯ and a family {Qh } of finite-dimensional subspaces of Q . For numerical stability, we will restrict our attention to the choices of V¯ h and Qh for which a discrete analogue of Babushka–Brezzi condition (13) holds. We set Vh = V¯ h × Qh . We define a projection operator Ph : V → Vh , applied componentwise. That is, Ph = (P¯ h , Pˆ h ) with lim ∥P¯ h v¯ − v¯ ∥ = 0,
for every v¯ ∈ V¯
lim ∥Pˆ h q − q∥ = 0,
for every q ∈ Q .
h→0
h→0
(18a) (18b)
Analogously, we assume that Bh is a family of finite-dimensional subspaces of B. We set Ah = Bh A and assume that ∞ h=1 Ah ̸= ∅. We assume that Ah is uniformly bounded. Furthermore, we assume that for every µ ∈ A there exists a sequence {µ ˜ h } with µ ˜ h ∈ Ah such that µ ˜ h ⇀ µ in H 2 (Ω ) and hence µ ˜ h → µ in L∞ (Ω ). For any (µ, uh ) ∈ A × Vh , the element eh (µ, uh ) ∈ Vh is such that for every vh = (¯vh , qh ) ∈ Vh , we have
⟨eh (µh , uh ), vh ⟩V =
Ω
2µh ϵ(¯uh ) · ϵ(¯vh ) +
− Γ2
Ω
(div v¯ h )ph −
Ω
(div u¯ h )qh +
1 Ω
λ
ph qh −
Ω
f v¯ h
v¯ h f1 .
(19)
We consider the following discrete minimization problem: find µh ∈ Ah by solving min JEh (µ) =
µ∈Ah
1 2
∥eh (µ, z )∥2V +
κ 2
∥µ∥2H 2 (Ω ) .
(20)
The following result ensures that the continuous problem can be approximated by its discrete analog. Theorem 2.2. The discrete minimization problem (20) is uniquely solvable. If {µ ˜ h }h>0 is a sequence of minimizers of (20), then there is a subsequence which converges to a minimizer of the continuous problem (17). Proof. The existence of minimizers of (20) can be proved using the same arguments as employed in the proof of Theorem 2.1. Let {µ ˜ h } be a sequence of minimizers of JEh . Then {µ ˜ h } remains bounded in B = H 2 (Ω ) norm. This further ensures the existence of a subsequence, still denoted by {µ ˜ h }, which converges to some µ ˜ ∈ A in the L∞ (Ω ) norm. We claim that eh (µ ˜ h , z ) → e(µ, ˜ z ) in V . In fact, for v = (¯v , q) ∈ V we define (¯vh , qh ) = (P¯h (¯v ), Pˆh (q)). Then, after a rearrangements of terms, we have
⟨eh (µ ˜ h , z ) − e(µ, ˜ z ), v⟩ =
Ω
2(µh − µ)ϵ(¯z ) · ϵ(¯vh ) +
− Ω
(div (¯z ))(qh − q) +
Ω
Ω
2µϵ(¯z ) · ϵ(¯vh − v¯ ) +
zˆ (qh − q) +
Ω
(div (¯vh − v¯ ))ˆz f (¯vh − v¯ ) + (¯vh − v¯ )f1 , Ω
Γ2
and by appropriately choosing v , the above expression confirms that eh (µ ˜ h , z ) → e(µ, ˜ z ) as h → 0.
128
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
Let µ ∈ A be arbitrary. Then, there exists a sequence {µ ˆ h } with µ ˆ h ∈ Ah such that µ ˆ h → µ in L∞ (Ω ). Then, J (µ) ˜ =
1 2
∥e(µ, ˜ z )∥2H 2 (Ω ) +
κ 2
∥µ∥ ˜ 2H 2 (Ω )
1 κ ≤ lim inf ∥eh (µ ˜ h , z )∥2H 2 (Ω ) + lim inf ∥µ ˜ h ∥2H 2 (Ω ) h→0 2 h→0 2 1 κ 2 2 ≤ lim inf ∥eh (µ ˜ h , z )∥H 2 (Ω ) + ∥µ ˜ h ∥H 2 ( Ω ) h→0 2 2 1 κ ≤ lim inf ∥eh (µ ˆ h , z )∥2H 2 (Ω ) + ∥µ ˆ h ∥2H 2 (Ω ) 2
h→0
=
1 2
∥e(µ, z )∥2H 2 (Ω ) +
2
κ 2
∥µ∥2H 2 (Ω ) .
Since µ ∈ A was chosen arbitrarily, we have shown that µ ˜ ∈ A is a minimizer.
2.2. Equation error approach for the identification of the discontinuous coefficient µ We now focus on the case when the coefficient µ is either discontinuous or sharply-varying. In recent years, the total variation regularization has been used extensively to reconstruct sharply-varying images (see [34–36]) and has been applied to identify discontinuous coefficients in elliptic inverse problems (see [37,12,14]). Before any continuation, we recall that the total variation of f ∈ L1 (Ω ) is defined by TV (f ) = sup
f (divg ) : g ∈ C01 (Ω )
Ω
2
, |g (x)| ≤ 1 for every x ∈ Ω ,
where | · | represents the Euclidean norm of a vector. If f ∈ L1 (Ω ) satisfies TV (f ) < ∞, then f is said to have bounded variation, and the space of functions of bounded variations BV (Ω ) is defined by BV (Ω ) = f ∈ L1 (Ω ) : TV (f ) < ∞ ,
which is equipped with the norm ∥f ∥BV (Ω ) = ∥f ∥L1 (Ω ) + TV (f ). The functional TV (·) is the BV-seminorm on the space BV (Ω ). It is known that BV (Ω ) is compactly embedded in L1 (Ω ). We define the set of admissible parameters:
A = µ ∈ L1 (Ω )|TV (µ) < ∞ and α1 ≤ µ(x) ≤ α2 almost everywhere in Ω ,
(21)
where α1 and α2 are two positive constants. For sufficiently smooth data z = (¯z , zˆ ) and µ ∈ A, we define the following regularized equation error functional
JE (µ) =
1 2
∥e(µ, z )∥2V + κ TV (µ),
where κ > 0 is the regularization parameter and TV (µ) is the regularizing term. We consider the following minimization problem: find µ ∈ A by solving min JE (µ).
(22)
µ∈A
The existence of the above problem can be obtained by using the arguments given in Theorem 2.1 and in [12,14]. Let {Th } be a family of triangulations of the domain Ω . Let V˘ h be the space of all continuous linear polynomials relative to Th , and let Vh be the subspace of Vˇ h of functions subject to the constraints that the homogeneous Dirichlet boundary conditions on Γ1 are satisfied. We set V¯ h = Vh × Vh . Let the discrete constraint set Ah of admissible coefficients is given by
Ah = {µh ∈ V˘ h : 0 < α1 ≤ µh (x) ≤ α2 , ∀x ∈ Ω } or equivalently,
Ah = {µh ∈ V˘ h : 0 < α1 ≤ µh (xi ) ≤ α2 , for i = 1, 2, . . . , n}. Finally, we formulate the following discrete minimization problem: find µh ∈ Ah by solving min JhE (µ) =
µ∈Ah
1 2
∥eh (µ, z )∥2V + κϕh (µh ),
where ϕh (·) is given by ϕh (µh ) =
(23)
|∇µh |2 + δ(h), and limh→0 δ(h) = δ(0) = 0. We would like to point out the involvement of the space Vˇ h in the constraint set Ah and that in (23), the TV-regularizer is being replaced by ϕh (·). Note that, the main difficulty associated with the TV regularization is that the closure of the C ∞ (Ω ) Ω
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
129
functions with respect to the BV-norm is the space H 1,1 (Ω ). Therefore, the BV functions, which do not belong to H 1,1 (Ω ) cannot be approximated by a family of smooth functions with respect to the bounded variation norm (see [34]). However, the bounds imposed in the definition of the constraint sets A and its discrete analogue allow a specific construction (see the proof of Theorem 2.3) by involving the space V˘ h . Furthermore, the definition of Ah allows the replacement of the TV regularizer by its smooth analogue as given in (23) above. Therefore, the discontinuous coefficients are being approximated by discrete smooth elements. The following is the main convergence result. Theorem 2.3. Assume that µ∗h is a sequence of solutions to the discrete optimization problem (23). Then µ∗h has a convergent subsequence with limit point as a solution of the minimization problem (22).
Proof. The following proof is based on the arguments given in the proof of Theorem 2.2 above and the techniques developed in [12,14,34,36]. By choosing µ = α1 , we readily obtain that {JhE (µ∗h )} is bounded above by a constant which is independent of h, and hence we deduce that {µ∗h } is a bounded sequence. Recalling that the space BV (Ω ) is compact in L1 (Ω ), we extract a subsequence {µ∗h } converging to µ∗ ∈ A. We choose an arbitrary µ ∈ A. Then, for every constant τ > 0, µτ ∈ C ∞ (Ω ) (see Evans and there exists a function Gariepy [38, p. 127 and p. 172]) satisfying ∥µτ − µ∥L1 (Ω ) ≤ τ and Ω |∇µτ | − TV (µ) ≤ τ . By using the ideas of Zou [36], we define
µτ (x), µ ˆ τ (x) = α1 , α2 ,
if α1 ≤ µτ (x) ≤ α2 ; if µτ (x) < α1 ; if µτ (x) > α2 ,
so that µ ˆ τ ∈ W 1,∞ (Ω ). Furthermore, ∥µ ˆ τ − µ∥L1 (Ω ) ≤ ∥µτ − µ∥L1 (Ω ) ≤ τ and
Ω
|∇ µ ˆ τ| =
{µ ˆ τ =µτ }
|∇ µ ˆτ| =
{µ ˆ τ =µτ }
|∇µτ | ≤
Ω
|∇µτ | ≤ TV (µ) + τ .
We define µ ˆ h = Ih µ ˆ τ , where Ih (·) is the nodal interpolant. Consequently, µ ˆ h ∈ Ah , which further implies that for every h, we have JhE (µ∗h ) ≤ Jh (µ ˆ h ). Notice that µ ˆ h ∈ H 1 (Ω ) which implies that TV (µ ˆ h ) = Ω |∇µh |. By using the convergence properties of the interpolation operator, we have limh→0 ∥µ ˆτ −µ ˆ h ∥W 1,1 (Ω ) = 0, confirming that µ ˆh → µ ˆ τ in L1 (Ω ) and TV (µ ˆ h ) → TV (µ ˆ τ ). Finally, by using the lower-semicontinuity of the BV seminorm, and the above observations, we get
JE (µ∗ ) =
1 2
∥e(µ∗ , z )∥2V + κ TV (µ∗ ) 1
≤ lim ∥e(µ∗h , z )∥2V + κ lim inf TV (µ∗h ) h→0
2
h→0
≤ lim inf JhE (µ∗h ) ≤ lim inf JhE (µ ˆ h ) = Jκ ( µ ˆτ) h→0 h→0 1 ∇ µ ˆτ ˆ τ , z )∥2 + κ = ∥e(µ ≤
2 1 2
V
Ω
∥e(µ ˆ τ , z )∥2V + κ TV (µ) + κτ ,
which, when passed to the limit τ → 0, using the previous results, and noting that µ was chosen arbitrarily, yields
JE (µ∗ ) ≤ JE (µ),
for every µ ∈ A,
ensuring that µ ∈ A solves (22). The proof is complete. ∗
3. Computational framework In this section, we collect general discrete formulae for the proposed equation error formulation. We assume, therefore, that Th is a triangulation on Ω , Lh is the space of all piecewise continuous polynomials of degree dµ relative to Th , U¯ h is the space of all piecewise continuous polynomials of degree du relative to Th , and Qh is the space of all piecewise continuous polynomials of degree dq relative to Th . In order to represent the discrete saddle point problem in a computable form we proceed as follows. We represent bases for Lh , U¯ h and Qh by {ϕ1 , ϕ2 , . . . , ϕm } , {ψ1 , ψ2 , . . . , ψn }, and {χ1 , χ2 , . . . , χk }, respectively. The space Lh is then isomorphic to Rm and for any µ ∈ Lh , we define L ∈ Rm by Li = µ(xi ),
i = 1, 2, . . . , m,
130
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
where the nodal basis {ϕ1 , ϕ2 , . . . , ϕm } corresponds to the nodes {x1 , x2 , . . . , xm }. Conversely, each L ∈ Rm corresponds to µ ∈ Lh defined by
µ=
m
Li ϕi .
i =1
Similarly, u¯ ∈ U¯ h will correspond to U¯ ∈ Rn , where U¯ i = u(yi ), i = 1, 2, . . . , n, and u¯ =
n
U¯ i ψi ,
i =1
where y1 , y2 , . . . , yn are the nodes of the mesh defining U¯ h . Finally, q ∈ Qh will correspond to Q ∈ Rk , where Qi = q(zi ), i = 1, 2, . . . , k, and k
q=
Qi χi ,
i=1
where z1 , z2 , . . . , zk are the nodes of the mesh defining Qh . (The spaces Lh , U¯ h , and Qh are defined relative to the same elements, but the nodes will be different if dµ ̸= du ̸= dq .) Recall that the discrete saddle point problem seeks, for each µh , the unique (¯uh , ph ) ∈ U¯ h × Qh such that
2µh ϵ(¯uh ) · ϵ(¯vh ) +
Ω
Ω
(div u¯ h )qh −
1 Ω
λ
Ω
(div v¯ h )ph =
ph qh = 0,
Ω
f v¯ h +
Γ2
v¯ h f1 ,
for every v¯ ∈ U¯ h ,
for every q ∈ Qh .
(24a)
(24b)
We define S : Rm → Rn+k to be the finite element solution operator that assigns to each coefficient µh ∈ Ah , the unique approximate solution uh = (¯uh , ph ) ∈ U¯ h × Qh . Then S (L) = U, where U is defined by K (L)U = F .
(25)
The stiffness matrix K (L) ∈ R
Kn×n (L)
K ( L) =
BTn×k
(n+k)×(n+k)
n +k
and the load vector F ∈ R
are given by
−Ck×k
Bk×n
with
K (L)i,j = 2µϵ(ψj ) · ϵ(ψi ), i, j = 1, 2, . . . , n, Ω Bi,j = (div ψj )χi , i = 1, 2, . . . , k, n = 1, 2, . . . , n Ω Ci,j = χj χi , i, j = 1, 2, . . . , k, Ω Fi = f ψj + ψj f1 , i = 1, 2, . . . , n, Ω
Γ2
Fj = 0,
j = n + 1, n + 2, . . . , n + k.
For future reference, it will be useful to notice that
K (L)ij = Tijk Lk , where the summation convention is used and T is the tensor defined by
Tijk =
Ω
2ϕk ϵ(ψj ) · ϵ(ψi ),
for every i, j = 1, . . . , n, k = 1, . . . , m.
It is convenient to approximate the components of U¯ h in a single finite element space Uh where U¯ h = Uh × Uh . Therefore, if {ψ1 , . . . , ψl } are the basic of Uh then the vector-valued basis of U¯ h can be constructed by choosing n i i=1
{ψ }
ψ1 ψ2 ψl 0 0 0 = , ,..., , , ,..., 0 0 0 ψ1 ψ2 ψl
we obtain
K¯ K =
O
O . K¯
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
131
In the following U¯ and P are the discrete solutions which get replaced by the data (Z¯ , Zˆ ). From (15) and (16), we have
(K + M )E1 = K (µ)U¯ + BT P − F ¯ MQ E2 = BU − CP and consequently E1 = (K + M )−1 K (µ)U¯ + BT P − F
E2 = MQ−1 BU¯ − CP .
Using the above formulas, we have JE (µ) =
= =
1
⟨E1 , (M + K )E1 ⟩ +
2 1
1 2
E2 , MQ E2
1 K (µ)U¯ + BT P − F , (K + M )−1 K (µ)U¯ + BT P − F + BU¯ − CP , MQ−1 BU¯ − CP
2 1 2
2
1 L(U¯ )µ + BT P − F , (K + M )−1 L(U¯ )µ + BT P − F + BU¯ − CP , MQ−1 BU¯ − CP . 2
The derivative is given by DJE (µ)(δµ) = L(U¯ )δµ, (K + M )−1 L(U¯ )µ + BT P − F
and consequently,
∇ JE = L(U¯ )T (K + M )−1 L(U¯ )µ + BT P − F . 4. Numerical results In this section, we use the proposed equation error approach to identify the variable coefficient µ(x, y) in two representative examples of mixed boundary value problems in linear elasticity. In all examples, we consider an isotropic membrane of unit square dimensions, that is, Ω = (0, 1) × (0, 1), with boundary ∂ Ω = Γ1 ∪ Γ2 , where Γ1 and Γ2 are the sub-boundaries with Dirichlet and Neumann conditions, respectively. Since the main focus of this paper is the recovery of parameters in nearly incompressible materials, λ is taken as a large constant, typically λ = 106 . The examples are of a synthetic nature and the examples use the data from computed solutions using mixed finite element methods rather than from either a direct measurement or an analytic solution. In addition, although the method for solving the general saddle point problem also provides values for the pressure term Z , this data is generally unavailable in practical situations. Thus some other method must be employed to obtain Z data. Assuming that Z¯ is relatively close to the actual solution U¯ , Z can be estimated using Z = C −1 BZ¯ . This estimation, rather than an arbitrary value, was used throughout the numerical experiments. It should be noted that no systematic methods of optimal regularization parameter selection were used due to the lack of such theoretically-grounded automatic methods for problems of this type. As a result, the choice of regularization parameter was largely heuristic. 4.1. Example 1 In this first example, we consider a smooth µ (x, y) on a uniform 62 × 62 quadrangular mesh with 7688 total solution degrees of freedom. The top and bottom of the membrane (Γ1 ) are fixed with constant Dirichlet condition g (x, y) and the left and right sides (Γ2 ) have applied traction h (x, y). The functions defining the coefficient, load, and boundary conditions are as follows:
−1 1 + 0.1x2 , f (x, y) = , µ (x, y) = 1 − 0.12 cos(3π x2 + y2 ) 0.1(1 + y) 0.1x 0 g ( x, y ) = on Γ1 , h (x, y) = on Γ2 . 0.1 0.5 + y2 The numerical results presented in Fig. 1 for a smooth coefficient on the domain use an appropriate Tikhonov-type regularization scheme. The data was generated by solving the direct problem numerically. The results are encouraging, however, the method appears to be highly sensitive to any noise in Pˆ (data for the pressure term). No random noise was added. Taking this example as representative we also consider the effects of mesh refinement and the choice of regularization parameter on the maximum error. This is outlined in detail in Table 1.
132
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
(a) Exact coefficient.
(b) Early stage of reconstruction.
(c) Progressive reconstruction.
(d) Progressive reconstruction.
(e) Computed µ.
(f) Coefficient error.
Fig. 1. Numerical results for example 4.1 using smooth regularization with parameter κ = 10−6 . Table 1 Comparison of maximum error (∥µcomp − µexact ∥∞ ) for Example 1. Mesh size
κ = 10−4
κ = 10−5
κ = 10−6
κ = 10−7
κ = 10−8
32 62 92 122 152
7.4373E−03 9.4410E−04 2.8512E−04 1.1917E−04 7.6370E−05
7.2089E−03 1.0511E−03 2.7059E−04 8.5532E−05 3.0538E−05
7.1667E−03 1.1360E−03 3.1352E−04 8.3987E−05 2.2674E−05
7.1900E−03 1.1672E−03 3.4203E−04 8.6359E−05 3.0461E−05
7.1818E−03 1.1817E−03 3.5654E−04 9.6161E−05 2.5921E−05
4.2. Example 2 In this example, we consider a discontinuous µ (x, y) on a uniform 62 × 62 quadrangular mesh with 7688 total solution degrees of freedom. The top of the membrane (Γ1 ) is fixed with constant Dirichlet condition g (x, y) and the remaining
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
(a) Exact coefficient.
(b) Early stage of reconstruction.
(c) Progressive reconstruction.
(d) Progressive reconstruction.
(e) Computed µ.
133
(f) Coefficient error.
Fig. 2. Numerical results for example 4.2 using smooth regularization with parameter κ = 10−8 .
boundaries (Γ2 ) have applied traction h (x, y). The functions defining the coefficient, load, and boundary conditions are as follows:
1 + 0.2 sin(π x) 1 + 0.3 sin(π x) µ (x, y) = 1 + 0.1 sin(π x) 1 + 0.5 sin(π x) 1
f ( x, y ) =
1 + 0.1x2 , 0.1(1 + y)
for {(x, y) : for {(x, y) : for {(x, y) : for {(x, y) : otherwise, g (x, y) =
0.2 0.2 0.6 0.6
0 0
≤ x ≤ 0.4, ≤ x ≤ 0.4, ≤ x ≤ 0.8, ≤ x ≤ 0.8, on Γ1 ,
0.2 0.6 0.2 0.6
≤ y ≤ 0.4} ≤ y ≤ 0.8} ≤ y ≤ 0.4} ≤ y ≤ 0.8}
h ( x, y ) =
0.5 + x2 0 10 1
on Γ2 .
For the numerical experiment outlined in Fig. 2 we recover a discontinuous coefficient by using the BV regularization method as described above with a maximum error ∥µcomp − µexact ∥∞ = 0.001368. Again, the data was generated by
134
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
solving the direct problem numerically. In this example too, the method appears to be highly sensitive to any noise in Pˆ (data for the pressure term). No random noise was added. 5. Concluding remarks We have presented a theoretical and numerical treatment of the inverse problem of identifying variable parameters in a saddle point problem. We applied our results to identifying parameters in a system of nearly incompressible linear elasticity with an emphasis on the elasticity imaging inverse problem. We proposed a new equation error formulation and give existence results and a convergence analysis. Our results are quite encouraging. However, it has been noticed that the equation error approach is quite sensitive to the noise, an observation that has been made for the use of equation error for other inverse problems as well. Another limitation of the method is that, due to the mixed formulation, the mixed data consist of not only the displacement but also the pressure. We have given a method of computing data for P from data for U. However, our numerical experiments indicate that for good reconstruction, the data for U should be sufficiently accurate, a noisy measurement has an adverse impact on the data for P (obtained from U). This also signifies the importance of the data smoothing. It should be noticed that our approach requires full interior data (see [25]). In the future, we plan to investigate the effect of noise on parameter recovery and the subsequent issues of data-smoothing and preconditioning. In recent years, preconditioners based on Sobolev gradients have attracted a great deal of attention for the differentiation of such noisy functions and we look to apply these methods here in the context of the tumor identification problem (see [39]). Acknowledgments The authors are grateful to the referees for their careful reading and insightful remarks and suggestions that brought substantial improvement to the paper. The work of A.A. Khan is partially supported by RIT’s COS D-RIG Acceleration Research Funding Program 2012–2013 and a grant from the Simons Foundation (#210443 to Akhtar Khan). The work of B. Jadamba is partially supported by RIT’s COS Faculty Development Grant (FEAD) for 2012. References [1] K.R. Raghavan, A.E. Yagle, Forward and inverse problems in elasticity imaging of soft tissues, IEEE Trans. Nucl. Sci. 41 (1994) 1639–1648. [2] F. Kallel, M. Bertrand, Tissue elasticity reconstruction using linear perturbation method, IEEE Trans. Med. Imaging 15 (1996) 299–313. [3] M.M. Doyley, P.M. Meaney, J.C. Bamber, Evaluation of an iterative reconstruction method for quantitative elastography, Phys. Med. Biol. 45 (2000) 1521–1540. [4] M.A. Aguilo, W. Aquino, J.C. Brigham, M. Fatemi, An inverse problem approach for elasticity imaging through vibroacoustics, IEEE Trans. Med. Imaging 29 (2010) 1012–1021. [5] H. Ammari, P. Garapon, F. Jouve, Separation of scales in elasticity imgaing: a numerical study, J. Comput. Math. 28 (2010) 354–370. [6] A. Arnold, S. Reichling, O. Bruhns, J. Mosler, Efficient computation of the elastography inverse problem by combining variational mesh adaption and clustering technique, Phys. Med. Biol. 55 (2010) 2035–2056. [7] J. Fehrenbach, M. Masmoudi, R. Souchon, P. Trompette, Detection of small inclusions by elastography, Inverse Problems 22 (2006) 1055–1069. [8] J. McLaughlin, J.R. Yoon, Unique identifiability of elastic parameters from time-dependent interior displacement measurement, Inverse Problems 20 (2004) 25–45. [9] A.A. Oberai, N.H. Gokhale, G.R. Feijóo, Solution of inverse problems in elasticity imaging using the adjoint method, Inverse Problems 19 (2003) 297–313. [10] R. Acar, Identification of the coefficient in elliptic equations, SIAM J. Control Optim. 31 (1993) 1221–1244. [11] T. Kärkkäinen, An equation error method to recover diffusion from the distributed observation, Inverse Problems 13 (1997) 1033–1051. [12] M.S. Gockenbach, B. Jadamba, A.A. Khan, Numerical estimation of discontinuous coefficients by the method of equation error, Int. J. Math. Comput. Sci. 1 (2006) 343–359. [13] M.F. Al-Jamal, M.S. Gockenbach, Stability and error estimates for an equation error method for elliptic equations, Inverse Problems 28 (2012) 095006. p. 15. [14] M.S. Gockenbach, B. Jadamba, A.A. Khan, Equation error approach for elliptic inverse problems with an application to the identification of Lamé parameters, Inverse Probl. Sci. Eng. 16 (2008) 349–367. [15] N. Cahill, B. Jadamba, A.A. Khan, M. Sama, B. Winkler, A first-order adjoint and a second-order hybrid method for an energy output least squares elastography inverse problem of identifying tumor location, Bound. Value Probl. (2013) in press. [16] M.M. Doyley, B. Jadamba, A.A. Khan, M. Sama, B. Winkler, A new energy inversion for parameter identification in saddle point problems with an application to the elasticity imaging inverse problem of predicting tumor location, 2013, submitted for publication. [17] R.S. Falk, Error estimates for the numerical identification of a variable coefficient, Math. Comp. 40 (1983) 537–546. [18] T. Harrigan, E.E. Konofagou, Estimation of material elastic moduli in elastography: a local method, and an investigation of Poisson ratio sensitivity, J. Biomech. 37 (2004) 1215–1221. [19] B. Jadamba, A.A. Khan, G. Rus, M. Sama, B. Winkler, A new convex inversion framework for parameter identification in saddle point problems with an application to the elasticity imaging inverse problem of predicting tumor location, 2013, submitted for publication. [20] B. Jadamba, A.A. Khan, M. Sama, Inverse problems of parameter identification in partial differential equations, in: Mathematics in Science and Technology, World Sci. Publ., Hackensack, NJ, 2011, pp. 228–258. [21] T. Kärkkäinen, Error estimates for distributed parameter identification in parabolic problems with output least squares and Crank–Nicolson method, Appl. Math. 42 (1997) 259–277. [22] A.A. Khan, B.D. Rouhani, Iterative regularization for elliptic inverse problems, Comput. Math. Appl. 54 (2007) 850–860. [23] U. Albocher, A.A. Oberai, P.E. Barbone, I. Harari, Adjoint-weighted equation for inverse problems of incompressible plane-stress elasticity, Comput. Methods Appl. Mech. Engrg. 198 (2009) 2412–2420. [24] P.E. Barbone, A.A. Oberai, I. Harari, Adjoint-weighted variational foprmulation for direct solution of inverse heat conduction problem, Inverse Problems 23 (2007) 2325–2342. [25] P.E. Barbone, C.E. Rivas, I. Harari, U. Albocher, A.A. Oberai, Y. Zhang, Adjoint-weighted variational formulation for the direct solution of inverse problems of general linear elasticity with full interior data, Internat. J. Numer. Methods Engrg. 81 (2010) 1713–1736.
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
135
[26] B. Banerjee, T.F. Walsh, W. Aquino, M. Bonnet, Large scale parameter estimation problems in frequency-domain elastodynamics using an error in constitutive equation functional, Comput. Methods Appl. Mech. Engrg. 253 (2013) 60–72. [27] F. Pierron, S. Avril, V.T. Tran, Extension of the virtual fields method to elasto-plastic material identfication with cyclic loads and kinematic hardening, Int. J. Solids Struct. 47 (2010) 2993–3010. [28] M.S. Gockenbach, A.A. Khan, Identification of Lamé parameters in linear elasticity: a fixed point approach, J. Ind. Manag. Optim. 1 (2005) 487–497. [29] B. Jadamba, A.A. Khan, F. Raciti, On the inverse problem of identifying Lamé coefficients in linear elasticity, Comput. Math. Appl. 56 (2008) 431–443. [30] D. Braess, Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics, third ed., Cambridge University Press, 2007. [31] Z. Cai, G. Starke, Least-squares methods for linear elasticity, SIAM J. Numer. Anal. 42 (2004) 826–842. [32] P. Houston, D. Schotzau, T.P. Wihler, An hp-adaptive mixed discontinuous Galerkin FEM for nearly incompressible linear elasticity, Comput. Methods Appl. Mech. Engrg. 195 (2006) 25–28. 3224–3246. [33] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, New York, 1991. [34] M.Z. Nashed, O. Scherzer, Stable approximation of nondifferentiable optimization problems with variational inequalities, in: Recent Developments in Optimization Theory and Nonlinear Analysis (Jerusalem, 1995), in: Contemp. Math., vol. 204, Amer. Math. Soc., Providence, RI, 1997, pp. 155–170. [35] S. Osher, A. Solé, L. Vese, Image decomposition and restoration using total variation minimization and the H 1 norm, Multiscale Model. Simul. 1 (2003) 349–370. [36] J. Zou, Numerical methods for elliptic inverse problems, Int. J. Comput. Math. 70 (1998) 211–232. [37] M.S. Gockenbach, A.A. Khan, An abstract framework for elliptic inverse problems. Part 1: an output least-squares approach, Math. Mech. Solids 12 (2007) 259–276. [38] L. Evans, R. Gariepy, Measure Theory and Fine Properties of Functions, CRC Press, Boca Raton, 1992. [39] P. Kazemi, R. Renka, A Levenberg–Marquardt method based on Sobolev gradients, Nonlinear Anal. 75 (2012) 6170–6179.
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
An equation error approach for the elasticity imaging inverse problem for predicting tumor location E. Crossen a , M.S. Gockenbach b , B. Jadamba a , A.A. Khan a,∗ , B. Winkler a a
Center for Applied and Computational Mathematics, School of Mathematical Sciences, Rochester Institute of Technology, 85 Lomb Memorial Drive, Rochester, NY, 14623, USA b Department of Mathematical Sciences, Michigan Technological University, Houghton, MI, USA
article
info
Article history: Received 1 April 2013 Received in revised form 20 August 2013 Accepted 21 October 2013 Keywords: Inverse problems Parameter identification Tumor identification Equation error approach Regularization Total variation
abstract The primary objective of this work is to study the elasticity imaging inverse problem of identifying cancerous tumors in the human body. This nonlinear inverse problem not only represents an important and interesting application, it also brings forth noteworthy mathematical challenges since the underlying model is a system of elasticity equations involving incompressibility. Due to the locking effect, classical finite element methods are not effective for incompressible elasticity equations. Therefore, special treatment is necessary for both the direct and inverse problems. To study the inverse problem in an optimization framework, we propose an extension of the equation error approach. We focus on two cases, namely when the material parameter is sufficiently smooth and when it is may be discontinuous. For the latter case, we extend the total variation regularization method to the elasticity imaging inverse problem. We give the existence results for the proposed equation error approach and give the convergence analysis for the discretized problem. We give sufficient details on the discrete formulas as well as on the implementation issues. Numerical examples for smooth and discontinuous coefficients are given. © 2013 Elsevier Ltd. All rights reserved.
1. Introduction In living soft tissue, differences in molecular makeup as well as in microscopic and macroscopic structure result in significant differences in tissue stiffness. Moreover, changes in the tissue stiffness generally correlate with pathological changes in the tissue’s state. Many cancers appear as hard nodules within the surrounding softer tissue. Such diseases are often detectable by the standard medical practice of palpation, but palpation remains a subjective technique and is mostly limited to the detection of large, stiff tissue abnormalities that lie near the skin’s surface. In contrast, a tool like ultrasound can be used to diagnose tumors deeper within the body, but it can also fail to find lesions that lack certain acoustical properties. The potential for using varying elastic properties to differentiate between healthy and diseased tissue in a more quantitative manner is clear and has been recognized by many authors. Recently, several researchers have proposed elastic imaging methods that operate on just such principals. The elasticity imaging inverse problem is an innovative approach that makes use of the varying elastic features of healthy and diseased tissue to detect lesions. The basic principle is to exert a small, external, and quasistatic compression force to the tissue and to measure its axial displacement field or, more indirectly, the tissue’s overall motion. From this measurement, a tumor can be identified by solving the inverse problem of determining the tissue’s underlying elasticity. Although several
∗
Corresponding author. E-mail addresses: [email protected] (E. Crossen), [email protected] (M.S. Gockenbach), [email protected] (B. Jadamba), [email protected] (A.A. Khan), [email protected] (B. Winkler). 0898-1221/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.camwa.2013.10.006
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
123
researchers have proposed elastic imaging methods, Raghavan and Yagle [1] were among the first authors to realize that this study can be best done in a framework of an inverse problem of reconstructing elasticity from measured strains and the equations of equilibrium. In their contribution, the equations of equilibrium were solved using finite difference schemes. Kallel and Bertrand [2] used the Newton–Raphson method to vary a finite-element-based model of the elasticity equations to fit, in an output least-squares sense, a set of axial tissue displacement fields estimated using a correlation technique applied to ultrasound signals. The ill-conditioning of the Hessian matrix in this method was eliminated employing the Tikhonov regularization technique. In the same vein, Doyley et al. [3], used a modified Newton–Raphson based iterative scheme for solving the inverse problem of computing the spatial distribution of the Young’s modulus. Computer simulations were conducted to determine the relative merits of reconstructing tissue elasticity using the knowledge of known displacement and stress boundary conditions. Some other important developments can be found in [4–9]. The most popular optimization approach for parameter identification in PDE models is output least-squares (OLS). In the context of the following elliptic problem
− ∇ · ( a∇ u) = f
in Ω ,
(1)
augmented with suitable boundary conditions, the OLS approach consists of minimizing the functional JOLS (a) = ∥u(a) − z ∥2 ,
(2)
where z is the data (the measurement of u), ∥ · ∥ is a suitable norm and u(a) solves the variational problem associated with (1). A drawback of OLS-based approaches is that each iteration of the optimization algorithm requires a solution of the weak form of (1), making these methods relatively expensive. The equation error approach is a general approach in which the unknown parameter is chosen to minimize the residual error. In the context of (1), the equation error approach chooses a to minimize JEE (a) =
1 2
∥∇ · (a∇ z ) + f ∥2H −1 (Ω ) .
(3)
The functional JEE is quadratic in a, so that (after discretization) minimizing JEE reduces to solving a positive (semi)definite linear system. This method has been analyzed by Acar [10], Kärkkäinen [11], Gockenbach and Khan [12], Gockenbach, Jadamba and Khan [12] and recently by Al-Jamal and Gockenbach [13] for (1), and by Gockenbach, Jadamba and Khan [14] for general elliptic inverse problems. See [15–22] for further details. At this juncture, we would like to discuss the so-called adjoint-weighted equation (AWE) approach that has been employed recently by Albocher, Oberai, Barbone and Harari [23] to solve the inverse problem of identifying the shear modulus in incompressible plane-stress elasticity. The AWE approach, introduced by Barbone, Oberai and Harari [24] to study the heat conductivity problem, is a variational framework that is weighted by the adjoint operator and is designed to obtain the direct solution of the inverse problem. The AWE approach has the advantage that it produces a well-posed problem under conditions which are typically milder than the ones required for the well-posedness of the strong form. It is known that under a suitable constraint, the AWE approach is identical to a least-squares formulation of the direct problem. In the above, (3) corresponds to the least-squares formulation for (1), where the dual norm is being used. As compared to the AWE approach, in the equation error approach, a variational formulation of the underlying direct problem is used to define the residual error which is then to be minimized. The uniqueness in the equation error approach is obtained via a regularization whereas an additional constraint leads to the uniqueness in the AWE approach. The AWE for the general linear elasticity problem is given in [25]. Recently, the so-called error in constitutive equation approach in studied by Banerjee et al. [26] to solve an inverse problem in frequency-domain elastodynamics. Pierron et al. [27] presented an extension of the virtual fields method (VFM) to elasto-plastic material identification. The VFM has a similar flavor and builds upon minimizing the error in the constitutive equations. We remark that the equation error approach has two distinct advantages over the OLS approach. Firstly, as mentioned above, it leads to a convex optimization problem and hence it only possesses global minimizers. Secondly, the equation approach is computationally inexpensive as there is no underlying variational problem to be solved. To be more specific, the optimization problem resulting from the equation error approach is an unconstrained optimization problem whereas the optimization problem resulting from the OLS approach is a constrained optimization problem where the PDE is the constraint. (Of course, in both formulations, we must take into account some simple constraints on the parameter.) On the other hand, a deficiency of the equation error approach is that it relies on differentiating the data and hence it is quite sensitive to noise. Furthermore, the equation error approach requires that the data be available on the whole domain. 1.1. Elasticity imaging inverse problem. Associated difficulties and our approach The elasticity imaging inverse problem builds upon a combination of mathematical, computational, and experimental techniques based on linear elasticity models. The underlying mathematical model for this nonlinear inverse problem is the following system of partial differential equations that describes the response of an isotropic elastic object to certain body forces and traction applied to its boundary:
−∇ · σ = f
in Ω ,
(4a)
124
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
σ = 2µϵ(u) + λdiv u I , u = f0 on Γ1 , σ n = f1 on Γ2 .
(4b) (4c) (4d)
Here the domain Ω is a subset of R or R and ∂ Ω = Γ1 ∪ Γ2 is its boundary. In (4), the vector-valued function u = u(x) represents the displacement of the elastic object, f is the applied body force, n is the unit outward normal, and ϵ(u) = 21 (∇ u + ∇ uT ) is the linearized strain tensor. The resulting stress tensor σ and the stress–strain law (4b) are obtained under the assumption that the elastic object is isotropic and the displacement is small enough so that a linear relationship holds. The coefficients µ and λ are the Lamé parameters that quantify the elastic properties of the elastic material. In this setting, the direct problem for (4) is to find the displacement u when functions f0 , f1 , the variable coefficients µ and λ, and the body force f are known. In the elasticity imaging inverse problem, the focus, however, is on the inverse problem of identifying the parameter µ when a certain measurement z of the displacement u is available. As mentioned above, using ultrasound, one can measure displacements in human tissue. Since cancerous tumors differ markedly in their elastic properties from the surrounding healthy tissue, these tumors are then located by solving the inverse problem of identifying the variable parameters that describe the elastic properties of the tissue. We remark that in many engineering applications the corresponding inverse problem is to find both µ and λ which typically vary in a small range (see [28,29] and the cited references therein). The primary objective of this work is to extend the equation error approach to the elasticity imaging inverse problem. The main technical challenge in the study of elasticity imaging inverse problem stems from the fact that human tissue is nearly incompressible, that is, the elasticity modulus λ is significantly large and the standard finite element methods become inefficient for the direct as well as the inverse problem. To have a better understanding of the difficulties associated with a large λ, we introduce some notation. The dot product of two tensors Υ1 and Υ2 will be denoted by Υ1 · Υ2 . Given a sufficiently smooth domain Ω ⊂ R2 , the L2 -norm of a tensor-valued function Υ = Υ (x) is given by 2
∥Υ ∥2L2 = ∥Υ ∥2L2 (Ω ) =
Υ ·Υ =
Ω
3
Ω
2 2 2 2 Υ11 + Υ12 + Υ21 + Υ22 .
On the other hand, for a vector-valued function u(x) = (u1 (x), u2 (x))T , the L2 -norm is given by
∥u∥2L2 = ∥u∥2L2 (Ω ) =
u21 + u22 ,
Ω
1
whereas the H -norm by
∥u∥2H 1 = ∥u∥2H 1 (Ω ) = ∥u∥2L2 + ∥∇ u∥2L2 . In the following discussion, for the sake of simplicity, in (4) we set f0 = 0. For this case, the space of test functions, denoted by V¯ , is given by V¯ = {¯v ∈ H 1 (Ω ) × H 1 (Ω ) : v¯ = 0 on Γ1 }.
(5)
By using the Green’s identity and the boundary conditions (4c) and (4d), we obtain the following weak form of the elasticity system (4): find u¯ ∈ V¯ such that
Ω
2µϵ(¯u) · ϵ(¯v ) +
Ω
λ(div u¯ )(div v¯ ) =
Ω
f v¯ +
Γ2
v¯ f1 ,
for every v¯ ∈ V¯ .
(6)
We define Ψ : V¯ × V¯ → R by
Ψ (¯u, v¯ ) =
Ω
2µϵ(¯u) · ϵ(¯v ) +
Ω
λ(div u¯ )(div v¯ ).
Provided that µ and µ + λ are both away from zero, it can be shown that there are two positive constants ψ1 > 0 and ψ2 > 0 with ψ1 ≤ µ and ψ2 ≥ λ + µ such that ψ1 ∥¯v ∥2¯ ≤ Ψ (¯v , v¯ ) ≤ ψ2 ∥¯v ∥2¯ , for every v¯ ∈ V¯ . V
V
Since λ ≫ µ, the ratio ψ3 = ψ2 /ψ1 is very large. However, since the constant ψ3 determines the error estimates defined by Cea’s lemma, it is natural to suspect that the error estimates could be much larger than the actual approximation error. In fact, this phenomenon is often observed and is one of the reasons that is attributed to the so-called locking effect (see Braess [30]). In the literature, several approaches have been proposed to overcome the locking effect. This includes the least-squares finite element methods (see [31]), the discontinuous Galerkin methods (see [32]), and the mixed finite element methods (see [33]), among others. One of the most popular of these methods is the mixed finite elements, which, in the present context, consists of introducing a pressure term p ∈ Q = L2 (Ω ) by p = λdiv u¯ .
(7)
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
125
The weak formulation of (7) reads:
Ω
(div u¯ )q −
1 Ω
λ
pq = 0,
for every q ∈ Q .
(8)
By using relation (7), the weak form (6) can be expressed as follows: find u¯ ∈ V¯ such that
Ω
2µϵ(¯u) · ϵ(¯v ) +
Ω
p(div v¯ ) =
Ω
f v¯ +
Γ2
v¯ f1 ,
for every v¯ ∈ V¯ ,
(9)
where the pressure p is also an unknown. Therefore, the problem of finding u¯ ∈ V¯ satisfying (6) has now been converted into the problem of finding (¯u, p) ∈ V¯ × Q satisfying (8) and (9). It turns out that (8) and (9) can conveniently be studied in the framework of saddle point problems, as is done in this work. We conclude this section, by briefly discussing an aspect of combining (8) and (9). Let us now define a bilinear form : V × V → R by Ψ
(u, v) = Ψ
Ω
2µϵ(¯u) · ϵ(¯v ) +
Ω
(div v¯ )p −
Ω
(div u¯ )q +
1 Ω
λ
pq,
where u = (¯u, p) and v = (¯v , q). By imposing suitable conditions on the elasticity parameters µ, it can be shown that there exists a constant ψ4 > 0 such that
(u, v) ≤ ψ4 ∥u∥V ∥v∥V , Ψ
for every u, v ∈ V .
Furthermore, assuming uniform positivity of µ and using the Korn’s inequality, we ensure that for every v = (¯v , q) ∈ V , we have
(v, v) = Ψ
Ω
2µϵ(¯v ) · ϵ(¯v ) +
1 Ω
λ
qq ≥ ψ5 (∥¯v ∥2V¯ + ∥q∥2Q )
(10)
is coercive and continuous. Therefore, the Lax–Milgram where ψ5 > 0 is a positive constant. That is, the bilinear map Ψ lemma ensures the existence of a unique solution u = (¯u, p) ∈ V of the saddle point problem (8) and (9). However it should (·, ·), although sufficient to give the unique solvability, is be noted that, unfortunately, the coercivity of the bilinear form Ψ influenced by the value Ψ5 ≈ λ1 , which for a large λ, is very small and serves as an indicator of the numerical instability. In conclusion, although (8) and (9) can be combined into a single variational problem, such a strategy does not alone alleviate the locking effect. In this work, we propose a new equation error approach for the elasticity imaging inverse problem of identifying the elasticity parameter µ in the saddle point problem (8) and (9). The proposed framework leads to a convex optimization problem with no PDEs constraints. Besides giving theoretical details, we also present encouraging numerical results to show the efficiency of the proposed approach. The contents of this paper are organized into five sections. Section 2 introduces the equation error formulation for the identification of the smooth parameter. The existence theorem and convergence results for the discrete problem are given. Discrete formulas are given in Section 3 and numerical results are given in Section 4. The paper concludes with some remarks and open questions. 2. Equation error approach for the elasticity imaging inverse problem We consider the following saddle point problem: given a suitable µ, find (¯u, p) ∈ V¯ × Q such that
Ω
Ω
2µϵ(¯u) · ϵ(¯v ) +
(div u¯ )q −
1 Ω
λ
Ω
p(div v¯ ) =
pq = 0,
Ω
f v¯ +
Γ2
v¯ f1 ,
for every v¯ ∈ V¯ ,
for every q ∈ Q ,
(11a)
(11b)
where Q = L2 (Ω ) and the space V¯ is given by (5). Given all the data, the direct problem in the context of the above saddle point problem is to find (¯u, p). However, our interest is in the inverse problem of finding a parameter µ that makes (11) nearly true for a measurement (¯z , zˆ ) of (¯u, p). For u = (¯u, p) ∈ V and v = (¯v , q) ∈ V , we define the following trilinear form: T (µ, u, v) =
Ω
2µϵ(¯u) · ϵ(¯v ) +
Ω
p(div v¯ ) −
Ω
(div u¯ )q +
1 Ω
λ
pq.
(12)
126
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
By (10), the above trilinear form is coercive. However, since λ1 is very small, for stable numerical approximation, it is necessary to obtain estimates independent of λ. Such estimates are achieved by using the fact that Q and V satisfy the following Babushka–Brezzi condition (see [33]): there exists a constant β > 0 such that
β ≤ inf sup q∈Q
(div v¯ ) q . ∥q∥ ∥¯v ∥
Ω
v¯ ∈V¯
(13)
The above condition is an abstract condition of the angle between the spaces V¯ and Q . Among other things, it provides a guidance for the choice of the finite element spaces used to discretize the saddle point problem. For every (µ, u) ∈ L∞ × V and m(v) := Ω f v¯ + Γ v¯ f1 , we have T (µ, u, ·) − m(·) ∈ V ∗ . We define e(µ, u) to be the 2 preimage of this element under the Riesz map. That is,
⟨e(µ, u), v⟩ =
Ω
2µϵ(¯u) · ϵ(¯v ) +
Ω
p(div v¯ ) −
Ω
(div u¯ )q +
1 Ω
λ
pq − Ω
f v¯ −
Γ2
v¯ f1 .
(14)
We have e(µ, u) = (e1 (µ, u), e2 (µ, u)) ∈ V¯ × Q , where e1 (µ, u) and e2 (µ, u) are defined by
⟨e1 (µ, u), v¯ ⟩ =
⟨e2 (µ, u), q⟩ =
Ω
Ω
2µϵ(¯u) · ϵ(¯v ) +
(div u¯ )q −
1 Ω
λ
Ω
pq,
p(div v¯ ) −
Ω
f v¯ −
Γ2
v¯ f1 ,
for every v¯ ∈ V¯ ,
for every q ∈ Q .
(15)
(16)
2.1. Equation error approach for the identification of smooth coefficient µ To define the equation error approach, we choose B = H 2 (Ω ) to be the parameter space and by A we denote the set of admissible coefficients. We assume that the set A is closed and convex in B. For a given z = (¯z , zˆ ) ∈ V and µ ∈ A, we define the following regularized equation error functional JE (µ) =
1 2
∥e(µ, z )∥2V +
κ 2
∥µ∥2H 2 (Ω ) ,
where κ > 0 is the regularization parameter and ∥µ∥2H 2 (Ω ) is the regularizing term. We now pose the following minimization problem: find µ ∈ A by solving: min JE (µ).
(17)
µ∈A
Notice that if µ = µ∗ and u = u∗ = (¯u∗ , p∗ ) satisfy (11), then e(µ∗ , u∗ ) = 0. The idea of the equation error method is that if z is a measurement of u∗ , then we estimate µ∗ by choosing µ to make the ‘‘residual’’ e(µ, z ) as small as possible in V ∗. We begin with the following existence result. Theorem 2.1. The convex minimization problem (17) is uniquely solvable. Proof. The convexity of the functional JE (·) is a direct consequence of the linearity of e(µ, z ) in µ. Due to the fact that JE (µ) ≥ 0 for every µ ∈ A, we can find a minimizing sequence {µn } ⊂ A such that lim JE (µn ) = inf J (µ).
n→∞
µ∈A
The inequality
κ 2
∥µn ∥2H 2 (Ω ) ≤
1 2
∥e(µn , z )∥2V +
κ 2
∥µn ∥2H 2 (Ω ) ,
confirms that the sequence {µn } is bounded in ∥ · ∥H 2 (Ω ) . By the compact embedding of H 2 (Ω ) into L∞ (Ω ), there exists a subsequence that converges weakly in H 2 (Ω ) and strongly in L∞ (Ω ). By using the same notation for the subsequences as well, we obtain that µn ⇀ µ ˜ ∈ A in H 2 Ω and µn → µ ˜ in L∞ (Ω ). By using the definition of e(·, ·), we have the following two equations
1 ⟨e(µn , z ), v⟩V = 2µn ϵ(¯z ) · ϵ(¯v ) + zˆ (div v¯ ) − (div z¯ )q + zˆ q − f v¯ − v¯ f1 , ∀ v ∈ V , Ω Ω Ω Ω λ Ω Γ2 1 ⟨e(µ, ˜ z ), v⟩V = 2µϵ(¯ ˜ z ) · ϵ(¯v ) + zˆ (div v¯ ) − (div z¯ )q + zˆ q − f v¯ − v¯ f1 , ∀ v ∈ V . Ω Ω Ω Ω λ Ω Γ2
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
127
By subtracting the above two equations, we obtain
⟨e(µn , z ) − e(µ, ˜ z ), v⟩ =
2(µn − µ)ϵ(¯ ˜ z ) · ϵ(¯v ),
Ω
for every v ∈ V .
We set v = e(µn , z ) − e(µ, z ) and perform a simple manipulation to conclude that e(µn , z ) → e(µ, ˜ z ) in V . Finally, using the lower-semicontinuity of the norm ∥ · ∥H 2 (Ω ) , we get JE (µ) ˜ =
1 2
∥e(µ, ˜ z )∥2V +
κ 2
∥µ∥ ˜ 2H 2 (Ω )
κ ∥e(µn , z )∥2V + lim inf ∥µn ∥2H 2 (Ω ) n →∞ 2 2 1 κ 2 = lim inf ∥e(µn , z )∥V + ∥µn ∥2H 2 (Ω ) ≤ lim
1
n→∞
n→∞
2
2
= inf JE (µ), µ∈A
which confirms that µ ˜ ∈ A is a solution of (17) and the proof is complete.
To develop a computational framework for the direct as well as for the inverse problem, discretization is necessary. In this work, we employ the finite element discretization for the direct and the inverse problems. To continue with an abstract treatment of the saddle point problem, we assume that we are given a parameter h converging to 0 and a family V¯ h of finite-dimensional subspaces of V¯ and a family {Qh } of finite-dimensional subspaces of Q . For numerical stability, we will restrict our attention to the choices of V¯ h and Qh for which a discrete analogue of Babushka–Brezzi condition (13) holds. We set Vh = V¯ h × Qh . We define a projection operator Ph : V → Vh , applied componentwise. That is, Ph = (P¯ h , Pˆ h ) with lim ∥P¯ h v¯ − v¯ ∥ = 0,
for every v¯ ∈ V¯
lim ∥Pˆ h q − q∥ = 0,
for every q ∈ Q .
h→0
h→0
(18a) (18b)
Analogously, we assume that Bh is a family of finite-dimensional subspaces of B. We set Ah = Bh A and assume that ∞ h=1 Ah ̸= ∅. We assume that Ah is uniformly bounded. Furthermore, we assume that for every µ ∈ A there exists a sequence {µ ˜ h } with µ ˜ h ∈ Ah such that µ ˜ h ⇀ µ in H 2 (Ω ) and hence µ ˜ h → µ in L∞ (Ω ). For any (µ, uh ) ∈ A × Vh , the element eh (µ, uh ) ∈ Vh is such that for every vh = (¯vh , qh ) ∈ Vh , we have
⟨eh (µh , uh ), vh ⟩V =
Ω
2µh ϵ(¯uh ) · ϵ(¯vh ) +
− Γ2
Ω
(div v¯ h )ph −
Ω
(div u¯ h )qh +
1 Ω
λ
ph qh −
Ω
f v¯ h
v¯ h f1 .
(19)
We consider the following discrete minimization problem: find µh ∈ Ah by solving min JEh (µ) =
µ∈Ah
1 2
∥eh (µ, z )∥2V +
κ 2
∥µ∥2H 2 (Ω ) .
(20)
The following result ensures that the continuous problem can be approximated by its discrete analog. Theorem 2.2. The discrete minimization problem (20) is uniquely solvable. If {µ ˜ h }h>0 is a sequence of minimizers of (20), then there is a subsequence which converges to a minimizer of the continuous problem (17). Proof. The existence of minimizers of (20) can be proved using the same arguments as employed in the proof of Theorem 2.1. Let {µ ˜ h } be a sequence of minimizers of JEh . Then {µ ˜ h } remains bounded in B = H 2 (Ω ) norm. This further ensures the existence of a subsequence, still denoted by {µ ˜ h }, which converges to some µ ˜ ∈ A in the L∞ (Ω ) norm. We claim that eh (µ ˜ h , z ) → e(µ, ˜ z ) in V . In fact, for v = (¯v , q) ∈ V we define (¯vh , qh ) = (P¯h (¯v ), Pˆh (q)). Then, after a rearrangements of terms, we have
⟨eh (µ ˜ h , z ) − e(µ, ˜ z ), v⟩ =
Ω
2(µh − µ)ϵ(¯z ) · ϵ(¯vh ) +
− Ω
(div (¯z ))(qh − q) +
Ω
Ω
2µϵ(¯z ) · ϵ(¯vh − v¯ ) +
zˆ (qh − q) +
Ω
(div (¯vh − v¯ ))ˆz f (¯vh − v¯ ) + (¯vh − v¯ )f1 , Ω
Γ2
and by appropriately choosing v , the above expression confirms that eh (µ ˜ h , z ) → e(µ, ˜ z ) as h → 0.
128
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
Let µ ∈ A be arbitrary. Then, there exists a sequence {µ ˆ h } with µ ˆ h ∈ Ah such that µ ˆ h → µ in L∞ (Ω ). Then, J (µ) ˜ =
1 2
∥e(µ, ˜ z )∥2H 2 (Ω ) +
κ 2
∥µ∥ ˜ 2H 2 (Ω )
1 κ ≤ lim inf ∥eh (µ ˜ h , z )∥2H 2 (Ω ) + lim inf ∥µ ˜ h ∥2H 2 (Ω ) h→0 2 h→0 2 1 κ 2 2 ≤ lim inf ∥eh (µ ˜ h , z )∥H 2 (Ω ) + ∥µ ˜ h ∥H 2 ( Ω ) h→0 2 2 1 κ ≤ lim inf ∥eh (µ ˆ h , z )∥2H 2 (Ω ) + ∥µ ˆ h ∥2H 2 (Ω ) 2
h→0
=
1 2
∥e(µ, z )∥2H 2 (Ω ) +
2
κ 2
∥µ∥2H 2 (Ω ) .
Since µ ∈ A was chosen arbitrarily, we have shown that µ ˜ ∈ A is a minimizer.
2.2. Equation error approach for the identification of the discontinuous coefficient µ We now focus on the case when the coefficient µ is either discontinuous or sharply-varying. In recent years, the total variation regularization has been used extensively to reconstruct sharply-varying images (see [34–36]) and has been applied to identify discontinuous coefficients in elliptic inverse problems (see [37,12,14]). Before any continuation, we recall that the total variation of f ∈ L1 (Ω ) is defined by TV (f ) = sup
f (divg ) : g ∈ C01 (Ω )
Ω
2
, |g (x)| ≤ 1 for every x ∈ Ω ,
where | · | represents the Euclidean norm of a vector. If f ∈ L1 (Ω ) satisfies TV (f ) < ∞, then f is said to have bounded variation, and the space of functions of bounded variations BV (Ω ) is defined by BV (Ω ) = f ∈ L1 (Ω ) : TV (f ) < ∞ ,
which is equipped with the norm ∥f ∥BV (Ω ) = ∥f ∥L1 (Ω ) + TV (f ). The functional TV (·) is the BV-seminorm on the space BV (Ω ). It is known that BV (Ω ) is compactly embedded in L1 (Ω ). We define the set of admissible parameters:
A = µ ∈ L1 (Ω )|TV (µ) < ∞ and α1 ≤ µ(x) ≤ α2 almost everywhere in Ω ,
(21)
where α1 and α2 are two positive constants. For sufficiently smooth data z = (¯z , zˆ ) and µ ∈ A, we define the following regularized equation error functional
JE (µ) =
1 2
∥e(µ, z )∥2V + κ TV (µ),
where κ > 0 is the regularization parameter and TV (µ) is the regularizing term. We consider the following minimization problem: find µ ∈ A by solving min JE (µ).
(22)
µ∈A
The existence of the above problem can be obtained by using the arguments given in Theorem 2.1 and in [12,14]. Let {Th } be a family of triangulations of the domain Ω . Let V˘ h be the space of all continuous linear polynomials relative to Th , and let Vh be the subspace of Vˇ h of functions subject to the constraints that the homogeneous Dirichlet boundary conditions on Γ1 are satisfied. We set V¯ h = Vh × Vh . Let the discrete constraint set Ah of admissible coefficients is given by
Ah = {µh ∈ V˘ h : 0 < α1 ≤ µh (x) ≤ α2 , ∀x ∈ Ω } or equivalently,
Ah = {µh ∈ V˘ h : 0 < α1 ≤ µh (xi ) ≤ α2 , for i = 1, 2, . . . , n}. Finally, we formulate the following discrete minimization problem: find µh ∈ Ah by solving min JhE (µ) =
µ∈Ah
1 2
∥eh (µ, z )∥2V + κϕh (µh ),
where ϕh (·) is given by ϕh (µh ) =
(23)
|∇µh |2 + δ(h), and limh→0 δ(h) = δ(0) = 0. We would like to point out the involvement of the space Vˇ h in the constraint set Ah and that in (23), the TV-regularizer is being replaced by ϕh (·). Note that, the main difficulty associated with the TV regularization is that the closure of the C ∞ (Ω ) Ω
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
129
functions with respect to the BV-norm is the space H 1,1 (Ω ). Therefore, the BV functions, which do not belong to H 1,1 (Ω ) cannot be approximated by a family of smooth functions with respect to the bounded variation norm (see [34]). However, the bounds imposed in the definition of the constraint sets A and its discrete analogue allow a specific construction (see the proof of Theorem 2.3) by involving the space V˘ h . Furthermore, the definition of Ah allows the replacement of the TV regularizer by its smooth analogue as given in (23) above. Therefore, the discontinuous coefficients are being approximated by discrete smooth elements. The following is the main convergence result. Theorem 2.3. Assume that µ∗h is a sequence of solutions to the discrete optimization problem (23). Then µ∗h has a convergent subsequence with limit point as a solution of the minimization problem (22).
Proof. The following proof is based on the arguments given in the proof of Theorem 2.2 above and the techniques developed in [12,14,34,36]. By choosing µ = α1 , we readily obtain that {JhE (µ∗h )} is bounded above by a constant which is independent of h, and hence we deduce that {µ∗h } is a bounded sequence. Recalling that the space BV (Ω ) is compact in L1 (Ω ), we extract a subsequence {µ∗h } converging to µ∗ ∈ A. We choose an arbitrary µ ∈ A. Then, for every constant τ > 0, µτ ∈ C ∞ (Ω ) (see Evans and there exists a function Gariepy [38, p. 127 and p. 172]) satisfying ∥µτ − µ∥L1 (Ω ) ≤ τ and Ω |∇µτ | − TV (µ) ≤ τ . By using the ideas of Zou [36], we define
µτ (x), µ ˆ τ (x) = α1 , α2 ,
if α1 ≤ µτ (x) ≤ α2 ; if µτ (x) < α1 ; if µτ (x) > α2 ,
so that µ ˆ τ ∈ W 1,∞ (Ω ). Furthermore, ∥µ ˆ τ − µ∥L1 (Ω ) ≤ ∥µτ − µ∥L1 (Ω ) ≤ τ and
Ω
|∇ µ ˆ τ| =
{µ ˆ τ =µτ }
|∇ µ ˆτ| =
{µ ˆ τ =µτ }
|∇µτ | ≤
Ω
|∇µτ | ≤ TV (µ) + τ .
We define µ ˆ h = Ih µ ˆ τ , where Ih (·) is the nodal interpolant. Consequently, µ ˆ h ∈ Ah , which further implies that for every h, we have JhE (µ∗h ) ≤ Jh (µ ˆ h ). Notice that µ ˆ h ∈ H 1 (Ω ) which implies that TV (µ ˆ h ) = Ω |∇µh |. By using the convergence properties of the interpolation operator, we have limh→0 ∥µ ˆτ −µ ˆ h ∥W 1,1 (Ω ) = 0, confirming that µ ˆh → µ ˆ τ in L1 (Ω ) and TV (µ ˆ h ) → TV (µ ˆ τ ). Finally, by using the lower-semicontinuity of the BV seminorm, and the above observations, we get
JE (µ∗ ) =
1 2
∥e(µ∗ , z )∥2V + κ TV (µ∗ ) 1
≤ lim ∥e(µ∗h , z )∥2V + κ lim inf TV (µ∗h ) h→0
2
h→0
≤ lim inf JhE (µ∗h ) ≤ lim inf JhE (µ ˆ h ) = Jκ ( µ ˆτ) h→0 h→0 1 ∇ µ ˆτ ˆ τ , z )∥2 + κ = ∥e(µ ≤
2 1 2
V
Ω
∥e(µ ˆ τ , z )∥2V + κ TV (µ) + κτ ,
which, when passed to the limit τ → 0, using the previous results, and noting that µ was chosen arbitrarily, yields
JE (µ∗ ) ≤ JE (µ),
for every µ ∈ A,
ensuring that µ ∈ A solves (22). The proof is complete. ∗
3. Computational framework In this section, we collect general discrete formulae for the proposed equation error formulation. We assume, therefore, that Th is a triangulation on Ω , Lh is the space of all piecewise continuous polynomials of degree dµ relative to Th , U¯ h is the space of all piecewise continuous polynomials of degree du relative to Th , and Qh is the space of all piecewise continuous polynomials of degree dq relative to Th . In order to represent the discrete saddle point problem in a computable form we proceed as follows. We represent bases for Lh , U¯ h and Qh by {ϕ1 , ϕ2 , . . . , ϕm } , {ψ1 , ψ2 , . . . , ψn }, and {χ1 , χ2 , . . . , χk }, respectively. The space Lh is then isomorphic to Rm and for any µ ∈ Lh , we define L ∈ Rm by Li = µ(xi ),
i = 1, 2, . . . , m,
130
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
where the nodal basis {ϕ1 , ϕ2 , . . . , ϕm } corresponds to the nodes {x1 , x2 , . . . , xm }. Conversely, each L ∈ Rm corresponds to µ ∈ Lh defined by
µ=
m
Li ϕi .
i =1
Similarly, u¯ ∈ U¯ h will correspond to U¯ ∈ Rn , where U¯ i = u(yi ), i = 1, 2, . . . , n, and u¯ =
n
U¯ i ψi ,
i =1
where y1 , y2 , . . . , yn are the nodes of the mesh defining U¯ h . Finally, q ∈ Qh will correspond to Q ∈ Rk , where Qi = q(zi ), i = 1, 2, . . . , k, and k
q=
Qi χi ,
i=1
where z1 , z2 , . . . , zk are the nodes of the mesh defining Qh . (The spaces Lh , U¯ h , and Qh are defined relative to the same elements, but the nodes will be different if dµ ̸= du ̸= dq .) Recall that the discrete saddle point problem seeks, for each µh , the unique (¯uh , ph ) ∈ U¯ h × Qh such that
2µh ϵ(¯uh ) · ϵ(¯vh ) +
Ω
Ω
(div u¯ h )qh −
1 Ω
λ
Ω
(div v¯ h )ph =
ph qh = 0,
Ω
f v¯ h +
Γ2
v¯ h f1 ,
for every v¯ ∈ U¯ h ,
for every q ∈ Qh .
(24a)
(24b)
We define S : Rm → Rn+k to be the finite element solution operator that assigns to each coefficient µh ∈ Ah , the unique approximate solution uh = (¯uh , ph ) ∈ U¯ h × Qh . Then S (L) = U, where U is defined by K (L)U = F .
(25)
The stiffness matrix K (L) ∈ R
Kn×n (L)
K ( L) =
BTn×k
(n+k)×(n+k)
n +k
and the load vector F ∈ R
are given by
−Ck×k
Bk×n
with
K (L)i,j = 2µϵ(ψj ) · ϵ(ψi ), i, j = 1, 2, . . . , n, Ω Bi,j = (div ψj )χi , i = 1, 2, . . . , k, n = 1, 2, . . . , n Ω Ci,j = χj χi , i, j = 1, 2, . . . , k, Ω Fi = f ψj + ψj f1 , i = 1, 2, . . . , n, Ω
Γ2
Fj = 0,
j = n + 1, n + 2, . . . , n + k.
For future reference, it will be useful to notice that
K (L)ij = Tijk Lk , where the summation convention is used and T is the tensor defined by
Tijk =
Ω
2ϕk ϵ(ψj ) · ϵ(ψi ),
for every i, j = 1, . . . , n, k = 1, . . . , m.
It is convenient to approximate the components of U¯ h in a single finite element space Uh where U¯ h = Uh × Uh . Therefore, if {ψ1 , . . . , ψl } are the basic of Uh then the vector-valued basis of U¯ h can be constructed by choosing n i i=1
{ψ }
ψ1 ψ2 ψl 0 0 0 = , ,..., , , ,..., 0 0 0 ψ1 ψ2 ψl
we obtain
K¯ K =
O
O . K¯
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
131
In the following U¯ and P are the discrete solutions which get replaced by the data (Z¯ , Zˆ ). From (15) and (16), we have
(K + M )E1 = K (µ)U¯ + BT P − F ¯ MQ E2 = BU − CP and consequently E1 = (K + M )−1 K (µ)U¯ + BT P − F
E2 = MQ−1 BU¯ − CP .
Using the above formulas, we have JE (µ) =
= =
1
⟨E1 , (M + K )E1 ⟩ +
2 1
1 2
E2 , MQ E2
1 K (µ)U¯ + BT P − F , (K + M )−1 K (µ)U¯ + BT P − F + BU¯ − CP , MQ−1 BU¯ − CP
2 1 2
2
1 L(U¯ )µ + BT P − F , (K + M )−1 L(U¯ )µ + BT P − F + BU¯ − CP , MQ−1 BU¯ − CP . 2
The derivative is given by DJE (µ)(δµ) = L(U¯ )δµ, (K + M )−1 L(U¯ )µ + BT P − F
and consequently,
∇ JE = L(U¯ )T (K + M )−1 L(U¯ )µ + BT P − F . 4. Numerical results In this section, we use the proposed equation error approach to identify the variable coefficient µ(x, y) in two representative examples of mixed boundary value problems in linear elasticity. In all examples, we consider an isotropic membrane of unit square dimensions, that is, Ω = (0, 1) × (0, 1), with boundary ∂ Ω = Γ1 ∪ Γ2 , where Γ1 and Γ2 are the sub-boundaries with Dirichlet and Neumann conditions, respectively. Since the main focus of this paper is the recovery of parameters in nearly incompressible materials, λ is taken as a large constant, typically λ = 106 . The examples are of a synthetic nature and the examples use the data from computed solutions using mixed finite element methods rather than from either a direct measurement or an analytic solution. In addition, although the method for solving the general saddle point problem also provides values for the pressure term Z , this data is generally unavailable in practical situations. Thus some other method must be employed to obtain Z data. Assuming that Z¯ is relatively close to the actual solution U¯ , Z can be estimated using Z = C −1 BZ¯ . This estimation, rather than an arbitrary value, was used throughout the numerical experiments. It should be noted that no systematic methods of optimal regularization parameter selection were used due to the lack of such theoretically-grounded automatic methods for problems of this type. As a result, the choice of regularization parameter was largely heuristic. 4.1. Example 1 In this first example, we consider a smooth µ (x, y) on a uniform 62 × 62 quadrangular mesh with 7688 total solution degrees of freedom. The top and bottom of the membrane (Γ1 ) are fixed with constant Dirichlet condition g (x, y) and the left and right sides (Γ2 ) have applied traction h (x, y). The functions defining the coefficient, load, and boundary conditions are as follows:
−1 1 + 0.1x2 , f (x, y) = , µ (x, y) = 1 − 0.12 cos(3π x2 + y2 ) 0.1(1 + y) 0.1x 0 g ( x, y ) = on Γ1 , h (x, y) = on Γ2 . 0.1 0.5 + y2 The numerical results presented in Fig. 1 for a smooth coefficient on the domain use an appropriate Tikhonov-type regularization scheme. The data was generated by solving the direct problem numerically. The results are encouraging, however, the method appears to be highly sensitive to any noise in Pˆ (data for the pressure term). No random noise was added. Taking this example as representative we also consider the effects of mesh refinement and the choice of regularization parameter on the maximum error. This is outlined in detail in Table 1.
132
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
(a) Exact coefficient.
(b) Early stage of reconstruction.
(c) Progressive reconstruction.
(d) Progressive reconstruction.
(e) Computed µ.
(f) Coefficient error.
Fig. 1. Numerical results for example 4.1 using smooth regularization with parameter κ = 10−6 . Table 1 Comparison of maximum error (∥µcomp − µexact ∥∞ ) for Example 1. Mesh size
κ = 10−4
κ = 10−5
κ = 10−6
κ = 10−7
κ = 10−8
32 62 92 122 152
7.4373E−03 9.4410E−04 2.8512E−04 1.1917E−04 7.6370E−05
7.2089E−03 1.0511E−03 2.7059E−04 8.5532E−05 3.0538E−05
7.1667E−03 1.1360E−03 3.1352E−04 8.3987E−05 2.2674E−05
7.1900E−03 1.1672E−03 3.4203E−04 8.6359E−05 3.0461E−05
7.1818E−03 1.1817E−03 3.5654E−04 9.6161E−05 2.5921E−05
4.2. Example 2 In this example, we consider a discontinuous µ (x, y) on a uniform 62 × 62 quadrangular mesh with 7688 total solution degrees of freedom. The top of the membrane (Γ1 ) is fixed with constant Dirichlet condition g (x, y) and the remaining
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
(a) Exact coefficient.
(b) Early stage of reconstruction.
(c) Progressive reconstruction.
(d) Progressive reconstruction.
(e) Computed µ.
133
(f) Coefficient error.
Fig. 2. Numerical results for example 4.2 using smooth regularization with parameter κ = 10−8 .
boundaries (Γ2 ) have applied traction h (x, y). The functions defining the coefficient, load, and boundary conditions are as follows:
1 + 0.2 sin(π x) 1 + 0.3 sin(π x) µ (x, y) = 1 + 0.1 sin(π x) 1 + 0.5 sin(π x) 1
f ( x, y ) =
1 + 0.1x2 , 0.1(1 + y)
for {(x, y) : for {(x, y) : for {(x, y) : for {(x, y) : otherwise, g (x, y) =
0.2 0.2 0.6 0.6
0 0
≤ x ≤ 0.4, ≤ x ≤ 0.4, ≤ x ≤ 0.8, ≤ x ≤ 0.8, on Γ1 ,
0.2 0.6 0.2 0.6
≤ y ≤ 0.4} ≤ y ≤ 0.8} ≤ y ≤ 0.4} ≤ y ≤ 0.8}
h ( x, y ) =
0.5 + x2 0 10 1
on Γ2 .
For the numerical experiment outlined in Fig. 2 we recover a discontinuous coefficient by using the BV regularization method as described above with a maximum error ∥µcomp − µexact ∥∞ = 0.001368. Again, the data was generated by
134
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
solving the direct problem numerically. In this example too, the method appears to be highly sensitive to any noise in Pˆ (data for the pressure term). No random noise was added. 5. Concluding remarks We have presented a theoretical and numerical treatment of the inverse problem of identifying variable parameters in a saddle point problem. We applied our results to identifying parameters in a system of nearly incompressible linear elasticity with an emphasis on the elasticity imaging inverse problem. We proposed a new equation error formulation and give existence results and a convergence analysis. Our results are quite encouraging. However, it has been noticed that the equation error approach is quite sensitive to the noise, an observation that has been made for the use of equation error for other inverse problems as well. Another limitation of the method is that, due to the mixed formulation, the mixed data consist of not only the displacement but also the pressure. We have given a method of computing data for P from data for U. However, our numerical experiments indicate that for good reconstruction, the data for U should be sufficiently accurate, a noisy measurement has an adverse impact on the data for P (obtained from U). This also signifies the importance of the data smoothing. It should be noticed that our approach requires full interior data (see [25]). In the future, we plan to investigate the effect of noise on parameter recovery and the subsequent issues of data-smoothing and preconditioning. In recent years, preconditioners based on Sobolev gradients have attracted a great deal of attention for the differentiation of such noisy functions and we look to apply these methods here in the context of the tumor identification problem (see [39]). Acknowledgments The authors are grateful to the referees for their careful reading and insightful remarks and suggestions that brought substantial improvement to the paper. The work of A.A. Khan is partially supported by RIT’s COS D-RIG Acceleration Research Funding Program 2012–2013 and a grant from the Simons Foundation (#210443 to Akhtar Khan). The work of B. Jadamba is partially supported by RIT’s COS Faculty Development Grant (FEAD) for 2012. References [1] K.R. Raghavan, A.E. Yagle, Forward and inverse problems in elasticity imaging of soft tissues, IEEE Trans. Nucl. Sci. 41 (1994) 1639–1648. [2] F. Kallel, M. Bertrand, Tissue elasticity reconstruction using linear perturbation method, IEEE Trans. Med. Imaging 15 (1996) 299–313. [3] M.M. Doyley, P.M. Meaney, J.C. Bamber, Evaluation of an iterative reconstruction method for quantitative elastography, Phys. Med. Biol. 45 (2000) 1521–1540. [4] M.A. Aguilo, W. Aquino, J.C. Brigham, M. Fatemi, An inverse problem approach for elasticity imaging through vibroacoustics, IEEE Trans. Med. Imaging 29 (2010) 1012–1021. [5] H. Ammari, P. Garapon, F. Jouve, Separation of scales in elasticity imgaing: a numerical study, J. Comput. Math. 28 (2010) 354–370. [6] A. Arnold, S. Reichling, O. Bruhns, J. Mosler, Efficient computation of the elastography inverse problem by combining variational mesh adaption and clustering technique, Phys. Med. Biol. 55 (2010) 2035–2056. [7] J. Fehrenbach, M. Masmoudi, R. Souchon, P. Trompette, Detection of small inclusions by elastography, Inverse Problems 22 (2006) 1055–1069. [8] J. McLaughlin, J.R. Yoon, Unique identifiability of elastic parameters from time-dependent interior displacement measurement, Inverse Problems 20 (2004) 25–45. [9] A.A. Oberai, N.H. Gokhale, G.R. Feijóo, Solution of inverse problems in elasticity imaging using the adjoint method, Inverse Problems 19 (2003) 297–313. [10] R. Acar, Identification of the coefficient in elliptic equations, SIAM J. Control Optim. 31 (1993) 1221–1244. [11] T. Kärkkäinen, An equation error method to recover diffusion from the distributed observation, Inverse Problems 13 (1997) 1033–1051. [12] M.S. Gockenbach, B. Jadamba, A.A. Khan, Numerical estimation of discontinuous coefficients by the method of equation error, Int. J. Math. Comput. Sci. 1 (2006) 343–359. [13] M.F. Al-Jamal, M.S. Gockenbach, Stability and error estimates for an equation error method for elliptic equations, Inverse Problems 28 (2012) 095006. p. 15. [14] M.S. Gockenbach, B. Jadamba, A.A. Khan, Equation error approach for elliptic inverse problems with an application to the identification of Lamé parameters, Inverse Probl. Sci. Eng. 16 (2008) 349–367. [15] N. Cahill, B. Jadamba, A.A. Khan, M. Sama, B. Winkler, A first-order adjoint and a second-order hybrid method for an energy output least squares elastography inverse problem of identifying tumor location, Bound. Value Probl. (2013) in press. [16] M.M. Doyley, B. Jadamba, A.A. Khan, M. Sama, B. Winkler, A new energy inversion for parameter identification in saddle point problems with an application to the elasticity imaging inverse problem of predicting tumor location, 2013, submitted for publication. [17] R.S. Falk, Error estimates for the numerical identification of a variable coefficient, Math. Comp. 40 (1983) 537–546. [18] T. Harrigan, E.E. Konofagou, Estimation of material elastic moduli in elastography: a local method, and an investigation of Poisson ratio sensitivity, J. Biomech. 37 (2004) 1215–1221. [19] B. Jadamba, A.A. Khan, G. Rus, M. Sama, B. Winkler, A new convex inversion framework for parameter identification in saddle point problems with an application to the elasticity imaging inverse problem of predicting tumor location, 2013, submitted for publication. [20] B. Jadamba, A.A. Khan, M. Sama, Inverse problems of parameter identification in partial differential equations, in: Mathematics in Science and Technology, World Sci. Publ., Hackensack, NJ, 2011, pp. 228–258. [21] T. Kärkkäinen, Error estimates for distributed parameter identification in parabolic problems with output least squares and Crank–Nicolson method, Appl. Math. 42 (1997) 259–277. [22] A.A. Khan, B.D. Rouhani, Iterative regularization for elliptic inverse problems, Comput. Math. Appl. 54 (2007) 850–860. [23] U. Albocher, A.A. Oberai, P.E. Barbone, I. Harari, Adjoint-weighted equation for inverse problems of incompressible plane-stress elasticity, Comput. Methods Appl. Mech. Engrg. 198 (2009) 2412–2420. [24] P.E. Barbone, A.A. Oberai, I. Harari, Adjoint-weighted variational foprmulation for direct solution of inverse heat conduction problem, Inverse Problems 23 (2007) 2325–2342. [25] P.E. Barbone, C.E. Rivas, I. Harari, U. Albocher, A.A. Oberai, Y. Zhang, Adjoint-weighted variational formulation for the direct solution of inverse problems of general linear elasticity with full interior data, Internat. J. Numer. Methods Engrg. 81 (2010) 1713–1736.
E. Crossen et al. / Computers and Mathematics with Applications 67 (2014) 122–135
135
[26] B. Banerjee, T.F. Walsh, W. Aquino, M. Bonnet, Large scale parameter estimation problems in frequency-domain elastodynamics using an error in constitutive equation functional, Comput. Methods Appl. Mech. Engrg. 253 (2013) 60–72. [27] F. Pierron, S. Avril, V.T. Tran, Extension of the virtual fields method to elasto-plastic material identfication with cyclic loads and kinematic hardening, Int. J. Solids Struct. 47 (2010) 2993–3010. [28] M.S. Gockenbach, A.A. Khan, Identification of Lamé parameters in linear elasticity: a fixed point approach, J. Ind. Manag. Optim. 1 (2005) 487–497. [29] B. Jadamba, A.A. Khan, F. Raciti, On the inverse problem of identifying Lamé coefficients in linear elasticity, Comput. Math. Appl. 56 (2008) 431–443. [30] D. Braess, Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics, third ed., Cambridge University Press, 2007. [31] Z. Cai, G. Starke, Least-squares methods for linear elasticity, SIAM J. Numer. Anal. 42 (2004) 826–842. [32] P. Houston, D. Schotzau, T.P. Wihler, An hp-adaptive mixed discontinuous Galerkin FEM for nearly incompressible linear elasticity, Comput. Methods Appl. Mech. Engrg. 195 (2006) 25–28. 3224–3246. [33] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, New York, 1991. [34] M.Z. Nashed, O. Scherzer, Stable approximation of nondifferentiable optimization problems with variational inequalities, in: Recent Developments in Optimization Theory and Nonlinear Analysis (Jerusalem, 1995), in: Contemp. Math., vol. 204, Amer. Math. Soc., Providence, RI, 1997, pp. 155–170. [35] S. Osher, A. Solé, L. Vese, Image decomposition and restoration using total variation minimization and the H 1 norm, Multiscale Model. Simul. 1 (2003) 349–370. [36] J. Zou, Numerical methods for elliptic inverse problems, Int. J. Comput. Math. 70 (1998) 211–232. [37] M.S. Gockenbach, A.A. Khan, An abstract framework for elliptic inverse problems. Part 1: an output least-squares approach, Math. Mech. Solids 12 (2007) 259–276. [38] L. Evans, R. Gariepy, Measure Theory and Fine Properties of Functions, CRC Press, Boca Raton, 1992. [39] P. Kazemi, R. Renka, A Levenberg–Marquardt method based on Sobolev gradients, Nonlinear Anal. 75 (2012) 6170–6179.