Residual-based stabilized formulation for the solution of inverse elliptic partial differential equations

Abstract

We consider an inverse problem where the forward problem is that of linear plane stress elasticity, or equivalently, that of linear heat/hydraulic conduction. We demonstrate that the linearized version of the saddle point problem obtained from the minimization problem inherits some stability from the forward elliptic problem. In particular, it is stable for the response variable and the Lagrange multiplier, but not for the material property field. This lack of stability implies that we are unable to prove optimal convergence with mesh refinement for the overall problem. We overcome this difficulty by adding to the saddle point problem a residual-based term that provides sufficient stability, and prove optimal convergence in an energy-like norm. We verify these estimates through simple numerical examples. We note that while we have considered a specific model for an inverse elliptic problem in this manuscript, similar ideas could be developed for a broad class of inverse elliptic problems.

Introduction

The problem of determining the material parameters from a given set of interior measurements is an inverse problem. The determination of shear modulus from displacement measurements [1], thermal conductivity from temperature measurements [2], and the aquifer permeability from water pressure head measurements [3], are examples of these types of problems.

In the inverse two-dimensional plane stress problem, given the measured displacement components $\tilde{\mathit{u}}=[{\tilde{u}}_1,{\tilde{u}}_2]$, the aim is to determine the shear modulus distribution, $\mu$, such that,

$$\nabla \cdot (\mu ℂ\nabla \tilde{\mathit{u}})=\mathbf{0}.$$

Here $ℂ$ is a 4th order tensor given by ${ℂ}_{ijkl}={\delta }_{ij}{\delta }_{kl}+\frac{1}{2}({\delta }_{ik}{\delta }_{jl}+{\delta }_{jk}{\delta }_{il})$.

The appropriate boundary condition for the above system to be well posed is that $\mu$ is prescribed at a point of the domain, or the average value of $\mu$ is prescribed. That is,

$$\frac{{\int }_{\Omega }\mu d\Omega }{{\int }_{\Omega }d\Omega }={\mu }_0.$$

We note that the system (1), along with the condition (2) can represent other inverse problems, such as the inverse heat conduction problem, and the inverse hydraulic problem. In the two-dimensional heat conduction (hydraulic) problems ${\tilde{u}}_1$ and ${\tilde{u}}_2$ represent two independent temperature (pressure) field measurements, $\mu$ represents the thermal (hydraulic) conductivity, and the tensor $ℂ$ is given by ${ℂ}_{ijkl}={\delta }_{ijkl}$.

Several methods have been proposed to solve the inverse problem defined in (1), (2). They can be broadly classified as “direct” or minimization methods. In direct methods, the measured displacement field is directly inserted in (1), and this equation is interpreted as an equation for the shear modulus $\mu$. This system yields two hyperbolic PDEs for $\mu$ that must be solved simultaneously with data (2). The simplest approach to solve these equations is to consider a least-squares formulation [4]. However, it has been shown that this results in numerical results that are overly dissipative, especially for rough measured data [2]. An alternate is to consider the adjoint-weighted variational equation, where the original equation for $\mu$ is weighted by its adjoint operating on the weighting functions [1], [5]. This is in contrast to the least squares approach, where the original equation is weighted by the same operator acting on the weighting function. In results presented in [2] the adjoint-weighted equations have demonstrated superior performance (particularly for sharp spatial changes) when compared with the least-squares method.

Direct methods are appealing in that they are easily implemented and computationally inexpensive. On the other hand, since they use the measured displacement data and its derivatives directly in the equation of equilibrium, they have difficulty in handling noisy data obtained from real experiments. Recently, methods that smooth this data under minimal assumptions on the underlying modulus distributions and address this issue to some extent have been proposed [6].

By minimizing a data mismatch, optimization methods of solving the inverse problem circumvent the issue of noisy data in pde (1) problem entirely. In these methods the displacement field that appears in (1) is a smoother, predicted displacement field, which is computed by solving this equation for a given distribution of $\mu$. The difference between the predicted and measured displacement field is minimized by varying the spatial distribution of the shear modulus. The shear modulus is itself often smoothed through the use of a regularization term. This approach leads to a regularized constrained minimization problem, where the difference between the predicted and measured displacement fields is minimized under the constraint of the equation of equilibrium. When the difference between the predicted and measured fields is measured in the $L_2$ norm, we are lead to the problem considered in [7], [8]. Another choice is to measure the difference between the predicted and measured displacements in the “energy” norm induced by the shear modulus distribution, as described in [9], [10].

The solution of the constrained minimization problem leads to variational equations that are most naturally approximated using the finite element method. A question then arises about the convergence with mesh refinement of the finite dimensional solution to the exact solution. In [3], for the $L_2$ minimization approach, in the context of the inverse hydraulic problem, with a single measurement it was shown that for perfect data, the hydraulic conductivity converges sub-optimally in the $L_2$ norm to the exact modulus distribution. Further, in [9], [11], it was shown that the finite element discretization of the variational formulation associated with a mismatch measured in the energy norm, converges at optimal rates with the mesh refinement.

In this paper, following [3], we consider the $L_2$ minimization problem. However we consider the two-dimensional plane stress problem, or equivalently the heat conduction problem with two measurements. Our Lagrangian is comprised of an $L_2$ displacement mismatch term, and a constraint equation. We observe that the linearized version of the variational equations obtained from the Lagrangian lacks stability. We address this problem by adding to the variational equations terms that are driven by the residual of the constraint. We note that these additional terms add stability to the linearized version of the variational equation. This approach is similar to an augmented Lagrangian method [11], [12], [13], [14], [15]. In the augmented Lagrangian approach, a penalty term that is driven by the constraint is added to the Lagrangian. In contrast to this, we append a term driven by the constraint directly to the variational equation derived from the stationarity of the original Lagrangian. This allows us to consider a larger range of the possible terms in the variational equations than the augmented Lagrangian method, which only allows for least squares type terms.

We note that setting the first variation of the Lagrangian to zero leads to a nonlinear saddle point problem, the stability of whose linear counterpart can be analyzed using techniques that are typically used for analyzing mixed variational problems [16], [17]. This includes probing the so-called inf–sup stability of the variational problem. A common approach to developing viable numerical methods for mixed problems involves selecting function spaces for trial solutions and weighting functions that guarantee inf–sup stability [18], [19]. Another approach involves adding terms that are consistent, in that they are proportional to residual of the original problem, and enhance the stability of the original variational problem. In this case often the goal is to derive a variational formulation that is coercive in a problem-dependent norm for a large class of weighting functions and trial solutions [20], [21], [22]. The approach described in this paper first makes the connection between the solution of an inverse problem and a mixed variational problem, and thereafter it relies on adding residual-based terms to derive a provably stabilized formulation.

The layout of the remainder of this manuscript is as follows. In Section 2, we pose the inverse problem as $L_2$ constrained minimization problem. We construct the Lagrangian to incorporate the constraint using Lagrange multipliers, and then derive the non-linear saddle point problem (SPP) by setting its first variation to zero. We discretize the non-linear SPP using standard finite element shape functions. We show that the linearized version of the saddle point problem lacks stability. We account for this by adding a residual based stabilization term to the saddle point problem. In Section 3, we consider the convergence of the stabilized formulation with mesh refinement and prove optimal convergence rates in an “energy norm”. In Section 4, we perform numerical verification of the stabilized formulation, and demonstrate the optimal convergence rates using two sample problems.