1 Introduction

Variational method for multiple parameter identification in elliptic PDEs

Tran Nhan Tam Quyen

University of Hamburg, Bundesstrasse 55, 20146 Hamburg, Germany
Email: quyen.tran@uni-hamburg.de

Abstract: In the present paper we investigate the inverse problem of identifying simultaneously the diffusion matrix, source term and boundary condition as well as the state in the Neumann boundary value problem for an elliptic partial differential equation (PDE) from a measurement data, which is weaker than required of the exact state. A variational method based on energy functions with Tikhonov regularization is here proposed to treat the identification problem. We discretize the PDE with the finite element method and prove the convergence as well as analyse error bounds of this approach.

Key words and phrases: Parameter identification, energy functional, finite element methods, diffusion matrix, source term, boundary condition

AMS Subject Classifications: 65N21, 65N12, 35J25, 35R30

## 1 Introduction

Let be an open bounded connected domain of with polygonal boundary . In this paper we study the problem of identifying simultaneously the diffusion matrix , source term and boundary condition as well as the state in the Neumann boundary value problem for the elliptic PDE

 −∇⋅(Q∇Φ) = f~{}in~{}Ω, (1.1) Q∇Φ⋅→n = g~{}on~{}∂Ω (1.2)

from a measurement of the solution , where is the unit outward normal on .

To formulate precisely our problem, let us first denote by the set of all symmetric, real -matrices equipped with the inner product and the corresponding norm , where . Furthermore, for we set

 Lp\tiny sym(Ω):={H:=(hij)i,j=¯¯¯¯¯¯1,d∈Lp(Ω)d×d | H(x):=(hij(x))i,j=¯¯¯¯¯¯1,d∈Sd~{}a.e. in~{}Ω}.

In we use the scalar product and the corresponding norm , while the space is endowed with the norm .

Let us denote by

with

and , being given constants satisfying . Let

 γ:H1(Ω)→H1/2(∂Ω)

be the continuous Dirichlet trace operator and be the closed subspace of consisting all functions with zero-mean on the boundary, i.e.,

 H1⋄(Ω):={u∈H1(Ω) ∣∣ ∫∂Ωγudx=0}

while stands for the positive constant appearing in the Poincaré-Friedrichs inequality (cf. [38])

 CΩ∫Ωφ2dx≤∫Ω|∇φ|2dx~{}for all~{}φ∈H1⋄(Ω). (1.4)

Then, due to the coervicity condition

 ∥φ∥2H1(Ω)≤1+CΩCΩ∫Ω|∇φ|2dx≤1+CΩCΩq–∫ΩQ∇φ⋅∇φdx (1.5)

holding for all and the Lax-Milgram lemma, we conclude for each , there exists a unique weak solution of (1.1)–(1.2) in the sense that and satisfies the identity

 ∫ΩQ∇Φ⋅∇φdx=(f,φ)+⟨g,γφ⟩ (1.6)

for all . Here the expressions and stand for the scalar product on space and , respectively. Furthermore, there holds the priori estimate

 ∥Φ∥H1(Ω) ≤ (1.7) ≤ CN(∥g∥L2(∂Ω)+∥f∥L2(Ω))

with

Then we can define the non-linear coefficient-to-solution operator

which maps each to the unique weak solution of the problem (1.1)–(1.2). Here, for convenience in computing numerical solutions of the pure Neumann problem (cf., e.g., [23, Subsection 5.2.1], [28, Section 2]) we normalize the solution with vanishing mean on the boundary; however, all results performed in the present paper are still valid for the normalization of solutions of the Neumann problem with zero-mean over the domain, i.e.,

The identification problem is now stated as follows:

Given , find an element

such that (1.6) is satisfied with and .

This inverse problem may have more than one solution. Following the general convergence theory for ill-posed problems (see, e.g., [9, Chapter 5] and [42, Subsection 3.2.1]), we shall use the concept of the unique minimum norm solution which is defined as

 (Q†,f†,g†):=argmin(Q,f,g)∈I(Φ†)R(Q,f,g), (1.8)

where and

 R(Q,f,g):=∥Q∥2L2\tiny sym(Ω)+∥f∥2L2(Ω)+∥g∥2L2(∂Ω).

In problems (1.10) and (1.11) below the function appears as the penalty term of the Tikhonov regularization method, and is the solution to the identification problem with the smallest penalty value.

We note that the admissible set of the problem (1.8) is non-empty, convex and weakly closed in , so that the minimizer is defined uniquely. Furthermore, the exact data may not be known in practice, thus we assume instead of to have a measurement such that

 ∥∥Φ†−zδ∥∥L2(Ω)≤δ,δ>0. (1.9)

Our identification problem is now to reconstruct from .

We also mention that the condition is weaker than required of the exact state . In practice, the observation is measured at some points and an interpolation process of these point observations is required to get the distributed observation . In the numerical implementation of §6 the observation is only assumed to be given at nodes of the coarsest triangulation grid of the domain , and the interpolation with piecewise linear, continuous elements is used to fully obtain the data . After that, the interpolation of the computed numerical state (corresponded to the coarsest grid and followed by an algorithm presented in §5) on the next finer grid is considered as an observation of the exact state on this finer grid, and so on.

Let denote a family of triangulations of the domain with the mesh size and be the approximation of the operator on the piecewise linear, continuous finite element space associated with . Furthermore, let be the Clément’s mollification interpolation operator (cf. §2). The standard method for solving the above mentioned identification problem is the output least squares one with Tikhonov regularization, i.e., one considers a minimizer of the problem

as a discrete approximation of the identified coefficient , here is the regularization parameter. However, due to the non-linearity of the coefficient-to-solution operator, we are faced with certain difficulties in holding the non-convex minimization problem (1.10). Thus, instead of working with the above least squares functional and following the use of energy functions (cf. [37, 35, 47]), in the present work the convex cost function (cf. §2)

will be taken into account. We then consider a unique minimizer of the strictly convex problem

as a discrete regularized solution of the identification problem. Note that, by using variational discretization concept introduced in [22], every solution of the minimization problem (1.11) is proved to automatically belong to finite dimensional spaces. Thus, a discretization of the admissible set can be avoided. Furthermore, for simplicity of exposition we here restrict ourselves to the case of one set of data . In case with several sets of data being available, we can replace the misfit term in the problem (1.11) by the sum .

In §3 we will show the convergence of these approximation solutions to the identification in the -norm as well as the convergence of corresponding approximation states to the exact in the -norm. Under the structural source condition — but without the smallness requirement — of the general convergence theory for non-linear, ill-posed problems (cf. [15, 16]), we prove in §4 error bounds for these discrete approximations. For the numerical solution of the minimization problem (1.11) we in §5 employ a gradient projection algorithm with Armijo steplength rule. Finally, a numerical implementation will be performed to illustrate the theoretical findings.

The coefficient identification problem in PDEs arises from different contexts of applied sciences, e.g., from aquifer analysis, geophysical prospecting and pollutant detection, and attracted great attention from many scientists in the last 30 years or so. For surveys on the subject one may consult in [3, 9, 29, 42, 44, 45]. So far there is no paper devoted to such a simultaneous identification problem. The problem of identifying the scalar diffusion coefficient has been extensively studied for both theoretical research and numerical implementation, see e.g., [7, 8, 10, 11, 17, 18, 19, 27, 30, 32, 33, 36, 40, 47]. Some contributions for the case of the simultaneous identification can be found in [2, 20, 21, 34] while some works treated the diffusion matrix case have been obtained in [14, 24, 25, 26, 39].

We conclude this introduction with the following mention. This work is a continuance of [24] to the multiple parameter identification case. On the other hand, we here used a different analysis technique which is based on the convexity of the cost functional. Surprisingly, by using the H-convergent concept, the convergence analysis presented in [14, 24] can not be applied directly to the problem of identifying scalar diffusion coefficients. There are two main difficulties for the scalar coefficient identification. First, the set

 D:={qId ∣∣ q∈L∞(Ω)~{}with~{}% q–≤q(x)≤¯¯¯q~{}a.e. in ~{}Ω~{}and% ~{}Id~{}is the unit~{}d×d-matrix}

is in general not a closed subset of under the topology of the H-convergence (cf. [46]), i.e., if the sequence is H-convergent to , then is not necessarily proportional to in dimension or . Second, the forward operator is not weakly sequentially closed in , i.e., if weakly in , it is not guaranteed that (see [14] and the references therein for counterexamples). However, it is wroth to note that above is a weakly closed subset of (cf. Remark 2.1) and so that the convergence analysis performed in the present paper covers the scalar diffusion identification case. Furthermore, in [14, 24] the source term was assumed to be given. The present situation is different, the source term is an unknown which has to be found simultaneously together with the diffusion matrix and the boundary condition from observations.

Throughout the paper we write instead of for the convenience of relevant notations. We use the standard notion of Sobolev spaces , , , etc from, e.g., [1].

## 2 Finite element discretization

### 2.1 Preliminaries

In product spaces and we use respectively the norm

 ∥(H,l,s)∥L2\tiny sym(Ω)×L2(Ω)×L2(∂Ω) = ⎛⎝∥H∥2L2\tiny sym(Ω)+∥l∥2L2(Ω)+∥s∥2L2(∂Ω)⎞⎠1/2~{}and ∥(H,l,s)∥L∞\tiny sym(Ω)×L2(Ω)×L2(∂Ω) = ∥H∥L∞\tiny sym(Ω)+∥l∥L2(Ω)+∥s∥L2(∂Ω).

We note that the coefficient-to-solution operator

with

is Fréchet differentiable on . For each the action of its Fréchet derivative in direction denoted by is the unique weak solution in to the equation

 ∫ΩQ∇ξλ⋅∇φ=−∫ΩH∇UΓ⋅∇φ+(l,φ)+⟨s,γφ⟩ (2.1)

for all .

In we introduce the convex subset

 K:={M∈Sd | q–≤Mξ⋅ξ≤¯¯¯q~{}for all~{}ξ∈Rd}

together with the orthogonal projection that is characterised by

 (A−PK(A))⋅(B−PK(A))≤0

for all and . Furthermore, let and be two arbitrary vectors in , we use the notation

 (ξ⊗η)1≤i,j≤d∈Sd with (ξ⊗η)ij:=12(ξiηj+ξjηi) for % all i,j=1,⋯,d.

We close this subsection by the following note.

###### Remark 2.1.

Let

 D:={q∈L∞(Ω) ∣∣ q–≤q(x)≤¯¯¯q~{}a.e. in ~{}Ω}.

Then is a weakly compact subset of , i.e., for any sequence a subsequence and an element exist such that is weakly convergent in to . In other words, for all there holds the limit

 limm→∞∫Ωqnmθ1=∫Ωξ∞θ1.

We also remark that any can be considered as an element in by

 ⟨Ψ,ψ⟩(L∞(Ω)∗,L∞(Ω)):=∫ΩΨψ (2.2)

for all and . Therefore, due to (2.2), the assertion of Remark 2.1 is a direct consequence of the Banach-Alaoglu theorem.

### 2.2 Discretization

Let be a family of regular and quasi-uniform triangulations of the domain with the mesh size such that each vertex of the polygonal boundary is a node of . For the definition of the discretization space of the state functions let us denote

 Vh1:={φh∈C(¯¯¯¯Ω)∩H1⋄(Ω) | φh|T∈P1(T)~{}% for all~{}T∈Th} (2.3)

with consisting all polynomial functions of degree at most . Similar to the continuous case, we have the following result.

###### Lemma 2.2.

Let be in . Then the variational equation

 ∫ΩQ∇Φh⋅∇φh=(f,φh)+⟨g,γφh⟩ (2.4)

for all admits a unique solution . Furthermore, the priori estimate

 ∥Φh∥H1(Ω)≤CN(∥f∥L2(Ω)+∥g∥L2(∂Ω)) (2.5)

is satisfied.

The map from each to the unique solution of (2.4) is called the discrete coefficient-to-solution operator. This operator is also Fréchet differentiable on the set . For each and the Fréchet differential is an element of and satisfies for all in the equation

 ∫ΩQ∇ξhλ⋅∇φh=−∫ΩH∇UhΓ⋅∇φh+(l,φh)+⟨s,γφh⟩. (2.6)

Due to the standard theory of the finite element method for elliptic problems (cf. [6, 12]), for any fixed there holds the limit

 limh→0∥∥UΓ−UhΓ∥∥H1(Ω)=0. (2.7)

Let

 Πh:L1(Ω)→{φh∈C(¯¯¯¯Ω) | φh|T∈P1(T)~{}for all~{}T∈Th}

be the Clément’s mollification interpolation operator with properties

 limh→0∥∥ϕ−Πhϕ∥∥Hk(Ω)=0\enskipfor% all\enskipk∈{0,1} (2.8)

and

 ∥∥ϕ−Πhϕ∥∥Hk(Ω)≤Chl−k∥ϕ∥Hl(Ω) (2.9)

for , where is independent of and (cf. [13, 4, 5, 43]). Then, using the discrete operator and the interpolation operator , we can now introduce the discrete cost functional

 Jhδ(Q,f,g):=∫ΩQ∇(UhQ,f,g−Πhzδ)⋅∇(UhQ,f,g−Πhzδ), (2.10)

where .

###### Lemma 2.3.

Assume that the sequence weakly converges to in . Then for any fixed the sequence converges to in the -norm.

###### Proof.

Due to Remark 2.1, has a subsequence denoted by the same symbol which is weakly convergent in to . Furthermore, by (2.5), the corresponding state sequence is bounded in the finite dimensional space . A subsequence which is not relabelled and an element then exist such that converges to in the -norm. It follows from the equation (2.4) that

 ∫ΩQn∇(UhΓn−UhΓ)⋅∇φh = ∫Ω(Q−Qn)∇UhΓ⋅∇φh (2.11) +(fn−f,φh)+⟨gn−g,γφh⟩

for all . Taking , by (1.5), we obtain that

 CΩq–1+CΩ∥∥UhΓn−UhΓ∥∥2H1(Ω) ≤ ∫Ω(Q−Qn)∇UhΓ⋅∇(UhΓn−Θh+Θh−UhΓ) +(fn−f,UhΓn−Θh+Θh−UhΓ) +⟨gn−g,γ(UhΓn−Θh+Θh−UhΓ)⟩ ≤ C∥∥UhΓn−Θh∥∥H1(Ω)+∫Ω(Q−Qn)∇UhΓ⋅∇(Θh−UhΓ)

Since weakly in , we get Sending to , we thus obtain from the last inequality that , which finishes the proof. ∎

We now state the following useful result on the convexity of the cost functional.

###### Lemma 2.4.

is convex and continuous on with respect to the -norm.

###### Proof.

The continuity of follows directly from Lemma 2.3. We show that is convex.

Let and . We have that

 UhΓ′(λ)=∂UhΓ∂QH+∂UhΓ∂fl+∂UhΓ∂gs~{}and~{}Jhδ′(Γ)(λ)=∂Jhδ(Γ)∂QH+∂Jhδ(Γ)∂fl+∂Jhδ(Γ)∂gs.

We compute for each term in the right hand side of the last equation. First we get

 ∂Jhδ(Γ)∂QH = ∫ΩH∇(UhΓ−Πhzδ)⋅∇(UhΓ−Πhzδ) +2∫ΩQ∇(∂UhΓ∂QH)⋅∇(UhΓ−Πhzδ).

For the second term we have

 ∂Jhδ(Γ)∂fl=2∫ΩQ∇(∂UhΓ∂fl)⋅∇(UhΓ−Πhzδ).

Finally, we have

 ∂Jhδ(Γ)∂gs=2∫ΩQ∇(∂UhΓ∂gs)⋅∇(UhΓ−Πhzδ).

Therefore,

 Jhδ′(Γ)(λ) = +∫ΩH∇(UhΓ−Πhzδ)⋅∇(UhΓ−Πhzδ) = 2∫ΩQ∇UhΓ′(λ)⋅∇(UhΓ−Πhzδ)+∫ΩH∇(UhΓ−Πhzδ)⋅∇(UhΓ−Πhzδ) = 2∫ΩQ∇UhΓ′(λ)⋅∇(UhΓ−¯Πhzδ)+∫ΩH∇(UhΓ−¯Πhzδ)⋅∇(UhΓ−¯Πhzδ),

where

 (2.13)

By (2.6), we infer that

 Jhδ′(Γ)(λ) = −2∫ΩH∇UhΓ⋅∇(UhΓ−¯Πhzδ)+2(l,UhΓ−¯Πhzδ)+2⟨s,γ(UhΓ−¯Πhzδ)⟩ (2.14) +∫ΩH∇(UhΓ−¯Πhzδ)⋅∇(UhΓ−¯Πhzδ) = −∫ΩH∇UhΓ⋅∇UhΓ+∫ΩH∇¯Πhzδ⋅∇¯Πhzδ +2(l,UhΓ−¯Πhzδ)+2⟨s,γ(UhΓ−¯Πhzδ)⟩.

Therefore, by (2.6) again, we arrive at

 Jhδ′′(Γ)(λ,λ) = −2∫ΩH∇UhΓ⋅∇UhΓ′(λ)+2(l,UhΓ′(λ))+2⟨s,γUhΓ′(λ)⟩ = 2∫ΩQ∇UhΓ′(λ)⋅∇UhΓ′(λ)≥2CΩq–1+CΩ∥∥UhΓ′(λ)∥∥2H1(Ω)≥0,

by (1.5), which completes the proof. ∎

Now we are in position to prove the main result of this section.

###### Theorem 2.5.

The strictly convex minimization problem

attains a unique minimizer. Furthermore, an element is the unique minimizer to if and only if the system

 Q(x) = PK(12ρ(∇UhΓ(x)⊗∇UhΓ(x)−∇¯Πhzδ(x)⊗∇¯Πhzδ(x))), (2.15) f(x) = 1ρ(¯Πhzδ(x)−UhΓ(x)), (2.16) g(x) = 1ργ(¯Πhzδ(x)−UhΓ(x)) (2.17)

holds for a.e. in , where was generated from according to (2.13).

###### Proof.

Let be a minimizing sequence of , i.e.,