Data assimilation for the heat equation using stabilized finite element methods

We consider data assimilation for the heat equation using a finite element space semi-discretization. The approach is optimization based, but the design of regularization operators and parameters rely on techniques from the theory of stabilized finite elements. The space semi-discretized system is shown to admit a unique solution. Combining sharp estimates of the numerical stability of the discrete scheme and conditional stability estimates of the ill-posed continuous pde-model we then derive error estimates that reflect the approximation order of the finite element space and the stability of the continuous model. Two different data assimilation situations with different stability properties are considered to illustrate the framework. Full detail on how to adapt known stability estimates for the continuous model to work with the numerical analysis framework is given in “Appendix”.


Introduction
We consider two data assimilation problems for the heat equation where T > 0 and ⊂ R n is open. Let ω, B ⊂ be open and let 0 < T 1 < T 2 ≤ T . Both the data assimilation problems are of the following general form: determine the restriction u| (T 1 ,T 2 )×B of a solution to the heat equation (1) given f and the restriction u| (0,T )×ω . The crucial difference between the two problems is that in the first problem we do not make any assumptions on the boundary condition satisfied by u, whereas in the second problem we assume that it satisfies the lateral Dirichlet boundary condition u| (0,T )×∂ = 0. We emphasize that, in both the problems, no information on the initial condition satisfied by u is assumed to be known.
For the first problem we assume the following geometry, (A1) B is connected, ω ⊂ B, and the closure B is compact and contained in , and in the second case we may choose B = while assuming that (A2) is a compact, convex, polyhedral domain.
The first problem is conditionally Hölder stable, and this stability is optimal. More precisely, the following continuum stability estimate holds.
In Theorem 1, the factor with the power κ can be viewed as the size of the data q = u| (0,T )×ω , f = Lu, (3) and the norm u L 2 ((0,T )× ) as an apriori bound for the unknown u. Let us also emphasize that the assumption B ⊂ is essential, indeed, if B ∩ ∂ = ∅ then the optimal stability is conditionally logarithmic. The second problem is stable.
The estimates in Theorems 1 and 2 are tailored to work well with the finite element discretizations that are the main topic of the present paper. We review the literature on continuum stability estimates and give the proofs bridging the gap between the existing estimates and those in Theorems 1 and 2 in an "Appendix".
In what follows, we call the two data assimilation problems the unstable and the stable problem, respectively. We use the shorthand notation a(u, z) = (∇u, ∇z), G f (u, z) = (∂ t u, z) + a(u, z) − f, z , G = G 0 , where (·, ·) is the inner product of L 2 ((0, T ) × ) and ·, · is the dual pairing between L 2 (0, T ; H −1 ( )) and L 2 (0, T ; H 1 0 ( )). Note that for u ∈ H 1 ((0, T ) × ), the equations give the weak formulation of ∂ t u− u = f . Our approach to solve the data assimilation problem, in both the cases, is based on minimizing the Lagrangian functional where the data q and f are fixed. Here · ω is the norm of L 2 ((0, T ) × ω), and s and s * are symmetric bilinear forms, that are chosen differently in the two cases, and that we call the primal and dual stabilizers, respectively. Note that minimizing L q, f can be seen as fitting u| (0,T )×ω to the data q under the constraint (5), z can be interpreted as a Lagrange multiplier, and s/2 and s * /2 as regularizing penalty terms. In this paper we consider only semi-discretizations, that is, we minimize L q, f on a scale of spaces that are discrete in the spatial variable but not in the time variable. The spatial mesh size h > 0 will be the only parameter controlling the convergence of the approximation, and we use piecewise affine finite elements.
An important feature of the present work is that the choice of the regularizing terms is driven by the analysis and designed to give error estimates that reflect the approximation properties of the finite element space and the stability of the continuous model. In particular, when the continuous model is Lipschitz stable, we obtain optimal error estimates in the usual sense in numerical analysis, that is, the estimates are optimal compared to interpolation of the exact solution. This is the case for the second model problem introduced above. The regularization is constructed on the discrete level, that is, s and s * are defined on the semi-discrete spaces, and in some cases, they may not even make sense on the continuous level.
We show how the different stability properties of the two model problems lead to different regularization operators. The choice for a given problem is not unique, and the design of regularizing terms leading to optimal estimates can also be driven by computational considerations, such as computational cost or couplings in the discrete formulations. Therefore we present an abstract framework for the design of regularization operators. When measurement errors are present we also show how to quantify the damping of perturbations.
For the unstable problem, under suitable choices of the semi-discrete spaces and the regularization, we show that (6) has a unique minimizer (u h , z h ), h > 0, satisfying the following convergence rate where u is the unique solution of the continuum data assimilation problem, κ ∈ (0, 1) is the constant in Theorem 1, and the norms on the right-hand side capture the assumed apriori smoothness, see Theorem 3 below for the precise formulation. For the stable problem, we establish an analogous result but with κ = 1, that is, the convergence rate reflects now the stability estimate in Theorem 2. Moreover, in the stable case, we can replace the norm on the left-hand side of (7) with the stronger norm in Theorem 2, see Theorem 6 below for the precise formulation. We consider the time discretization of the stable problem in a separate paper [10]. The approach there is an extension of the strategy in the present paper: a Lagrangian of the form (6) is first discretized and then the stabilizers are chosen in the discrete spaces. The main difference is that an additional stabilization in time is required. In [10] we describe also two computational implementations of the method: one using an off-the-self linear solver, and the other using an iterative scheme that is based on the interpretation of the Euler-Lagrange equation for (6) as a system of two coupled heat equations. Let us also remark that the convergence analysis of our method is unchanged if the stabilizers in (6) are rescaled. We provide a computational study of such rescaling in [10].

Previous literature
The classical approach to data assimilation and inverse problems is to introduce Tikhonov regularization on the continuous level [2,15,26] and then discretize the well-posed regularized model. The use of Carleman estimates for uniqueness of inverse problems and the quantification of stability was first proposed in [7]. The use of conditional stability to choose a Tikhonov regularization parameter has been explored in [13]. Nevertheless a common feature of all methods where regularization takes place on the continuous level is that the accuracy of the reconstruction will be determined by the regularization error, leaving no room to exploit the interplay between discretization and regularization.
The method we propose is an instance of the so-called 4DVAR method [21]. The analysis in the present paper, which draws on previous ideas in the stationary case [8,9,11], is to our best knowledge the first complete numerical analysis of a 4DVAR type method. The main characteristic of our approach is to separate the numerical stability and the stability of the continuous model, and use the continuous stability estimate in the perturbation analysis.
In contrast, another recent approach to parabolic data assimilation problems is to derive Carleman estimates directly on the discrete scheme [3,4], which may then be used for convergence analysis. Other methods for data assimilation include the so called back and forth nudging [1] and methods using null-controllability [18,23]. In the former case, forward and backward solves are combined with filtering techniques (or penalty) to drive the approximate solution close to data, and the latter case leads to an approximation algorithm which uses auxiliary optimal control problems. For another example of a finite element method for identification of the initial data using Tikhonov regularization see [12] and for a discussion of different approaches to numerical approximation of control and inverse problems we refer to the monograph [17, Chapter 5].

Spatial discretization
We begin by observing that we may assume without loss of generality that ⊂ R n is a compact and connected polyhedral domain also in the case (A1). Indeed, we can choose a compact and connected polygonal domain˜ ⊂ such that B ⊂˜ and then replace by˜ .
Consider a family T = {T h ; h > 0} of triangulations of consisting of simplices such that the intersection of any two distinct simplices is either a common vertex, a common edge or a common face. Moreover, assume that the family T is quasi uniform, see e.g. [16,Def. 1.140], and indexed by The family T satisfies the following trace inequality, see e.g. [6,Eq. 10.3.9], where the constant C is independent of K ∈ T h and h > 0. Let V h be the H 1 -conformal approximation space based on the P 1 finite element, that is, where P 1 (K ) denotes the set of polynomials of degree less than or equal to 1 on K . We recall that the family T satisfies the following discrete inverse inequality, see e.g. [16,Lem. 1.138], where the constant C is independent of K ∈ T h and h > 0. We choose an interpolator that satisfies the following stability and approximation properties where k = 1, 2 and m = 0, k − 1, and that preserves vanishing Dirichlet boundary conditions, that is, using the notation A possible choice is the Scott-Zhang interpolator [25].

Spatial jump stabilizer
In the unstable case, the regularization will be based on the spatial jump stabilizer that we will introduce next. We denote by F h the set of internal faces of T h , and define for where K 1 , K 2 ∈ T h are the two simplices satisfying K 1 ∩ K 2 = F, and n j is the outward unit normal vector of K j , j = 1, 2. We define the jump stabilizer where ds is the Euclidean surface measure on F.
Proof The original results on Poincaré inequalities for piecewise H 1 -functions may be found in [5]. For a detailed proof of (15) see [11].
Proof Towards (16) we integrate by parts and recall that u h is an affine function on each element to obtain Thus by the Cauchy-Schwarz inequality , and the inequality (16) follows by applying (8) to the last factor. Let us now turn to (17). As n · ∇w F = 0, we have using (8) Summing over all internal faces and observing that The first term on the right-hand side satisfies by (12).

The unstable problem
There are several possible choices for the primal and dual stabilizers s and s * in (6), and different choices lead to numerical methods that may differ in terms of practical performance. To illustrate the main ideas of our approach, we begin by considering a concrete choice of the stabilizers in Sect. 3.1, and give an abstract framework in Sect. 4.

A model case
We choose the following primal and dual stabilizers defined on the respective semi-discrete spaces We define also the semi-norm and norm Recall that W h = V h ∩ H 1 0 ( ) and whence the Poincaré inequality implies that · W h is indeed a norm on W h .

Lemma 3
The semi-norm | · | V h and norm · W h satisfy the following inequalities with a constant C > 0 that is independent from h > 0: lower bounds upper bounds and finally, for the semi-norm | · | V h only, the Poincaré type inequality Proof We have where the bound for a(u h , z) follows from (18), analogously to (16), after integration in time. This implies that (23) holds. The lower bound (24) follows from the Cauchy-Schwarz inequality together with the Poincaré inequality (22). Towards the upper bound (25), we have where the bound for the first term follows from (17). The bound for the second term holds by the H 1 -stability (11) of the interpolator π h . The upper bound (26) follows immediately from the H 1 -stability of π h , and the lower bound (27) follows from (15).
Observe that, by the lower bound (27), where A h is the symmetric bilinear form Note that and therefore A h is weakly coercive in the following sense The Babuska-Lax-Milgram theorem implies that the equation (28) has a unique solu-

Error estimates
In this section we show that the solution (u h , z h ) of (28) satisfies u h → u in (T 1 , T 2 )×B as h → 0, with the convergence rate (7). Here u is a smooth enough solution of the unstable data assimilation problem in the continuum. We use the shorthand notation (28). Then there exists C > 0 such that for all h ∈ (0, 1), Proof The equations ∂ t u − u = f and u| ω = q are equivalent with and the Eqs. (28) and (31)

By (30) it is enough to show that
We use (12) to bound the first term in (32), for the second term we use (24) and (12), and for the third term we use (25), (28). Then there exists C > 0 such that for all h ∈ (0, 1) We arrive to the estimate after we bound the first term in (33) by using (23) and (11)-(12) the second term by using (12) 1) , and the last term by using (26) We bound the first term in (34) by using Lemma 4 and (25) and observe that the analogous bound for the second term follows immediately from Lemma 4. Thus By applying Theorem 1 to u h − u, and using Lemma 4 and (12) to bound the term 1) we see that It remains to show the bound The Poincaré type inequality (27) implies where we used again Lemma 4. The inequality (36) follows since π h u − u can be bounded by u * , in fact, even by Ch 2 u * .

The effect of perturbations in data
In practice the data f, q in (3) are typically polluted by measurement errors. Such effects can easily be included in the analysis and we show in this section how the above results must be modified to account for this case. We consider the case where f, q in (28) are replaced by the perturbed counterparts where u is the solution to the underlying unpolluted problem, that is, ∂ t u − u = f . We assume that the perturbations δq, δ f satisfy δ f ∈ H −1 ( ) and δq ∈ L 2 (ω), and use the following measure of the size of the perturbation δ(q,f ) = δq L 2 (0,T ;L 2 (ω)) + δ f (0,−1) .

The critical points
We consider a perturbation of q in L 2 since this is used in the data fitting term in the Lagrangian. It should be noted that the formulation (38) remains well defined for rougher perturbations δq, but they will lead to worse estimates in h. First we consider again the residual quantities as in Lemma 4.

Proof Proceeding as in the proof of Lemma 4, the Eqs. (38) and (31) imply for all
The only terms that are not covered by the previous analysis are those with the perturbations, and we have The choice (19) of the dual stabilizer s * guarantees that w (0,1) ≤ w W h .
The error estimate of Theorem 3 may now be modified to include the perturbations.

Theorem 4
Let ω, B ⊂ , 0 < T 1 < T 2 < T and κ ∈ (0, 1) be as in Theorem 1. Let u ∈ H (1,1) ∩ H (0,2) and definef andq by (37). Let (u h , z h ) ∈ V h × W h be the solution of (38). Then there exists C > 0 such that for all h ∈ (0, 1) Proof We proceed as in the proof of Theorem 3, but due to the perturbation, the Galerkin orthogonality (33) no longer holds exactly and we obtain for w ∈ H (0,1) 0 . Analogously to the proof of Theorem 3, we obtain We bound the first term in (40) by using Lemma 5 and (25) and the second term by using again Lemma 5. Thus As before, applying Theorem 1 to u h − u, and using Lemma 5 and (12) to bound the term u h − u ω leads to It remains to bound u h − u . The Poincaré type inequality (27) implies where we used again Lemma 5.
According to Theorem 4, the error of the numerical approximation can diverge if no functionũ exists such that ∂ tũ − ũ =f andũ| (0,T )×ω =q. On the other hand, in the context of a specific application, an estimate on the noise level might be available, and this can be used to choose a suitable mesh size. That is, assuming that there exists h 0 > 0 such that δ(q,f ) ≤ h 0 ( u * + f ),

A framework for stabilization
Before proceeding to the stable problem, we introduce an abstract stabilization framework based on the essential features of the model case in Sect. 3.1. Let s h and s * h be bilinear forms on the spaces V h and W h , respectively. Let | · | V h be a semi-norm on V h and let · W h be a norm on W h . We relax (21) by requiring only that s h and s * h are continuous with respect to | · | V h and · W h , that is, Let · * be the norm of a continuously embedded subspace H * of the energy space (2). The space H * encodes the apriori smoothness. We relax the lower bounds (23) and (24) as follows where π h is an interpolator satisfying (11)- (13). We replace the upper bound (25) by its abstract analogue and require that (26) and (27) hold as before. Finally, in the abstract setting, we assume that the weak coercivity (30) holds where A h and | · | h are defined as above.
It can be verified that the following analogue of Theorem 3 holds under these assumptions.

The stable problem
Let us now consider the stable problem with the additional lateral boundary condition u| (0,T )×∂ = 0. We may impose the boundary condition on the discrete level by changing the space V h , defined previously by (20), to In the stable case we do not need an inequality analogous to (36) since the estimate in Theorem 2 does not contain an apriori term, that is, a term analogous to u L 2 ((0,T )× ) in Theorem 1. In the previous section, the Poincaré type inequality (27) was used only to obtain (36), and whence we can relax the conditions there by dropping (27). However, we still need to require that | · | h is a norm on V h × W h .
We will now proceed to a concrete case. The choices in Sect. 4 work also in the stable case, however, the additional structure allows us to choose a weaker stabilization. The use of a weaker stabilization is motivated by the fact that it leads to a weaker coupling of the two heat equations in (28), which again can be exploited when solving (28) in practice. See [10] for a computational implementation based on this.
The stabilizers and semi-norms are chosen as follows, and we define H * = H (1,1) 0 . To counter the lack of primal stabilization on most of the cylinder (0, T ) × , we choose π h to be the orthogonal projection π h : H 1 0 ( ) → W h with respect to the inner product (∇u, ∇v) L 2 ( ) . That is, π h is defined by As is a convex polyhedron, this choice satisfies the estimates (11)- (12), see e.g. [16,.

Lemma 6
The choices (45)-(47) satisfy the properties (41)-(44), (26) and (30). More- Proof It is clear that the continuities (41) hold. We begin with the lower bound (42). By the orthogonality (47), Towards the lower bound (43), we use the orthogonality (47) as above, The bound (43) follows from the Poincaré inequality (22). The bound (44) follows from the continuity of the trace and the continuity of the projection π h . The bound (26) follows immediately from the continuity of π h . We turn to the weak coercivity (30). The essential difference between the unstable and the stable case is that in the latter case ∂ t u h ∈ W h when u h ∈ V h . We have and thus using bilinearity of A h , where α > 0. We will establish (30) by showing that there is α ∈ (0, 1) such that where w h = αh 2 ∂ t u h . Towards (50) we observe that We use the discrete inverse inequality (10) to bound the second term It remains to show (51). Towards bounding the first cross term in (49) we observe that Hence We use the arithmetic-geometric inequality, and the discrete inverse inequality (10) to bound the second cross term in (49), and therefore (51) holds with small enough α > 0. Finally, using the Poincaré inequality (22), we see that | (u h , z h ) | h = 0 implies that z h = 0 and u h (0, ·) = 0. As also ∂ t u h = 0, we have u h = 0, and therefore | · | h is a norm.
In the stable case we have the following convergence rate.

Theorem 6
Let ω ⊂ be open and non-empty and let 0 < T 1 < T . Suppose that (A2) holds. Let u ∈ H * and define f = ∂ t u − u and q = u| ω . Suppose that the primal and dual stabilizers satisfy (41)-(44), (26) and (30). Then (28) has a unique solution (u h , z h ) ∈ V h × W h , and there exists C > 0 such that for all h ∈ (0, 1) The proof is analogous to that of Theorem 3, however, we give it here for the sake of completeness.
Proof We begin again by showing the estimate The weak coercivity (30) implies that in order to show (52) it is enough bound the three terms on the right hand side of (32). For the first of them, that is, (u − π h u, v h ) ω , we use (12) as in the proof of Lemma 4. The lower bound (43) applies to the second term G(u − π h u, w h ), and for the third one we use the continuity (41) and the upper bound (44), We define the residual r as in the proof of Theorem 3, and show next that (34) holds. It is enough to bound the three terms on the right hand side of (33). The lower bound (42) applies to the first term G(u h , w − π h w), for the second term ( f, w − π h w) we use (12) as in the proof of Theorem 3, and for the third term we use the continuity (41) and the upper bound (26) The inequalities (34), (52) and (44) imply Finally, Theorem 2 implies that The claim follows by using (52) and (12), Here we used also the assumption that H * is a continuously embedded subspace of the energy space (2), namely, this implies that the embedding H * ⊂ H (0,1) is continuous.

Remark 1
If the data q, f is perturbed in the stable case, the data assimilation problem behaves like a typical well posed problem, that is, the term δ(q,f ) needs to be added on the right-hand side of the estimate in Theorem 6, but this time without any negative power of h.

Conclusion
In the present paper our aim was to show how methods known from the theory of stabilized finite element methods can be applied to the design of computational methods for non-stationary data assimilation problems by first discretizing and then regularizing the corresponding 4DVAR optimization system.
A key feature of this framework is that the error analysis is based on numerical stability of residual quantities that are independent of the stability properties of the continuous model. The residual quantities are then bounded by using conditional stability estimates, based on Carleman estimates, to derive error estimates that reflect the approximation properties of the finite element space and the stability of the continuous problem in an optimal way. An upshot of this approach is that it gives a clear lead on how to design regularization for a given problem in order to obtain the best accuracy of the approximation with the least effect of perturbations in data.
Observe that the stabilization operators proposed herein are not unique, for instance it is straightforward to show that the dual stabilizer in Eq. (19) may be chosen as the first part of the primal stabilizer, leading to similar error estimates for unperturbed data. The error estimates also gives an indication on what form Carleman estimates should take to make them immediately applicable for error analysis.

Appendix: Continuum estimates
Theorem 1 is a consequence of the so-called three cylinders inequality that goes back to [19], see [27,28] for an overview of the related literature. The variant of the inequality, needed in the above convergence analysis, seems not to appear in the literature, and we prove it here by using the Carleman estimate [20].
Theorem 2 is a variant of [14], the difference being that in [14] the domain is assumed to have C 2 -smooth boundary. We outline below the modifications needed in the case that is a convex polyhedron.
Proof of Theorem 1 We use notation from the proof of Theorem 7. By replacing ω with a smaller set we may assume without loss of generality that it is a ball of the above form B 1 . We will show a local version of the claimed estimate where B is replaced by a ball of the form B 2 . The general case follows by covering B by finite chains of balls starting from ω, and by iterating the local result. Let us first consider the case u ∈ C ∞ (R × ). We choose 0 < r 0 < r 1 and r 2 < r 3 < d(x 0 , ) and write B j = B(x 0 , r j ), j = 0, 3. Let 0 < 0 < and choose η ∈ C ∞ (R) such that η = 0 near the origin and η(t) = 1 for t > 0 . Let w be the solution of and we have shown the local estimate in the case that u is smooth. To conclude we observe that smooth functions are dense in the space (2).
Let us now turn to Theorem 2. In the case of convex, polygonal , the following simple lemma gives a weight function satisfying Condition 1.1 of [14].
The claim follows by applying (57) once more.