a priori error estimates for the finite element discretization of elliptic ...
A PRIORI ERROR ESTIMATES FOR THE FINITE ELEMENT DISCRETIZATION OF ELLIPTIC PARAMETER IDENTIFICATION PROBLEMS WITH POINTWISE MEASUREMENTS R. RANNACHER† ‡ AND B. VEXLER† ‡ Abstract. We develop an a priori error analysis for the finite element Galerkin discretization of parameter identification problems. The state equation is given by an elliptic partial differential equation of second order with a finite number of unknown parameters, which are estimated using point-wise measurements of the state variable. Key words. Parameter identification, finite elements, pointwise measurements, L∞ -error estimates. AMS subject classifications. 65K10, 65N30, 65N21, 49M25, 49K20.
1. Introduction. We consider parameter identification problems governed by an elliptic partial differential equation of second order. The finitely many unknown parameters are estimated using the measurements of point values of the state variable. Let Ω ⊂ R2 be a convex polygonal domain, L2 (Ω) the corresponding Lebesgue space with inner product and norm denoted by (·, ·) and k · k2 , respectively, and H m (Ω) the Sobolev space of order m ∈ N. With this notation, we set V := {v ∈ H 1 (Ω) ∩ C(Ω) | v = 0 on ∂Ω}. The state variable u ∈ V is determined by an elliptic partial differential equation (the state equation) −∇ · (A(q)∇u) = f u=0
in Ω on ∂Ω,
(1.1)
for given H¨older continuous f ∈ C α (Ω), α ∈ (0, 1) . Here, q ∈ Q0 ⊂ Q = Rnp denotes the set of unknown parameters and A(q) = (Aij (q)) is a symmetric and positive definite 2×2 matrix with twice continuously differentiable entries Aij : Q → C 1+α (Ω). The above conditions guarantee that for any admissible value of the parameter q , the corresponding solution u of the state equation (1.1) is in H 2 (Ω) . At the corner points of ∂Ω , the second derivatives of the solution may become singular. However, u has H¨older continuous second derivatives, u ∈ C 2+α (Ωd ), for each subdomain Ωd ⊂ Ω with positive distance to the corner points. The usual weak formulation of (1.1) is: a(q)(u, φ) = (f, φ) ∀φ ∈ V,
(1.2)
where the bilinear form a(q)(·, ·) is defined by a(q)(u, φ) := (A(q)∇u, ∇φ). † Institute
(1.3)
of Applied Mathematics, University of Heidelberg, INF 294, 69120 Heidelberg, Germany. ‡ This work has been supported by the German Research Foundation (DFG), through the Sonderforschungsbereich 359 ’Reactive Flows, Diffusion and Transport’, and the Graduiertenkolleg ’Complex Processes: Modeling, Simulation and Optimization’ at the Interdisciplinary Center of Scientific Computing (IWR), University of Heidelberg. 1
2
R. RANNACHER AND B. VEXLER
Further, the observation operator C(·) describing the mapping of the state variable u to the space of measurements Z = Rnm is given by Ci (v) = v(ξi ),
i = 1, 2, . . . , nm ,
(1.4)
where {ξi } ⊂ Ω is a finite set of measurement points. We assume that nm ≥ np . The euclidian product and norm on Q and Z are denoted by h·, ·i and k · k , respectively, and the same notation is also used for the corresponding natural norms of matrices. The values of the parameters are estimated from a given set of measurements Cˆ ∈ Z using a least squares approach such that we obtain a constrained optimization problem with the cost functional J : V → R: Minimize
ˆ 2, J(u) := 21 kC(u) − Ck
(1.5)
under the constraint (1.2). Throughout we assume the existence of a solution (u, q) ∈ V ×Q0 of the problem (1.2,1.5). For an analysis of existence of solutions for parameter identification problems see, e.g., Banks & Kunisch [2], Kravaris & Seinfeld [14] and Litvinov [15]. The state equation is discretized by a conforming finite element Galerkin method defined on a family {Th }h>0 of shape regular quasi-uniform meshes Th = {K} consisting of closed cells K which are either triangles or quadrilaterals. The straight parts which make up the boundary ∂K of a cell K are called faces. The mesh parameter h is defined as a cell-wise constant function by setting h|K = hK and hK is the diameter of K. Usually we use the symbol h also for the maximal cell size, i.e. h = max hK . K∈Th
(1.6)
For convenience, we assume that 0 < h < 1 . On the mesh Th we define finite element spaces Vh ⊂ V consisting of linear or bilinear shape functions, see e.g. Brenner & Scott [5] or Johnson [12]. The corresponding discrete state uh ∈ Vh and parameter qh ∈ Q0 are determined by: Minimize
J(uh ),
(1.7)
under the constraint a(qh )(uh , φh ) = (f, φh ) ∀φh ∈ Vh .
(1.8)
Since Q is finite dimensional, the parameter qh is determined in the same space Q. The main purpose of this paper is to analyze the behavior of the error in parameters kq − qh k for h tending to zero. There is a number of publications, where a priori error estimates are derived for optimal control problems governed by partial differential equations, see Falk [7], Arada, Casas & Tr¨oltzsch [1], Deckelnic & Hinze [6] and Gunzburger & Hou [11]. However, there are only few published results on this topic in the context of parameter identification problems, see Falk [8], Neittaanm¨aki & Tai [16] and K¨ arkk¨ ainen [13]. In Vexler [21], an a priori error analysis for the case of V -stable observation operators Ci (·) is developed and optimal-order convergence is shown, kq − qh k = O(h2 ),
(1.9)
FINITE ELEMENTS IN PARAMETER IDENTIFICATION
3
essentially under the assumption that |Ci (u) − Ci (uh )| ≤ ch2 kukH 2 . However, the generalization of this result to point-wise observations is not straightforward. For this case, we prove in this paper that, under certain regularity conditions, kq − qh k = O(h2 | log(h)|2 ).
(1.10)
The proof uses the technique for estimating discrete Green functions developed in Frehse & Rannacher [9]. A complementary result of a posteriori error analysis for parameter identification problems is given in Becker & Vexler [4]. To the authors knowledge, this is the first a priori error analysis for parameter identification problems with point-wise observations. The consideration of point-wise observations in determining discrete parameters seems very natural in view of practical measurement techniques; see [3] for applications in reactive flow analysis. The paper is organized as follows: In the next section, we describe an algorithm for solving problem (1.7, 1.8). In Section 3, we present a paradigm for a priori error analysis for discretization of a class of optimization problems. Thereafter, in Section 4, we derive the announced error estimate using an L∞ -stability result, which is proven in Section 5. In Section 6, we present a numerical example conforming the asymptotic sharpness of our error estimate. Possible extensions are addressed in the last section. 2. Optimization algorithm. In this section, we reformulate the problem under consideration as an unconstrained optimization problem and describe a solution algorithm for it. Throughout, we assume the coefficient matrix A(q) to be positive definite at the solution of the problem. By virtue of the implicit function theorem for Banach spaces, this implies the existence of an open (bounded) set Q0 ⊂ Q containing the optimal parameter q , and of a continuously differentiable solution operator S : Q0 → V satisfying: a(q)(S(q), φ) = (f, φ) ∀φ ∈ V
∀q ∈ Q0 .
(2.1)
Without loss of generality, we assume that A(q) is positive definite for q ∈ Q0 , i.e. there exists γ ∈ R+ , such that p∗ A(q)p ≥ γkpk2 ,
∀p ∈ Q ∀q ∈ Q0 .
(2.2)
The convexity of the domain Ω implies the following regularity property (see e.g. Grisvard [10]): S(q) ∈ V ∩ H 2 (Ω) ∀q ∈ Q0 .
(2.3)
We introduce the reduced observation operator c : Q0 → Z by c(q) := C(S(q)).
(2.4)
This allows us to reformulate the problem under consideration as an unconstrained optimization problem with the reduced cost functional j : Q0 → R : Minimize
ˆ 2, j(q) := 21 kc(q) − Ck
q ∈ Q0 .
(2.5)
Denoting by G = c0 (q) ∈ Rnp ×nm the Jacobian matrix of the reduced observation operator c(·), the first-order necessary optimality condition j 0 (q) = 0 for (2.5) reads ˆ = 0, G∗ (c(q) − C)
(2.6)
4
R. RANNACHER AND B. VEXLER
where G∗ denotes the transpose of G. The positive semi-definiteness of the Hessian matrix H := ∇2 j(q) is the second-order necessary optimality condition. A solution q of problem (2.5) is called stable, if the sufficient optimality condition holds, i.e. if the Hessian H is positive definite. Throughout, we will assume the solution q to be stable. The stability of the solution is given, for instance, if the value of the ˆ is small enough and the matrix G has full rank np , see cost functional kC(u) − Ck e.g. [21] for details. Since by assumption the matrix coefficient A(·) is twice continuously differentiable, there holds sup |||A(ξ)|||1,∞ + sup |||A0qj (ξ)|||1,∞ + sup |||A00qj qk (ξ)|||1,∞ < ∞,
ξ∈Q0
ξ∈Q0
(2.7)
ξ∈Q0
where |||B|||1,∞ := maxi,j=1,2 kBij k1,∞ for a matrix function B = (Bij ) ∈ C 1 (Ω)2×2 . In the following propositions, we give representations of the Jacobian G of c(·) , the Hessian H of j(·) , and the Hessian of ci (·) . Theorem 2.1. Let the reduced observation operator c(·) and the reduced functional j(·) be defined as in (2.4) and (2.5), respectively. Then, the elements of the Jacobian of c(·) at some q ∈ Q0 are given by Gij =
∂ci (q) = Ci (wj ), ∂qj
i = 1, . . . , nm , j = 1, . . . , np ,
(2.8)
where wj ∈ V are the solutions of the problems a(q)(wj , φ) = −(A0qj (q)∇u, ∇φ)
∀φ ∈ V,
(2.9)
with u = S(q) . The Hessian of j(·) can be expressed by H = G∗ G + M,
(2.10)
where the matrix M ∈ Rnp ×np is given by M=
nm X
c00i (q)(ci (q) − Cˆi ).
(2.11)
i=1
The Hessian of ci (q) is given by ∂2 ci (q) = Ci (vjk ), ∂qj ∂qk
(2.12)
where the vjk ∈ V are the solutions of the problems a(q)(vjk , φ) = −(A0qj (q)∇wk , ∇φ) − (A0qk (q)∇wj , ∇φ) − (A00qj qk (q)∇u, ∇φ)
∀φ ∈ V,
(2.13)
with wj as defined in (2.9). Proof. The derivation of the derivatives of c(·) is straightforward using the implicit function theorem and the chain rule. In practice the Hessian H of j(·) is computed using the representation Mjk = −(A0qj (q)∇wk , ∇z) − (A0qk (q)∇wj , ∇z) − (A00qj qk (q)∇u, ∇z),
(2.14)
FINITE ELEMENTS IN PARAMETER IDENTIFICATION
5
with the function z determined by the dual equation ¯ C(φ)i. (A(q)∇φ, ∇z) = hC(u) − C,
(2.15)
Similar to the continuous case, we introduce a discrete solution operator Sh : Q0 → Vh by the equation a(qh )(Sh (qh ), φh ) = (f, φh ) ∀φh ∈ Vh , qh ∈ Q0 .
(2.16)
As before, we turn the discrete problem (1.7, 1.8) into an unconstrained minimization problem: Minimize
ˆ 2, jh (qh ) := 21 kch (qh ) − Ck
qh ∈ Q0 ,
(2.17)
where the discrete reduced observation operator ch is defined by ch (qh ) = C(Sh (q)).
(2.18)
Denoting the corresponding Jacobian by Gh = c0h (qh ), the necessary optimality condition jh0 (qh ) = 0 reads ˆ = 0. G∗h (ch (qh ) − C)
(2.19)
The derivatives of the discrete observation operator ch can be computed in an analogous way as in Theorem 2.1. Problem (2.17) is solved iteratively starting with an initial guess qh0 and using the recursive setting qhk+1 = qhk + δqh . The update δqh is obtained as the solution of the system of linear equations Hk δqh = G∗h (Cˆ − ch (qhk )),
(2.20)
where Gh = c0h (qhk ) , and Hk is an appropriate symmetric approximation of the Hessian ∇2 jh (qhk ) . The most widely used choice of the matrix Hk = G∗h Gh leads us to the Gauß-Newton algorithm, see, e.g., Nocedal & Wright [19]. For one step of the Gauß-Newton algorithm the state equation and np tangent problems (2.9) have to be solved which involves the same linear operator but with different right-hand sides. Due to the small dimension np of the parameter space Q the solution of (2.20) is uncritical. For discussing other Newton type methods and trust-region techniques for globalization of the convergence in this context, see [21]. 3. A paradigm for a priori error analysis. In this section we present a general approach to the error analysis of a class of optimization problems such as considered in this paper. The main result is stated in the following theorem. It is a variant of well-known perturbation theorems for differentiable mappings which is particularly tailored to the present situation. However, it seems easier to include the elementary proof than to search for the precise reference. Theorem 3.1. Let F, Fh : Rn → Rn , for a discretization parameter h ∈ R+ , be continuously differentiable operators , and x ∈ Rn be a solution of F (x) = 0 . Let the following conditions be fulfilled: (i) The derivative F 0 (x) is positive definite, i.e., there is a constant γ > 0 , such that p∗ F 0 (x)p ≥ γkpk2 ,
p ∈ Rn .
(3.1)
6
R. RANNACHER AND B. VEXLER
(ii) There is a neighborhood U of x and a positive number L(h) ∈ R+ , such that kFh0 (ξ) − Fh0 (η)k ≤ L(h)kξ − ηk
∀ξ, η ∈ U.
(3.2)
(iii) With the h-dependent constant L(h) , there holds lim L(h)kF (x) − Fh (x)k = 0.
(3.3)
lim kF 0 (x) − Fh0 (x)k = 0.
(3.4)
h→0
(iv) There holds h→0
Then, for h small enough, there exists xh ∈ U , such that Fh (xh ) = 0 , and Fh0 (xh ) is positive definite uniformly in h . Further, there holds the a priori error estimate kx − xh k ≤
2 kF (x) − Fh (x)k. γ
(3.5)
Proof. Due to condition (iv), we can choose a positive number h1 ∈ R+ , such that for h ≤ h1 there holds kF 0 (x) − Fh0 (x)k ≤ 41 γ . Moreover, for ρ = ρ(h) =
γ kL(h)
(3.6)
, with some k ≥ 4 sufficiently large, there holds
Bρ (x) = {ξ ∈ Rn , kx − ξk ≤ ρ} ⊂ U.
(3.7)
For this choice, we obtain that, for h ≤ h1 , Fh0 (·) is positive definite on Bρ (x) : p∗ Fh0 (ξ)p = p∗ F 0 (x)p + p∗ (Fh0 (x) − F 0 (x))p + p∗ (Fh0 (ξ) − Fh0 (x))p ≥ γkpk2 − kFh0 (x) − F 0 (x)k kpk2 − kFh0 (ξ) − Fh0 (x)k kpk2 ≥ (γ − 41 γ − L(h)ρ)kpk2 ≥ 12 γkpk2 . In a similar way, we conclude that, for h ≤ h1 , Fh0 (·) is also bounded on Bρ (x): kFh0 (ξ)k ≤ β := kF 0 (x)k + 21 γ . Next, we prove, that there exists a unique xh ∈ Bρ (x) with Fh (xh ) = 0 . To this end, we define an operator Ds : Rn → Rn , for s ∈ R+ , by Ds (ξ) = ξ − sFh (ξ). For a certain choice of s , we show that Ds is a contraction on Bρ (x) and use the Banach fixed point theorem. For ξ ∈ Bρ (x), h ≤ h1 , and an arbitrary p ∈ Rn , there holds kDs0 (ξ)pk2 = kp − sFh0 (ξ)pk2 = kpk2 − 2sp∗ Fh0 (ξ)p + s2 kFh0 (ξ)pk2 ≤ (1 − sγ + s2 β 2 )kpk2 . For the choice s = γ(2β 2 )−1 , we obtain γ2 kDs0 (ξ)pk2 ≤ 1 − 2 kpk2 , 4β
FINITE ELEMENTS IN PARAMETER IDENTIFICATION
7
and consequently, γ 2 1/2 kDs0 (ξ)k ≤ 1 − 2 < 1. 4β Moreover, for arbitrary ξ ∈ Bρ (x), there holds: kx − Ds (ξ)k = kDs (x) − Ds (ξ) + sFh (x)k ≤ kDs (x) − Ds (ξ)k + skFh (x) − F (x)k ≤ kDs0 (η)k kx − ξk + skFh (x) − F (x)k, for a certain η ∈ Bρ . Hence, the above estimate implies γ 2 1/2 kx − Ds (ξ)k ≤ 1 − 2 ρ + s kFh (x) − F (x)k 4β γ 2 1/2 k =ρ 1− 2 + s L(h)kFh (x) − F (x)k . 4β γ Due to condition (iii), there is a number h2 ∈ R+ , such that, for h ≤ h2 , there holds γ 2 1/2 γ 1− 1− 2 L(h)kFh (x) − F (x)k ≤ . ks 4β Hence, for h < h0 := min{h1 , h2 }, kx − Ds (ξ)k ≤ ρ, and consequently Ds (ξ) ∈ Bρ (x) . For h ≤ h0 , by the Banach fixed point theorem, we obtain the existence of xh ∈ Bρ (x) with Fh (xh ) = 0. By construction of Bρ (x) the derivative Fh0 (xh ) is positive definite with the h-independent constant 21 γ . This implies that, for a certain ξ ∈ Bρ (x), (x − xh )∗ (Fh (x) − Fh (xh )) = (x − xh )∗ Fh0 (ξ)(x − xh ) ≥
γ kx − xh k2 . 2
Hence, using F (x) = Fh (xh ) = 0, kx − xh k2 ≤
2 2 (x − xh )∗ (Fh (x) − Fh (xh )) = (x − xh )∗ (Fh (x) − F (x)) γ γ 2 ≤ kFh (x) − F (x)k kx − xh k. γ
This completes the proof. 4. A priori error estimation. In this section we apply the paradigm presented in Section 3 to the problem under consideration. We prove the following theorem. Theorem 4.1. Let q ∈ Q be a stable solution of (2.5). Then, for h small enough, there exist a stable solution qh ∈ Q of (2.17), and there holds the following a priori error estimate: kq − qh k = O(h2 | log(h)|2 ).
(4.1)
8
R. RANNACHER AND B. VEXLER
On the basis of the estimate (4.1), we can also derive optimal-order estimates for the error u − uh in the corresponding states. However, since this would be a simple exercise using the arguments developed below, and since the optimal states are of only minor practical interest in parameter estimation problems, we do not state these estimates. The proof of Theorem 4.1 is given by checking the conditions from Theorem 3.1 for the operators F (ξ) := ∇j(ξ)
Fh (ξ) := ∇jh (ξ).
The constant in (4.1) turns out to depend in a reciprocal way on the distance δ :=
min
i=1,...,nm
dist(ξi , Σ)
of the set of measurement points ξi to the set Σ of corner points of ∂Ω. Therefore, we will use generic constants c and cδ , where c only depends on the domain Ω, the force f , and the characteristics of the mesh family {Th }h , while cδ may additionally depend on the distance δ like cδ ≈ δ −1 . Further, by Lp (Ω) and W m,p (Ω) , for m ∈ N and 1 ≤ p ≤ ∞ , we denote the standard Lebesgue and Sobolev spaces and by k · kp and k · km,p the corresponding norms. The restriction of such a norm to a subset Ω0 ⊂ Ω is indicated by k · km,p;Ω0 . An important ingredient of the proof of Theorem 4.1 is the following L∞ -stability theory. For d > 0 , we define the subset Ωd ⊂ Ω by Ωd := {x ∈ Ω, dist(x, Σ) > d}. ¯ , Theorem 4.2 (Stability theorem). Let be given q ∈ Q0 , ψ ∈ H01 (Ω) ∩ C(Ω) 1,∞ 2×2 and a matrix B = B(x) ∈ W (Ω) . Moreover, let vh ∈ Vh be a solution of a(q)(vh , φh ) = (B∇ψ, ∇φh )
∀φh ∈ Vh .
(4.2)
Then, there hold the L2 -stability estimate kvh k2 + hk∇vh k2 ≤ c |||B|||1,∞ kψk2 + hk∇ψk2 , and the local L∞ -stability estimate kvh k∞;Ωd ≤ cd |||B|||1,∞ | log(h)| kψk∞;Ωd/2 + kψk2 + hk∇ψk2 ,
(4.3)
(4.4)
with a constant cd ≈ d−1 . The L2 estimate (4.3) is a standard result from finite element analysis, while the L∞ estimate (4.4) follows by estimates of discrete Green functions such as developed in Frehse & Rannacher [9]. A similar L∞ -stability result has been proven in Rannacher [20] in the time-dependent parabolic case. The proof is given in Section 5. For the solution q of problem (2.5), we introduce u ¯h ∈ Vh determined by a(q)(¯ uh , φh ) = (f, φh ) ∀φh ∈ Vh .
(4.5)
Further, we define wj,h ∈ Vh and vjk,h ∈ Vh , for j, k = 1, 2, . . . , np , as the solutions of the problems a(q)(wj,h , φh ) = −(A0qj (q)∇¯ uh , ∇φh ) ∀φh ∈ Vh ,
(4.6)
9
FINITE ELEMENTS IN PARAMETER IDENTIFICATION
and a(q)(vjk,h , φh ) = −(A0qj (q)∇wk,h , ∇φh ) − (A0qk (q)∇wj,h , ∇φh )− (A00qj qk (q)∇¯ uh , ∇φh ) ∀φh ∈ Vh ,
(4.7)
respectively. The next lemma provides necessary estimates for the errors u − u ¯h , wj − wj,h and vjk − vjk,h . We recall the notation δ := mini=1,...,nm dist(ξi , Σ) . Lemma 4.3. Under the above assumptions the following estimates hold: kC(u − u ¯h )k ≤ cδ h2 | log(h)|, 2
2
kC(wj − wj,h )k ≤ cδ h | log(h)| , kC(vjk − vjk,h )k ≤ cδ h2 | log(h)|3 ,
(4.8) j = 1, 2 . . . np , j, k = 1, 2 . . . np .
(4.9) (4.10)
Proof. (i) The proof uses some a priori bounds for u and the auxiliary functions wj , vjk , j, k = 1, . . . , np , corresponding to arbitrary q ∈ Q0 . By our regularity assumptions, we have the following a priori estimates: kuk2,2 ≤ c, kwj k2,2 ≤ ckuk2,2 ≤ c, kvjk k2,2 ≤ c max kwj k2,2 , kwk k2,2 + ckuk2,2 ≤ c,
(4.11) (4.12) (4.13)
where c is a generic constant depending only on the data of the problem. Further, by elliptic regularity theory, for d > 0 , there hold the local a priori estimates (4.14) kukC 2+α (Ωd ) ≤ cd kf kC α (Ωd/2 ) + kuk2,2 ≤ cd , kwj kC 2+α (Ωd ) ≤ cd kukC 2+α (Ωd/2 ) + kuk2,2 ≤ cd , (4.15) (4.16) kvjk kC 2+α (Ωd ) ≤ cd kwj kC 2+α (Ωd/2 ) + kwj k2,2 ≤ cd , with a generic constant cd ≈ d−1 . Here, k · kC 2+α (Ωd ) denotes the norm of the space C 2+α (Ωd ) of functions possessing H¨ older continuous second derivatives on Ωd . (ii) By definition, u ¯h is the Ritz projection of u corresponding to the energy form a(q)(·, ·) , i.e., a(q)(¯ uh , φh ) = a(q)(u, φh ) ∀φh ∈ Vh . By the standard L2 -error estimate for finite elements, there holds ku − u ¯h k2 + hk∇(u − u ¯h )k2 ≤ ch2 .
(4.17)
Further, let ih u ∈ Vh be the usual nodal interpolant of u . Then, applying the L∞ -stability estimate (4.4) of Theorem 4.2 for the equation a(q)(ih u − u ¯h , φh ) = a(q)(ih u − u, φh ) ∀φh ∈ Vh , yields the estimate kih u − u ¯h k∞;Ωδ ≤ cδ | log(h)| kih u − uk∞;Ωδ/2 + kih u − uk2 + hk∇(ih u − u)k2 .
10
R. RANNACHER AND B. VEXLER
In view of the approximation properties of ih and the a priori bounds (4.11) and (4.14), we conclude the error estimate ku − u ¯h k∞;Ωδ ≤ cδ h2 | log(h)|.
(4.18)
Since ξi ∈ Ωδ , we obtain the estimate (4.8). (iii) For proving (4.9), we introduce an additional discrete variable w ¯j,h determined by the equation a(q)(w ¯j,h , φh ) = −(A0qj (q)∇u, ∇φh ) ∀φh ∈ Vh . The error e = wj − wj,h is split like e = e1 + e2 , with e1 = wj − w ¯j,h and e2 = w ¯j,h − wj,h . For the Ritz-projection error e1 , similar as before, there holds the L2 error estimate ke1 k2 + hk∇e1 k2 ≤ ch2 kwj k2,2 ≤ ch2 , and the pointwise error estimate ke1 k∞;Ωδ ≤ cδ h2 | log(h)| k∇2 wj k∞;Ωδ/2 + kwj k2,2 ≤ cδ h2 | log(h)|. For e2 ∈ Vh , we have a(q)(e2 , φh ) = −(A0qj (q)∇(u − u ¯h ), ∇φh ) ∀φh ∈ Vh . Hence, the L2 -stability estimate (4.3) of Theorem 4.2 and the estimate (4.17) imply ke2 k2 + hk∇e2 k2 ≤ c ku − u ¯h k2 + hk∇(u − u ¯h )k2 ≤ ch2 . This shows that, for j = 1, . . . , np , kwj − wj,h k2 + hk∇(wj − wj,h )k2 ≤ ch2 .
(4.19)
Further, the L∞ -stability estimate (4.4) of Theorem 4.2 yields ke2 k∞;Ωδ ≤ cδ | log(h)| ku − u ¯h k∞;Ωδ/2 + ku − u ¯h k2 + k∇(u − u ¯h )k2 . which, by (4.17) and (4.18), implies ke2 k∞;Ωδ ≤ cδ | log(h)|2 h2 . We obtain kwj − wj,h k∞;Ωδ ≤ cδ | log(h)|2 h2 ,
j = 1, . . . , np ,
(4.20)
which implies the desired estimate (4.9). (iv) The proof of (4.10) uses the same line of argument as before. Using the additional discrete variable v¯jk,h determined by the equation a(q)(¯ vjk,h , φh ) = −(A0qj (q)∇wk , ∇φh ) − (A0qk (q)∇wj , ∇φh ) − (A00qj qk (q)∇u, ∇φh ) ∀φh ∈ Vh , the error e = vjk − vjk,h is split like e = e1 + e2 , with e1 = vjk − v¯jk,h and e2 = v¯jk,h − vjk,h . For the Ritz-projection error e1 , similar as before, we conclude the the pointwise error estimate ke1 k∞;Ωδ ≤ cδ h2 | log(h)| k∇2 vjk k∞;Ωδ/2 + kvjk k2,2 ≤ cδ h2 | log(h)|.
11
FINITE ELEMENTS IN PARAMETER IDENTIFICATION
For e2 ∈ Vh , we have a(q)(e2 , φh ) = −(A0qj (q)∇(wk − wk,h ), ∇φh ) − (A0qk (q)∇(wj − wj,h ), ∇φh ) − (A00qj qk (q)∇(u − u ¯h ), ∇φh ) ∀φh ∈ Vh , and therefore, again by the L∞ -stability estimate (4.4) of Theorem 4.2, ke2 k∞;Ωδ ≤ cδ | log(h)| max kwj − wj,h k∞;Ωδ/2 + ku − u ¯h k∞;Ωδ/2 j=1,...,np kwj − wj,h k2 + hk∇(wj − wj,h )k2 + c max j=1,...,np + c ku − u ¯h k2 + hk∇(u − u ¯h )k2 . Then, by the foregoing error estimates, we obtain ke2 k∞;Ωδ ≤ cδ h2 | log(h)|3 , and consequently, kvjk − vjk,h k∞;Ωδ ≤ cδ | log(h)|3 h2 ,
j, k = 1, . . . , np .
(4.21)
This eventually yields the desired estimated (4.10). A direct application of Lemma 4.3 leads to the following result. Lemma 4.4. Under the above assumptions, there holds ∂ (j − jh )(q) ≤ cδ h2 | log(h)|2 , j = 1, 2, . . . , np , ∂qj ∂2 (j − jh )(q) ≤ cδ h2 | log(h)|3 , j, k = 1, 2, . . . , np . ∂qj ∂qk
(4.22) (4.23)
Proof. We have the representation ∂ ˆ C(wj )i − hC(¯ ˆ C(wj,h )i (j − jh )(q) = hC(u) − C, uh ) − C, ∂qj ˆ C(wj − wj,h )i + hC(u − u = hC(u) − C, ¯h ), C(wj,h )i, from which we obtain ∂ ˆ kC(wj − wj,h )k + kC(u − u (j − jh )(q) ≤ kC(u) − Ck ¯h )k kC(wj,h )k. ∂qj By the a priori bounds (4.11)-(4.13) and the Sobolev imbedding theorem, we see that kC(u)k + kC(wj )k + kC(vjk )k ≤ c.
(4.24)
Combining this with the error estimate (4.9) implies kC(wj,h )k ≤ c . Then, we can conclude the first estimate (4.22) from the error estimates of Lemma 4.3. To prove (4.23), we write ∂2 ˆ C(vjk )i (j − jh )(q) = hC(wj ), C(wk )i + hC(u) − C, ∂qj qk ˆ C(vjk,h )i − hC(wj,h ), C(wk,h )i − hC(¯ uh ) − C, = hC(wj − wj,h ), C(wk )i + hC(wj,h ), C(wk − wk,h )i ˆ C(vjk − vjk,h )i. + hC(u − u ¯h ), C(vjk )i + hC(¯ uh ) − C,
12
R. RANNACHER AND B. VEXLER
Using as before the bounds (4.24) and the error estimates of Lemma 4.3 completes the proof. For the application of Theorem 3.1 it remains to check the Lipschitz condition (3.2). For two arbitrary parameter sets ξ, η ∈ Q0 , we set uξ = Sh (ξ) and uη = Sh (η) . Correspondingly, we define wj,ξ , wj,η ∈ Vh and vjk,ξ , vjk,η ∈ Vh similarly to wj,h and vj,h for q = ξ and q = η , respectively. Lemma 4.5. For ξ, η ∈ Q0 , there holds kC(uξ − uη )k ≤ cδ | log(h)| kξ − ηk,
(4.25)
2
(4.26)
3
(4.27)
kC(wj,ξ − wj,η )k ≤ cδ | log(h)| kξ − ηk, kC(vjk,ξ − vjk,η )k ≤ cδ | log(h)| kξ − ηk. Proof. Due to the definition of uξ and uη , we have (A(ξ)∇(uξ − uη ), ∇φh ) = −((A(ξ) − A(η))∇uη , ∇φh ) ∀φh ∈ Vh .
Using Theorem 4.2, with d = δ , we obtain: kC(uξ − uη )k ≤ c |||(A(ξ) − A(η)|||1,∞ | log(h)| kuη k∞;Ωδ/2 + kuη k2 + hk∇uη k2 . Since uη is the Ritz projection of an H 2 function, all its norms occurring on the right-hand side are bounded independent of h and η ∈ Q0 by standard estimates from finite element analysis. This implies (4.25) since |||A(ξ) − A(η)|||1,∞ ≤ ckξ − ηk. The estimates (4.26) and (4.27) are obtained in the similar way. Lemma 4.6. For ξ, η ∈ Q0 , there holds: ∂2 ∂2 jh (ξ) − jh (η) ≤ L(h)kξ − ηk, ∂qj qk ∂qj qk where L(h) = cδ | log(h)|3 . Proof. We have ∂2 ∂2 ˆ C(vjk,ξ )i jh (ξ) − jh (η) = hC(wj,ξ ), C(wk,ξ )i + hC(uξ ) − C, ∂qj qk ∂qj qk ˆ C(vjk,η )i − hC(wj,η ), C(wk,η )i − hC(uη ) − C, = hC(wj,ξ − wj,η ), C(wk,ξ )i + hC(wj,η ), C(wk,ξ − wk,η )i ˆ C(vjk,ξ − vjk,η )i, + hC(uξ − uη ), C(vjk,ξ )i + hC(uη ) − C, and, consequently, ∂2 ∂2 jh (ξ) − jh (η) ≤ kC(wj,ξ − wj,η )k kC(wk,ξ )k ∂qj qk ∂qj qk + kC(wj,η )k kC(wk,ξ − wk,η )k + kC(uξ − uη )k kC(vjk,ξ )k ˆ kC(vjk,ξ − vjk,η )k. + kC(uη ) − Ck
(4.28)
13
FINITE ELEMENTS IN PARAMETER IDENTIFICATION
Now, the assertion follows by the estimates of Lemma 4.5 if we can bound the terms C(uη ), C(wj,η ), and C(vjk,ξ ) . This is achieved by using the bounds for C(u), C(wj ), and C(vjk ) in (4.24) together with the error estimates of Lemma 4.3. To complete the proof of Theorem 4.1 we check the conditions of Theorem 3.1. Condition (3.1) is fulfilled due to the stability of the solution q of the problem (2.5). Condition (3.2) is shown in Lemma 4.6. Condition (3.3) is obtained by Lemma 4.4 and Lemma 4.6 using limh→0 h2 | log(h)|5 = 0 . Finally, condition (3.4) holds due to Lemma 4.4. Hence, the estimate (3.5) of Theorem 3.1 completes the proof. 5. Proo f of Theorem 4.2. (i) We begin with the L2 -stability estimate. Taking φh := vh in (4.2), we obtain k∇vh k2 ≤ c|||B|||1,∞ k∇ψk2 .
(5.1)
To estimate kvh k2 , we use the solution z ∈ V ∩ H 2 (Ω) of the auxiliary equation a(q)(φ, z) = (φ, vh )kvh k−1 2
∀φ ∈ V.
Then, taking φ := vh and using the approximation properties of the interpolant ih z ∈ Vh , we conclude that kvh k2 = a(q)(vh , z) = a(q)(vh , z − ih z) + a(q)(vh , ih z) = a(q)(vh , z − ih z) + (B∇ψ, ∇ih z) = a(q)(vh , z − ih z) + (B∇ψ, ∇(ih z − z)) + (B∇ψ, ∇z) ≤ ck∇vh k2 k∇(z − ih z)k2 + |||B|||1,∞ k∇ψk2 k∇(z − ih z)k2 + |||B|||1,∞ kψk2 kzk2,2 ≤ c hk∇vh k + |||B|||1,∞ hk∇ψk2 + |||B|||1,∞ kψk2 kzk2,2 . Hence, observing (5.1) and the bound kzk2,2 ≤ c , we obtain kvh k2 ≤ c|||B|||1,∞ kψk2 + hk∇ψk2 .
(5.2)
(ii) Next, we prove the L∞ -stability estimate. Let a ∈ Ωδ be an arbitrary point lying in a cell K . For any fixed h , there exists a cell-wise polynomial function δh with supp(δh ) ⊂ K, such that (φh , δh ) = φh (a) ∀φh ∈ Vh . The function δh plays the role of an approximate Dirac function. Correspondingly, we introduce a regularized Green function g ∈ V ∩ H 2 (Ω) by (A(q)∇φ, ∇g) = (δh , φ) ∀φ ∈ V, and the corresponding Ritz projection gh ∈ Vh by (A(q)∇φh , ∇gh ) = (δh , φh ) ∀φh ∈ Vh . By ih : C(Ω) → Vh , we denote the usual (linear) operator of nodal interpolation for which the following cell-wise estimate is well known: −1 2 0 2 h−2 T kv − ih vkp;T + hT k∇(v − ih v)kp:T + k∇ ih vkp;T ≤ ck∇ vkp;T ,
(5.3)
14
R. RANNACHER AND B. VEXLER
for 1 ≤ p ≤ ∞ , with constants c independent P of h . For functions which are only cell-wise defined, we will use the norm kvk0p := K∈Th kvkp;T . The following three lemmas provide the key estimates for the proof of the theorem. Lemma 5.1. The following global L2 estimates hold: kgk2 + | log(h)|−1/2 k∇gk2 + hk∇2 gk2 ≤ c, h
−1
kg − gh k2 + k∇(g − gh )k2 + hk∇
2
gh k02
≤ c.
(5.4) (5.5)
Proof. The assertion follows by standard L2 a priori and error estimates for g and g − gh , respectively. We skip the details and refer to [9]. Note that k∇2 gh k02 vanishes for linear finite elements. In the case of bilinear elements, we estimate as follows: k∇2 gh k02 ≤ k∇2 (gh − ih g)k02 + k∇2 (g − ih g)k02 + k∇2 gk2 . For the first term, we obtain, using an inverse inequality, k∇2 (gh − ih g)k02 ≤ ch−1 k∇(gh − ih g)k2 ≤ ch−1 k∇(g − ih g)k2 + k∇(g − gh )k2 and obtain by the interpolation estimate (5.3), with p = 2 , and by the other estimates derived before: k∇2 gh k02 ≤ c h−1 k∇(g − gh )k2 + k∇2 gk2 ≤ ch−1 . This completes the proof. Lemma 5.2. For sufficiently small h δ , the following local L2 estimate holds: k∇(g − gh )k2;Ω\Ωδ/2 + hk∇2 gk2;Ω\Ωδ/2 ≤ cδ h,
(5.6)
with a constant cδ ≈ δ −1 but independent of h . Proof. The assertion follows by standard local elliptic a priori estimates and by arguments from the local L2 error analysis for finite elements as provided in Nitsche & Schatz [18]: k∇2 gk2;Ω\Ωδ/2 ≤ ck∆gk2;Ω\Ω3δ/4 + cδ kgk2 , k∇(g − gh )k2;Ω\Ωδ/2 ≤ ck∇(g − ih g)k2;Ω\Ω3δ/4 + cδ kg − gh k2 , with constants cδ ≈ δ −1 . Now, the assertion follows by the interpolation estimate (5.3) and the other estimates already proven. Lemma 5.3. The following L1 a priori and error estimates hold k∇gk1 + k∇2 gk01 ≤ c| log(h)|,
(5.7)
gh )k01
(5.8)
2
k∇(g − gh )k1 + hk∇ (g − with a constant c independent of h and δ . Proof. The proof can be found in [9].
≤ ch| log(h)|,
15
FINITE ELEMENTS IN PARAMETER IDENTIFICATION
For the point a ∈ Ωd , there holds vh (a) = (vh , δh ) = (A(q)∇vh , ∇gh ) = (B∇ψ, ∇gh ). We employ a standard localization argument. Let ω ∈ C0∞ (Ω) be a a smooth function with the properties 0 ≤ ω ≤ 1,
ω|Ωδ/2 ≡ 1,
ω|Ω\Ωδ/4 ≡ 0.
With this notation, we have (B∇ψ, ∇gh ) = (B∇(ωψ), ∇gh ) + (B∇((1 − ω)ψ), ∇gh ) =: Σ1 + Σ2 . First, we estimate the term Σ1 . By integration by parts and observing that ψ|∂Ω = 0 , we obtain X (B∇(ωψ), ∇gh ) = (ωψ, −∇ · (B∇gh ))T + (ωψ, n · B∇gh )∂K\∂Ω , K∈Th
where n is the outward unit normal vector to ∂K . Let [∇gh ] denote the jump of the gradient across the interior faces Γ ⊂ ∂K . Using this notation, we obtain X Σ1 ≤ kψk∞;Ωδ/2 k∇ · (B∇gh )k1;K + 21 kn · [B∇gh ]k1;∂K\∂Ω . K∈Th
First, the estimates of Lemma 5.3 yield X k∇ · (B∇gh )k1;K ≤ c |||B|||1,∞ k∇gh k1 + k∇2 gh k01 ≤ c |||B|||1,∞ | log(h)|. K∈Th
Next, observing that g ∈ H 2 (Ω) and therefore [B∇g] = 0 , we obtain by a trace theorem: X X kn · [B∇gh ]k1;K\∂Ω = kn · [B∇(gh − g)]k1;∂K\∂Ω K∈Th
K∈Th
≤ c |||B|||1,∞
X
h
−1
k∇(g − gh )k1;K + k∇2 (g − gh )k01;K .
K∈Th
Hence, collecting the foregoing estimates, Σ1 ≤ c |||B|||1,∞ kψk∞;Ωδ/2 | log(h)| + h−1 k∇(g − gh )k1 + k∇2 (g − gh )k01 . Again the estimates of Lemma 5.3, we obtain Σ1 ≤ c |||B|||1,∞ kψk∞;Ωδ/2 | log(h)|.
(5.9)
For the term Σ2 , we estimate as follows: Σ2 = (B∇((1 − ω)ψ), ∇(gh − g)) + (B∇((1 − ω)ψ), ∇g) ≤ |||B|||1,∞ k∇ψk2 k∇(gh − g)k2;Ω\Ωδ/2 + cδ kψk2 k∇(gh − g)k2 + c |||B|||1,∞ kψk2 k∇2 gk2;Ω\Ωδ/2 + k∇gk2;Ω\Ωδ/2 . Then, by the L2 estimates of Lemmas 5.1 and 5.2, it follows that Σ2 ≤ cδ |||B|||1,∞ kψk2 + hk∇ψk2 . This completes the proof of the theorem.
(5.10)
16
R. RANNACHER AND B. VEXLER
6. Numerical results. In this section, we discuss a sample problem confirming the a priori error estimate of Theorem 4.1. The state equation is given by −∇ · (A(q)∇u) = 2 in Ω, u = 0 on ∂Ω,
(6.1)
where Ω is the unit square. The matrix A(q) is a function of the parameter q = (q1 , q2 ) ∈ Q = R2 , given by 2 q1 q1 q2 . A(q) = q1 q2 exp(q2 ) The parameters are estimated from the measurements of the state variable at nine different points ξi ∈ Ω , see Figure 6.1.
Fig. 6.1. The computational domain with measurement points marked by circles.
The vector of measurements Cˆ is given by Cˆi = Ci (S(ˆ q ))(1 + εi ) i = 1, . . . , 9, where the reference parameter is qˆ = (5, 6) and ε = (εi ) describes the data perturbation. We consider two cases: a) ε = 0,
b) ε ≈ (0.12, −0.26, 0.29, −0.37, −0.49, 0.13, −0.04, −0.45, 0.20).
For case (a) the cost functional J(u) vanishes in the optimal point. Case (b) is more realistic because of the “measurement errors” modelled by a randomly chosen ε . Moreover, in this case, in contrast to case (a) the solution q of the corresponding parameter identification problem and the reference parameter qˆ differ. The parameter identification problem is discretized using bilinear finite elements on uniformly refined meshes. The results are listed in Table 6.1 and Table 6.2. For both cases the theoretically predicted orders of convergence are achieved. 7. Conclusions and Extensions. In this paper we have derived an a priori error estimate for the finite element discretization of an elliptic discrete parameter identification problem with pointwise measurements. The crucial point in our argument is the stability estimate of Theorem 4.2. The result of Theorem 4.1 can be extended as far as such a stability estimate is available. We list some possible directions of generalization. 1. More general meshes: For simplicity, we have assumed a quasi uniform mesh family {Th }h . The analysis can be extended to locally refined meshes, provided that
FINITE ELEMENTS IN PARAMETER IDENTIFICATION
17
Table 6.1 The error and the order of convergence with respect to the components of q without data perturbation.
N 81 289 1089 4225 16641 66049 order
q1 − q1,h 5.955e-1 1.407e-1 3.436e-2 8.509e-3 2.098e-3 4.993e-4 2.05
q2 − q2,h 9.902e-4 1.731e-4 4.343e-5 1.080e-5 2.668e-6 6.352e-7 2.04
Table 6.2 The error and the order of convergence with respect to the components of q with data perturbation.
N 81 289 1089 4225 16641 66049 order
q1 − q1,h 2.059e-0 5.172e-1 1.467e-1 3.771e-2 9.832e-3 2.350e-3 2.01
q2 − q2,h 1.874e-2 2.999e-3 8.341e-4 2.111e-4 5.640e-5 1.348e-5 1.98
the ratio of hmin and h = hmax is polynomial, hmin ≈ hp , with some p ≥ 1. For such meshes the stability result of Theorem 4.2 holds true with | log(hmin )| ≈ p| log(h)|. 2. More general domains: Our argument uses that the solution operator S(·) maps Q into H01 (Ω) ∩ H 2 (Ω) which is guaranteed on smoothly bounded or convex domains. In the case of a domain with reentrant corners or edges this regularity property is lost. This lack of regularity of the solution can be compensated by an appropriate refinement of the mesh near the critical corner points or edges. The stability estimate of Theorem 4.2 also holds in this situation. This is true in two as well as in tree dimensions. This will be shown in a forthcoming paper. 3. Higher order approximation: The result of Theorem 4.1 can be also extended to the case of higher-order finite elements, similar to the analysis of Nitsche [17]. In this case the logarithmic factor | log(h)| can be dropped in the stability Theorem 4.2. 4. More general equations: Theorem 4.1 can also be extended to more general elliptic equations or systems of the form −∇ · (A(q)∇u + b1 (q)u) + b2 (q)u = f, with parameter depended coefficients A(q), b1 (q), b2 (q). REFERENCES ¨ ltzsch Error estimates for the numerical approximation of a [1] N. Arada, E. Casas and F. Tro semilinear elliptic control problem, J. Comput. Optim. Appl., 23 (2002), pp. 201–229.
18
R. RANNACHER AND B. VEXLER
[2] H. T. Banks and K. Kunisch Estimation Techniques for Distributed Parameter Systems, Boston: Birkh¨ auser, 1989. [3] R. Becker, M. Braack and B. Vexler Numerical parameter estimation for chemical models in multidimensional reactive flows, Combustion Theory and Modelling, to appear, 2004. [4] R. Becker and B. Vexler A posteriori error estimation for finite element discretization of parameter identification problems, Numer. Math. 96 (2004), pp. 435–459. [5] S. Brenner and R. L. Scott The mathematical Theory of Finite Element Methods, Springer, Berlin Heidelberg New York, 1994. [6] K. Deckelnick and M. Hinze Error estimates in space and time for tracking-type control of the instationary Stokes system, in Proc. 8th Conference on Control of Distributed Parameter Systems Graz, July 15-21, 2001, Intern. Series of Numerical Mathematics 143, 2002, pp. 87–104. [7] R. S. Falk Approximation of a class of optimal control problems with order of convergence estimates, J. Math. Anal. Appl., 44, (1973), pp. 28–47. [8] R. S. Falk Error estimates of the numerical identification of a variable coefficient, Math. Comp. 40 (1983), pp. 537–546. [9] J. Frehse and R. Rannacher Eine L1 -Fehlerabsch¨ atzung f¨ ur diskrete Grundl¨ osungen in der Methode der finiten Elemente, in Tagungsband Finite Elemente, Bonn. Math. Schr. 89 (1976), pp. 92–114. [10] P. Grisvard Singularities in Boundary Value Problems, Masson, Paris, and Springer, Berlin, 1992. [11] M. Gunzburger and L. Hou Finite-dimensional approximation of a class of constrained nonlinear optimal control problems, SIAM J. Control Optim., 34 (1996), pp. 1001–1043. [12] C. Johnson Numerical Solution of Partial Differential Equations by Finite Element Method, Cambridge University Press, Cambridge, 1987. ¨rkka ¨inen Error estimates for distributed parameter identification in linear elliptic equa[13] T. Ka tions, J. Math. Syst. Estim. Control, 6 (1996), pp. 117–120. [14] C. Kravaris and J. H. Seinfeld Identification of parameters in distributed parameter systems by regularization, SIAM J. Control Optim., 23 (1985), pp. 217–241. [15] V. G. Litvinov Optimization in elliptic problems with applications to mechanics of deformable bodies and fluid mechanics, Operator Theory: Advances and Applications. 119. Basel: Birkh¨ auser, 2000. ¨ki and X.-C. Tai Error estimates for numerical identification of distributed [16] P. Neittaanma parameters, J. Comput. Math., 10 (1992), pp. 66–78. ¨ [17] J. Nitsche Uber L∞ -Absch¨ atzungen von Projektoren auf Finite Elemente, Tagungsband Finite Elemente, Bonn. Math. Schr. 89 (1976), pp. 13–30 . [18] J. Nitsche and A. Schatz Interior estimates for Ritz-Galerkin methods, Math. Comput., 28 (1976), pp. 937–958. [19] J. Nocedal and S. J. Wright Numerical Optimization, Springer Series in Operations Research, Springer New York, 1999. [20] R. Rannacher L∞ -Stability estimates and asymptotic error expansion for parabolic finite element equations, in Tagungsband Extrapolation and Defect Correction, Bonn. Math. Schr.,228 (1991), pp. 74–94. [21] B. Vexler Adaptive Finite Element Methods for Parameter Identification Problems, PhD thesis, University of Heidelberg, 2004.
1. Introduction. We consider parameter identification problems governed by an elliptic partial differential equation of second order. The finitely many unknown parameters are estimated using the measurements of point values of the state variable. Let Ω ⊂ R2 be a convex polygonal domain, L2 (Ω) the corresponding Lebesgue space with inner product and norm denoted by (·, ·) and k · k2 , respectively, and H m (Ω) the Sobolev space of order m ∈ N. With this notation, we set V := {v ∈ H 1 (Ω) ∩ C(Ω) | v = 0 on ∂Ω}. The state variable u ∈ V is determined by an elliptic partial differential equation (the state equation) −∇ · (A(q)∇u) = f u=0
in Ω on ∂Ω,
(1.1)
for given H¨older continuous f ∈ C α (Ω), α ∈ (0, 1) . Here, q ∈ Q0 ⊂ Q = Rnp denotes the set of unknown parameters and A(q) = (Aij (q)) is a symmetric and positive definite 2×2 matrix with twice continuously differentiable entries Aij : Q → C 1+α (Ω). The above conditions guarantee that for any admissible value of the parameter q , the corresponding solution u of the state equation (1.1) is in H 2 (Ω) . At the corner points of ∂Ω , the second derivatives of the solution may become singular. However, u has H¨older continuous second derivatives, u ∈ C 2+α (Ωd ), for each subdomain Ωd ⊂ Ω with positive distance to the corner points. The usual weak formulation of (1.1) is: a(q)(u, φ) = (f, φ) ∀φ ∈ V,
(1.2)
where the bilinear form a(q)(·, ·) is defined by a(q)(u, φ) := (A(q)∇u, ∇φ). † Institute
(1.3)
of Applied Mathematics, University of Heidelberg, INF 294, 69120 Heidelberg, Germany. ‡ This work has been supported by the German Research Foundation (DFG), through the Sonderforschungsbereich 359 ’Reactive Flows, Diffusion and Transport’, and the Graduiertenkolleg ’Complex Processes: Modeling, Simulation and Optimization’ at the Interdisciplinary Center of Scientific Computing (IWR), University of Heidelberg. 1
2
R. RANNACHER AND B. VEXLER
Further, the observation operator C(·) describing the mapping of the state variable u to the space of measurements Z = Rnm is given by Ci (v) = v(ξi ),
i = 1, 2, . . . , nm ,
(1.4)
where {ξi } ⊂ Ω is a finite set of measurement points. We assume that nm ≥ np . The euclidian product and norm on Q and Z are denoted by h·, ·i and k · k , respectively, and the same notation is also used for the corresponding natural norms of matrices. The values of the parameters are estimated from a given set of measurements Cˆ ∈ Z using a least squares approach such that we obtain a constrained optimization problem with the cost functional J : V → R: Minimize
ˆ 2, J(u) := 21 kC(u) − Ck
(1.5)
under the constraint (1.2). Throughout we assume the existence of a solution (u, q) ∈ V ×Q0 of the problem (1.2,1.5). For an analysis of existence of solutions for parameter identification problems see, e.g., Banks & Kunisch [2], Kravaris & Seinfeld [14] and Litvinov [15]. The state equation is discretized by a conforming finite element Galerkin method defined on a family {Th }h>0 of shape regular quasi-uniform meshes Th = {K} consisting of closed cells K which are either triangles or quadrilaterals. The straight parts which make up the boundary ∂K of a cell K are called faces. The mesh parameter h is defined as a cell-wise constant function by setting h|K = hK and hK is the diameter of K. Usually we use the symbol h also for the maximal cell size, i.e. h = max hK . K∈Th
(1.6)
For convenience, we assume that 0 < h < 1 . On the mesh Th we define finite element spaces Vh ⊂ V consisting of linear or bilinear shape functions, see e.g. Brenner & Scott [5] or Johnson [12]. The corresponding discrete state uh ∈ Vh and parameter qh ∈ Q0 are determined by: Minimize
J(uh ),
(1.7)
under the constraint a(qh )(uh , φh ) = (f, φh ) ∀φh ∈ Vh .
(1.8)
Since Q is finite dimensional, the parameter qh is determined in the same space Q. The main purpose of this paper is to analyze the behavior of the error in parameters kq − qh k for h tending to zero. There is a number of publications, where a priori error estimates are derived for optimal control problems governed by partial differential equations, see Falk [7], Arada, Casas & Tr¨oltzsch [1], Deckelnic & Hinze [6] and Gunzburger & Hou [11]. However, there are only few published results on this topic in the context of parameter identification problems, see Falk [8], Neittaanm¨aki & Tai [16] and K¨ arkk¨ ainen [13]. In Vexler [21], an a priori error analysis for the case of V -stable observation operators Ci (·) is developed and optimal-order convergence is shown, kq − qh k = O(h2 ),
(1.9)
FINITE ELEMENTS IN PARAMETER IDENTIFICATION
3
essentially under the assumption that |Ci (u) − Ci (uh )| ≤ ch2 kukH 2 . However, the generalization of this result to point-wise observations is not straightforward. For this case, we prove in this paper that, under certain regularity conditions, kq − qh k = O(h2 | log(h)|2 ).
(1.10)
The proof uses the technique for estimating discrete Green functions developed in Frehse & Rannacher [9]. A complementary result of a posteriori error analysis for parameter identification problems is given in Becker & Vexler [4]. To the authors knowledge, this is the first a priori error analysis for parameter identification problems with point-wise observations. The consideration of point-wise observations in determining discrete parameters seems very natural in view of practical measurement techniques; see [3] for applications in reactive flow analysis. The paper is organized as follows: In the next section, we describe an algorithm for solving problem (1.7, 1.8). In Section 3, we present a paradigm for a priori error analysis for discretization of a class of optimization problems. Thereafter, in Section 4, we derive the announced error estimate using an L∞ -stability result, which is proven in Section 5. In Section 6, we present a numerical example conforming the asymptotic sharpness of our error estimate. Possible extensions are addressed in the last section. 2. Optimization algorithm. In this section, we reformulate the problem under consideration as an unconstrained optimization problem and describe a solution algorithm for it. Throughout, we assume the coefficient matrix A(q) to be positive definite at the solution of the problem. By virtue of the implicit function theorem for Banach spaces, this implies the existence of an open (bounded) set Q0 ⊂ Q containing the optimal parameter q , and of a continuously differentiable solution operator S : Q0 → V satisfying: a(q)(S(q), φ) = (f, φ) ∀φ ∈ V
∀q ∈ Q0 .
(2.1)
Without loss of generality, we assume that A(q) is positive definite for q ∈ Q0 , i.e. there exists γ ∈ R+ , such that p∗ A(q)p ≥ γkpk2 ,
∀p ∈ Q ∀q ∈ Q0 .
(2.2)
The convexity of the domain Ω implies the following regularity property (see e.g. Grisvard [10]): S(q) ∈ V ∩ H 2 (Ω) ∀q ∈ Q0 .
(2.3)
We introduce the reduced observation operator c : Q0 → Z by c(q) := C(S(q)).
(2.4)
This allows us to reformulate the problem under consideration as an unconstrained optimization problem with the reduced cost functional j : Q0 → R : Minimize
ˆ 2, j(q) := 21 kc(q) − Ck
q ∈ Q0 .
(2.5)
Denoting by G = c0 (q) ∈ Rnp ×nm the Jacobian matrix of the reduced observation operator c(·), the first-order necessary optimality condition j 0 (q) = 0 for (2.5) reads ˆ = 0, G∗ (c(q) − C)
(2.6)
4
R. RANNACHER AND B. VEXLER
where G∗ denotes the transpose of G. The positive semi-definiteness of the Hessian matrix H := ∇2 j(q) is the second-order necessary optimality condition. A solution q of problem (2.5) is called stable, if the sufficient optimality condition holds, i.e. if the Hessian H is positive definite. Throughout, we will assume the solution q to be stable. The stability of the solution is given, for instance, if the value of the ˆ is small enough and the matrix G has full rank np , see cost functional kC(u) − Ck e.g. [21] for details. Since by assumption the matrix coefficient A(·) is twice continuously differentiable, there holds sup |||A(ξ)|||1,∞ + sup |||A0qj (ξ)|||1,∞ + sup |||A00qj qk (ξ)|||1,∞ < ∞,
ξ∈Q0
ξ∈Q0
(2.7)
ξ∈Q0
where |||B|||1,∞ := maxi,j=1,2 kBij k1,∞ for a matrix function B = (Bij ) ∈ C 1 (Ω)2×2 . In the following propositions, we give representations of the Jacobian G of c(·) , the Hessian H of j(·) , and the Hessian of ci (·) . Theorem 2.1. Let the reduced observation operator c(·) and the reduced functional j(·) be defined as in (2.4) and (2.5), respectively. Then, the elements of the Jacobian of c(·) at some q ∈ Q0 are given by Gij =
∂ci (q) = Ci (wj ), ∂qj
i = 1, . . . , nm , j = 1, . . . , np ,
(2.8)
where wj ∈ V are the solutions of the problems a(q)(wj , φ) = −(A0qj (q)∇u, ∇φ)
∀φ ∈ V,
(2.9)
with u = S(q) . The Hessian of j(·) can be expressed by H = G∗ G + M,
(2.10)
where the matrix M ∈ Rnp ×np is given by M=
nm X
c00i (q)(ci (q) − Cˆi ).
(2.11)
i=1
The Hessian of ci (q) is given by ∂2 ci (q) = Ci (vjk ), ∂qj ∂qk
(2.12)
where the vjk ∈ V are the solutions of the problems a(q)(vjk , φ) = −(A0qj (q)∇wk , ∇φ) − (A0qk (q)∇wj , ∇φ) − (A00qj qk (q)∇u, ∇φ)
∀φ ∈ V,
(2.13)
with wj as defined in (2.9). Proof. The derivation of the derivatives of c(·) is straightforward using the implicit function theorem and the chain rule. In practice the Hessian H of j(·) is computed using the representation Mjk = −(A0qj (q)∇wk , ∇z) − (A0qk (q)∇wj , ∇z) − (A00qj qk (q)∇u, ∇z),
(2.14)
FINITE ELEMENTS IN PARAMETER IDENTIFICATION
5
with the function z determined by the dual equation ¯ C(φ)i. (A(q)∇φ, ∇z) = hC(u) − C,
(2.15)
Similar to the continuous case, we introduce a discrete solution operator Sh : Q0 → Vh by the equation a(qh )(Sh (qh ), φh ) = (f, φh ) ∀φh ∈ Vh , qh ∈ Q0 .
(2.16)
As before, we turn the discrete problem (1.7, 1.8) into an unconstrained minimization problem: Minimize
ˆ 2, jh (qh ) := 21 kch (qh ) − Ck
qh ∈ Q0 ,
(2.17)
where the discrete reduced observation operator ch is defined by ch (qh ) = C(Sh (q)).
(2.18)
Denoting the corresponding Jacobian by Gh = c0h (qh ), the necessary optimality condition jh0 (qh ) = 0 reads ˆ = 0. G∗h (ch (qh ) − C)
(2.19)
The derivatives of the discrete observation operator ch can be computed in an analogous way as in Theorem 2.1. Problem (2.17) is solved iteratively starting with an initial guess qh0 and using the recursive setting qhk+1 = qhk + δqh . The update δqh is obtained as the solution of the system of linear equations Hk δqh = G∗h (Cˆ − ch (qhk )),
(2.20)
where Gh = c0h (qhk ) , and Hk is an appropriate symmetric approximation of the Hessian ∇2 jh (qhk ) . The most widely used choice of the matrix Hk = G∗h Gh leads us to the Gauß-Newton algorithm, see, e.g., Nocedal & Wright [19]. For one step of the Gauß-Newton algorithm the state equation and np tangent problems (2.9) have to be solved which involves the same linear operator but with different right-hand sides. Due to the small dimension np of the parameter space Q the solution of (2.20) is uncritical. For discussing other Newton type methods and trust-region techniques for globalization of the convergence in this context, see [21]. 3. A paradigm for a priori error analysis. In this section we present a general approach to the error analysis of a class of optimization problems such as considered in this paper. The main result is stated in the following theorem. It is a variant of well-known perturbation theorems for differentiable mappings which is particularly tailored to the present situation. However, it seems easier to include the elementary proof than to search for the precise reference. Theorem 3.1. Let F, Fh : Rn → Rn , for a discretization parameter h ∈ R+ , be continuously differentiable operators , and x ∈ Rn be a solution of F (x) = 0 . Let the following conditions be fulfilled: (i) The derivative F 0 (x) is positive definite, i.e., there is a constant γ > 0 , such that p∗ F 0 (x)p ≥ γkpk2 ,
p ∈ Rn .
(3.1)
6
R. RANNACHER AND B. VEXLER
(ii) There is a neighborhood U of x and a positive number L(h) ∈ R+ , such that kFh0 (ξ) − Fh0 (η)k ≤ L(h)kξ − ηk
∀ξ, η ∈ U.
(3.2)
(iii) With the h-dependent constant L(h) , there holds lim L(h)kF (x) − Fh (x)k = 0.
(3.3)
lim kF 0 (x) − Fh0 (x)k = 0.
(3.4)
h→0
(iv) There holds h→0
Then, for h small enough, there exists xh ∈ U , such that Fh (xh ) = 0 , and Fh0 (xh ) is positive definite uniformly in h . Further, there holds the a priori error estimate kx − xh k ≤
2 kF (x) − Fh (x)k. γ
(3.5)
Proof. Due to condition (iv), we can choose a positive number h1 ∈ R+ , such that for h ≤ h1 there holds kF 0 (x) − Fh0 (x)k ≤ 41 γ . Moreover, for ρ = ρ(h) =
γ kL(h)
(3.6)
, with some k ≥ 4 sufficiently large, there holds
Bρ (x) = {ξ ∈ Rn , kx − ξk ≤ ρ} ⊂ U.
(3.7)
For this choice, we obtain that, for h ≤ h1 , Fh0 (·) is positive definite on Bρ (x) : p∗ Fh0 (ξ)p = p∗ F 0 (x)p + p∗ (Fh0 (x) − F 0 (x))p + p∗ (Fh0 (ξ) − Fh0 (x))p ≥ γkpk2 − kFh0 (x) − F 0 (x)k kpk2 − kFh0 (ξ) − Fh0 (x)k kpk2 ≥ (γ − 41 γ − L(h)ρ)kpk2 ≥ 12 γkpk2 . In a similar way, we conclude that, for h ≤ h1 , Fh0 (·) is also bounded on Bρ (x): kFh0 (ξ)k ≤ β := kF 0 (x)k + 21 γ . Next, we prove, that there exists a unique xh ∈ Bρ (x) with Fh (xh ) = 0 . To this end, we define an operator Ds : Rn → Rn , for s ∈ R+ , by Ds (ξ) = ξ − sFh (ξ). For a certain choice of s , we show that Ds is a contraction on Bρ (x) and use the Banach fixed point theorem. For ξ ∈ Bρ (x), h ≤ h1 , and an arbitrary p ∈ Rn , there holds kDs0 (ξ)pk2 = kp − sFh0 (ξ)pk2 = kpk2 − 2sp∗ Fh0 (ξ)p + s2 kFh0 (ξ)pk2 ≤ (1 − sγ + s2 β 2 )kpk2 . For the choice s = γ(2β 2 )−1 , we obtain γ2 kDs0 (ξ)pk2 ≤ 1 − 2 kpk2 , 4β
FINITE ELEMENTS IN PARAMETER IDENTIFICATION
7
and consequently, γ 2 1/2 kDs0 (ξ)k ≤ 1 − 2 < 1. 4β Moreover, for arbitrary ξ ∈ Bρ (x), there holds: kx − Ds (ξ)k = kDs (x) − Ds (ξ) + sFh (x)k ≤ kDs (x) − Ds (ξ)k + skFh (x) − F (x)k ≤ kDs0 (η)k kx − ξk + skFh (x) − F (x)k, for a certain η ∈ Bρ . Hence, the above estimate implies γ 2 1/2 kx − Ds (ξ)k ≤ 1 − 2 ρ + s kFh (x) − F (x)k 4β γ 2 1/2 k =ρ 1− 2 + s L(h)kFh (x) − F (x)k . 4β γ Due to condition (iii), there is a number h2 ∈ R+ , such that, for h ≤ h2 , there holds γ 2 1/2 γ 1− 1− 2 L(h)kFh (x) − F (x)k ≤ . ks 4β Hence, for h < h0 := min{h1 , h2 }, kx − Ds (ξ)k ≤ ρ, and consequently Ds (ξ) ∈ Bρ (x) . For h ≤ h0 , by the Banach fixed point theorem, we obtain the existence of xh ∈ Bρ (x) with Fh (xh ) = 0. By construction of Bρ (x) the derivative Fh0 (xh ) is positive definite with the h-independent constant 21 γ . This implies that, for a certain ξ ∈ Bρ (x), (x − xh )∗ (Fh (x) − Fh (xh )) = (x − xh )∗ Fh0 (ξ)(x − xh ) ≥
γ kx − xh k2 . 2
Hence, using F (x) = Fh (xh ) = 0, kx − xh k2 ≤
2 2 (x − xh )∗ (Fh (x) − Fh (xh )) = (x − xh )∗ (Fh (x) − F (x)) γ γ 2 ≤ kFh (x) − F (x)k kx − xh k. γ
This completes the proof. 4. A priori error estimation. In this section we apply the paradigm presented in Section 3 to the problem under consideration. We prove the following theorem. Theorem 4.1. Let q ∈ Q be a stable solution of (2.5). Then, for h small enough, there exist a stable solution qh ∈ Q of (2.17), and there holds the following a priori error estimate: kq − qh k = O(h2 | log(h)|2 ).
(4.1)
8
R. RANNACHER AND B. VEXLER
On the basis of the estimate (4.1), we can also derive optimal-order estimates for the error u − uh in the corresponding states. However, since this would be a simple exercise using the arguments developed below, and since the optimal states are of only minor practical interest in parameter estimation problems, we do not state these estimates. The proof of Theorem 4.1 is given by checking the conditions from Theorem 3.1 for the operators F (ξ) := ∇j(ξ)
Fh (ξ) := ∇jh (ξ).
The constant in (4.1) turns out to depend in a reciprocal way on the distance δ :=
min
i=1,...,nm
dist(ξi , Σ)
of the set of measurement points ξi to the set Σ of corner points of ∂Ω. Therefore, we will use generic constants c and cδ , where c only depends on the domain Ω, the force f , and the characteristics of the mesh family {Th }h , while cδ may additionally depend on the distance δ like cδ ≈ δ −1 . Further, by Lp (Ω) and W m,p (Ω) , for m ∈ N and 1 ≤ p ≤ ∞ , we denote the standard Lebesgue and Sobolev spaces and by k · kp and k · km,p the corresponding norms. The restriction of such a norm to a subset Ω0 ⊂ Ω is indicated by k · km,p;Ω0 . An important ingredient of the proof of Theorem 4.1 is the following L∞ -stability theory. For d > 0 , we define the subset Ωd ⊂ Ω by Ωd := {x ∈ Ω, dist(x, Σ) > d}. ¯ , Theorem 4.2 (Stability theorem). Let be given q ∈ Q0 , ψ ∈ H01 (Ω) ∩ C(Ω) 1,∞ 2×2 and a matrix B = B(x) ∈ W (Ω) . Moreover, let vh ∈ Vh be a solution of a(q)(vh , φh ) = (B∇ψ, ∇φh )
∀φh ∈ Vh .
(4.2)
Then, there hold the L2 -stability estimate kvh k2 + hk∇vh k2 ≤ c |||B|||1,∞ kψk2 + hk∇ψk2 , and the local L∞ -stability estimate kvh k∞;Ωd ≤ cd |||B|||1,∞ | log(h)| kψk∞;Ωd/2 + kψk2 + hk∇ψk2 ,
(4.3)
(4.4)
with a constant cd ≈ d−1 . The L2 estimate (4.3) is a standard result from finite element analysis, while the L∞ estimate (4.4) follows by estimates of discrete Green functions such as developed in Frehse & Rannacher [9]. A similar L∞ -stability result has been proven in Rannacher [20] in the time-dependent parabolic case. The proof is given in Section 5. For the solution q of problem (2.5), we introduce u ¯h ∈ Vh determined by a(q)(¯ uh , φh ) = (f, φh ) ∀φh ∈ Vh .
(4.5)
Further, we define wj,h ∈ Vh and vjk,h ∈ Vh , for j, k = 1, 2, . . . , np , as the solutions of the problems a(q)(wj,h , φh ) = −(A0qj (q)∇¯ uh , ∇φh ) ∀φh ∈ Vh ,
(4.6)
9
FINITE ELEMENTS IN PARAMETER IDENTIFICATION
and a(q)(vjk,h , φh ) = −(A0qj (q)∇wk,h , ∇φh ) − (A0qk (q)∇wj,h , ∇φh )− (A00qj qk (q)∇¯ uh , ∇φh ) ∀φh ∈ Vh ,
(4.7)
respectively. The next lemma provides necessary estimates for the errors u − u ¯h , wj − wj,h and vjk − vjk,h . We recall the notation δ := mini=1,...,nm dist(ξi , Σ) . Lemma 4.3. Under the above assumptions the following estimates hold: kC(u − u ¯h )k ≤ cδ h2 | log(h)|, 2
2
kC(wj − wj,h )k ≤ cδ h | log(h)| , kC(vjk − vjk,h )k ≤ cδ h2 | log(h)|3 ,
(4.8) j = 1, 2 . . . np , j, k = 1, 2 . . . np .
(4.9) (4.10)
Proof. (i) The proof uses some a priori bounds for u and the auxiliary functions wj , vjk , j, k = 1, . . . , np , corresponding to arbitrary q ∈ Q0 . By our regularity assumptions, we have the following a priori estimates: kuk2,2 ≤ c, kwj k2,2 ≤ ckuk2,2 ≤ c, kvjk k2,2 ≤ c max kwj k2,2 , kwk k2,2 + ckuk2,2 ≤ c,
(4.11) (4.12) (4.13)
where c is a generic constant depending only on the data of the problem. Further, by elliptic regularity theory, for d > 0 , there hold the local a priori estimates (4.14) kukC 2+α (Ωd ) ≤ cd kf kC α (Ωd/2 ) + kuk2,2 ≤ cd , kwj kC 2+α (Ωd ) ≤ cd kukC 2+α (Ωd/2 ) + kuk2,2 ≤ cd , (4.15) (4.16) kvjk kC 2+α (Ωd ) ≤ cd kwj kC 2+α (Ωd/2 ) + kwj k2,2 ≤ cd , with a generic constant cd ≈ d−1 . Here, k · kC 2+α (Ωd ) denotes the norm of the space C 2+α (Ωd ) of functions possessing H¨ older continuous second derivatives on Ωd . (ii) By definition, u ¯h is the Ritz projection of u corresponding to the energy form a(q)(·, ·) , i.e., a(q)(¯ uh , φh ) = a(q)(u, φh ) ∀φh ∈ Vh . By the standard L2 -error estimate for finite elements, there holds ku − u ¯h k2 + hk∇(u − u ¯h )k2 ≤ ch2 .
(4.17)
Further, let ih u ∈ Vh be the usual nodal interpolant of u . Then, applying the L∞ -stability estimate (4.4) of Theorem 4.2 for the equation a(q)(ih u − u ¯h , φh ) = a(q)(ih u − u, φh ) ∀φh ∈ Vh , yields the estimate kih u − u ¯h k∞;Ωδ ≤ cδ | log(h)| kih u − uk∞;Ωδ/2 + kih u − uk2 + hk∇(ih u − u)k2 .
10
R. RANNACHER AND B. VEXLER
In view of the approximation properties of ih and the a priori bounds (4.11) and (4.14), we conclude the error estimate ku − u ¯h k∞;Ωδ ≤ cδ h2 | log(h)|.
(4.18)
Since ξi ∈ Ωδ , we obtain the estimate (4.8). (iii) For proving (4.9), we introduce an additional discrete variable w ¯j,h determined by the equation a(q)(w ¯j,h , φh ) = −(A0qj (q)∇u, ∇φh ) ∀φh ∈ Vh . The error e = wj − wj,h is split like e = e1 + e2 , with e1 = wj − w ¯j,h and e2 = w ¯j,h − wj,h . For the Ritz-projection error e1 , similar as before, there holds the L2 error estimate ke1 k2 + hk∇e1 k2 ≤ ch2 kwj k2,2 ≤ ch2 , and the pointwise error estimate ke1 k∞;Ωδ ≤ cδ h2 | log(h)| k∇2 wj k∞;Ωδ/2 + kwj k2,2 ≤ cδ h2 | log(h)|. For e2 ∈ Vh , we have a(q)(e2 , φh ) = −(A0qj (q)∇(u − u ¯h ), ∇φh ) ∀φh ∈ Vh . Hence, the L2 -stability estimate (4.3) of Theorem 4.2 and the estimate (4.17) imply ke2 k2 + hk∇e2 k2 ≤ c ku − u ¯h k2 + hk∇(u − u ¯h )k2 ≤ ch2 . This shows that, for j = 1, . . . , np , kwj − wj,h k2 + hk∇(wj − wj,h )k2 ≤ ch2 .
(4.19)
Further, the L∞ -stability estimate (4.4) of Theorem 4.2 yields ke2 k∞;Ωδ ≤ cδ | log(h)| ku − u ¯h k∞;Ωδ/2 + ku − u ¯h k2 + k∇(u − u ¯h )k2 . which, by (4.17) and (4.18), implies ke2 k∞;Ωδ ≤ cδ | log(h)|2 h2 . We obtain kwj − wj,h k∞;Ωδ ≤ cδ | log(h)|2 h2 ,
j = 1, . . . , np ,
(4.20)
which implies the desired estimate (4.9). (iv) The proof of (4.10) uses the same line of argument as before. Using the additional discrete variable v¯jk,h determined by the equation a(q)(¯ vjk,h , φh ) = −(A0qj (q)∇wk , ∇φh ) − (A0qk (q)∇wj , ∇φh ) − (A00qj qk (q)∇u, ∇φh ) ∀φh ∈ Vh , the error e = vjk − vjk,h is split like e = e1 + e2 , with e1 = vjk − v¯jk,h and e2 = v¯jk,h − vjk,h . For the Ritz-projection error e1 , similar as before, we conclude the the pointwise error estimate ke1 k∞;Ωδ ≤ cδ h2 | log(h)| k∇2 vjk k∞;Ωδ/2 + kvjk k2,2 ≤ cδ h2 | log(h)|.
11
FINITE ELEMENTS IN PARAMETER IDENTIFICATION
For e2 ∈ Vh , we have a(q)(e2 , φh ) = −(A0qj (q)∇(wk − wk,h ), ∇φh ) − (A0qk (q)∇(wj − wj,h ), ∇φh ) − (A00qj qk (q)∇(u − u ¯h ), ∇φh ) ∀φh ∈ Vh , and therefore, again by the L∞ -stability estimate (4.4) of Theorem 4.2, ke2 k∞;Ωδ ≤ cδ | log(h)| max kwj − wj,h k∞;Ωδ/2 + ku − u ¯h k∞;Ωδ/2 j=1,...,np kwj − wj,h k2 + hk∇(wj − wj,h )k2 + c max j=1,...,np + c ku − u ¯h k2 + hk∇(u − u ¯h )k2 . Then, by the foregoing error estimates, we obtain ke2 k∞;Ωδ ≤ cδ h2 | log(h)|3 , and consequently, kvjk − vjk,h k∞;Ωδ ≤ cδ | log(h)|3 h2 ,
j, k = 1, . . . , np .
(4.21)
This eventually yields the desired estimated (4.10). A direct application of Lemma 4.3 leads to the following result. Lemma 4.4. Under the above assumptions, there holds ∂ (j − jh )(q) ≤ cδ h2 | log(h)|2 , j = 1, 2, . . . , np , ∂qj ∂2 (j − jh )(q) ≤ cδ h2 | log(h)|3 , j, k = 1, 2, . . . , np . ∂qj ∂qk
(4.22) (4.23)
Proof. We have the representation ∂ ˆ C(wj )i − hC(¯ ˆ C(wj,h )i (j − jh )(q) = hC(u) − C, uh ) − C, ∂qj ˆ C(wj − wj,h )i + hC(u − u = hC(u) − C, ¯h ), C(wj,h )i, from which we obtain ∂ ˆ kC(wj − wj,h )k + kC(u − u (j − jh )(q) ≤ kC(u) − Ck ¯h )k kC(wj,h )k. ∂qj By the a priori bounds (4.11)-(4.13) and the Sobolev imbedding theorem, we see that kC(u)k + kC(wj )k + kC(vjk )k ≤ c.
(4.24)
Combining this with the error estimate (4.9) implies kC(wj,h )k ≤ c . Then, we can conclude the first estimate (4.22) from the error estimates of Lemma 4.3. To prove (4.23), we write ∂2 ˆ C(vjk )i (j − jh )(q) = hC(wj ), C(wk )i + hC(u) − C, ∂qj qk ˆ C(vjk,h )i − hC(wj,h ), C(wk,h )i − hC(¯ uh ) − C, = hC(wj − wj,h ), C(wk )i + hC(wj,h ), C(wk − wk,h )i ˆ C(vjk − vjk,h )i. + hC(u − u ¯h ), C(vjk )i + hC(¯ uh ) − C,
12
R. RANNACHER AND B. VEXLER
Using as before the bounds (4.24) and the error estimates of Lemma 4.3 completes the proof. For the application of Theorem 3.1 it remains to check the Lipschitz condition (3.2). For two arbitrary parameter sets ξ, η ∈ Q0 , we set uξ = Sh (ξ) and uη = Sh (η) . Correspondingly, we define wj,ξ , wj,η ∈ Vh and vjk,ξ , vjk,η ∈ Vh similarly to wj,h and vj,h for q = ξ and q = η , respectively. Lemma 4.5. For ξ, η ∈ Q0 , there holds kC(uξ − uη )k ≤ cδ | log(h)| kξ − ηk,
(4.25)
2
(4.26)
3
(4.27)
kC(wj,ξ − wj,η )k ≤ cδ | log(h)| kξ − ηk, kC(vjk,ξ − vjk,η )k ≤ cδ | log(h)| kξ − ηk. Proof. Due to the definition of uξ and uη , we have (A(ξ)∇(uξ − uη ), ∇φh ) = −((A(ξ) − A(η))∇uη , ∇φh ) ∀φh ∈ Vh .
Using Theorem 4.2, with d = δ , we obtain: kC(uξ − uη )k ≤ c |||(A(ξ) − A(η)|||1,∞ | log(h)| kuη k∞;Ωδ/2 + kuη k2 + hk∇uη k2 . Since uη is the Ritz projection of an H 2 function, all its norms occurring on the right-hand side are bounded independent of h and η ∈ Q0 by standard estimates from finite element analysis. This implies (4.25) since |||A(ξ) − A(η)|||1,∞ ≤ ckξ − ηk. The estimates (4.26) and (4.27) are obtained in the similar way. Lemma 4.6. For ξ, η ∈ Q0 , there holds: ∂2 ∂2 jh (ξ) − jh (η) ≤ L(h)kξ − ηk, ∂qj qk ∂qj qk where L(h) = cδ | log(h)|3 . Proof. We have ∂2 ∂2 ˆ C(vjk,ξ )i jh (ξ) − jh (η) = hC(wj,ξ ), C(wk,ξ )i + hC(uξ ) − C, ∂qj qk ∂qj qk ˆ C(vjk,η )i − hC(wj,η ), C(wk,η )i − hC(uη ) − C, = hC(wj,ξ − wj,η ), C(wk,ξ )i + hC(wj,η ), C(wk,ξ − wk,η )i ˆ C(vjk,ξ − vjk,η )i, + hC(uξ − uη ), C(vjk,ξ )i + hC(uη ) − C, and, consequently, ∂2 ∂2 jh (ξ) − jh (η) ≤ kC(wj,ξ − wj,η )k kC(wk,ξ )k ∂qj qk ∂qj qk + kC(wj,η )k kC(wk,ξ − wk,η )k + kC(uξ − uη )k kC(vjk,ξ )k ˆ kC(vjk,ξ − vjk,η )k. + kC(uη ) − Ck
(4.28)
13
FINITE ELEMENTS IN PARAMETER IDENTIFICATION
Now, the assertion follows by the estimates of Lemma 4.5 if we can bound the terms C(uη ), C(wj,η ), and C(vjk,ξ ) . This is achieved by using the bounds for C(u), C(wj ), and C(vjk ) in (4.24) together with the error estimates of Lemma 4.3. To complete the proof of Theorem 4.1 we check the conditions of Theorem 3.1. Condition (3.1) is fulfilled due to the stability of the solution q of the problem (2.5). Condition (3.2) is shown in Lemma 4.6. Condition (3.3) is obtained by Lemma 4.4 and Lemma 4.6 using limh→0 h2 | log(h)|5 = 0 . Finally, condition (3.4) holds due to Lemma 4.4. Hence, the estimate (3.5) of Theorem 3.1 completes the proof. 5. Proo f of Theorem 4.2. (i) We begin with the L2 -stability estimate. Taking φh := vh in (4.2), we obtain k∇vh k2 ≤ c|||B|||1,∞ k∇ψk2 .
(5.1)
To estimate kvh k2 , we use the solution z ∈ V ∩ H 2 (Ω) of the auxiliary equation a(q)(φ, z) = (φ, vh )kvh k−1 2
∀φ ∈ V.
Then, taking φ := vh and using the approximation properties of the interpolant ih z ∈ Vh , we conclude that kvh k2 = a(q)(vh , z) = a(q)(vh , z − ih z) + a(q)(vh , ih z) = a(q)(vh , z − ih z) + (B∇ψ, ∇ih z) = a(q)(vh , z − ih z) + (B∇ψ, ∇(ih z − z)) + (B∇ψ, ∇z) ≤ ck∇vh k2 k∇(z − ih z)k2 + |||B|||1,∞ k∇ψk2 k∇(z − ih z)k2 + |||B|||1,∞ kψk2 kzk2,2 ≤ c hk∇vh k + |||B|||1,∞ hk∇ψk2 + |||B|||1,∞ kψk2 kzk2,2 . Hence, observing (5.1) and the bound kzk2,2 ≤ c , we obtain kvh k2 ≤ c|||B|||1,∞ kψk2 + hk∇ψk2 .
(5.2)
(ii) Next, we prove the L∞ -stability estimate. Let a ∈ Ωδ be an arbitrary point lying in a cell K . For any fixed h , there exists a cell-wise polynomial function δh with supp(δh ) ⊂ K, such that (φh , δh ) = φh (a) ∀φh ∈ Vh . The function δh plays the role of an approximate Dirac function. Correspondingly, we introduce a regularized Green function g ∈ V ∩ H 2 (Ω) by (A(q)∇φ, ∇g) = (δh , φ) ∀φ ∈ V, and the corresponding Ritz projection gh ∈ Vh by (A(q)∇φh , ∇gh ) = (δh , φh ) ∀φh ∈ Vh . By ih : C(Ω) → Vh , we denote the usual (linear) operator of nodal interpolation for which the following cell-wise estimate is well known: −1 2 0 2 h−2 T kv − ih vkp;T + hT k∇(v − ih v)kp:T + k∇ ih vkp;T ≤ ck∇ vkp;T ,
(5.3)
14
R. RANNACHER AND B. VEXLER
for 1 ≤ p ≤ ∞ , with constants c independent P of h . For functions which are only cell-wise defined, we will use the norm kvk0p := K∈Th kvkp;T . The following three lemmas provide the key estimates for the proof of the theorem. Lemma 5.1. The following global L2 estimates hold: kgk2 + | log(h)|−1/2 k∇gk2 + hk∇2 gk2 ≤ c, h
−1
kg − gh k2 + k∇(g − gh )k2 + hk∇
2
gh k02
≤ c.
(5.4) (5.5)
Proof. The assertion follows by standard L2 a priori and error estimates for g and g − gh , respectively. We skip the details and refer to [9]. Note that k∇2 gh k02 vanishes for linear finite elements. In the case of bilinear elements, we estimate as follows: k∇2 gh k02 ≤ k∇2 (gh − ih g)k02 + k∇2 (g − ih g)k02 + k∇2 gk2 . For the first term, we obtain, using an inverse inequality, k∇2 (gh − ih g)k02 ≤ ch−1 k∇(gh − ih g)k2 ≤ ch−1 k∇(g − ih g)k2 + k∇(g − gh )k2 and obtain by the interpolation estimate (5.3), with p = 2 , and by the other estimates derived before: k∇2 gh k02 ≤ c h−1 k∇(g − gh )k2 + k∇2 gk2 ≤ ch−1 . This completes the proof. Lemma 5.2. For sufficiently small h δ , the following local L2 estimate holds: k∇(g − gh )k2;Ω\Ωδ/2 + hk∇2 gk2;Ω\Ωδ/2 ≤ cδ h,
(5.6)
with a constant cδ ≈ δ −1 but independent of h . Proof. The assertion follows by standard local elliptic a priori estimates and by arguments from the local L2 error analysis for finite elements as provided in Nitsche & Schatz [18]: k∇2 gk2;Ω\Ωδ/2 ≤ ck∆gk2;Ω\Ω3δ/4 + cδ kgk2 , k∇(g − gh )k2;Ω\Ωδ/2 ≤ ck∇(g − ih g)k2;Ω\Ω3δ/4 + cδ kg − gh k2 , with constants cδ ≈ δ −1 . Now, the assertion follows by the interpolation estimate (5.3) and the other estimates already proven. Lemma 5.3. The following L1 a priori and error estimates hold k∇gk1 + k∇2 gk01 ≤ c| log(h)|,
(5.7)
gh )k01
(5.8)
2
k∇(g − gh )k1 + hk∇ (g − with a constant c independent of h and δ . Proof. The proof can be found in [9].
≤ ch| log(h)|,
15
FINITE ELEMENTS IN PARAMETER IDENTIFICATION
For the point a ∈ Ωd , there holds vh (a) = (vh , δh ) = (A(q)∇vh , ∇gh ) = (B∇ψ, ∇gh ). We employ a standard localization argument. Let ω ∈ C0∞ (Ω) be a a smooth function with the properties 0 ≤ ω ≤ 1,
ω|Ωδ/2 ≡ 1,
ω|Ω\Ωδ/4 ≡ 0.
With this notation, we have (B∇ψ, ∇gh ) = (B∇(ωψ), ∇gh ) + (B∇((1 − ω)ψ), ∇gh ) =: Σ1 + Σ2 . First, we estimate the term Σ1 . By integration by parts and observing that ψ|∂Ω = 0 , we obtain X (B∇(ωψ), ∇gh ) = (ωψ, −∇ · (B∇gh ))T + (ωψ, n · B∇gh )∂K\∂Ω , K∈Th
where n is the outward unit normal vector to ∂K . Let [∇gh ] denote the jump of the gradient across the interior faces Γ ⊂ ∂K . Using this notation, we obtain X Σ1 ≤ kψk∞;Ωδ/2 k∇ · (B∇gh )k1;K + 21 kn · [B∇gh ]k1;∂K\∂Ω . K∈Th
First, the estimates of Lemma 5.3 yield X k∇ · (B∇gh )k1;K ≤ c |||B|||1,∞ k∇gh k1 + k∇2 gh k01 ≤ c |||B|||1,∞ | log(h)|. K∈Th
Next, observing that g ∈ H 2 (Ω) and therefore [B∇g] = 0 , we obtain by a trace theorem: X X kn · [B∇gh ]k1;K\∂Ω = kn · [B∇(gh − g)]k1;∂K\∂Ω K∈Th
K∈Th
≤ c |||B|||1,∞
X
h
−1
k∇(g − gh )k1;K + k∇2 (g − gh )k01;K .
K∈Th
Hence, collecting the foregoing estimates, Σ1 ≤ c |||B|||1,∞ kψk∞;Ωδ/2 | log(h)| + h−1 k∇(g − gh )k1 + k∇2 (g − gh )k01 . Again the estimates of Lemma 5.3, we obtain Σ1 ≤ c |||B|||1,∞ kψk∞;Ωδ/2 | log(h)|.
(5.9)
For the term Σ2 , we estimate as follows: Σ2 = (B∇((1 − ω)ψ), ∇(gh − g)) + (B∇((1 − ω)ψ), ∇g) ≤ |||B|||1,∞ k∇ψk2 k∇(gh − g)k2;Ω\Ωδ/2 + cδ kψk2 k∇(gh − g)k2 + c |||B|||1,∞ kψk2 k∇2 gk2;Ω\Ωδ/2 + k∇gk2;Ω\Ωδ/2 . Then, by the L2 estimates of Lemmas 5.1 and 5.2, it follows that Σ2 ≤ cδ |||B|||1,∞ kψk2 + hk∇ψk2 . This completes the proof of the theorem.
(5.10)
16
R. RANNACHER AND B. VEXLER
6. Numerical results. In this section, we discuss a sample problem confirming the a priori error estimate of Theorem 4.1. The state equation is given by −∇ · (A(q)∇u) = 2 in Ω, u = 0 on ∂Ω,
(6.1)
where Ω is the unit square. The matrix A(q) is a function of the parameter q = (q1 , q2 ) ∈ Q = R2 , given by 2 q1 q1 q2 . A(q) = q1 q2 exp(q2 ) The parameters are estimated from the measurements of the state variable at nine different points ξi ∈ Ω , see Figure 6.1.
Fig. 6.1. The computational domain with measurement points marked by circles.
The vector of measurements Cˆ is given by Cˆi = Ci (S(ˆ q ))(1 + εi ) i = 1, . . . , 9, where the reference parameter is qˆ = (5, 6) and ε = (εi ) describes the data perturbation. We consider two cases: a) ε = 0,
b) ε ≈ (0.12, −0.26, 0.29, −0.37, −0.49, 0.13, −0.04, −0.45, 0.20).
For case (a) the cost functional J(u) vanishes in the optimal point. Case (b) is more realistic because of the “measurement errors” modelled by a randomly chosen ε . Moreover, in this case, in contrast to case (a) the solution q of the corresponding parameter identification problem and the reference parameter qˆ differ. The parameter identification problem is discretized using bilinear finite elements on uniformly refined meshes. The results are listed in Table 6.1 and Table 6.2. For both cases the theoretically predicted orders of convergence are achieved. 7. Conclusions and Extensions. In this paper we have derived an a priori error estimate for the finite element discretization of an elliptic discrete parameter identification problem with pointwise measurements. The crucial point in our argument is the stability estimate of Theorem 4.2. The result of Theorem 4.1 can be extended as far as such a stability estimate is available. We list some possible directions of generalization. 1. More general meshes: For simplicity, we have assumed a quasi uniform mesh family {Th }h . The analysis can be extended to locally refined meshes, provided that
FINITE ELEMENTS IN PARAMETER IDENTIFICATION
17
Table 6.1 The error and the order of convergence with respect to the components of q without data perturbation.
N 81 289 1089 4225 16641 66049 order
q1 − q1,h 5.955e-1 1.407e-1 3.436e-2 8.509e-3 2.098e-3 4.993e-4 2.05
q2 − q2,h 9.902e-4 1.731e-4 4.343e-5 1.080e-5 2.668e-6 6.352e-7 2.04
Table 6.2 The error and the order of convergence with respect to the components of q with data perturbation.
N 81 289 1089 4225 16641 66049 order
q1 − q1,h 2.059e-0 5.172e-1 1.467e-1 3.771e-2 9.832e-3 2.350e-3 2.01
q2 − q2,h 1.874e-2 2.999e-3 8.341e-4 2.111e-4 5.640e-5 1.348e-5 1.98
the ratio of hmin and h = hmax is polynomial, hmin ≈ hp , with some p ≥ 1. For such meshes the stability result of Theorem 4.2 holds true with | log(hmin )| ≈ p| log(h)|. 2. More general domains: Our argument uses that the solution operator S(·) maps Q into H01 (Ω) ∩ H 2 (Ω) which is guaranteed on smoothly bounded or convex domains. In the case of a domain with reentrant corners or edges this regularity property is lost. This lack of regularity of the solution can be compensated by an appropriate refinement of the mesh near the critical corner points or edges. The stability estimate of Theorem 4.2 also holds in this situation. This is true in two as well as in tree dimensions. This will be shown in a forthcoming paper. 3. Higher order approximation: The result of Theorem 4.1 can be also extended to the case of higher-order finite elements, similar to the analysis of Nitsche [17]. In this case the logarithmic factor | log(h)| can be dropped in the stability Theorem 4.2. 4. More general equations: Theorem 4.1 can also be extended to more general elliptic equations or systems of the form −∇ · (A(q)∇u + b1 (q)u) + b2 (q)u = f, with parameter depended coefficients A(q), b1 (q), b2 (q). REFERENCES ¨ ltzsch Error estimates for the numerical approximation of a [1] N. Arada, E. Casas and F. Tro semilinear elliptic control problem, J. Comput. Optim. Appl., 23 (2002), pp. 201–229.
18
R. RANNACHER AND B. VEXLER
[2] H. T. Banks and K. Kunisch Estimation Techniques for Distributed Parameter Systems, Boston: Birkh¨ auser, 1989. [3] R. Becker, M. Braack and B. Vexler Numerical parameter estimation for chemical models in multidimensional reactive flows, Combustion Theory and Modelling, to appear, 2004. [4] R. Becker and B. Vexler A posteriori error estimation for finite element discretization of parameter identification problems, Numer. Math. 96 (2004), pp. 435–459. [5] S. Brenner and R. L. Scott The mathematical Theory of Finite Element Methods, Springer, Berlin Heidelberg New York, 1994. [6] K. Deckelnick and M. Hinze Error estimates in space and time for tracking-type control of the instationary Stokes system, in Proc. 8th Conference on Control of Distributed Parameter Systems Graz, July 15-21, 2001, Intern. Series of Numerical Mathematics 143, 2002, pp. 87–104. [7] R. S. Falk Approximation of a class of optimal control problems with order of convergence estimates, J. Math. Anal. Appl., 44, (1973), pp. 28–47. [8] R. S. Falk Error estimates of the numerical identification of a variable coefficient, Math. Comp. 40 (1983), pp. 537–546. [9] J. Frehse and R. Rannacher Eine L1 -Fehlerabsch¨ atzung f¨ ur diskrete Grundl¨ osungen in der Methode der finiten Elemente, in Tagungsband Finite Elemente, Bonn. Math. Schr. 89 (1976), pp. 92–114. [10] P. Grisvard Singularities in Boundary Value Problems, Masson, Paris, and Springer, Berlin, 1992. [11] M. Gunzburger and L. Hou Finite-dimensional approximation of a class of constrained nonlinear optimal control problems, SIAM J. Control Optim., 34 (1996), pp. 1001–1043. [12] C. Johnson Numerical Solution of Partial Differential Equations by Finite Element Method, Cambridge University Press, Cambridge, 1987. ¨rkka ¨inen Error estimates for distributed parameter identification in linear elliptic equa[13] T. Ka tions, J. Math. Syst. Estim. Control, 6 (1996), pp. 117–120. [14] C. Kravaris and J. H. Seinfeld Identification of parameters in distributed parameter systems by regularization, SIAM J. Control Optim., 23 (1985), pp. 217–241. [15] V. G. Litvinov Optimization in elliptic problems with applications to mechanics of deformable bodies and fluid mechanics, Operator Theory: Advances and Applications. 119. Basel: Birkh¨ auser, 2000. ¨ki and X.-C. Tai Error estimates for numerical identification of distributed [16] P. Neittaanma parameters, J. Comput. Math., 10 (1992), pp. 66–78. ¨ [17] J. Nitsche Uber L∞ -Absch¨ atzungen von Projektoren auf Finite Elemente, Tagungsband Finite Elemente, Bonn. Math. Schr. 89 (1976), pp. 13–30 . [18] J. Nitsche and A. Schatz Interior estimates for Ritz-Galerkin methods, Math. Comput., 28 (1976), pp. 937–958. [19] J. Nocedal and S. J. Wright Numerical Optimization, Springer Series in Operations Research, Springer New York, 1999. [20] R. Rannacher L∞ -Stability estimates and asymptotic error expansion for parabolic finite element equations, in Tagungsband Extrapolation and Defect Correction, Bonn. Math. Schr.,228 (1991), pp. 74–94. [21] B. Vexler Adaptive Finite Element Methods for Parameter Identification Problems, PhD thesis, University of Heidelberg, 2004.
