A mixed formulation for the direct approximation of L2-weighted controls for the linear heat equation

A mixed formulation for the direct approximation of L2-weighted controls for the linear heat equation Arnaud M¨ unch∗Diego A. Souza† June 2, 2014 Abstract This paper deals with the numerical computation of null controls for the linear heat equation. The goal is to compute approximations of controls that drive the solution from a prescribed initial state to zero at a given positive time. In [Fernandez-Cara & M¨unch, Strong convergence approximations of null controls for the 1D heat equation, 2013], a so-called primal method is described leading to a strongly convergent approximation of distributed control: the controls minimize quadratic weighted functionals involving both the control and the state and are obtained by solving the corresponding optimality conditions. In this work, we adapt the method to approximate the control of minimal square integrable-weighted norm. The optimality conditions of the problem are reformulated as a mixed formulation involving both the state and its adjoint. We prove the well-posedeness of the mixed formulation (in particular the inf-sup condition) then discuss several numerical experiments. The approach covers both the boundary and the inner situation and is valid in any dimension. Keywords: Linear heat equation; Null controllability; Finite element methods; Mixed formulation. Mathematics Subject Classification (2010)- 35K35, 65M12, 93B40, 65K10 Contents 1 Introduction. The null controllability problem 2 2 Control of minimal weighted L2-norm: mixed reformulations 5 2.1 The penalized case: Mixed formulation I . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1.1 Mixed formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1.2 Dual problem of the extremal problem (15) . . . . . . . . . . . . . . . . . . 8 2.2 The penalized case : Mixed formulation II (relaxing the condition L⋆ϕε= 0 in L2(QT)).......................................... 10 2.3 Third mixed formulation of the controllability problem : the limit case ε= 0. . . . 13 ∗Laboratoire de Math´ematiques, Universit´e Blaise Pascal (Clermont-Ferrand 2), UMR CNRS 6620, Campus de C´ezeaux, 63177, Aubi`ere, France. E-mails: [email protected]. †Dpto. EDAN, University of Sevilla, 41080 Sevilla, Spain and Departamento de Matem´atica, Universidade Federal da Para´ıba, 58051-900, Jo˜ao Pessoa–PB, Brazil. E-mail: [email protected]. Partially supported by CAPES (Brazil) and grant MTM2010-15592 (DGI-MICINN, Spain). 1 1 INTRODUCTION. THE NULL CONTROLLABILITY PROBLEM 2 3 Numerical approximation and experiments 16 3.1 Discretization of the mixed formulation (16) . . . . . . . . . . . . . . . . . . . . . . 16 3.2 Normalization and discretization of the mixed formulation (32) . . . . . . . . . . . 18 3.3 The discrete inf-sup test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.4 Numerical experiments for the mixed formulation (16) . . . . . . . . . . . . . . . . 21 3.5 Conjugate gradient for J⋆⋆ ε,r ............................... 24 3.6 Numerical experiments for the mixed formulation (32) - limit case ε= 0. . . . . . . 27 4 Concluding remarks and Perspectives 31 A Appendix 32 A.1 Appendix : Fourier expansion of the control of minimal L2(ρ0, qT) norm. . . . . . . 32 A.2 Appendix: Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 1 Introduction. The null controllability problem Let Ω ⊂RNbe a bounded connected open set whose boundary ∂Ω is regular enough (for instance of class C2). Let ω⊂Ω be a (small) nonempty open subset and assume that T > 0. In the sequel, for any τ > 0 we denote by Qτ,qτand Στthe sets Ω×(0, τ), ω×(0, τ) and ∂Ω×(0, τ), respectively. This work is concerned with the null controllability problem for the heat equation    yt−∇·(c(x)∇y) + d(x, t)y=v1ω,in QT, y= 0,on ΣT, y(x, 0) = y0(x),in Ω. (1) Here, we assume that c:= (ci,j)∈C1(Ω; MN(R)) with (c(x)ξ, ξ)≥c0|ξ|2in Ω (c0>0), d∈ L∞(QT) and y0∈L2(Ω); v=v(x, t) is the control (a function in L2(qT)) and y=y(x, t) is the associated state. Moreover, 1ωis the characteristic function associated to the set ω. In the sequel, we shall use the following notation : L y := yt−∇·(c(x)∇y) + d(x, t)y, L⋆ϕ:= −ϕt−∇·(c(x)∇ϕ) + d(x, t)ϕ. For any y0∈L2(Ω) and v∈L2(qT), there exists exactly one solution yto (1), with the regularity y∈C0([0, T]; L2(Ω)) ∩L2(0, T;H1 0(Ω)) (see [29, 7]). Accordingly, for any final time T > 0, the associated null controllability problem at time Tis the following : for each y0∈L2(Ω), find v∈L2(qT) such that the corresponding solution to (1) satisfies y(·, T) = 0 in Ω.(2) The controllability of PDEs is an important area of research and has been the subject of many papers in recent years. Some relevant references are [26, 28, 34] and [12]. In particular, we refer to [20] and [27] where the null controllability of (1) is proved. The numerical approximation is also a fundamental issue, since it is not in general possible to get explicit expression of controls. Due to the strong regularization property of the heat kernel, numerical approximation of controls is a rather delicate issue. The same holds in inverse problems theory when parabolic equations and systems are involved (see [15]). This have been exhibited numerically in [5] who made use of duality argument and focused on the control of minimal square integrable norm: the problem reads      Minimize J1(y, v) := 1 2ZZqT|v(x, t)|2dx dt Subject to (y, v)∈ C(y0, T) (3) 1 INTRODUCTION. THE NULL CONTROLLABILITY PROBLEM 3 where C(y0;T) denotes the linear manifold C(y0;T) := {(y, v) : v∈L2(qT), y solves (1) and satisfies (2) }. The earlier contribution is due to Glowinski and Lions in [22] (updated in [23]) and relies on duality arguments. Duality allows to replace the original constrained minimization problem by an unconstrained and a priori easier minimization (dual) problem. The dual problem associated with (3) is : min ϕT∈H J⋆ 1(ϕT) := 1 2ZZqT|ϕ(x, t)|2dxdt +ZΩ y0(x)ϕ(x, 0)dx (4) where the variable ϕsolves the backward heat equation : L⋆ϕ= 0 in QT, ϕ = 0 on ΣT;ϕ(·, T) = ϕTin Ω,(5) and the Hilbert space His defined as the completion of D(Ω) with respect to the norm kϕTkH:= kϕkL2(qT). In view of the unique continuation property to (5), the mapping ϕT7→ kϕTkHis a Hilbertian norm in D(Ω). Hence, we can certainly consider the completion of D(Ω) for this norm. The coercivity of the functional J⋆ 1in His a consequence of the so-called observability inequality kϕ(·,0)k2 L2(Ω) ≤CZZqT|ϕ(x, t)|2dx dt ∀ϕT∈ H,(6) where ϕsolves (5). This inequality holds for some constant C=C(ω, T) and, in turn, is a consequence of some appropriate global Carleman inequalities; see [20]. The minimization of J⋆ 1is numerically ill-posed, essentially because of the hugeness of the completed space H. The control of minimal square integrable norm highly oscillates near the final time T, property which is hard to capture numerically. We refer to [1, 24, 30, 33] where this phenomenon is highlighted under several perspectives. Moreover, at the level of the approximation, the minimization of J⋆ 1requires to find a finite dimensional and conformal approximation of Hsuch that the corresponding discrete adjoint solution satisfies (5), which is in general impossible for polynomial piecewise approximation. In practice, the trick initially described in [22], consists first to introduce a discrete and consistent approximation of (1) and then to minimize the corresponding discrete conjugate functional. However, this requires to get some uniform discrete observability inequalities; for the heat equation with non constant coefficient, this remains an open questions. This fact and the hugeness of Hhas raised many authors to relax the controllability problem: precisely, the constraint (2). We mention the references [5, 3, 35] and notably [2, 18, 25] for some numerical realizations. In [17] (see also [16] in a semi-linear case), a different - so-called primal approach - allowing more general results has been used and consists to solve directly optimality conditions : specifically, the following general extremal problem (initially introduced by Fursikov and Imanuvilov in [20]) is considered :      Minimize J(y, v) := 1 2ZZQT ρ2|y|2dx dt +1 2ZZqT ρ2 0|v|2dx dt Subject to (y, v)∈ C(y0, T). (7) The weights ρ=ρ(x, t) and ρ0=ρ0(x, t) are continuous, uniformly positive and are assumed to belong to L∞(QT−δ) for any δ > 0 (hence, they can blow up as t→T−). Under those conditions, the extremal problem (7) is well-posed (see [17]). Moreover, the explicit occurrence of the term yin the functional allow to solve directly the optimality conditions associated with (7): defining the Hilbert space Pas the completion of the linear space P0={q∈C∞(QT) : q= 0 on ΣT}with respect to the scalar product (p, q)P:= ZZQT ρ−2L⋆p L⋆q dx dt +ZZqT ρ−2 0p q dx dt, (8) 1 INTRODUCTION. THE NULL CONTROLLABILITY PROBLEM 4 the optimal pair (y, v) for Jis characterized as follows y=ρ−2L⋆pin QT, v =−ρ−2 0p1ωin QT(9) in term of an additional variable p∈Punique solution to the following variational equality : (p, q)P=ZΩ y0(x)q(x, 0) dx, ∀q∈P. (10) The well-posedeness of this formulation is ensured as soon as the weights ρ0, ρ are of Carleman type (in particular ρand ρ0blow up exponentially as t→T−); this specific behavior near Treinforces the null controllability requirement and prevents the control of any oscillations near the final time. The search of a control vin the manifold C(y0, T) is reduced to solve the (elliptic) variational formulation (10). In [17], the approximation of (10) is performed in the framework of the finite element theory through a discretization of the space-time domain QT. In practice, an approximation phof pis obtained in a direct way by inverting a symmetric positive definite matrix, in contrast with the iterative (and possibly divergent) methods used within dual methods. Moreover, a major advantage of this approach is that a conformal approximation, say Phof P, leads to the strong convergence of phtoward pin P, and consequently from (9), to a strong convergence in L2(qT) of vh:= −ρ−2 0ph1ωtoward v, a null control for (1). It is worth to mention that, for any h > 0, vh is not a priori an exact control for any finite dimensional system (which is not necessary at all in practice) but an approximation for the L2-norm of the control v. The variational formulation (10) derived from the optimality conditions (9) is obtained assuming that the weights ρand ρ0are both strictly positive in QTand qTrespectively. In particular, this approach does not apply for the control of minimal L2-norm, for which simply ρ:= 0 and ρ0:= 1. The main reason of the present work is to adapt this approach to cover the case ρ:= 0 and therefore obtain directly an approximation vhof the control of some minimal weighted L2-norm. To do so, we adapt the idea developed in [11] devoted to the wave equation. We also mention [32] where a different space-time variational approach (based on Least-squares principles) is used to approximate null controls for the heat equation. The paper is organized as follows. In Section 2, we associate to the dual problem (4) an equivalent mixed formulation which relies on the optimality conditions associated to the problem (7) with ρ= 0. In Section 2.1, we first address the penalization case and write the constraint L⋆ϕ= 0 as an equality in L2(QT). We then show the well-posedness of this mixed formulation, in particular we check the inf-sup condition (Theorem 2.1). The mixed formulation allows to approximate simultaneously the dual variable and the primal one, controlled solution of (1). Interestingly, we also derive an equivalent extremal problem in the primal variable yonly (see Prop 2.2, Section 2.1.2). In Section 2.2, we reproduce the analysis relaxing the condition L⋆ϕ= 0 in the weaker space L2(0, T, H−1(Ω)). Then, in Section 2.3, by using the Global Carleman estimate (34), we show that a well-posed mixed formulation is also available for the limit and singular case for which ε= 0 leading to Theorem 2.3. Section 3 is devoted to the numerical approximation of the mixed formulation (16) in the case ε > 0 (Section 3.1) and of the mixed formulation (32) in the case ε= 0 (section 3.2). Conformal approximations based on space-time finite elements are employed. In Section 3.3, we numerically check that the approximations used lead to discrete inf-sup properties, uniformly w.r.t. the discretization parameter h. Then the remaining of Section 3 is devoted to some experiments which emphasize the remarkable robuteness of the method. Section 4 concludes with some perspectives. 2 CONTROL OF MINIMAL WEIGHTED L2-NORM: MIXED REFORMULATIONS 5 2 Control of minimal weighted L2-norm: mixed reformulations In order to avoid the minimization of the conjugate functional J⋆with respect to the final state ϕT by an iterative process, we now present a direct way to approximate the control of minimal square integrable norm in the spirit of the primal approach recalled in the introduction and developed in [17]. We adapt the case of the wave equation studied in [11]. 2.1 The penalized case: Mixed formulation I Let ρ⋆∈R+ ⋆and let ρ0∈ R defined by R:= {w:w∈C(QT); w≥ρ⋆>0 in QT;w∈L∞(QT−δ)∀δ > 0}(11) so that in particular, the weight ρ0may blow up as t→T−. We first consider the approximate controllability case. For any ε > 0, the problem reads as follows:      Minimize Jε(y, v) := 1 2ZZqT ρ2 0|v|2dt +1 2εky(·, T)k2 L2(Ω) Subject to (y, v)∈ A(y0;T) (12) where A(y0;T) denotes the linear manifold A(y0;T) := {(y, v) : v∈L2(qT), y solves (1) }and where εdenotes a penalty parameter (see [2, 5, 18]). The corresponding conjugate and well-posed problem is given by      Minimize J⋆ ε(ϕT) := 1 2ZZqT ρ−2 0|ϕ(x, t)|2dx dt +ε 2kϕTk2 L2(Ω) + (y0, ϕ(·,0))L2(Ω) Subject to ϕT∈L2(Ω). (13) where ϕsolves (5). We recall that the penalized problem (12) is a consistent approximation of the original null controllability problem, in the sense that its unique solution converges to the solution of (7) with ρ= 0 as ε→0. We refer for instance to [18], Prop. 3.3 for a proof of the following result, consequence of the null controllability for the heat equation. Proposition 2.1 Let (yε, vε)be the solution of Problem (12) and let (y, v)be the solution of Problem (7) with ρ= 0. Then, one has yε→ystrongly in L2(QT), vε→vstrongly in L2(qT) as ε→0+. 2.1.1 Mixed formulation Since the variable ϕ, solution of (5), is completely and uniquely determined by the data ϕT, the main idea of the reformulation is to keep ϕas main variable. We introduce the linear space Φ0:= {ϕ∈C2(QT), ϕ = 0 on ΣT}. For any η > 0, we define the bilinear form (ϕ, ϕ)Φ0:= ZZqT ρ−2 0ϕϕ dx dt +ε(ϕ(·, T), ϕ(·, T))L2(Ω) +ηZZQT L⋆ϕ L⋆ϕ dx dt, ∀ϕ, ϕ ∈Φ0. 2 CONTROL OF MINIMAL WEIGHTED L2-NORM: MIXED REFORMULATIONS 6 From the unique continuation property for the heat equation, this bilinear form defines for any ε≥0 a scalar product. For any ε > 0, let Φεbe the completion of Φ0for this scalar product. We denote the norm over Φεby k·kΦεsuch that kϕk2 Φε:= kρ−1 0ϕk2 L2(qT)+εkϕ(·, T)k2 L2(Ω) +ηkL⋆ϕk2 L2(QT),∀ϕ∈Φε.(14) Finally, we defined the closed subset Wεof Φεby Wε={ϕ∈Φε:L⋆ϕ= 0 in L2(QT)} and we endow Wεwith the same norm than Φε. Then, we define the following extremal problem : min ϕ∈Wε ˆ J⋆ ε(ϕ) := 1 2ZZqT ρ−2 0|ϕ(x, t)|2dx dt +ε 2kϕ(·, T)k2 L2(Ω) + (y0, ϕ(·,0))L2(Ω).(15) Standard energy estimates for the heat equation imply that, for any ϕ∈Wε,ϕ(·,0) ∈L2(Ω) so that the functional ˆ J⋆ εis well-defined over Wε. Moreover, since for any ϕ∈Wε,ϕ(·, T) belongs to L2(Ω), Problem (15) is equivalent to the minimization problem (13). As announced, the main variable is now ϕsubmitted to the constraint equality (in L2(QT)) L⋆ϕ= 0. This constraint equality is addressed by introducing a Lagrangian multiplier. We consider the following mixed formulation : find (ϕε, λε)∈Φε×L2(QT) solution of (aε(ϕε,ϕ) + b(ϕ, λε) = l(ϕ),∀ϕ∈Φε b(ϕε,λ) = 0,∀λ∈L2(QT),(16) where aε: Φε×Φε→R, aε(ϕ, ϕ) := ZZqT ρ−2 0ϕϕ dx dt +ε(ϕ(·, T), ϕ(·, T))L2(Ω) b: Φε×L2(QT)→R, b(ϕ, λ) := −ZZQT L⋆ϕ λ dx dt l: Φε→R, l(ϕ) := −(y0, ϕ(·,0))L2(Ω). We have the following result : Theorem 2.1 (i) The mixed formulation (16) is well-posed. (ii) The unique solution (ϕε, λε)∈Φε×L2(QT)is the unique saddle-point of the Lagrangian Lε: Φε×L2(QT)→Rdefined by Lε(ϕ, λ) := 1 2aε(ϕ, ϕ) + b(ϕ, λ)−l(ϕ).(17) (iii) The optimal function ϕεis the minimizer of ˆ J⋆ εover Wεwhile the optimal multiplier λε∈ L2(QT)is the state of the heat equation (1) in the weak sense. Proof - We easily check that the bilinear form aεis continuous over Φε×Φε, symmetric and positive and that the bilinear form bεis continuous over Φε×L2(QT). Furthermore, for any fixed ε, the continuity of the linear form lover Φεcan be viewed from the energy estimate : kϕ(·,0)k2 L2(Ω) ≤CZZQT|L⋆ϕ|2dx dt +kϕ(·, T )k2 L2(Ω),∀ϕ∈Φε, for some C > 0, so that kϕ(·,0)k2 L2(Ω) ≤Cmax(η−1, ε−1)kϕk2 Φε. Therefore, the well-posedness of the mixed formulation is a consequence of the following two properties (see [4]): 2 CONTROL OF MINIMAL WEIGHTED L2-NORM: MIXED REFORMULATIONS 7 •aεis coercive on N(b), where N(b) denotes the kernel of b: N(b) := {ϕ∈Φε:b(ϕ, λ) = 0 for every λ∈L2(QT)}; •bsatisfies the usual “inf-sup” condition over Φε×L2(QT): there exists δ > 0 such that inf λ∈L2(QT)sup ϕ∈Φε b(ϕ, λ) kϕkΦεkλkL2(QT)≥δ. (18) From the definition of aε, the first point is clear : for all ϕ∈ N(b) = Wε,aε(ϕ, ϕ) = kϕk2 Φε. Let us check the inf-sup condition. For any fixed λ0∈L2(QT), we define the (unique) element ϕ0of L⋆ϕ0=−λ0in QT, ϕ0= 0 on ΣT;ϕ0(·, T) = 0 in Ω, so that ϕ0solves the backward heat equation with source term −λ0∈L2(QT), null Dirichlet boundary condition and zero initial state. Since −λ0∈L2(QT), then using energy estimates, there exists a constant CΩ,T >0 such that the solution ϕ0of the backward heat equation with source term λ0satisfies the inequality ZZqT ρ−2 0|ϕ0|2dx dt ≤ρ−2 ⋆ZZqT|ϕ0|2dx dt ≤ρ2 ⋆CΩ,T kλ0k2 L2(QT). Consequently, ϕ0∈Φε. In particular, we have b(ϕ0, λ0) = kλ0k2 L2(QT)and sup ϕ∈Φε b(ϕ, λ0) kϕkΦεkλ0kL2(QT)≥b(ϕ0, λ0) kϕ0kΦεkλ0kL2(QT) =kλ0k2 L2(QT) kρ−1 0ϕ0k2 L2(qT)+ηkλ0k2 L2(QT)1 2kλ0kL2(QT) . Combining the above two inequalities, we obtain sup ϕ0∈Φε b(ϕ0, λ0) kϕ0kΦεkλ0kL2(QT)≥1 pρ2 ⋆CΩ,T +η and, hence, (18) holds with δ=ρ2 ⋆CΩ,T +η−1/2. The point (ii) is due to the symmetry and to the positivity of the bilinear form aε. (iii) Concerning the third point, the equality b(ϕε,λ) = 0 for all λ∈L2(QT) implies that L⋆ϕε= 0 as an L2(QT) function, so that if (ϕε, λε)∈Φε×L2(QT) solves the mixed formulation, then ϕε∈Wε and Lε(ϕε, λε) = ˆ J⋆ ε(ϕε). Finally, the first equation of the mixed formulation reads as follows: ZZqT ρ−2 0ϕεϕ dx dt +ε(ϕε(·, T), ϕ(·, T)) −ZZQT L⋆ϕ(x, t)λε(x, t)dx dt =l(ϕ),∀ϕ∈Φε, or equivalently, since the control is given by vε:= ρ−2 0ϕε1ω, ZZqT vεϕ dx dt + (εϕε(·, T), ϕ(·, T)) −ZZQT L⋆ϕ(x, t)λε(x, t)dx dt =l(ϕ),∀ϕ∈Φε. But this means that λε∈L2(QT) is solution of the heat equation in the transposition sense. Since y0∈L2(Ω) and vε∈L2(qT), λεmust coincide with the unique weak solution to (1) (yε=λε) such that λε(·, T) = −εϕε(·, T). As a conclusion, the optimal pair (yε, vε) to (12) is characterized in term of the adjoint variable ϕεsolution of (16) by vε=ρ−2 0ϕε1ωand yε(·, T) = −εϕε(·, T). ✷ Theorem 2.1 reduces the search of the approximated control to the resolution of the mixed formulation (16), or equivalently the search of the saddle point for Lε. In general, it is convenient 2 CONTROL OF MINIMAL WEIGHTED L2-NORM: MIXED REFORMULATIONS 8 to “augment” the Lagrangian (see [19]), and consider instead the Lagrangian Lε,r defined for any r > 0 by        Lε,r(ϕ, λ) := 1 2aε,r(ϕ, ϕ) + b(ϕ, λ)−l(ϕ), aε,r(ϕ, ϕ) := aε(ϕ, ϕ) + rZZQT|L⋆ϕ|2dx dt. Since aε(ϕ, ϕ) = aε,r(ϕ, ϕ) on Wεand since the function ϕεsuch that (ϕε, λε) is the saddle point of Lεverifies ϕε∈Wε, the lagrangian Lεand Lε,r share the same saddle-point. 2.1.2 Dual problem of the extremal problem (15) The mixed formulation allows to solve simultaneously the dual variable ϕε, argument of the conjugate functional (15), and the Lagrange multiplier λε. Since λεturns out to be the (approximate) controlled state of (1), we may qualify λεas the primal variable of the problem. We derive in this section the corresponding extremal problem involving only that variable λε. For any r > 0, let us define the linear operator Aε,r from L2(QT) into L2(QT) by Aε,rλ:= L⋆ϕ, ∀λ∈L2(QT) where ϕ∈Φεis the unique solution to aε,r(ϕ, ϕ) = −b(ϕ, λ),∀ϕ∈Φε.(19) Note that the assumption r > 0 is necessary here in order to guarantee the well-posedness of (19). Precisely, for any r > 0, the form aε,r defines a norm equivalent to the norm on Φε(see (14)). We have the following crucial lemma : Lemma 2.1 For any r > 0, the operator Aε,r is a strongly elliptic, symmetric isomorphism from L2(QT)into L2(QT). Proof- From the definition of aε,r, we easily get that kAε,rλkL2(QT)≤r−1kλkL2(QT)and the continuity of Aε,r. Next, consider any λ′∈L2(QT) and denote by ϕ′the corresponding unique solution of (19) so that Aε,rλ′:= L⋆ϕ′. Relation (19) with ϕ=ϕ′then implies that ZZQT (Aε,rλ′)λ dx dt =aε,r(ϕ, ϕ′) (20) and therefore the symmetry and positivity of Aε,r. The last relation with λ′=λimplies that Aε,r is also positive definite. Finally, let us check the strong ellipticity of Aε,r, equivalently that the bilinear functional (λ, λ′)7→ ZZQT (Aε,rλ)λ′dx dt is L2(QT)-elliptic. Thus we want to show that ZZQT (Aε,rλ)λ dx dt ≥Ckλk2 L2(QT),∀λ∈L2(QT) (21) for some positive constant C. Suppose that (21) does not hold; there exists then a sequence {λn}n≥0of L2(QT) such that kλnkL2(QT)= 1,∀n≥0,lim n→∞ ZZQT (Aε,rλn)λndx dt = 0. 2 CONTROL OF MINIMAL WEIGHTED L2-NORM: MIXED REFORMULATIONS 9 Let us denote by ϕnthe solution of (19) corresponding to λn. From (20), we then obtain that lim n→∞ kL⋆ϕnkL2(QT)= 0,lim n→∞ kρ−1 0ϕnkL2(qT)= 0,lim n→∞ kϕn(·, T)kL2(Ω) = 0.(22) From (19) with λ=λnand ϕ=ϕn, we have ZZqT ρ−2 0ϕnϕ dx dt +εZ1 0 ϕn(·, T)ϕ(·, T)dx +ZZQT (rL⋆ϕn−λn)L⋆ϕ dx dt = 0,∀ϕ∈Φε.(23) We define the sequence {ϕn}n≥0as follows : L⋆ϕn=r L⋆ϕn−λnin QT, ϕn= 0 on ΣT;ϕ(·, T) = 0 in Ω, so that, for all n≥0, ϕnis the solution of the backward heat equation with zero initial datum and source term r L⋆ϕn−λnin L2(QT). Using again energy type estimates, we get kρ−1 0ϕnkL2(qT)≤ρ−1 ⋆kϕnkL2(qT)≤ρ−1 ⋆CΩ,T krL⋆ϕn−λnkL2(QT), so that ϕn∈Φε. Then, using (23) with ϕ=ϕn, we get krL⋆ϕn−λnkL2(QT)≤ρ−1 ⋆CΩ,T kρ−1 0ϕnkL2(qT). Then, from (22), we conclude that limn→+∞kλnkL2(QT)= 0 leading to a contradiction and to the strong ellipticity of the operator Aε,r.✷ The introduction of the operator Aε,r is motivated by the following proposition: Proposition 2.2 For any r > 0, let ϕ0∈Φεbe the unique solution of aε,r(ϕ0,ϕ) = l(ϕ),∀ϕ∈Φε and let J⋆⋆ ε,r :L2(QT)→L2(QT)be the functional defined by J⋆⋆ ε,r(λ) := 1 2ZZQT (Aε,rλ)λ dx dt −b(ϕ0, λ). The following equality holds : sup λ∈L2(QT) inf ϕ∈ΦεLε,r(ϕ, λ) = −inf λ∈L2(QT)J⋆⋆ ε,r(λ) + Lε,r(ϕ0,0). Proof- For any λ∈L2(QT), let us denote by ϕλ∈Φεthe minimizer of ϕ7→ Lε,r(ϕ, λ); ϕλ satisfies the equation aε,r(ϕλ,ϕ) + b(ϕ, λ) = l(ϕ),∀ϕ∈Φε and can be decomposed as follows : ϕλ=ψλ+ϕ0where ψλ∈Φεsolves aε,r(ψλ,ϕ) + b(ϕ, λ) = 0,∀ϕ∈Φε. We then have inf ϕ∈ΦεLε,r(ϕ, λ) = Lε,r(ϕλ, λ) = Lε,r(ψλ+ϕ0, λ) =1 2aε,r(ψλ+ϕ0, ψλ+ϕ0) + b(ψλ+ϕ0, λ)−l(ψλ+ϕ0) := X1+X2+X3 with      X1:= 1 2aε,r(ψλ, ψλ) + b(ψλ, λ) + b(ϕ0, λ) X2:= aε,r(ϕ0, ψλ)−l(ψλ), X3:= 1 2aε,r(ϕ0, ϕ0)−l(ϕ0). 3 NUMERICAL APPROXIMATION AND EXPERIMENTS 16 formulation (29). Therefore, this weight may be seen as a parameter to improve some specific properties of the mixed formulation, specifically at the numerical level. Precisely, in the limit case ε= 0, we recall that the trace ϕ|t=Tof the solution does not belong to L2(Ω) but to a much larger and singular space. Very likely, a similar behavior occurs for the function L⋆ϕnear Ω×{T}so that the constraint L⋆ϕ= 0 in L2(QT) introduced in Section 2.1 is too “strong” and must be replaced at the limit in εby the relaxed one ρ−1L⋆ϕ= 0 in L2(QT) with ρ−1“small” near Ω×{T}. Remark that this is actually the effect and the role of the Carleman type weights ρcdefined by (33) and initially introduced in [20]. As in Section 2.1, it is convenient to “augment” the Lagrangian and consider instead the Lagrangian Lrdefined for any r > 0 by        Lr(ϕ, λ) := 1 2˜ar(ϕ, ϕ) + ˜ b(ϕ, λ)−˜ l(ϕ), ˜ar(ϕ, ϕ) := ˜a(ϕ, ϕ) + rZZQT|ρ−1L⋆ϕ|2dx dt. Finally, similarly to Lemma 2.1 and Proposition 2.2, we have the following result. Let ρ0∈ R and ρ∈ R∩L∞(QT) Proposition 2.4 For any r > 0, let ρ0∈ R and ρ∈ R∩L∞(QT)verifying (35). Let us define the linear operator Arfrom L2(QT)into L2(QT)by Arλ:= ρ−1L⋆ϕ, ∀λ∈L2(QT), where ϕ∈e Φρ0,ρ is the unique solution to ar(ϕ, ϕ) = −b(ϕ, λ),∀ϕ∈e Φρ0,ρ. Aris a strongly elliptic, symmetric isomorphism from L2(QT)into L2(QT). Let ˆ J⋆⋆ r:L2(QT)→ L2(QT)be the functional defined by ˆ J⋆⋆ r(λ) := 1 2ZZQT (Arλ)λ dx dt −˜ b(ϕ0, λ). where ϕ0∈e Φρ0,ρ is the unique solution of ˜ar(ϕ0,ϕ) = ˜ l(ϕ),∀ϕ∈e Φρ0,ρ. The following equality holds : sup λ∈L2(QT) inf ϕ∈ e Φρ0,ρ Lr(ϕ, λ) = −inf λ∈L2(QT) ˆ J⋆⋆ r(λ) + Lr(ϕ0,0). 3 Numerical approximation and experiments 3.1 Discretization of the mixed formulation (16) We now turn to the discretization of the mixed formulation (16) assuming r > 0. Let then Φε,h and Mε,h be two finite dimensional spaces parametrized by the variable hsuch that, for any ε > 0, Φε,h ⊂Φε, Mε,h ⊂L2(QT),∀h > 0. Then, we can introduce the following approximated problems : find (ϕh, λh)∈Φε,h ×Mε,h solution of (aε,r(ϕh,ϕh) + b(ϕh, λh) = l(ϕh),∀ϕh∈Φε,h b(ϕh,λh) = 0,∀λh∈Mε,h.(40) 3 NUMERICAL APPROXIMATION AND EXPERIMENTS 17 The well-posedness of this mixed formulation is again a consequence of two properties : the coercivity of the bilinear form aε,r on the subset Nh(b) = {ϕh∈Φε,h;b(ϕh, λh) = 0 ∀λh∈Mε,h}. Actually, from the relation aε,r(ϕ, ϕ)≥Cr,ηkϕk2 Φε,∀ϕ∈Φε, where Cr,η = min{1, r/η}, the form aε,r is coercive on the full space Φε, and so a fortiori on Nh(b)⊂Φε,h ⊂Φε. The second property is a discrete inf-sup condition : there exists δh>0 such that inf λh∈Mε,h sup ϕh∈Φε,h b(ϕh, λh) kϕhkΦε,h kλhkMε,h ≥δh.(41) For any fixed h, the spaces Mε,h and Φε,h are of finite dimension so that the infimum and supremum in (41) are reached: moreover, from the property of the bilinear form aε,r, it is standard to prove that δhis strictly positive (see Section 3.3). Consequently, for any fixed h > 0, there exists a unique couple (ϕh, λh) solution of (40). On the other hand, the property infhδh>0 is in general difficult to prove and depends strongly on the choice made for the approximated spaces Mε,h and Φε,h. We shall analyze numerically this property in Section 3.3. Remark 2 For r= 0, the discrete mixed formulation (40) is not well-posed over Φε,h ×Mε,h because the bilinear form aε,r=0 is not coercive over the discrete kernel of b: the equality b(λh, ϕh) = 0for all λh∈Mε,h does not imply that L⋆ϕhvanishes. Therefore, the term rkL⋆ϕhk2 L2(QT)may be understood as a numerical stabilization term: for any h > 0, it ensures the uniform coercivity of the form aε,r (and so the well-posedness) and vanishes at the limit in h. We also emphasize that this term is not a regularization term as it does not add any regularity to the solution ϕh. As in [10], the finite dimensional and conformal space Φε,h must be chosen such that L⋆ϕh belongs to L2(QT) for any ϕh∈Φε,h. This is guaranteed as soon as ϕhpossesses second-order derivatives in L2 loc(QT). Any conformal approximation based on standard triangulation of QT achieves this sufficient property as soon as it is generated by spaces of functions continuously differentiable with respect to the variable xand spaces of continuous functions with respect to the variable t. We introduce a triangulation Thsuch that QT=∪K∈ThKand we assume that {Th}h>0is a regular family. Then, we introduce the space Φε,h as follows : Φε,h ={ϕh∈C1(QT) : ϕh|K∈P(K)∀K∈ Th, ϕh= 0 on ΣT} where P(K) denotes an appropriate space of polynomial functions in xand t. In this work, we consider for P(K) the so-called Bogner-Fox-Schmit (BFS for short) C1-element defined for rectangles. In the one dimensional setting considered in the sequel, it involves 16 degrees of freedom, namely the values of ϕh, ϕh,x, ϕh,t, ϕh,xt on the four vertices of each rectangle K. Therefore P(K) = P3,x ⊗P3,t where Pr,ξ is by definition the space of polynomial functions of order rin the variable ξ. We refer to [9] page 76. We also define the finite dimensional space Mε,h ={λh∈C0(QT) : λh|K∈Q(K)∀K∈ Th}, where Q(K) denotes the space of affine functions both in xand ton the element K. Again, in the one dimensional setting, for rectangle, we simply have Q(K) = P1,x ⊗P1,t. We also mention that the approximation is conformal : for any h > 0, we have Φε,h ⊂Φεand Mε,h ⊂L2(QT). 3 NUMERICAL APPROXIMATION AND EXPERIMENTS 18 Let nh= dim Φε,h, mh= dim Mε,h and let the real matrices Aε,r,h ∈Rnh,nh,Bh∈Rmh,nh, Jh∈Rmh,mhand Lh∈Rnhbe defined by                aε,r(ϕh,ϕh) =< Aε,r,h{ϕh},{ϕh}>Rnh,Rnh,∀ϕh,ϕh∈Φε,h, b(ϕh, λh) =< Bh{ϕh},{λh}>Rmh,Rmh,∀ϕh∈Φε,h,∀λh∈Mε,h, ZZQT λhλhdx dt =< Jh{λh},{λh}>Rmh,Rmh,∀λh,λh∈Mε,h, l(ϕh) =< Lh,{ϕh}>, ∀ϕh∈Φε,h where {ϕh} ∈ Rnhdenotes the vector associated to ϕhand <·,·>Rnh,Rnhthe usual scalar product over Rnh. With these notations, Problem (40) reads as follows : find {ϕh} ∈ Rnhand {λh} ∈ Rmh such that Aε,r,h BT h Bh0Rnh+mh,nh+mh{ϕh} {λh}Rnh+mh =Lh 0Rnh+mh .(42) The matrix Aε,r,h as well as the mass matrix Jhare symmetric and positive definite for any h > 0 and any r > 0. On the other hand, the matrix of order mh+nhin (42) is symmetric but not positive definite. We use exact integration methods developed in [14] for the evaluation of the coefficients of the matrices. The system (42) is solved using the direct LU decomposition method. Let us also mention that for r= 0, although the formulation (16) is well-posed, numerically, the corresponding matrix Aε,0,h is not invertible in agreement with Remark 2. In the sequel, we shall consider strictly positive values for r. Once an approximation ϕhis obtained, an approximation vε,h of the control vεis given by vε,h =ρ−2 0ϕε,h 1ω. The corresponding controlled state yε,h may be obtained by solving (1) with standard forward approximation (we refer to [10], Section 4 where this is detailed). Here, since the controlled state is directly given by the multiplier λ, we simply use λhas an approximation of yand we do not report here the computation of yh. In the sequel, we only report numerical experiments in the one dimensional setting. We use uniform rectangular meshes. Each element is a rectangle of lengths ∆xand ∆t; ∆x > 0 and ∆t > 0 denote as usual the discretization parameters in space and time, respectively. We note h:= max{diam(K), K ∈ Th} where diam(K) denotes the diameter of K. 3.2 Normalization and discretization of the mixed formulation (32) The same approximation may be used for the mixed formulation (32). In particular, we easily check that the finite dimensional spaces Mε,h and Φε,h (which actually do not depend on ε) are conformal approximation of L2(QT) and e Φρ0,ρ respectively. However, in the limit case ε= 0, a normalization of the variable ϕ, which is singular and takes arbitrarily large amplitude in the neighborhood of Ω ×{T}is very convenient and suitable in practice. Following [17], we introduce the variable ψ:= ρ−1 0ϕ∈ρ−1 0e Φρ0,ρ and replace the mixed formulation (32) by the equivalent one: find (ψ, λ)∈ρ−1 0e Φρ0,ρ ×L2(QT) solution of (ˆa(ψ, ψ) + ˆ b(ψ, λ) = ˆ l(ψ),∀ψ∈ρ−1 0e Φρ0,ρ ˆ b(ψ, λ) = 0,∀λ∈L2(QT), (43) 3 NUMERICAL APPROXIMATION AND EXPERIMENTS 19 where ˆa:ρ−1 0e Φρ0,ρ ×ρ−1 0e Φρ0,ρ →R,ˆa(ψ, ψ) = ZZqT ψψ dx dt ˆ b:ρ−1 0e Φρ0,ρ ×L2(QT)→R,ˆ b(ψ, λ) = −ZZQT ρ−1L⋆(ρ0ψ)λ dx dt ˆ l:ρ−1 0e Φρ0,ρ →R,ˆ l(ϕ) = −(y0, ρ0(·,0)ψ(·,0))L2(Ω). Well-posedness of this formulation is the consequence of Theorem 2.3. Moreover, the optimal controlled state is still given by ρ−1λwhile the optimal control is expressed in term of the variable ψas v=ρ−1 0ψ1ω. The corresponding discretization approximation (augmented with the term rkρ−1L⋆(ρ0ψ)kL2(QT)) reads as follows: find (ψh, λh)∈Φh×Mhsolution of (ˆar(ψh,ψh) + ˆ b(ψh, λh) = ˆ l(ϕh),∀ψh∈Φh ˆ b(ψh,λh) = 0,∀λh∈Mh. (44) with ar(ψh,ψh) :=a(ψh, ψh) + r(ρ−1L⋆(ρ0ψh), ρ−1L⋆(ρ0ψh))L2(QT) =(ψh,ψh)L2(qT)+r(ρ−1L⋆(ρ0ψh), ρ−1L⋆(ρ0ψh))L2(QT). for any r > 0. Remark 3 When the weights ρ0and ρare chosen in such a way that they are compensated each other in the term ρ−1L⋆(ρ0ψ), the change of variable has the effect to reduced the amplitude (with respect to the time variable) of the coefficients in the integrals of ˆarand ˆ b, and therefore, at the discrete level, to improve significantly the condition number of square matrix ˆ Ar,h so that ˆar(ψh,ψh) =<ˆ Ar,h{ψh},{ψh}>Rnh,Rnh. In this respect, the change of variable, can be seen as a preconditioner for the mixed formulation (32). Similarly to (42), we note the matrix form of (44) as follows : ˆ Ar,h ˆ BT h ˆ Bh0!Rnh+mh,nh+mh{ψh} {λh}Rnh+mh =ˆ Lh 0Rnh+mh ,(45) where ˆ Bhis the matrix so that ˆ b(ψh, λh) =<ˆ Bh{ψh},{λh}>Rmh,Rmhand ˆ Lhis the matrix so that ˆ l(ψh) =<ˆ Lh,{ψh}>. 3.3 The discrete inf-sup test Before to give and discuss some numerical experiments, we first test numerically the discrete infsup condition (41). Taking η=r > 0 so that aε,r(ϕ, ϕ) = (ϕ, ϕ)Φεexactly for all ϕ, ϕ∈Φε, it is readily seen (see for instance [8]) that the discrete inf-sup constant satisfies δε,r,h = inf√δ:BhA−1 ε,r,hBT h{λh}=δ Jh{λh},∀{λh} ∈ Rmh\{0}.(46) The matrix BhA−1 ε,r,hBT henjoys the same properties than the matrix Aε,r,h: it is symmetric and positive definite so that the scalar δε,h defined in term of the (generalized) eigenvalue problem (46) is strictly positive. This eigenvalue problem is solved using the power iteration algorithm (assuming that the lowest eigenvalue is simple): for any {v0 h} ∈ Rnhsuch that k{v0 h}k2= 1, compute for any n≥0, {ϕn h} ∈ Rnh,{λn h} ∈ Rmhand {vn+1 h} ∈ Rmhiteratively as follows : (Aε,r,h{ϕn h}+BT h{λn h}= 0 Bh{ϕn h}=−Jh{vn h},{vn+1 h}={λn h} k{λn h}k2 . 3 NUMERICAL APPROXIMATION AND EXPERIMENTS 20 The scalar δε,r,h defined by (46) is then given by : δε,r,h = limn→∞(k{λn h}k2)−1/2. We now give some numerical values of δε,r,h with respect to hfor the C1-finite element introduced in Section 3.1. We consider the one dimensional case for which Ω = (0,1) and take for simplicity c:= 1/10 and d:= 0. Values of the diffusion cand of the potential ddo not affect qualitatively the results. In the spirit of the previous work [17], we consider the following choice for the weight ρ0∈ R: ρ0(x, t) := (T−t)3/2expK1 (T−t),(x, t)∈QT, K1:= 3 4(47) so that ρ0blows exponentially as t→T−. This allows a smooth behavior of the corresponding control v:= ρ−2 0ϕ1ω. Let us insist however that the mixed formulation is well-posed for any weight ρ0∈ R, in particular ρ0:= 1 (leading to the control of minimal L2-norm and for which we refer to [33]). ρ0is independent of the variable xfor simplicity. We consider the following data ω= (0.2,0.5), T= 1/2, and Ω = (0,1). Tables 1, 2 and 3 provides the values of δε,r,h with respect to hand εfor r= 10−2,1 and r= 102, respectively. In view of the definition, we check that δε,r,h increases as r→0 and ε→0. We also observe, that for rlarge enough (see Tables 2 and 3), the value of the inf-sup constant is almost constant with respect to εand behaves like δε,r,h ≈Cε,r,h ×r−1/2(48) for some constant Cε,r,h ∈(0,1). More importantly, we observe that for any rand ε, the value of δε,r,h is bounded by below uniformly with respect to the discretization parameter h. The same behavior is observed for other discretizations such that ∆t6= ∆x, other supports ωand other choices for the weight ρ0(in particular ρ0:= 1). Consequently, we may conclude that the finite approximation we have used do ”pass” the discrete inf-sup test. It is interesting to note that this is in contrast with the situation for the wave equation for which the parameter rhave to be adjusted carefully with respect to h; we refer to [11]. Moreover, as it is usual in mixed finite element theory, such a property together with the uniform coercivity of form aε,r then implies the convergence of the approximation sequence (ϕh, λh) solution of (40). h7.07 ×10−23.53 ×10−21.76 ×10−28.83 ×10−3 ε= 10−28.358 8.373 8.381 8.386 ε= 10−49.183 9.213 9.229 9.237 ε= 10−89.263 9.318 9.354 9.383 Table 1: δε,r,h w.r.t. εand h;r= 10−2; Ω = (0,1), ω= (0.2,0.5), T= 1/2. h7.07 ×10−23.53 ×10−21.76 ×10−28.83 ×10−3 ε= 10−29.933 ×10−19.938 ×10−19.940 ×10−19.941 ×10−1 ε= 10−49.933 ×10−19.938 ×10−19.941 ×10−19.942 ×10−1 ε= 10−89.933 ×10−19.938 ×10−19.941 ×10−19.942 ×10−1 Table 2: δε,r,h w.r.t. εand h;r= 1.; Ω = (0,1), ω= (0.2,0.5), T= 1/2. Similarly, Table 4 displays the discrete inf-sup constant corresponding to the limit case of the 3 NUMERICAL APPROXIMATION AND EXPERIMENTS 21 h7.07 ×10−23.53 ×10−21.76 ×10−28.83 ×10−3 ε= 10−29.933 ×10−29.939 ×10−29.940 ×10−29.941 ×10−2 ε= 10−49.933 ×10−29.939 ×10−29.941 ×10−29.942 ×10−2 ε= 10−89.933 ×10−29.939 ×10−29.941 ×10−29.942 ×10−2 Table 3: δε,r,h w.r.t. εand h;r= 102; Ω = (0,1), ω= (0.2,0.5), T= 1/2. mixed formulation (43): inf λh∈ c Mh sup ψh∈ b Φh ˆ b(ψh, λh) kλhkL2(QT)kψhkρ−1 0 e Φρ0,ρ = inf λh∈ c Mh sup ψh∈ b Φh ZZqT λhρ−1L⋆(ρ0ψh)dx dt kλhkL2(QT)(kψhk2 L2(qT)+rkρ−1L⋆(ρ0ψh)k2 L2(QT))1/2. We take here a weight ρindependent of the variable xgiven by ρ(x, t) := expK1 (T−t),(x, t)∈QT, K1:= 3 4.(49) Again, for the limit case, the value given in the Table suggest a similar behavior observed for ε > 0: the constant is uniformly bounded by below with respect to the parameter hand behaves like r−1/2for rlarge enough (up to 1). Remark that, due to the introduction of the weight ρ6= 1, the inf-sup constants given by Table 4 are not the limit (as ε→0) of the previous Tables. h7.07 ×10−23.53 ×10−21.76 ×10−28.83 ×10−34.41 ×10−3 r= 1026.9×10−26.91 ×10−27.06 ×10−28.08 ×10−29.52 ×10−2 r= 1 6.89 ×10−16.91 ×10−16.96 ×10−17.94 ×10−18.66 ×10−1 r= 10−21.944 1.922 1.845 1.775 1.731 Table 4: ε= 0; δr,h w.r.t. rand h; Ω = (0,1), ω= (0.2,0.5), T= 1/2. 3.4 Numerical experiments for the mixed formulation (16) We report in this section experiments for the mixed formulation (16) and for simplicity we consider only the one dimensional case: Ω = (0,1) and T= 1/2. Let us first remark that in general explicit solutions (ϕε, λε) of (16) are not available. However, when the coefficient cand dare constant, we may obtain a semi-explicit representation (using Fourier decomposition) of the minimizer ϕε,T of the conjugate functional J⋆ ε(see (13)), and consequently of the corresponding adjoint variable ϕε, the control of weighted minimal square integrable norm v=ρ−2 0ϕε1ωand finally the controlled state yεsolution of (1). In practice, the obtention of the Fourier representation amounts to solve a symmetric linear system. We refer to the Appendix for the details. Such representation allows to evaluate precisely the distance of the exact solution (ϕε, λε) from the approximation (ϕh, λh) with respect to hand validate the convergence of the approximation with respect to h. As for the initial data, we first simply take the first mode of the Laplacian, that is, y0(x) = sin(πx), x∈(0,1). In view of the regularization property of the heat equation, the regularity of the initial data has a very restricted effect on the optimal control and the robustness of the method. We take c(x) := 10−1,d(x, t) := 0 and recall that in the uncontrolled case (ω=∅), these data leads to ky(·, T)kL2(0,1) =p1/2e−π2cT ≈4.31 ×10−1. Finally, we take ω= (0.2,0.5). 3 NUMERICAL APPROXIMATION AND EXPERIMENTS 22 For r= 1, Tables 5, 6 and 7 report some norms with respect to hfor ε= 10−2, 10−4and ε= 10−8respectively. The cases r= 102and r= 10−2are reported in the Appendix, in Tables 16, 17, 18 and 19, 20, 21 respectively. In the Tables, ϕεand yεdenotes the unique solution of (16) given by (51) and (53). In the Tables, κε,h denotes the condition number associated to (42), independent of the initial data y01. We first check that the L2-norm kλε,h(·, T)kL2(0,1) of the final state is of the order of √εand that the condition number κε,h behave polynomially with respect to h; on the other hand, we observe a low variation of κε,h with respect to ε;κε,h ≈O(h5.9) for ε= 10−2and κε,h ≈O(h7.3) for ε= 10−8. Then, we check the convergence as htends to zero of the approximations (vε,h, λε,h) toward the optimal pair (vε, yε) in L2(qT)×L2(QT) for any values of εand r. More precisely, for large enough value of ε(here ε= 10−2), we observe a quasi linear rate of convergence for both kρ0(vε−vε,h)kL2(qT) kρ0vεkL2(qT)and kyε−λε,hkL2(QT) kyεkL2(QT)with respect to h, independent of the value of the parameter r. We refer to Figure 1. For small values of ε, we observe a reduced convergence both for the control and the state (see Figure 2 for ε= 10−4and Figure 3 for ε= 10−8). We recall that as εtends to zero, the space Φεdegenerates into a much larger space and ϕεhighly oscillates near T. Remark also that for ε= 10−8, the constraint L⋆ϕε= 0 as an L2(QT) function is badly represented: this is due to the loss of regularity on the variable ϕε(in the neighborhood of T) as ε→0+. This does not prevent the convergence of the variable ϕε,h for the norm Φε, in particular the control vε,h =ρ−2ϕε,h 1ω, and of the variable λε,h. We will come back to this situation in detail in the section devoted to the limit case ε= 0. Moreover, for small value of ε, the parameter rdoes have an influence; precisely, a low value of r(here r= 10−2) leads to better relative errors : this is in agreement with the behavior of the inf-sup constant δε,r,h which increases with r−1. h1.41 ×10−17.07 ×10−23.53 ×10−21.76 ×10−28.83 ×10−3 mh+nh330 1 155 4 305 16 605 65 205 kL⋆ϕε,hkL2(QT)1.32 ×10−13.75 ×10−29.66 ×10−32.42 ×10−37.82 ×10−4 kρ0(vε−vε,h)kL2(qT) kρ0vεkL2(qT)1.10 ×10−16.21 ×10−23.29 ×10−21.68 ×10−28.57 ×10−3 kyε−λε,hkL2(QT) kyεkL2(QT)5.13 ×10−22.84 ×10−21.48 ×10−27.60 ×10−33.89 ×10−3 kλε,h(·, T )kL2(0,1) 1.54 ×10−11.61 ×10−11.65 ×10−11.67 ×10−11.68 ×10−1 κε,h 1.52 ×1091.10 ×1011 6.80 ×1012 3.83 ×1014 1.96 ×1016 Table 5: Mixed formulation (16) - r= 1 and ε= 10−2with ω= (0.2,0.5). h1.41 ×10−17.07 ×10−23.53 ×10−21.76 ×10−28.83 ×10−3 kL⋆ϕε,hkL2(QT)1.383 1.471 9.05 ×10−12.56 ×10−16.54 ×10−2 kρ0(vε−vε,h)kL2(qT) kρ0vεkL2(qT)6.72 ×10−13.22 ×10−11.15 ×10−15.49 ×10−22.74 ×10−2 kyε−λε,hkL2(QT) kyεkL2(QT)2.73 ×10−11.86 ×10−15.89 ×10−22.51 ×10−21.26 ×10−2 kλε,h(·, T )kL2(0,1) 8.50 ×10−25.74 ×10−23.39 ×10−23.11 ×10−23.13 ×10−2 κε,h 3.02 ×1093.91 ×1011 3.86 ×1013 3.25 ×1015 2.46 ×1017 Table 6: Mixed formulation (16) - r= 1 and ε= 10−4with ω= (0.2,0.5). 1The condition number κ(Mh) of any square matrix Mhis defined by κ(Mh) = |||Mh|||2|||M−1 h|||2where the norm |||Mh|||2stands for the largest singular value of Mh. 3 NUMERICAL APPROXIMATION AND EXPERIMENTS 23 h1.41 ×10−17.07 ×10−23.53 ×10−21.76 ×10−28.83 ×10−3 kL⋆ϕε,hkL2(QT)1.48 2.03 2.50 2.52 2.61 kρ0(vε−vε,h)kL2(qT) kρ0vεkL2(qT)1.44 1.01 7.92 ×10−16.65 ×10−14.89 ×10−1 kyε−λε,hkL2(QT) kyεkL2(QT)8.42 ×10−18.27 ×10−15.73 ×10−14.35 ×10−12.89 ×10−1 kλε,h(·, T )kL2(0,1) 8.63 ×10−26.65 ×10−22.39 ×10−21.23 ×10−24.43 ×10−3 κε,h 3.12 ×1094.30 ×1011 6.05 ×1013 1.13 ×1016 1.90 ×1018 Table 7: Mixed formulation (16) - r= 1 and ε= 10−8with ω= (0.2,0.5). 10ï210ï1 10ï2 10ï1 h 10ï210ï1 10ï3 10ï2 10ï1 h Figure 1: ω= (0.2,0.5); y0(x) = sin(πx): ε= 10−2. ; kρ0(vε−vε,h)kL2(qT) kρ0vεkL2(qT)(Left) and kyε−λε,hkL2(QT) kyεkL2(QT) (Right) vs. hfor r= 102(◦), r= 1.(⋆) and r= 10−2(). 10ï210ï1 10ï2 10ï1 100 h 10ï210ï1 10ï2 10ï1 100 h Figure 2: ω= (0.2,0.5); y0(x) = sin(πx): ε= 10−4. ; kρ0(vε−vε,h)kL2(qT) kρ0vεkL2(qT)(Left) and kyε−λε,hkL2(QT) kyεkL2(QT) (Right) vs. hfor r= 102(◦), r= 1.(⋆) and r= 10−2(). 3 NUMERICAL APPROXIMATION AND EXPERIMENTS 24 10ï210ï1 10ï0.8 10ï0.6 10ï0.4 10ï0.2 100 h 10ï210ï1 10ï1 100 h Figure 3: ω= (0.2,0.5); y0(x) = sin(πx): ε= 10−8. ; kρ0(vε−vε,h)kL2(qT) kρ0vεkL2(qT)(Left) and kyε−λε,hkL2(QT) kyεkL2(QT) (Right) vs. hfor r= 102(◦), r= 1.(⋆) and r= 10−2(). Remarkably, we highlight that the variational approach developed here allows, for any ε, a direct and robust approximation of one control for the heat equation. As discussed at length in [18, 33], the minimization of the conjugate functional J⋆ εusing conjugate gradient algorithm requires a number of iterates for small ε(typically ε= 10−8and ω= (0.2,0.5)) and diverge for small values of h. 3.5 Conjugate gradient for J⋆⋆ ε,r We illustrate here the Section 2.1.2 and minimize the functional J⋆⋆ ε,r :L2(QT)→Rdefined in Proposition 2.2 with respect to the variable λε. From the ellipticity of the operator Aε,r, we use a conjugate gradient method which in the present case reads as follows : (i) Initialization Let λ0 ε∈L2(QT) be a given function; Solve ϕ0 ε∈Φε aε,r(ϕ0 ε, ϕ) + bε(ϕ, λ0 ε) = lε(ϕ),∀ϕ∈Φε and set g0 ε=L⋆ϕ0 εand set w0 ε=g0 ε. For n≥0, assuming that λn ε, gn εand wn εare known with gn ε6= 0 and wn ε6= 0, compute λn+1 ε, gn+1 εand wn+1 εas follows (ii) Steepest descent Compute ϕn∈Φεsolution to aε,r(ϕn ε, ϕ) = −bε(ϕ, wn ε),∀ϕ∈Φε and wn ε=L⋆ϕn εand then compute ρn=kgn εk2 L2(QT)/(wn ε, wn ε)L2(QT). and set λn+1 ε=λn ε−ρnwn ε. 3 NUMERICAL APPROXIMATION AND EXPERIMENTS 25 (iii) Testing the convergence and construction of the new descent direction Update gn εby gn+1 ε=gn ε−ρnwn ε. If kgn+1 εkL2(QT)/kg0 εkL2(QT)≤γ, take λε=λn+1 ε. Else, compute γn=kgn+1 εk2 L2(QT)/kgn εk2 L2(QT) and update wn εvia wn+1 ε=gn+1 ε+γnwn ε. Do n=n+ 1 and return to step (ii). As mentioned in [21] where this approach is discussed at length for Stokes and Navier-Stokes systems, this algorithm can be viewed as a sophisticated version of Uzawa type algorithm to solve the mixed formulation (16). Concerning the speed of convergence of this algorithm, it follows, for instance, from [13] that kλn ε−λεkL2(QT)≤2qν(Aε,r)pν(Aε,r)−1 pν(Aε,r) + 1n kλ0 ε−λεkL2(QT),∀n≥1 where λεminimizes J⋆⋆ ε,r.ν(Aε,r) = kAε,rkkA−1 ε,rkdenotes the condition number of the operator Aε,r. Eventually, once the above algorithm has converged we can compute ϕε∈Φεas solution of aε,r(ϕε,ϕ) + bε(ϕ, λε) = lε(ϕ),∀ϕ∈Φε. We use the same spaces Φε,h and Mε,h as described in Section 3.1. In practice, each iteration amounts to solve a linear system involving the matrix Aε,r,h of size nh= 4mh(see (42)) which is sparse, symmetric and positive definite. We use the Cholesky method. From the previous estimate, the performances of the algorithm are related to the condition number of the operator Aε,r restricted to Mε,h ⊂L2(QT), which coincides here (see [4]) with the condition number of the symmetric and positive definite matrix BhA−1 ε,r,hBT hintroduced in (46). Using again the power iteration algorithm, we obtain that, for any h, the largest eigenvalue of BhA−1 ε,r,hBT his very closed to r−1(and bounded by r−1). This is in agreement with the estimate kAε,rλkL2(QT)≤r−1kλkL2(QT)for all λ∈L2(QT). Consequently, the condition number is expressed in term of rand of the discrete inf-sup constant δε,h as follows : ν(BhA−1 ε,r,hBT h)≈r−1δ−2 ε,r,h. Since, from our observation in Section 3.3, the discrete inf-sup constant δε,r,h is uniformly bounded by below with respect to h, we deduce that the condition number is uniformly bounded by above with respect to the discretization parameter. This implies that the convergence of the sequence {λn ε,h}(n>0), minimizing for J⋆⋆ ε,r over Mε,h is independent of h. This is exactly what we observe from our numerical experiments. Morever, from (48), we get that the number ν(BhA−1 ε,r,hBT h)≈pCε,r,h is very closed to one. We refer to Tables 8 and 9 for the values. We consider the same data as in Section 3.4, that is, ω= (0.2,0.5), y0(x) = sin(πx) and T= 1/2. We take γ= 10−10 as a stopping threshold for the algorithm (that is the algorithm is stopped as soon as the norm of the residue gnat the iterate nsatisfies kgn εkL2(QT)≤10−10kg0 εkL2(QT)). The algorithm is initiated with λ0 ε,h = 0 in QT. We check that the method provides, for the same value of r,εand h, exactly the same approximation λε,h than the previous direct method (see Tables 5, etc). Table 10, 11 and 12, we simply give the number of iterates of the conjugate gradient algorithm for r= 102,r= 1 and A APPENDIX 32 where the support of the inner control evolves in time and takes the form qT:= {(x, t)∈QT;a(t)< x < b(t)t∈(0, T)} with any a, b in C0([0, T],]0,1[). We refer to [6] which examines this case for the wave equation. Eventually, we also mention that this approach which consists in solving directly the optimality conditions of a controllability problem may be employed to solve inverse problems where, for instance, the solution of the heat equation has to be recovered from a partial observation, typically localized on a sub-domain qTof the working domain: actually, the optimality conditions associated to a least-square type functional can be expressed as a mixed formulation very closed to (2.3). This issue will be analyzed in a future work. A Appendix A.1 Appendix : Fourier expansion of the control of minimal L2(ρ0, qT) norm. We expand in term of Fourier series the control of minimal L2(ρ0, qT) norm vfor the (1) and the corresponding controlled solution y. We use these expansions in Section 3.4 to evaluate with respect to hthe error kyε−λε,hkL2(QT)and kρ0(vε−vε,h)kL2(qT)where the sequence (ϕε,h, λε,h) solves the discrete mixed formulation (40). We use the characterization of the couple (yε, vε) in term of the adjoint solution ϕε(see (5)), unique minimizer in L2(Ω) of J⋆ εdefined by (13). We first note (aε,p)(p>0) the Fourier coefficients in l2(N) of the minimizer ϕT,ε ∈L2(0,1) of (13) such that ϕε,T (x) = X p>0 aε,p sin(pπx), x ∈(0,1).(51) The adjoint state takes the form ϕε(x, t) = P p≥1 aε,pecπ2p2(t−T)sin(pπx) in QT. The optimality equation associated to the functional J⋆ εthen reads, DJ⋆ ε(ϕε,T )·ϕT=ZZqT ρ−2 0ϕεϕ dx dt +εZ1 0 ϕε,T ϕT+ (y0, ϕ(·,0)) = 0,∀ϕT∈L2(0,1) and can be rewritten in terms of the (aε,p)p>0as follows : <{ap}p>0,MqT,ε{ap}p>0>=<{ap}p>0,Fy0>∀aε,p ∈l2(N) (52) where MqT,ε denotes a symmetric positive definite matrix and Fy0a vector obtained from the expansion (51). The resolution of the infinite dimensional system (reduced to a finite dimension one by truncation of the sums) allows an approximation of the minimizer ϕT,ε of J⋆ ε. Finally, we use that the control of minimal L2(ρ0, qT) norm is given by vε=ρ−2 0ϕε1ωand find that the corresponding controlled solution may be expanded as follows yε(x, t) = X q>0e−cπ2q2tb0 q+X p≥1 aε,pcq,p(ω)dq,p(t)]sin(pπx),(x, t)∈QT(53) with cp,q(ω) := 2 Zω sin(pπx) sin(qπx)dx;dp,q(t) := Zt 0 ρ−2 0(s)ecπ2(p2(s−T)+q2(s−t)) ds, t ∈(0, T). (b0 q)q>0denotes the Fourier coefficients of the initial data y0∈L2(0,1). A APPENDIX 33 h1.41 ×10−17.07 ×10−23.53 ×10−21.76 ×10−28.83 ×10−3 kL⋆ϕε,hkL2(QT)3.84 ×10−22.90 ×10−29.27 ×10−32.41 ×10−37.78 ×10−4 kρ0(vε−vε,h)kL2(qT) kρ0vεkL2(qT)1.32 ×10−15.90 ×10−23.24 ×10−21.68 ×10−28.57 ×10−3 kyε−λε,hkL2(QT) kyεkL2(QT)1.04 ×10−13.54 ×10−21.48 ×10−27.59 ×10−33.89 ×10−3 kλε,h(·, T )kL2(0,1) 2.02 ×10−11.68 ×10−11.65 ×10−11.67 ×10−11.68 ×10−1 κε4.44 ×1094.20 ×1011 3.84 ×1013 3.25 ×1015 5.72 ×1016 Table 16: Mixed formulation (16) - r= 102and ε= 10−2with ω= (0.2,0.5). h1.41 ×10−17.07 ×10−23.53 ×10−21.76 ×10−28.83 ×10−3 kL⋆ϕε,hkL2(QT)6.19 ×10−21.57 ×10−11.56 ×10−11.50 ×10−16.21 ×10−2 kρ0(vε−vε,h)kL2(qT) kρ0vεkL2(qT)1.02 7.36 ×10−13.65 ×10−11.52 ×10−13.01 ×10−2 kyε−λε,hkL2(QT) kyεkL2(QT)6.74 ×10−15.51 ×10−12.42 ×10−11.05 ×10−11.81 ×10−2 kλε,h(·, T )kL2(0,1) 2.23 ×10−11.76 ×10−17.86 ×10−24.87 ×10−23.28 ×10−2 κε5.31 ×1098.31 ×1011 9.64 ×1013 1.47 ×1016 1.50 ×1018 Table 17: Mixed formulation (16) - r= 102and ε= 10−4with ω= (0.2,0.5). h1.41 ×10−17.07 ×10−23.53 ×10−21.76 ×10−28.83 ×10−3 kL⋆ϕε,hkL2(QT)6.23 ×10−21.63 ×10−11.77 ×10−12.66 ×10−12.24 ×10−1 kρ0(vε−vε,h)kL2(qT) kρ0vεkL2(qT)1.50 1.11 9.53 ×10−18.33 ×10−17.19 ×10−1 kyε−λε,hkL2(QT) kyεkL2(QT)1.08 1.09 9.4×10−17.69 ×10−15.15 ×10−1 kλε,h(·, T )kL2(0,1) 2.24 ×10−11.79 ×10−18.10 ×10−25.67 ×10−21.71 ×10−2 κε5.32 ×1098.59 ×1011 9.86 ×1013 1.84 ×1016 3.07 ×1018 Table 18: Mixed formulation (16) - r= 102and ε= 10−8with ω= (0.2,0.5). h1.41 ×10−17.07 ×10−23.53 ×10−21.76 ×10−28.83 ×10−3 kL⋆ϕε,hkL2(QT)2.86 ×10−17.15 ×10−21.84 ×10−24.86 ×10−31.40 ×10−3 kρ0(vε−vε,h)kL2(qT) kρ0vεkL2(qT)1.11 ×10−16.21 ×10−23.29 ×10−21.68 ×10−28.57 ×10−3 kyε−λε,hkL2(QT) kyεkL2(QT)5.16 ×10−22.84 ×10−21.48 ×10−27.59 ×10−33.89 ×10−3 kλε,h(·, T )kL2(0,1) 1.53 ×10−11.61 ×10−11.65 ×10−11.67 ×10−11.68 ×10−1 κε9.15 ×1082.07 ×1010 8.05 ×1011 3.25 ×1013 1.45 ×1015 Table 19: Mixed formulation (16) - r= 10−2and ε= 10−2with ω= (0.2,0.5). h1.41 ×10−17.07 ×10−23.53 ×10−21.76 ×10−28.83 ×10−3 kL⋆ϕε,hkL2(QT)10.77 3.821 1.018 2.59 ×10−16.56 ×10−2 kρ0(vε−vε,h)kL2(qT) kρ0vεkL2(qT)4.63 ×10−12.23 ×10−11.10 ×10−15.52 ×10−22.74 ×10−2 kyε−λε,hkL2(QT) kyεkL2(QT)1.55 ×10−19.03 ×10−24.08 ×10−22.46 ×10−21.27 ×10−2 kλε,h(·, T )kL2(0,1) 3.22 ×10−22.85 ×10−22.99 ×10−23.08 ×10−23.12 ×10−2 κε3.04 ×1091.33 ×1011 7.55 ×1012 3.88 ×1014 1.96 ×1016 Table 20: Mixed formulation (16) - r= 10−2and ε= 10−4with ω= (0.2,0.5). REFERENCES 34 h1.41 ×10−17.07 ×10−23.53 ×10−21.76 ×10−28.83 ×10−3 kL⋆ϕε,hkL2(QT)21.872 19.388 26.098 28.310 21.249 kρ0(vε−vε,h)kL2(qT)14.989 9.459 6.606 4.175 1.556 kρ0(vε−vε,h)kL2(qT) kρ0vεkL2(qT)1.33 8.43 ×10−15.89 ×10−13.72 ×10−11.38 ×10−1 kyε−λε,hkL2(QT) kyεkL2(QT)5.73 ×10−14.71 ×10−3.51 ×10−12.11 ×10−16.82 ×10−2 kλε,h(·, T )kL2(0,1) 3.31 ×10−21.31 ×10−25.99 ×10−32.83 ×10−38.26 ×10−4 κε4.08 ×1093.04 ×1011 4.54 ×1013 6.79 ×1015 1.30 ×1018 Table 21: Mixed formulation (16) - r= 10−2and ε= 10−8with ω= (0.2,0.5). A.2 Appendix: Tables Acknowledgments: This work has been partially done while the second author was visiting the Blaise Pascal University (Clermont-Ferrand, France). He wishes to thank the members of the Laboratory of Mathematics for their kind hospitality. References [1] F. Ben Belgacem and S. M. Kaber,On the Dirichlet boundary controllability of the one-dimensional heat equation: semi-analytical calculations and ill-posedness degree, Inverse Problems, 27 (2011), pp. 055012, 19. [2] F. Boyer,On the penalised HUM approach and its applications to the numerical approximation of null-controls for parabolic problems, in CANUM 2012, Super-Besse, ESAIM Proc., EDP Sci., Les Ulis, 2013. [3] F. Boyer, F. Hubert, and J. Le Rousseau,Uniform controllability properties for space/time-discretized parabolic equations, Numer. Math., 118 (2011), pp. 601–661. [4] F. Brezzi and M. Fortin,Mixed and hybrid finite element methods, vol. 15 of Springer Series in Computational Mathematics, Springer-Verlag, New York, 1991. [5] C. Carthel, R. Glowinski, and J.-L. Lions,On exact and approximate boundary controllabilities for the heat equation: a numerical approach, J. Optim. Theory Appl., 82 (1994), pp. 429–484. [6] C. Castro, N. Cˆ ındea, and A. M¨ unch,Controllability of the linear wave equation with moving inner actions., (Preprint available at http://hal.archives-ouvertes.fr/hal-00927076). [7] T. Cazenave and A. Haraux,An introduction to semilinear evolution equations, vol. 13 of Oxford Lecture Series in Mathematics and its Applications, The Clarendon Press Oxford University Press, New York, 1998. Translated from the 1990 French original by Yvan Martel and revised by the authors. [8] D. Chapelle and K.-J. Bathe,The inf-sup test, Comput. & Structures, 47 (1993), pp. 537– 545. [9] P. G. Ciarlet,The finite element method for elliptic problems, vol. 40 of Classics in Applied Mathematics, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2002. Reprint of the 1978 original [North-Holland, Amsterdam; MR0520174 (58 #25001)]. REFERENCES 35 [10] N. Cˆ ındea, E. Fern´ andez-Cara, and A. M¨ unch,Numerical controllability of the wave equation through primal methods and carleman estimates, ESAIM Control Optim. Calc. Var., 19 (2013), pp. 1076–1108. [11] N. Cˆ ındea and A. M¨ unch,A mixed formulation for the direct approximation of the control of minimal L2-norm for linear type wave equations, To appear in Calcolo (http://hal.archivesouvertes.fr/hal-00853767). [12] J.-M. Coron,Control and nonlinearity, vol. 136 of Mathematical Surveys and Monographs, American Mathematical Society, Providence, RI, 2007. [13] J. W. Daniel,The approximate minimization of functionals, Prentice-Hall Inc., Englewood Cliffs, N.J., 1971. [14] D. A. Dunavant,High degree efficient symmetrical Gaussian quadrature rules for the triangle, Internat. J. Numer. Methods Engrg., 21 (1985), pp. 1129–1148. [15] H. W. Engl, M. Hanke, and A. Neubauer,Regularization of inverse problems, vol. 375 of Mathematics and its Applications, Kluwer Academic Publishers Group, Dordrecht, 1996. [16] E. Fern´ andez-Cara and A. M¨ unch,Numerical null controllability of semi-linear 1-D heat equations: fixed point, least squares and Newton methods, Math. Control Relat. Fields, 2 (2012), pp. 217–246. [17] E. Fern´ andez-Cara and A. M¨ unch,Strong convergence approximations of null controls for the 1D heat equation, S~eMA J., 61 (2013), pp. 49–78. [18] ,Numerical exact controllability of the 1d heat equation: Carleman weights and duality, Journal of Optimization Theory and Applications, (2014). [19] M. Fortin and R. Glowinski,Augmented Lagrangian methods, vol. 15 of Studies in Mathematics and its Applications, North-Holland Publishing Co., Amsterdam, 1983. Applications to the numerical solution of boundary value problems, Translated from the French by B. Hunt and D. C. Spicer. [20] A. V. Fursikov and O. Y. Imanuvilov,Controllability of evolution equations, vol. 34 of Lecture Notes Series, Seoul National University Research Institute of Mathematics Global Analysis Research Center, Seoul, 1996. [21] R. Glowinski,Handbook of numerical analysis. Vol. IX, Handbook of Numerical Analysis, IX, North-Holland, Amsterdam, 2003. Numerical methods for fluids. Part 3. [22] R. Glowinski and J.-L. Lions,Exact and approximate controllability for distributed parameter systems, in Acta numerica, 1995, Acta Numer., Cambridge Univ. Press, Cambridge, 1995, pp. 159–333. [23] R. Glowinski, J.-L. Lions, and J. He,Exact and approximate controllability for distributed parameter systems, vol. 117 of Encyclopedia of Mathematics and its Applications, Cambridge University Press, Cambridge, 2008. A numerical approach. [24] S. Kindermann,Convergence rates of the Hilbert uniqueness method via Tikhonov regularization, J. Optim. Theory Appl., 103 (1999), pp. 657–673. [25] S. Labb´ e and E. Tr´ elat,Uniform controllability of semidiscrete approximations of parabolic control systems, Systems Control Lett., 55 (2006), pp. 597–609. REFERENCES 36 [26] I. Lasiecka and R. Triggiani,Control theory for partial differential equations: continuous and approximation theories. I, vol. 74 of Encyclopedia of Mathematics and its Applications, Cambridge University Press, Cambridge, 2000. Abstract parabolic systems. [27] G. Lebeau and L. Robbiano,Contrˆole exact de l’´equation de la chaleur, Comm. Partial Differential Equations, 20 (1995), pp. 335–356. [28] J.-L. Lions,Contrˆolabilit´e exacte, perturbations et stabilisation de syst`emes distribu´es. Tome 1, vol. 8 of Recherches en Math´ematiques Appliqu´ees [Research in Applied Mathematics], Masson, Paris, 1988. Contrˆolabilit´e exacte. [Exact controllability], With appendices by E. Zuazua, C. Bardos, G. Lebeau and J. Rauch. [29] J.-L. Lions and E. Magenes,Non-homogeneous boundary value problems and applications. Vol. III, Springer-Verlag, New York-Heidelberg, 1973. Translated from the French by P. Kenneth, Die Grundlehren der mathematischen Wissenschaften, Band 183. [30] S. Micu and E. Zuazua,On the regularity of null-controls of the linear 1-d heat equation, C. R. Math. Acad. Sci. Paris, 349 (2011), pp. 673–677. [31] A. M¨ unch,A least-squares formulation for the approximation of controls for the stokes system, Math. Control Signals Systems, 27 (2014), pp. 1–27. [32] A. M¨ unch and P. Pedregal,Numerical null controllability of the heat equation through a least squares and variational approach, European J. Appl. Math., 25 (2014), pp. 277–306. [33] A. M¨ unch and E. Zuazua,Numerical approximation of null controls for the heat equation: ill-posedness and remedies, Inverse Problems, 26 (2010), pp. 085018, 39. [34] D. L. Russell,Controllability and stabilizability theory for linear partial differential equations: recent progress and open questions, SIAM Rev., 20 (1978), pp. 639–739. [35] E. Zuazua,Control and numerical approximation of the wave and heat equations, in International Congress of Mathematicians. Vol. III, Eur. Math. Soc., Z¨urich, 2006, pp. 1389–1417.