A stabilized finite element method for inverse problems subject to the convection–diffusion equation. I: diffusion-dominated regime

The numerical approximation of an inverse problem subject to the convection–diffusion equation when diffusion dominates is studied. We derive Carleman estimates that are of a form suitable for use in numerical analysis and with explicit dependence on the Péclet number. A stabilized finite element method is then proposed and analysed. An upper bound on the condition number is first derived. Combining the stability estimates on the continuous problem with the numerical stability of the method, we then obtain error estimates in local H1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H^1$$\end{document}- or L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document}-norms that are optimal with respect to the approximation order, the problem’s stability and perturbations in data. The convergence order is the same for both norms, but the H1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H^1$$\end{document}-estimate requires an additional divergence assumption for the convective field. The theory is illustrated in some computational examples.


Introduction
We consider the convection-diffusion equation (1) Lu where Ω ⊂ R n is open and bounded, µ > 0 is the diffusion coefficient and β ∈ [W 1,∞ (Ω)] n is the convective velocity field.For an open subset ω ⊂ Ω, let Ũω = u| ω + δ, where we assume u ∈ H 2 (Ω) is the solution to (1), and δ ∈ L 2 (ω) is some perturbation.The data assimilation (or unique continuation) problem consists in finding u given f and Ũω .
The aim is to design a finite element method for data assimilation with weakly consistent regularization applied to the convection-diffusion equation (1).In the present analysis we consider the regime where diffusion dominates and in the companion paper [7] we treat the one with dominating convective transport.To make this more precise we introduce the Péclet number associated to a given length scale l by P e(l) := |β|l µ , for a suitable norm | • | for β.If h denotes the characteristic length scale of the computation, we define the diffusive regime by P e(h) < 1 and the convective regime by P e(h) > 1.It is known that the character of the system changes drastically in the two regimes and we therefore need to apply different concepts of stability in the two cases.In the present paper we assume that the Péclet number is small and we use an approach similar to that employed for the Laplace equation in [5], for the Helmholtz equation in [8] and for the heat equation in [9], that is we combine conditional stability estimates for the physical problem with optimal numerical stability obtained using a bespoke weakly consistent stabilizing term.For high Péclet numbers on the other hand, we prove in [7] weighted estimates directly on the discrete solution, that reflect the anisotropic character of the convection-diffusion problem.
In the case of optimal control problems subject to convection-diffusion problems that are well-posed, there are several works in the literature on stabilized finite element methods.In [10] the authors considered stabilization using a Galerkin least squares approach in the Lagrangian.Symmetric stabilization in the form of local projection stabilization was proposed in [1] and using penalty on the gradient jumps in [19,15].The key difference between the well-posed case and the ill-posed case that we consider herein is that we can not use stability of neither the forward nor the backward equations.Crucial instead is the convergence of the weakly consistent stabilizing terms and the matching of the quantities in the discrete method and the available (best) stability of the continuous problem.Such considerations lead to results both in the case of high and low Péclet numbers, but the different stability properties in the two regimes lead to a different analysis for each case that will be considered in the two parts of this paper.
The main results of this current work are the convergence estimates with explicit dependence on the Péclet number in Theorem 1 and Theorem 2, that rely on the continuous three-ball inequalities in Lemma 2 and Corollary 2.

Stability estimates
We prove conditional stability estimates for the unique continuation problem subject to the convection-diffusion equation (1) in the form of three-ball inequalities, see e.g.[17] and the references therein.The novelty here is that we keep track of explicit dependence on the diffusion coefficient µ and the convective vector field β.The first such inequality is proven in Corollary 1, followed by Lemma 2 and Corollary 2, where the norms for measuring the size of the data are weakened to serve the purpose of devising a finite element method in Section 3.
First we prove an auxiliary logarithmic convexity inequality, which is a more explicit version of [16,Lemma 5.2].
Lemma 1. Suppose that a, b, c ≥ 0 and p, q > 0 satisfy c ≤ b and c ≤ e pλ a + e −qλ b for all λ > λ 0 ≥ 0. Then there are C > 0 and κ ∈ (0, 1) (depending only on p and q) such that Proof.We may assume that a, b > 0, since c = 0 if a = 0 or b = 0.The minimizer λ * of the function f (λ) = e pλ a + e −qλ b is given by and writing r = q/p, the minimum value is = r p/(p+q) + r −q/(p+q) a q/(p+q) b p/(p+q) .
That is, if where C 2 = r −q/(p+q) .As e qλ 0 ≥ 1 and The following Carleman inequality is well-known, see e.g.[16].For the convenience of the reader we have included an elementary proof in Appendix A.
Using the above Carleman estimate we prove a three-ball inequality that is explicit with respect to µ and β, i.e. the constants in the inequality are independent of the Péclet number.The corresponding inequality with constant depending implicitly on the Péclet number is proven for instance in [17].We denote by B(x, r) the open ball of radius r centred at x, and by d(x, ∂Ω) the distance from x to the boundary of Ω.
We now shift down the Sobolev indices in Corollary 1 by making a similar argument to that in Section 4 of [11] or Section 2.2 of [8], based on semiclassical pseudodifferential calculus. , where Proof.Let > 0 be the semiclassical parameter that satisfies = 1/τ , where τ is the parameter previously introduced in Proposition 1.The scale of semiclassical Bessel potentials is defined by , s ∈ R, and the semiclassical Sobolev spaces by We will make strong use of the following commutator and pseudolocal estimates for semiclassical pseudodifferential operators, see e.g.[20,Theorem 4.12].Suppose that η, ϑ ∈ C ∞ 0 (R n ) and that η = 1 near supp(ϑ), and let A, B be two semiclassical pseudodifferential operators of orders s, m, respectively.Then for all p, q, N ∈ R, there is C > 0, Let 0 < r j < r j+1 < d(x 0 , ∂Ω), j = 0, . . ., 4 and B j = B(x 0 , r j ), keeping B 1 , B 2 unchanged.Let rj ∈ (r j−1 , r j ) and Bj = B(x 0 , rj ), j = 0, . . ., 3, where As in Appendix A, by taking ℓ = φ/ and σ = ∆ℓ + 3αλφ/ , we obtain Scaling this with µ 2 4 , we insert the convective term and obtain that The last two terms in the right-hand side can be absorbed by the first one when , where we used (5) to absorb one term by the left-hand side.From ( 8) and (7) we have and the commutator estimate (4) gives Recalling the assumption (6), these terms can be absorbed by the left-hand side of ( 9), obtaining We now combine this estimate with the technique used to prove Corollary 1.Consider where we used the norm inequality Letting Φ(r) = e −αr and using a similar argument as in the proof of Corollary 1, we find that when satisfies (6) and is small enough.Absorbing the negative power of in the exponential, we then use Lemma 1 and conclude by absorbing the P e = 1 + |β|/µ factor into the exponential factor e C P e 2 .
Making the additional coercivity assumption ∇ • β ≤ 0, we can weaken the norms just in the right-hand side of Corollary 1 by using the stability estimate for a well-posed convectiondiffusion problem with homogeneous Dirichlet boundary conditions.Corollary 2. Let x 0 ∈ Ω and 0 < r 1 < r 2 < d(x 0 , ∂Ω).Define B j = B(x 0 , r j ), j = 1, 2. Then there are C > 0 and κ ∈ (0, 1) such that for µ > 0, , where Proof.Let the balls B 0 , B 3 ⊂ Ω such that B j ⊂ B j+1 , for j = 0, 2. Consider the well-posed problem Lw = Lu in B 3 , w = 0 on ∂B 3 .Since ess sup Ω ∇ • β ≤ 0, as a consequence of the divergence theorem we have Taking v = u − w, we have Lv = 0 in B 3 .The stability estimate in Corollary 1 used for B 0 , B 2 , B 3 reads as , and the following estimates hold and we obtain The same argument for B 3 ⊂ Ω gives thus leading to the conclusion after absorbing the P e = 1 + |β|/µ factor into the exponential factor e C P e 2 .

Finite element method
Let V h denote the space of piecewise affine finite element functions defined on a conforming computational mesh T h = {K}.T h consists of shape regular triangular elements K with diameter h K and is quasi-uniform.We define the global mesh size by h = max K∈T h h K .The interior faces of the triangulation will be denoted by F i , the jump of a quantity across a face F by • F , and the outward unit normal by n.
We will denote by C a generic positive constant independent of the mesh size and the Péclet number.Let π h : L 2 (Ω) → V h denote the standard L 2 -projection on V h , which for k = 1, 2 and m = 0, k − 1 satisfies We introduce the standard inner products with the induced norms and the following bilinear forms and The terms s Ω and s * are stabilizing terms, while the term s ω is aimed for data assimilation.
After scaling with the coefficients in the above forms, Lemma 2 in [6] writes as (12) (µ , ∀v h ∈ V h , and Lemma 2 in [9] gives the jump inequality (13) s , ∀u ∈ H 2 (Ω).The parameters γ and γ * in s Ω and s * , respectively, are fixed at the implementation level and, to alleviate notation, our analysis covers the choice γ = 1 = γ * .
We can then use the general framework in [4] to write the finite element method for unique continuation subject to (1) as follows.
where we recall that Ũω = u| ω + δ and u ∈ H 2 (Ω) is the solution to (1).We observe that by the ill-posed character of the problem, only the stabilization operators s Ω and s * provide some stability to the discrete system, and the corresponding system matrix is expected to be ill-conditioned.To quantify this effect we first prove an upper bound on the condition number.
Proposition 2. The finite element formulation (14) has a unique solution (u h , z h ) ∈ [V h ] 2 and the Euclidean condition number κ 2 of the system matrix satisfies

Proof. We write (14) as the linear system
Since A[(u h , z h ), (u h , −z h )] = s(u h , u h ) + s * (z h , z h ), using (12) the following inf-sup condition holds This provides the existence of a unique solution for the linear system.We use [13, Theorem 3.1] to estimate the condition number by ( 15) where .
We recall the following discrete inverse inequality, see for instance [12, Lemma 1.138], ( 16) We also recall the following continuous trace inequality, see for instance [18], and the discrete one Using the Cauchy-Schwarz inequality together with ( 18) and ( 16) we get s Ω (u h , v h ) = γµ(1 + P e(h)) Combining this with the Cauchy-Schwarz inequality and the inequalities ( 16) and ( 17), we obtain Again due to the Cauchy-Schwarz inequality, and trace and inverse inequalities, we have Collecting the above estimates we have Υ h ≤ Cµ(1+P e(h))h −2 , and we conclude by (15).

3.1.
Error estimates for the weakly consistent regularization.The error analysis proceeds in two main steps: (i) First we prove that the stabilizing terms and the data fitting term must vanish at an optimal rate for smooth solutions, with constant independent of the physical stability (Proposition 3).(ii) Then we show that the residual of the PDE is bounded by the stabilizing terms and the data fitting term.Using this result together with the first step and the continuous stability estimates in Section 2, we prove L 2 -and H 1 -convergence results (Theorems 1 and 2).To quantify stabilization and data fitting for We also define the "continuity norm" on H 3 2 +ǫ (Ω), for any ǫ > 0, Using standard approximation properties and the trace inequality ( 18), we have Using (13) and interpolation , where we used that s Ω (u, v h ) = 08, since u ∈ H 2 (Ω).Hence it follows that for u ∈ H 2 (Ω) (19) (u Observe that, when P e(h) < 1, the first term dominates and the estimate is O(h), whereas when P e(h) > 1 the bound is O(h 2 ).We note in passing that the same estimates hold for the nodal interpolant.
Lemma 3 (Consistency).Let u ∈ H 2 (Ω) be the solution to (1) Proof.The first claim follows from the definition of a h , since where in the last equality we integrated by parts.The second claim follows similarly from Lemma 4 (Continuity).Assume the low Péclet regime (11) Proof.Writing out the terms of a h and integrating by parts in the advective term leads to Using the Cauchy-Schwarz inequality and the trace inequality (17) for v, we see that By the Cauchy-Schwarz inequality and a discrete Poincaré inequality for w h , see e.g.[2], we bound Under the assumption |β| 1,∞ ≤ C|β|, we get We bound the remaining term by Finally, exploiting the low Péclet regime P e(h) < 1, we obtain the conclusion.
Proposition 3 (Convergence of regularization).Assume the low Péclet regime (11) and that |β| 1,∞ ≤ C|β|.Let u ∈ H 2 (Ω) be the solution to (1) Using both claims in Lemma 3 we may write The other terms are simply bounded using the Cauchy-Schwarz inequality as follows Collecting the above bounds we have and the claim follows by applying the approximation (19).
Proof.Let us consider the residual defined by r, w = a h (u h , w) − f, w , for w ∈ H 1 0 (Ω).Using (14) we obtain We split the first term in the right-hand side into convective and non-convective terms, and for the latter we integrate by parts on each element K and use Cauchy-Schwarz followed by the trace inequality (17) to get Using (21) and interpolation we obtain We bound the convective term in a h (u h , w − π h w) by Lemma 5, hence obtaining The next term in the residual is bounded by The last term left to bound from the residual is s * (z h , π h w) ≤ (0, z h ) s (0, π h w) s , and using (18) for the jump term, together with the H 1 -stability of π h , we see that where for the boundary term we used that w| ∂Ω = 0 together with interpolation and (17).Bounding (0, z h ) s by Proposition 3, we get Collecting the above estimates we bound the residual norm by We now use the stability estimate in Lemma 2 to write By Proposition 3 we have Using (12) and Proposition 3 again, we bound Hence we conclude by where we have absorbed the P e = 1 + |β|/µ factor into the exponential factor e C P e 2 .
Proof.Letting e h = u − u h , we combine the proof of Theorem 1 with the stability estimate in Corollary 2 to obtain

Numerical experiments
We illustrate the theoretical results with some numerical examples.The implementation of the stabilized FEM (14) has been carried out in FreeFem++ [14] on uniform triangulations with alternating left and right diagonals.The mesh size is taken as the inverse square root of the number of nodes.The parameters in s Ω and s * are set to γ = 10 −5 and γ * = 1.We also rescale the boundary term in s * by the factor 50, drawing on results from different numerical experiments.In this section we denote e h = π h u − u h .
We consider Ω to be the unit square and the exact solution with global unit L 2 -norm u(x, y) = 30x(1 − x)y(1 − y).
We take the diffusion coefficient µ = 1 and investigate two cases for the convection field: the coercive case of the constant field β c = (1, 0), and the case β nc = 100(x + y, y − x), plotted in Figure 1, for which ∇ • β = 200 and β 0, ∞ = 200.This makes the (well-posed) problem strongly non-coercive with a medium high Péclet number.The latter example was also considered in [4] for numerical experiments on a non-coercive convection-diffusion equation with Cauchy data.We consider the following domains for data assimilation, shown in Figure 2,  The condition number upper bound in Proposition 2 is illustrated for a particular case in Figure 1, where we plot the condition number κ 2 versus the mesh size h, together with reference dotted lines proportional to h −3 and h −4 .For five meshes with 2 N elements on each side, N = 3, . . ., 7, the approximate rates for κ 2 are -3.03,-3.16, -3.2, -3.34.Comparing the geometries in ( 22) and (23) we also expect to see different effects of the two convective fields β c and β nc .Notice that for both geometries the horizontal magnitude of β nc is greater than that of β c .In (22) the solution is continued in the crosswind direction for both β c and β nc , and a stronger convective field is not expected to improve the reconstruction.On the other side, in (23) information is propagated both downstream and with a j , b j > 0 depending only on α, φ and λ.Taking τ large enough and using the elementary inequality |∇v| 2 = e 2τ φ |τ w∇φ + ∇w| 2 ≥ e 2τ φ 1 2 |∇w| 2 − e 2τ φ |∇φ| 2 τ 2 w 2 , we conclude by integrating over K and using the divergence theorem.