https://doi.org/10.1365/s13291-021-00236-2 S U R V E Y A R T I C L E
An Introduction to Finite Element Methods for Inverse Coefficient Problems in Elliptic PDEs
Bastian Harrach1
Accepted: 3 May 2021 / Published online: 6 May 2021
© The Author(s) 2021
Abstract
Several novel imaging and non-destructive testing technologies are based on recon- structing the spatially dependent coefficient in an elliptic partial differential equation from measurements of its solution(s). In practical applications, the unknown coeffi- cient is often assumed to be piecewise constant on a given pixel partition (correspond- ing to the desired resolution), and only finitely many measurement can be made. This leads to the problem of inverting a finite-dimensional non-linear forward operator F: D(F)⊆Rn→Rm, where evaluatingFrequires one or several PDE solutions.
Numerical inversion methods require the implementation of this forward opera- tor and its Jacobian. We show how to efficiently implement both using a standard FEM package and prove convergence of the FEM approximations against their true- solution counterparts. We present simple example codes for Comsol with the Matlab Livelink package, and numerically demonstrate the challenges that arise from non- uniqueness, non-linearity and instability issues. We also discuss monotonicity and convexity properties of the forward operator that arise for symmetric measurement settings.
This text assumes the reader to have a basic knowledge on Finite Element Methods, including the variational formulation of elliptic PDEs, the Lax-Milgram- theorem, and the Céa-Lemma. Section3also assumes that the reader is familiar with the concept of Fréchet differentiability.
Keywords Finite element methods·Inverse problems·Finitely many measurements·Piecewise-constant coefficient
1 Introduction
Many practical reconstruction problems in the field of medical imaging and non- destructive testing lead to inverse coefficient problems in elliptic partial differential
B. Harrachharrach@math.uni-frankfurt.de
1 Institute for Mathematics, Goethe-University Frankfurt, Frankfurt am Main, Germany
equations. This text is meant to be an introductory tutorial for implementing such problems with Finite Element Methods (FEM).
We assume that the unknown coefficient is piecewise-constant on a given resolu- tion, and that finitely many linear measurements of one of several solutions are taken, where different solutions are generated by different linear excitation in the underly- ing physics model. This leads to the finite-dimensional non-linear inverse problem of determining
σ∈Rn from F(σ )∈Rm withn∈Nunknowns andm∈Nmeasurements.
Iterative numerical solution methods for this inverse problem require evaluating F and its derivatives at each iteration step, which means solving the underlying el- liptic PDE. In this work, we will demonstrate how FEM-based implementations for F and its Jacobian can be obtained very efficiently from standard FEM-solvers for the considered elliptic PDE. Roughly speaking, the sensitivity of a measurement with respect to changing the coefficient in one pixel can be simply calculated by multiply- ing FEM-solutions corresponding to the measurement and excitation patterns with so-called pixel stiffness matrices that are obtained from summing up all element stiff- ness matrices of elements belonging to the pixel where the change occurs. Hence, the FEM-based Jacobian can be obtained without any additional computational burden with just a few lines of extra code. Alternatively, for an even simpler implementa- tion, the pixel stiffness matrices can be easily obtained by subtracting global stiffness matrices without requiring any knowledge about the triangulation details.
This text is meant as a simple-to-read explanation of this approach in a sufficiently general but naturally arising setting. More precisely, we restrict ourselves to coercive and symmetric variational formulations that linearly depend on the unknown coeffi- cients, and to excitations and measurements that correspond to linear functionals. In this setting, we demonstrate how to obtain the Jacobian of the FEM-based forward map with the means of a standard FEM software package such as COMSOL. We also discuss monotonicity and convexity properties arising in symmetric measurement sit- uations that are the basis for recent research on rigorously justified reconstruction methods.
The purpose of this text is of introductory nature, but we proceed in a mathe- matically rigorous fashion to allow this text to also serve as a reference. We prove differentiability of the true-solution forward operator and its FEM-based approxima- tion, and show convergence of the FEM-approximated quantities to their true-solution counterparts.
Section2gives two examples to motivate our general setting: stationary diffusion and Elecrical Impedance Tomography. Section3introduces the forward operator us- ing the exact PDE solution and derives its properties. The FEM-approximation of the forward operator and its Jacobian is studied in Sect.4. In Sect.5we show numerical examples and demonstrate some of the major challenges that arise in solving inverse coefficient problems. The COMSOL/MATLAB source codes for all examples are given in theAppendix.
Fig. 1 Pixel partition and circular subdomains used for excitations and measurements
2 Motivation and Examples 2.1 Stationary Diffusion
We consider the stationary diffusion equation
−∇ ·(σ∇u)=g in (1) with homogeneous Dirichlet boundary conditionu|∂=0 in a Lipschitz bounded domain⊂Rd,d∈N. Foru∈H1(),σ∈L∞+()andg∈L2()the equation is equivalent to findingu∈H01()with
σ∇u· ∇vdx=
gvdx for allv∈H01(), (2) and unique solvability follows from the Lax-Milgram theorem.
We are interested in the inverse coefficient problem of determining the diffusivity coefficient σ in (1) from measurements of the solution for one or several source termsg, cf. [9] for an application in groundwater filtration. In practical applications with finitely many measurements, it is natural to only aim for a certain pixel-based resolution and therefore assume thatσis piecewise constant with respect to a partition =n
i=1Pi, i.e.
σ (x)= n
i=1
σiχPi(x) for allx∈,
where the pixelsPi ⊆are assumed to be measurable subsets. The left image in Fig.1shows a simple example where the unit square=(0,1)2is divided into 3×3 pixels. In the following, with a slight abuse of notation, we writeσ=(σ1, . . . , σn)T ∈ Rnfor the unknown diffusivity.
Fig. 2 PDE solution for source terms on circular subdomains
The source term g in the diffusion model (1) can be identified with the linear functional on the right hand side of the variational formulation (2)
l∈H−1(), l(v):=
gvdx,
which corresponds to identifyingL2()with a subset ofH−1(). Accordingly, we consider excitations in the form of linear functionals. Also, to emphasize that the solution depends on the diffusion coefficient and the excitation, we writeulσ in the following. The left image in Fig.2illustrates the concentration resulting from a con- stant source termg=χD, i.e.l(v)=
Dvdx, whereD=D2∪D4∪D5∪D7 is a union of four circular subdomains as sketched in the right image of Fig.1. The right image in Fig.2 shows the corresponding plot forD=D1∪D3∪D6∪D8. Both images show the solution of (1) with constant diffusion coefficientσ=1.
Natural models for measuring the solution of (1) also yield to linear functionals.
Measuring the total concentration in one of the circular subdomainsDj corresponds to measuringr(u):=
Djudx. Hence, the inverse problem of determining finitely many information about the diffusivity coefficient from finitely many measurements of the concentration (possibly but not necessarily resulting from different excitations) leads to the finite-dimensional inverse problem to
determine σ∈Rn+ from F(σ )∈Rm, where
F: Rn+→Rm, F(σ ):=(rj(ulσj))mj=1, andulσj∈H01()solves
n
i=1
σibi(ulσj, v)=lj(v) for allv∈H01(),
with givenlj, rj∈H−1(),j=1, . . . , m, and bi(u, v):=
Pi
∇u· ∇vdx, i=1, . . . , n.
2.2 Electrical Impedance Tomography (EIT)
We give another example for an application that leads to an inverse elliptic coefficient problem in a similar form as the diffusion example.
EIT aims to image the inner conductivity structure of a subject by current and voltage measurements through electrodes attached to the imaging subject. Let⊆ Rd,d ∈ {2,3}, be a smoothly bounded domain denoting the imaging subject. The electrodesEk,k=1, . . . , K, are assumed to be open connected subsets of∂with disjoint closures.
When currents with strengthJ =(J1, . . . , JK)∈RK are driven through the K electrodes (withK
k=1Jk=0), the resulting electric potentialu∈H1()inside, and the potentialU∈RKon the electrodes, solve
∇ ·(σ∇u)=0 in,
σ ∂νu=0 on∂\
K
k=1
Ek,
u+zσ ∂νu=const.=:Uk onEk, k=1, . . . , K,
Ek
σ ∂νu|Ekds=Jk onEk, k=1, . . . , K,
whereσ∈L∞+()is the conductivity inside, andz >0 is the contact impedance of the electrodes.
Under the gauge conditionU∈RK:= {V ∈RK: K
k=1Vk=0}, one can show (see [19]) that this so-called complete electrode model (CEM) for EIT is equivalent to the variational formulation that(u, U )∈H1()×RK solves
σ∇u· ∇wdx+ K
k=1
Ek
1
z(u−Uk)(w−Wk)ds= K
k=1
JkWk (3)
for all(w, W )∈H1()×RK, and unique solvability follows from the Lax-Milgram theorem.
We assume thatz >0 is known, and that σ (x)=n
i=1σiχPi(x)is piecewise constant with respect to a pixel partition=n
i=1Pi, and writeσ=(σ1, . . . , σn)∈ Rnfor the unknown conductivity values inside.
The applied current patternsJ =(J1, . . . , JK)∈RK can be identified with the functional
l∈H, l(w, W ):=
K
k=1
JkWk for all(w, W )∈H:=H1()×RK.
Likewise, measuring the voltage between the k1-th and the k2-th electrode corre- sponds to measuring the linear functional
r∈H, r(u, U ):=Uk1−Uk2, of the solution(u, U )∈H generated by some current pattern.
Hence, the problem of determining the interior conductivity with a fixed finite resolution from finitely many voltage-current measurements in EIT (with CEM) leads to the finite-dimensional inverse problem to
determine σ∈Rn+ from F(σ )∈Rm,
whereF: Rn+→Rm,F(σ ):=(rj(ulσj, Uσlj))mj=1, and(ulσj, Uσlj)∈H solves
b0((ulσj, Uσlj), (w, W ))+ n
i=1
σibi((ulσj, Uσlj), (w, W ))=lj(w, W )
for all(w, W )∈H, with givenlj, rj ∈H,j=1, . . . , m, and
b0((u, U ), (w, W )):=
K
k=1
Ek
1
z(u−Uk)(w−Wk)ds, bi((u, U ), (w, W )):=
Pi
∇u· ∇wdx.
Clearly, one could also extend this formulation to cover the case of unknown contact impedances.
3 The True-Solution Setting
The examples in Sect.2lead to inverse problems for a finite-dimensional non-linear forward operatorF: Rn+→Rm, where evaluations ofFrequire solving an infinite- dimensional linear problem (the PDE). In this section, we will first derive some prop- erties ofFfor the case that it is defined with the true infinite-dimensional PDE solu- tion. The properties of the operatorF≈F, that is defined with a FEM-approximation of the PDE solution, will be studied in Sect.4.
3.1 The True-Solution Forward Operator and Its Derivative
We will study problems that appear in the variational formulation of elliptic PDEs with piecewise constant coefficients on a fixed pixel partition, as in the examples in Sect.2.
The Variational Setting LetHbe a Hilbert space. We consider the problem of finding u∈Hthat solves
bσ(u, v)=l(v), (4)
wherebσ: H×H→Ris a bilinear form, andl∈H=L(H,R).bσ is assumed to linearly depend onnparametersσ=(σ1, . . . , σn)∈Rnin the following way
bσ(u, v)=b0(u, v)+ n
i=1
σibi(u, v),
whereb0, bi : H×H→Rare bounded, symmetric and positive semidefinite bi- linear forms. Writing1:=(1, . . . ,1)T ∈Rn, we also assume thatb1is bounded and coercive with constantsβ, C >0, i.e.,
Cv2≥b1(v, v)=b0(v, v)+ n
i=1
bi(v, v)≥βv2 ∀v∈H.
Clearly, this yields that for allσ∈Rn+
Cmax{1, σ1, . . . , σn}v2≥bσ(v, v)≥βmin{1, σ1, . . . , σn}v2 ∀v∈H, (5) so thatbσ is symmetric, bounded and coercive. Here and in the followingRn+denotes the set of allσ ∈Rn withσ >0 and “>” and “≥” are understood elementwise on Rn.
The True-Solution Forward Operator We now characterize the derivative of the solu- tion of (4) with respect toσ.
Lemma 1 Letl∈H. The solution operator
S: Rn+→H, S(σ ):=ulσ, whereulσ∈Hsolves (4), is infinitely often Fréchet differentiable. Its first derivative
S: Rn+→L(Rn, H )
fulfills that, for allσ∈Rn+andτ∈Rn,S(σ )τ∈His the unique solution of bσ(S(σ )τ, w)= −
n
i=1
τibi(ulσ, w) ∀w∈H.
Also, forr∈H,σ∈Rn+, andτ∈Rn, r(ulσ)=bσ(ulσ, urσ) and r
S(σ )τ = − n
i=1
τibi(ulσ, urσ).
Proof For σ ∈Rn+, the Riesz theorem yields that there exists a unique operator B(σ )∈L(H, H)associated to the bilinear formbσ(·,·), i.e.
B(σ )u, vH×H=bσ(u, v) for allu, v∈H.
Clearly,B(σ ) is symmetric and, by the Lax-Milgram theorem,B(σ ) is invertible with symmetric inverseB(σ )−1∈L(H, H ). Hence, (4) is uniquely solvable, and the solution operatorSis well-defined.
It is easily checked, thatB(σ )is Fréchet differentiable for everyσ∈Rn+, and that its derivativeB(σ )∈L(Rn,L(H, H))is given by
B(σ )τ= n
i=1
τiBi for allσ∈Rn+, τ∈Rn, whereBi∈L(H, H)is the unique operator fulfilling
(Biu, v)H×H=bi(u, v) for allu, v∈H.
Since B(σ ) does not depend on σ, this also shows that B(σ ) is infinitely often Fréchet differentiable with all second and higher derivatives being zero.
Using the derivative of operator inversion and the product and chain rule for the Fréchet derivative, we thus obtain thatS(σ )is infinitely often Fréchet differentiable with
S(σ )τ= −B(σ )−1(B(σ )(τ ))B(σ )−1l= − n
i=1
τiB(σ )−1Biulσ.
Hence,v=S(σ )τ∈H solves bσ(v, w)= −
n
i=1
τiBiulσ, wH×H= − n
i=1
τibi(ulσ, w) ∀w∈H.
Moreover, we obtain for allr∈H, by using the symmetry ofB(σ ), r(ulσ)= B(σ )B(σ )−1r, ulσH×H=bσ(ulσ, urσ), and
r
S(σ )τ = r,S(σ )τH×H= B(σ )S(σ )τ,B(σ )−1rH×H
= − n
i=1
τiBiulσ, urσ,H×H= − n
i=1
τibi(ulσ, urσ),
which finished the proof.
Corollary 1 Letl, r∈H. Then the mapping
Fl,r: Rn+→R, Fl,r(σ ):=r(ulσ)
fulfills
Fl,r(σ )=bσ(ulσ, urσ) for allσ∈Rn+.
Moreover, Fl,r : Rn+→R is infinitely often differentiable and its first derivatives fulfill
∂
∂σiFl,r(σ )= −bi(ulσ, urσ).
Proof This follows from Lemma1.
3.2 Convexity and Monotonicity for Symmetric Measurements
A special mathematical structure appears for measurementsFl,r, whenl andr are taken from the same subset ofH, and all combinations are used. In the stationary diffusion example this corresponds to using the same subsets of both for exci- tations and concentration measurements, in EIT this corresponds to using the same electrodes for voltage and current measurements.
Given a set ofm∈Nexcitations/measurements{l1, . . . , lm} ⊂H, we combine the measurements into a matrix-valued mapF: Rn+→Rm×m
F(σ )=(Fj,k(σ ))j,k=1,...,m∈Rm×m, Fj,k(σ )=Flj,lk(σ )=lk(ulσj).
As before, we write “≥” for the elementwise order on Rn. We also writeSm⊆ Rm×mfor the subset of symmetricm×m-matrices, and “” for the Loewner order onSm, i.e.BAdenotes thatB−Ais positive semi-definite.
Lemma 2 F: Rn+→Rm×mhas the following properties:
(a) Fis infinitely often differentiable.
(b) For all σ ∈ Rn+, F(σ ) ∈ Sm and F(σ ) 0. F(σ ) is positive definite if l1, . . . , lm∈Hare linearly independent.
(c) Fis monotonically non-increasing, i.e.
F(σ )τ0 for all σ ∈Rn+, 0≤τ ∈Rn, (6) and for allσ(1), σ(2)∈Rn+
σ(1)≤σ(2) implies F(σ(1))F(σ(2)). (7) (d) Fis convex, i.e., for allσ, σ(0)∈Rn+,
F(σ )−F(σ(0))F(σ(0))(σ−σ(0)), (8) and, for allt∈ [0,1],
F((1−t )σ(0)+t σ )(1−t )F(σ(0))+tF(σ ). (9)
Proof Corollary1shows that each component ofF is infinitely often differentiable so that (a) is proven.
For the rest of the proof letσ ∈Rn+,g∈Rm, and setl:=m
j=1gjlj. By Corol- lary1,
lk(ulσj)=bσ(ulσj, ulσk)=bσ(ulσk, ulσj)=lj(ulσk), so thatF(σ )is a symmetric matrix. Moreover,
gTF(σ )g= m
j,k=1
gjlk(ulσj)gk= m
j,k=1
gjgkbσ(ulσj, ulσk)=bσ(ulσ, ulσ)≥0, so thatF(σ )0. Ifg=0 andl1, . . . , lm∈H are linearly independent thenl=0, which impliesulσ =0 and thusgTF(σ )g >0. Hence, (b) is proven.
To prove (c) and (d), we start by using again Corollary1and obtain gT(F(σ )τ )g= −
n
i=1
τibi(ulσ, ulσ) for allτ∈Rn. Since the bilinear formsbi(·,·)are positive semi-definite, this implies (6).
To prove (8), letσ(0)∈Rn+. For brevity we writeul0:=ulσ
0. Using bσ(ulσ, ul0)=l(ul0)=bσ0(ul0, ul0),
we obtain that
0≤bσ(ulσ−ul0, ulσ −ul0)=bσ(ulσ, ulσ)−2bσ(ulσ, ul0)+bσ(ul0, ul0)
=bσ(ulσ, ulσ)−2bσ0(ul0, ul0)+bσ(ul0, ul0)
=gT(F(σ )−F(σ(0)))g+bσ(ul0, ul0)−bσ0(ul0, ul0)
=gT(F(σ )−F(σ(0))g+ n
i=1
(σi−σi(0))bi(ul0, ul0).
This shows that
gT(F(σ )−F(σ(0)))g≥ − n
i=1
(σi−σi(0))bi(ul0, ul0)=gTF(σ(0))(σ−σ(0))g, so that (8) holds. Together with (6) this also implies (7).
(9) follows from (8) by the following standard argument. Let σ, σ(0)∈Rn+,t∈ [0,1], and set
σ(t ):=t σ +(1−t )σ(0)∈Rn+.
Using (8) onF(σ )−F(σ(t ))andF(σ(0))−F(σ(t )), we then obtain that (1−t )F(σ(0))+tF(σ )−F(σ(t ))
=(1−t )(F(σ(0))−F(σ(t )))+t (F(σ )−F(σ(t ))) (1−t )F(σ(t ))(σ(0)−σ(t ))+tF(σ(t ))(σ−σ(t ))
=F(σ(t ))((1−t )σ(0)+t σ −σ(t ))=0,
which proves (9).
4 The FEM Setting
4.1 The FEM-Approximated Forward Operator and Its Derivative
The Finite Element Method The Finite Element Method numerically approximates the solution of (4) by solving it in a finite-dimensional subspaceV ⊂H, e.g. the subspace of continuous, piecewise linear functions on a fixed triangulation. Let
1, . . . , N denote a basis of V, e.g. the so-called hat functions for linear finite elements. Then the finite-dimensional variational problem
˜
ulσ ∈V solves bσ(u˜lσ, v)=l(v) for allv∈V (10) is equivalent to
˜ ulσ =
N
j=1
λj j, where Bσλ=yl, (11)
withλ=(λj)Nj=1∈RN, and the so-called stiffness matrix and load vector Bσ∈RN×N, with(j, k)-th entry given bybσ( j, k),
yl∈RN, withj-th entry given byl( j).
It follows from the Lax-Milgram theorem that (10) is uniquely solvable and thatBσ is a symmetric, positive definite (and thus invertible) matrix. Moreover, the Céa-Lemma yields that the FEM approximationu˜lσ ∈V is as good an approximation to the true solutionulσ∈Has elements of the finite-dimensional spaceV can be:
ulσ− ˜ulσ ≤Cσ βσ inf
v∈Vulσ−v, (12)
whereCσ:=Cmax{1, σ1, . . . , σn}, andβσ :=βmin{1, σ1, . . . , σn}are the continuity and coercivity constants ofbσ, cf. (5).
Pixel Stiffness Matrices Finite element software packages include triangulation algo- rithms, assembling routines for the global stiffness matrixBσ and the load vectoryl, and efficient solvers for the linear systemBσλ=yl. For our setting where
bσ(u, v)=b0(u, v)+ n
i=1
σibi(u, v),
Fig. 3 A coarser and a finer FEM-mesh for the diffusion example, both complying with the pixel partition and the measurement/excitation subdomains
we will also require the pixel stiffness matrices
Bi∈RN×N, with(j, k)-th entry given by bi( j, k).
The assembling ofBσ is usually done by writing it as a weighted sum of ele- ment stiffness matrices. In our setting, it is natural to assume that the pixel partition complies with the FEM triangulation, i.e., that each pixel is a union of triangulation elements. Figure3shows a coarser and a finer FEM mesh for the diffusion exam- ple, both complying with the pixel partition and with the subdomains that are used for measurements and excitations. Hence, during the assembly of the global stiffness matrixBσ, the pixel stiffness matrices can usually be obtained without any additional computational cost by the simple intermediate step of first summing up the element matrices for each pixel, and then summing up the pixel stiffness matrices to obtain Bσ. Alternatively, the pixel stiffness matrixBi can be conveniently obtained from global stiffness matrices by the simple identities
Bi=B1+ei−B1, and B0=B1− n
i=1
Bi,
whereB1+ei andB1denote the global stiffness matrixBσ forσ=1+eiandσ=1, respectively, andei∈Rn is thei-th unit vector. Note that this does not require any knowledge of the triangulation details.
The FEM-Approximated Forward Operator Givenl, r∈H, we approximate the true measurementFl,r(σ )=r(ulσ)by
Fl,r(σ ):=r(u˜lσ),
whereu˜lσ ∈V is the FEM-approximation to the true solutionulσ ∈H, i.e., the solu- tion of (10).
Algorithm 1 FEM-approximation ofFl,r(σ )and∂σ∂
iFl,r(σ ),i=1, . . . , n givenl, r∈H,σ∈Rn+
·use FEM package to calculate load vectorsylandyr
·use FEM package to calculate stiffness matricesB1, andB1+eifor alli=1, . . . , n
·setBi:=B1+ei−B1fori=1, . . . , n, andBσ:=B1+n
i=1(σi−1)Bi
·solveBσλl=ylandBσλr=yrforλlandλr return Fl,r(σ ):=(λl)Tyr and∂σ∂
iFl,r(σ ):= −(λl)TBiλr,i=1, . . . , n
Lemma 3 Letl, r∈H. Then
Fl,r(σ )=bσ(u˜lσ,u˜rσ) for allσ∈Rn+.
Moreover, Fl,r : Rn+→R is infinitely often differentiable and its first derivatives fulfill
∂
∂σi
Fl,r(σ )= −bi(u˜lσ,u˜rσ), i=1, . . . , n.
Proof This follows by applying Corollary1to the Hilbert spaceV. From Lemma3, we obtain a simple FEM-based implementation of the forward operator and its derivative.
Corollary 2 With
˜ ulσ=
N
j=1
λlj j, λl=(λlj)Nj=1∈RN,
˜ urσ=
N
j=1
λrj j, λr=(λrj)Nj=1∈RN,
we have that
Fl,r(σ )=(λl)TBσλr=(λl)Tyr, and ∂
∂σiFl,r(σ )= −(λl)TBiλr.
Proof This follows from Lemma3.
We summarize the consequences of Corollary2 in Algorithm 1. Using a FEM package that is capable of solving the considered PDE, and that allows access to the stiffness matrix and the load vector, one can simply implement the FEM- approximated forward operator and all its first derivatives by a few lines of extra code. This calculation merely requires solving two linear systems with the stiffness matrix (which is equivalent to two PDE solutions).
Convergence of the FEM-Approximated Forward Operator The following lemma shows that the FEM-approximated operator and its first derivatives agree with their true-solution counterparts as good as the FEM solution agrees with the true solu- tion. Hence, by the Céa-Lemma (12), Fl,r(σ ) and ∂σ∂
iFl,r(σ ) will be as good an approximation toFl,r(σ )and ∂σ∂
iFl,r(σ )as the true solutions can be approximated by elements of the finite-dimensional spaceV.
Lemma 4 For alll, r∈Handσ∈Rn+, we have that:
Fl,r(σ )−Fl,r(σ )=bσ(ulσ − ˜ulσ, urσ− ˜urσ), (13)
∂
∂σiFl,r(σ )− ∂
∂σiFl,r(σ )=bi(u˜lσ,u˜rσ−urσ)+bi(u˜lσ−ulσ, urσ). (14) Hence, by the Céa-Lemma (12),
0≤Fl,r(σ )−Fl,r(σ )≤Cσulσ − ˜ulσ urσ − ˜urσ
≤ Cσ3 βσ2 inf
v∈Vulσ−v inf
v∈Vurσ−v,
and
∂
∂σiFl,r(σ )− ∂
∂σiFl,r(σ )
≤Ci ˜ulσ ˜urσ−urσ +Ciurσ ˜ulσ −ulσ
≤CiCσ βσ
˜ulσ inf
v∈Vurσ−v + urσ inf
v∈V ulσ −v
, whereCi>0 is the continuity constant ofbi(·,·).
Proof Using
bσ(u˜lσ,u˜rσ)=l(u˜rσ)=bσ(ulσ,u˜rσ), and bσ(u˜lσ,u˜rσ)=r(u˜lσ)=bσ(u˜lσ, urσ), we obtain (13) from
Fl,r(σ )−Fl,r(σ )=bσ(ulσ, urσ)−bσ(u˜lσ,u˜rσ)=bσ(ulσ, urσ− ˜urσ)
=bσ(ulσ, urσ− ˜urσ)−bσ(u˜lσ, urσ− ˜urσ)
=bσ(ulσ− ˜ulσ, urσ− ˜urσ).
Also,
∂
∂σiFl,r(σ )− ∂
∂σiFl,r(σ )=bi(u˜lσ,u˜rσ)−bi(ulσ, urσ)
=bi(˜ulσ,u˜rσ−urσ)+bi(u˜lσ−ulσ, urσ),
which shows (14).
4.2 Convexity and Monotonicity for Symmetric Measurements
As in Sect.3.2we now consider the symmetric measurement case, wherelandrare taken from the same subset ofH (and all combinations are used). Given a set of m∈Nexcitations/measurements {l1, . . . , lm} ⊂H, we combine the measurements into a matrix-valued mapF : Rn+→Rm×m
F (σ )=(Fj,k(σ ))j,k=1,...,m∈Rm×m, Fj,k(σ )=Flj,lk(σ )=lk(ulσj).
The entries ofF (σ ) and its first derivatives ∂σ∂
iF (σ ), i=1, . . . , ncan be cal- culated as in Algorithm1. Let us stress that this approach is particularly efficient in this symmetric case as it requires onlymlinear system solutions with the stiff- ness matrix (i.e., the equivalent ofmPDE solutions) for calculating allm2entries of F (σ )∈Rm×mand allnm2entries of thenmatrices ∂σ∂
iF (σ )∈Rm×m.
As in Sect.3.2, the FEM-approximated forward operator is monotonically non- increasing and convex in the sense of the elementwise order “≥” on Rn, and the Loewner order “” on the set of symmetricm×m-matrices.
Lemma 5 F : Rn+→Rm×mhas the following properties:
(a) F is infinitely often differentiable.
(b) For allσ∈Rn+,F (σ )∈SmandF (σ )0.F (σ )is positive definite ifl1, . . . , lm∈ Hare linearly independent.
(c) F is monotonically non-increasing, i.e.
F(σ )τ0 for all σ∈Rn+, 0≤τ∈Rn, (15) and for allσ(1), σ(2)∈Rn+
σ(1)≤σ(2) implies F (σ(1))F (σ(2)). (16) (d) F is convex, i.e.
F (σ )−F (σ(0))F(σ(0))(σ−σ(0)) for all σ, σ(0)∈Rn+, (17) and, for allt∈ [0,1],
F ((1−t )σ(0)+t σ )(1−t )F (σ(0))+t F (σ ). (18) (e) F(σ )F (σ ).
Proof (a)–(d) follow from applying Lemma2on the Hilbert spaceV. (e) was proven
in Lemma4.
5 Numerical Examples and Inverse Problem Challenges
In this section, we will show some numerical results for the stationary diffusion ex- ample from Sect.2.1and demonstrate some major challenges that arise in solving
Fig. 4 Single measurement in the top right pixel for a source term in the lower left pixel as a function of changing the diffusivity in each of the 3×3 pixels
the inverse coefficient problem of recoveringσˆ ∈Rn fromF (σ )ˆ ∈Rm, or from a noisy versionYδ≈F (σ ). The source codes for the following examples (and also forˆ generating Fig.1and2) are given in the appendix for the reader’s reference.
5.1 Non-uniqueness
Even for m≥n, and a noise-free measurement Yˆ =F (σ )ˆ ∈Rm, it is not clear whether the measurements uniquely determine the unknownσˆ ∈Rn. To demonstrate this on a simple one-dimensional example, let us consider the stationary diffusion example with 3×3 pixels and circular excitation/measurement subdomains in each boundary pixel as in Fig.1. We apply a source term inD1in the lower left pixel, and measure the resulting total concentration in D8 in the top right pixel, so that l=χ1∈H−1()andr=χ8∈H−1(), where we writeχj:=χDj for the ease of notation. We chooseσ=1 in all pixels exceptPi, and onPi we vary the diffusivity in steps of 0.01 up to 3. Figure4showsFl,r(σ )for alli=1, . . . ,9, in the same or- der as the pixels, e.g., the lower left image showsFl,r(σ )forσ =(σ1,1, . . . ,1)for varyingσ1.
Intuitively speaking, one can see that rising the diffusivity in the middle pixel increases the measurement since particles can easier diffuse through the middle pixel
