APPROXIMATION OF CONTROLS FOR LINEAR WAVE EQUATIONS: A FIRST ORDER MIXED FORMULATION

. This paper deals with the numerical approximation of null controls for the wave equation posed in a bounded domain of R n . The goal is to compute approximations of controls that drive the solution from a prescribed initial state to zero at a large enough controllability time. In [Cindea & M¨unch, A mixed formulation for the direct approximation of the control of minimal L 2 -norm for linear type wave equations ], we have introduced a space-time variational approach ensuring strong convergent approximations with respect to the discretization parameter. The method, which relies on generalized observability inequality, requires H 2 -ﬁnite element approximation both in time and space. Following a similar approach, we present and analyze a variational method still leading to strong convergent results but using simpler H 1 -approximation. The main point is to preliminary restate the second order wave equation into a ﬁrst order system and then prove an appropriate observability inequality.


Introduction
This work is devoted to the numerical study of null-controllability for the wave equation by means of a first order equivalent formulation of the wave equation.In the work [11], a mixed formulation was introduced for the numerical approximation of the control of minimal L 2 -norm for the one dimensional wave equation with a potential.More precisely, in [11], under some regularity assumptions on the coefficients, the following kind of control systems were studied: u tt − (c(x)u x ) x + d(x, t)u = 0, (x, t) ∈ (0, 1) × (0, T ), u(0, t) = 0, u(1, t) = v(t), t ∈ (0, T ), u(x, 0) = u 0 (x), u t (x, 0) = u 1 (x), x ∈ (0, 1).
In this work we continue studying mixed formulations for the numerical approximation of controls for wave equations.We shall consider a formulation of the wave equation in terms of a first order hyperbolic system.The main observation is that, if u is a solution to the wave equation (1.1) u tt − ∆u = 0, and we define the new variables (1.2) v u t , p ∇u, then the variables (v, p) solve the following first order system (1.3) v t − div p = 0, Reciprocally, if we are able to establish well-posedness for (1.3) in some suitable functional spaces, by uniqueness we may recover solutions to (1.1); thus, it will be equivalent to solve (1.1) or (1.3).We will return to this point when discussing the well-posedness of initialboundary value problems associated to (1.3).There are two main reasons which motivate the introduction of the variables (v, p).The first reason is that the system (1.3)only involves first order differential operators, a property which will be handy for the finite element approximation of the control of minimal L 2 -norm, since it will allow to reduce the regularity required for the finite element spaces.The second reason is that it is a first step before considering other hyperbolic systems where variables analogous to p are sometime more important than the original variables, such as the system of linear elasticity.This kind of first order -or mixed-formulations of hyperbolic systems have been already considered in the finite element literature; see, for instance [3,4,12].
In this work, the first order formulation will appear through the conjugate functional J defined in (1.8) that one has to minimize in order to obtain the control of minimal L 2 -norm for the wave equation.In this sense, we will be mainly interested in understanding the adjoint state to (1.4) as a solution to a first order system.We will focus on the boundary control.
We recall some facts about the control of minimal L 2 -norm for the simple wave equation.
It is well-known ( [18,20]) that under some geometrical conditions on J the system (1.4) is null-controllable at any large time T > T for some T that depends on Ω and J .Moreover, as a consequence of the Hilbert Uniqueness Method of J.-L.Lions [20], the null-controllability of (1.4) is equivalent to an observability inequality for the associated adjoint problem.We recall that there is an infinite set of possible choices for the control at time T , thus, the problem of finding the control of minimal L 2 -norm arises naturally.In fact, the control of minimal L 2 -norm is unique and it can be found as the solution to a minimization problem involving the adjoint state [18,20].
Remark 1.We recall the reader the following existence, uniqueness and boundary regularity result [20]: given (ϕ 0 , ϕ 1 ) ∈ V and f ∈ L 1 (0, T ; L 2 (Ω)), there exists a unique solution ϕ ∈ C([0, T ]; Moreover, ∂ϕ ∂ν ∈ L 2 (Σ T ) and there exists a constant C = C(Ω, T ) such that A similar result is obtained for backward systems, i.e., when instead of initial conditions at time t = 0 we impose final conditions at t = T .In view of this regularity result, the boundary term which appears in (1.8) makes sense.
The coercivity and boundedness by below of J is the consequence of the following estimate: let ϕ solve (1.7), then there exists k T such that As it was pointed out in [11] (see also [2]), it is possible to reformulate the minimization problem (1.8) in terms of trajectories of the wave equation.This is possible because the wave equation is reversible in time, thus, there is a bijective relation between the trajectories and the initial data.For this purpose, in [11] the following spaces of trajectories was introduced which allows to restate (1.8) as At this point, it is possible to reformulate the minimization problem associated to J in terms of the solution to a first order system.Before doing so, we introduce the initialboundary value problem for the first order system which will be useful for us.
First, we must choose the adequate boundary conditions for (v, p); we observe that, since the solution ϕ to (1.7) satisfies the boundary condition ϕ = 0 on Σ T , then also ϕ t = 0 on Σ T , thus, we shall impose the boundary condition v = 0 on Σ T in order to be consistent with (1.2).In fact, this boundary condition will be enough for our purposes and it will not be necessary to impose boundary conditions on p. Regarding the conditions at time t = T , the natural choice is to set v(•, T ) = ϕ 1 and p(•, T ) = ∇ϕ 0 , i.e., the initial data for (v, p) should belong to where ∇H 1 ∇ϕ, ϕ ∈ H 1 0 (Ω) .The space ∇H 1 is a Hilbert space endowed with the inner product We define the Hilbert space H = L 2 (Ω) × ∇H 1 endowed with the inner product We also define V = H 1 0 (Ω) × H(div) and observe that H is the closure of V with respect to the norm In what follows, we will denote In order to state the suitable well-posedness for the first order system, we will consider the closed and densely defined operator A with domain D(A) = V defined by We show in Section 10 that A generates a C 0 -semigroup of contractions, which we denote e tA (see Lemma 9).
As a consequence of well-known facts concerning the theory of semigroups [23] and the definition of mild solutions, we have the following well-posedness statement for (1.16).
Remark 2. Since −A also generates a C 0 -semigroup, if (g, G) ∈ L 1 (0, T ; H) and (w T , q T ) ∈ H, then the backward problem is also well-posed in the sense of mild solutions.
Similarly to the boundary regularity result stated in Remark 1, we have a boundary regularity result for solutions to (1.18).In this situation we will prove in Lemma 3 that if (w, q) ∈ C([0, T ]; H) is a mild solution to (1.18), then q • ν ∈ L 2 (Σ T ) and for some constant will hold true.
We are now in position to restate the minimization problem (1.13) in terms of a solution to a first order system of the form (1.18).Since we want to identify w with ϕ t and q with ∇ϕ, the natural way to proceed is to consider a minimization problem analogous to (1.13) in which the quantities ∇ϕ and ϕ t are replaced by p and v respectively.
If we denote by U ((g, G), (w 0 , q 0 )) (w, q) the mild backward solution to (1.18) associated to a particular data (g, G) ∈ L 2 (0, T ; H) and (w 0 , q 0 ) ∈ H, then we define the following spaces of trajectories With this definition, the adequate minimization problem in terms of the first order system is (1.19) min As in the case of the wave equation, the coercivity and lower bound of (1.19) can be obtained as a consequence of an observability inequality which takes the form for any solution to (1.18) under suitable conditions on J and T .
The main objective of this work is to adapt the theoretical and numerical analysis done in [11], replacing the extremal problem (1.13) in ϕ, solution of a second order equation by the equivalent extremal problem (1.19) in (w, q), solution of a first order system.Following [11], (1.19) is addressed by means of a mixed formulation which involves (n + 1) Lagrange multipliers.The paper is organized as follows.In Section 2, we prove the generalized observability inequality (2.15) for mild solution of (w, q) of (1.18).This provides an appropriate functional setting and leads to the well-posedeness of the mixed formulation (3.1), which is the optimality system for (1.19).The dual problem, involving only the Lagrange multiplier variables is introduced and analyzed in Section 4. In Section 5, we exploit some properties of the multipliers to derive a stabilized mixed formulation à la Barbosa-Hughes [1] whose solution coincides with the initial mixed formulation.This allows to circumvent the Babuška-Brezzi inf-sup condition.Section 6 introduces a conformal H 1 -approximation of the mixed formulation leading to error estimates in the one dimensional case in Section 7. Section 8 is devoted to numerical experiments for a discontinuous initial condition.Eventually, Section 9 concludes with some perspectives.
As in [11], we use a conformal space-time finite element methods for the discretization of our infinite dimensional extremal problem (1.19), written in a variational form.This "optimize-then-discretize" approach guarantees the strong convergence of the approximation with respect to the discretization parameter.For other space-time finite element methods in similar contexts, we refer to [19] and the references therein.This contrasts with the well-known numerical pathologies (due to discrete high frequencies) appearing when the "discretize-then-optimize" approach is employed (exhibited in [16] and enhanced in [13]) and which requires very specific treatments ( [14,15,21]).

Boundary observability inequality
We will employ a Rellich-Necas type identity.Let β : Ω → R n and p : Ω → R n be two C 1 (Ω; R n ) vector fields.It is easy to check that the following identity holds 2div (( where Ju is the Jacobian matrix (Ju) ij = ∂ j u i of a vector field u = (u 1 , . . ., u n ).
Lemma 2. If p ∈ ∇H 1 is a smooth vector field, then identity (2.1) becomes Proof.If p ∈ V is smooth, then there exists a smooth function ϕ such that p = ∇ϕ.Therefore, the Hessian matrix D 2 ϕ = Jp is symmetric and we have We recall that, if p ∈ H(div), then its normal trace p • ν is well defined as a distribution in H − 1 2 (∂Ω) and Green's formula is satisfied: where •, • denotes the duality pairing between H − 1 2 (∂Ω) and H 1 2 (∂Ω).Next Lemma improves the regularity of p • ν when (v, p) is a solution to (1.3) belonging to C 1 (0, T ; H) ∩ C(0, T ; V ) and allows us to define the normal trace p•ν when (v, p) a mild solution of (1.16).
The regularity of the boundary implies that there exists a smooth vector field β satisfying β = ν on ∂Ω.With this choice of β we integrate identity (2.2) on Q T to obtain (2.7) Employing the equations solved by (v, p) and integrating by parts we find Consequently, the Cauchy-Schwarz inequality implies The inequality (2.8) together with (2.7) yields Finally, we apply (1.17). Let Lemma 4. Let (v, p) be a mild solution to (1.16) with (f, F) = 0.There is T = T (Ω) > 0 such that for any T > T there is a constant C(Ω, T ) > 0 such that: then there exists T > 0 such that for any T > T the following boundary observability estimate holds: (2.12) (u 0 , u 1 ) . We are going to use (2.12) to obtain a boundary observability estimate for (1.16) when (f, F) = 0.If (v 0 , p 0 ) ∈ H, then there exists ϕ 0 ∈ H 1 0 (Ω) such that p 0 = ∇ϕ 0 .We set û0 = ϕ 0 , û1 = v 0 , we let û be the solution to (2.11) with (u 0 , u 1 ) = (û 0 , û1 ) and define (v, p) = (û t , ∇û).If we apply the observability inequality (2.12) to û we obtain that for T > T the following estimate holds: . It is easy to check that (v, p) is a mild solution to (1.16) with (f, F) = 0; in fact it is the unique solution, so (v, p) = (v, p) and the proof is finished.
Lemma 6.The space W endowed with the inner product (•, •) W is a Hilbert space.
Proof.We first have to check that • W is indeed a norm.If (w, q) ∈ L 2 (0, T ; H) is such that (w, q) W = 0, then M (v, q) = 0, therefore, the observability inequality (2.16) implies (w, q) = 0; thus • W is a norm.We now check that W is closed with respect to the norm ⊆ W be a sequence converging to some (w, q) in the norm • W . Then M (w k , q k ) converges to some ( f , F) ∈ L 2 (0, T ; H) and the observability inequality (2.16) implies that (w k (T ), q k (T )) converges to some ( ŵ, q) ∈ H as k → +∞.By the mapping properties of the semigroup which defines the mild solutions, we have that (w k , q k ) converge to the mild solution of
Then, the unique solution ((w, q), (λ, µ)) iii) The solution (w, q) to (3.1) is the minimizer of J over W 0 .The Lagrange multiplier λ is the controlled state.
Since a r ((w, q), (w, q)) = a((w, q), (w, q)) in W 0 , L and L r share the same saddle point.
Proof.i) In virtue of [5, Theorem 4.2.1],we have to check: 1) a and b are continuous: this is obvious by the definition of • W . 2) l is continuous: this is a direct consequence of the generalized observability inequality (2.16) and the energy estimate (1.17).

Dual extremal problem
Here we derive an extremal problem which is dual to (1.19) and only involves the variable (λ, µ).
The condition r > 0 implies that the augmented bilinear form a r is coercive, hence (4.1) has a unique solution.
The ellipticity of A r allows us to introduce a coercive functional J in next Lemma.
From the ellipticity of the operator A r , the minimization of J in H is well posed.We observe that, in contrast with the initial problem of finding the control of minimal L 2 -norm, the minimization of J in H does not entail any constraint.
Due to the symmetry and ellipticity of the operator A r , the conjugate gradient method is well-suited for the numerical minimization of J .The Polak-Ribière version of the conjugate gradient method read as follows: for a given (w 0 , q 0 ) ∈ W (see Lemma 8) and > 0 (convergence criterion), we follow the steps 1), 2) and 3) below: 1) Initialization.An arbitrary initial candidate (λ 0 , µ 0 ) ∈ L 2 (Q T ) is chosen and then we set: 2) Iteration.For n ≥ 0, we compute 3) Convergence test.If (g n+1 , h n+1 ) H ≤ ε g 0 , h 0 H , stop and set ( λ, μ) as the approximation of the minimizer of J .Else, compute Do n = n + 1 and return to step 2).

Stabilized mixed formulation
In this section we introduce a stabilized mixed formulation equivalent to (3.1).It adds some uniform coercivity property with respect to the multiplier variables, property which allows to bypass the Babuška-Brezzi inf-sup condition.
Proof.This is a straightforward consequence of [5, Proposition 4.2.1], the boundedness of a r,α , b α and c α and the coercivity of a r,α and c α in W and Λ respectively.
Proposition 1.The solutions of (5.1) and of the augmented formulation associated to (3.1) coincide.

Finite dimensional conformal approximation
We now focus on the discretization of the mixed formulation (3.1) with the form a replaced by the augmented quadratic form a r introduced in (3.2) and assuming r > 0.
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 and a discrete inf-sup condition.Actually, from the relation (6.3) a r ((w, q), (w, q)) ≥ r η (w, q) 2 W , ∀(w, q) ∈ W the form a r is coercive on the whole space W , therefore, also on N h (b) ⊆ W h ⊆ W .We emphasize that for r = 0, the discrete formulation (6.2) may not be well-posed over W h × Z h because the form a r=0 may not be coercive over the discrete kernel of b: the equality b((w h , q h ), (λ h , µ h )) = 0 for all (λ h , µ h ) ∈ M h does not imply in general that M (w h , q h ) vanishes.Therefore, the term r M (w, q) L 2 (Q T ) which appears in the Lagrangien L r may be understood as a stabilization term: for any h > 0, it ensures the uniform coercivity of the form a r 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 on the variable (w, q).The discrete inf-sup condition reads as follows; for any h > 0, (6.4) Let us assume that this condition holds, so that for any fixed h > 0, there exists a unique couple (w h , q h ), (λ h , µ h ) solution of (6.2).Taking η = r in (2.17), we then have the following estimate which follows from the classical theory of approximations of saddle point problems (see [5,Theorem 5.2.2]).
Proposition 2. Let h > 0. Let (w, q), (λ, µ) and (w h , q h ), (λ h , µ h ) be the solution of (3.1) and (6.2).Let δ h be the discrete inf-sup constant defined by (6.4).Then, where d((w, q), W h ) := inf (w,q)∈W h (w, q)−(w h , q h ) W given by (2.17) and similarly d((λ, µ), where {w h , q h } ∈ R n h denotes the (column) vector associated to (w h , q h ) and the same notation holds for {λ h , µ h } ∈ R n h .Taking this into account, the discrete mixed formulation (6.2) reads as follows: find The matrix A r,h as well as the mass matrix J h are symmetric and positive definite for any h > 0 and any r > 0. On the other hand, the matrix of order m h + n h in (6.6) is symmetric but not necessarily positive definite.We recall (see [5,Theorem 3.2.1])that the inf-sup property (6.4) is equivalent to the injective character of the matrix B T h of size n h × m h , that is Ker(B T h ) = 0.If a necessary condition is given by m h ≤ n h , this property strongly depends on the structure of the space W h and M h .This issue will be numerically analyzed in Section 8.1 and will highlight definitively the role of the parameter r.

One-dimensional discretization
We observe that when n = 1, we have Similarly to what was done in [11], the finite dimensional and conformal space W h must be chosen so that M (w h , q h ) ∈ L 2 (0, T ; H) for any (w h , q h ) ∈ W h ; in the one-dimensional setting, this is guaranteed if both w h and q h have first order space and time derivatives in L 2 loc (Q T ), therefore, in order to have a conformal approximation based on standard triangulations of Q T it will be enough to consider H 1 (Q T ) functions, that is, functions having first order weak derivatives -in both space and time variables-in L 2 (Q T ).This is, at the practical level of the implementation, the main advantages with respect to [11].
We introduce a triangulation T h such that Q T = ∪ K∈T h K and we assume that {T h } h>0 is a regular family.Then we introduce the space W h as follows: with W 1 h and W 2 h are defined by: The space P 2 (K) denotes the space of second order polynomial functions in x and t defined in a triangle K; it involves 6 degrees of freedom, namely, the values of w h -or q h -on the vertices and on the midpoints of the edges of the triangle.
We also define the finite dimensional space (7.2) with M 1 h and M 2 h defined by The space P 1 (K) denotes the space of first order polynomial functions in x and t defined in a triangle K; it involves 3 degrees of freedom: the values of λ h -or µ h -on the vertices of the triangle.
Let us now describe the meshes we use in the 1D -in space-setting.We will use uniform and non-uniform regular meshes (see Figure 1): • Uniform triangular mesh: each element is a right isosceles triangle whose equal sides have length ∆t, ∆x > 0; • Non uniform but regular triangular mesh obtained by Delaunay triangulation.Each triangle of the mesh having a side on the boundaries (0, 1) × {0, T } has the side on the boundary of lenght ∆x; similarly, each triangle having a side on the boundaries {0, 1} × (0, T ) has the side on that boundary of length ∆t.In the numerical experiments we shall consider 5 levels of meshes, that is, we consider (∆x, ∆t) = (1/N, 1/N ) for N = 10 • 2 k for 0 ≤ k ≤ 4. We denote where diam(K) is the diameter of K, the parameter h measures the level of fineness of the mesh and decreases as N increases.We finally remark that other choices for the finite element spaces W h -M h could be considered, however, the numerical experiments suggest that choosing P 1 -based finite element spaces for both W h and M h is not a good choice.Nevertheless, when we consider stabilized formulations (see Section 8.4) we can bypass the inf-sup condition and good approximations of the control are obtained using P 1 -based finite element spaces for both W h and M h .

Numerical experiments
We now comment some computations performed with the FreeFem++ package developed at the University Paris 6 (see [17]), which is very well adapted to our space-time setting.Section 11 provides one of the code we have used.
8.1.The discrete inf-sup condition.Before showing the numerical experiments, we test numerically the discrete inf-sup condition for the choice of finite elements spaces W h and M h in (7.1)-(7.2).Taking η = r > 0 so that a r ((w, q), ( w, q)) = ((w, q), ( w, q)) W for any (w, q), ( w, q) ∈ W , it is readily seen ( [8]) that the discrete inf-sup constant satisfies (8.1) The matrix B h A −1 r,h B T h enjoys the same properties than the matrix A r,h : it is symmetric and positive definite so that the scalar δ h defined in terms of the (generalized) eigenvalue problem (8.1) is strictly positive.This eigenvalue problem is solved using the power iteration algorithm (assuming that the lower eigenvalue is simple): for any The scalar δ h defined by (8.1) is then given by We now give some numerical values of δ h with respect to N (equivalently with respect to h) for the chosen finite element spaces W h − M h with uniform and non uniform meshes.Tables 2 and 3 are concerned with the uniform and non uniform meshes respectively.We observe (see Tables 2 and 3 and Figure 2) that for some fixed values of r > 0, the quantity δ h is not bounded by below as h → 0, in that case, it can be concluded that the pair of finite element spaces considered do not pass the inf-sup test.However, if r is chosen to be equal h α with α ≥ 1, then we observe that for both the uniform and non uniform meshes, the discrete inf-sup constant δ h remains bounded by below.We remark that similar observations were done in [11,Section 4.2], but in that case it was necessary to choose r = h 2α with α ≥ 1 in order to have a lower bound for the discrete inf-sup constant, while in the present case a smaller power of h suffices.
With the non-uniform mesh we observe (Figure 2) that the inf-sup constant δ h behaves as: for some constant C r > 0 uniformly bounded with respect to r. Taking r = h, δ h is therefore of order one so that the general estimate of Proposition 2 leads to (w, q) − (w h , q h ) W + (λ, µ) − (λ h , µ h ) W ≤ h −1/2 d((w, q), W h ) + d((λ, µ), M h ) from which, using standard interpolation estimates for elliptic problems (given in [9, Chapter III]) and assuming regularity on the solution of (3.1), we can obtain the strong convergence of the approximation with respect to h (we refer to [10,Section 4]).Without additional regularity on the solution, precisely on the pair (w, q) ∈ C([0, T ], H), the convergence to zero of d((w, q), W h ) requires more care and will be detailed in a distinct work.3) for r = 1 ( ), r = 10 −1 ( ), r = 10 −2 ( ) , r = h ( ), r = h 2 (•).
On the other hand, when W h is replaced by W h based on P 1 approximation (see (8.5)), our simulation suggests that the discrete inf-sup constant is much smaller and not uniformly bounded by below with respect to h.As shown in Table 4, this property seems independent of r.This is a well-known stiff problem (considered in [21,11,7]) where the initial position belongs to L 2 (0, 1) but is discontinuous.The corresponding control of minimal L 2 -norm, discontinuous as well, is given by v(t) = 2(1 − t)1 (1/2,3/2) (t), t ∈ (0, T ), so that v L 2 (0,T ) = 1/ √ 3 ≈ 0.5773.We recall the reader that in the continuous problem the Lagrange multiplier λ coincides with the controlled solution of (8.3) and q = ϕ x , where ϕ is the solution to the minimization problem (1.13), therefore, both λ(1, t) and q(1, t) coincide with the theoretical control v and the restriction of the discrete solutions λ h and q h to x = 1 should approximate v.
Tables 5 and 6 collect some values of q h L 2 (0,T ) , v − q h L 2 (0,T ) , λ h L 2 (Q T ) and M (w h , q h ) L 2 (Q T ) for different values of r and h with uniform and non uniform meshes.Figure 3 depicts the dependence of v − q h (1, •) L 2 (0,T ) with respect to h for r = 1 and r = h.We observe similar rates of convergence for λ h and q h are similar, however, the error is smaller for r = h.We also observe (see Table 5) that the error is smaller for the non uniform mesh than for the uniform mesh.We observe the following rates of convergence (see Figure 3) for R 1 (h) = v − q h L 2 (0,T ) and R 2 (h) = M (w h , q h ) L 2 (Q T ) using the uniform mesh: r = 1 : R 1 (h) ≈ e −4.38 h 0.64 , R 2 (h) ≈ e −2.87 h 0.39 , r = 10 −2 : R 1 (h) ≈ e −5.35 h 0.7 , R 2 (h) ≈ e 0.89 h 0.51 , r = h : R 1 (h) ≈ e −4.74 h 0.62 , R 2 (h) ≈ e −1.57h 0.47 , r = h 2 : R 1 (h) ≈ e −5.19 h 0.66 , R 2 (h) ≈ e 0.21 h 0.51 .
For the non uniform mesh we observe the following rates of convergence:   6. r = h -non uniform mesh.
Figure 4 shows an initial non uniform but regular mesh and two adaptative refinements of it.We observe that the elements get concentrated along the jumps displayed by the primal variable λ, which are generated by the initial data u 0 , discontinuous at the point x = 1/2.The adaptative refinement together with the use of P 1 elements for the variable λ leads to an excellent approximation of the control (Figure 6) defined as the restriction of λ h at x = 1. Figure 5 depicts a level sel of λ h .We have used r = 10 −6 here.8.3.Conjugate gradient for J .We illustrate here Section 4 and minimize the functional J : L 2 (Q T ) → R with respect to the variable λ.This minimization corresponds to the resolution of the mixed formulation (3.1) by an iterative Uzawa type procedure.The conjugate gradient algorithm is given at the end of Section 4. In practice, each iteration amounts to solve a linear system involving the matrix A r,h (see (6.6)) which is sparse, symmetric and positive definite.We use the Cholesky method.The performances of the algorithm are related to the condition number of operator A r restricted to M h , which coincides here (see [5]) with the condition number, denoted by ν(B h A −1 r,h B T h ) of the symmetric and positive definite matrix B h A −1 r,h B T h introduced in (8.1).Arguing as in [11], Section 4.4, we obtain that ν(  better situation than the quadratic behavior observed in general for elliptic problems, in particular in [11]. We take ε = 10 −10 as a stopping threshold for the algorithm.The algorithm is initialized with (λ 0 , µ 0 ) = (0, 0) in Q T .Tables 8 display the results for r = 1, r = 10 −2 and r = h.We check that the minimization of J leads exactly to the same result.We recall that the norm of the control is v L 2 (0,T ) ≈ 0.5773.Moreover, we observe that the number of iterates is sub-linear with respect to h, with a low influence with respect to the value of r.This is in contrast with the behavior of the conjugate gradient algorithm when this latter is used to minimize J with respect to (ϕ 0 , ϕ 1 ) (see [7,21]).
Eventually, additional experiments not reported here suggest the divergence of the solution of (6.2) as h decreases when only P 1 approximation is used.This is in agreement with the behavior of δ h in that case, see Table 4. 8.4.Discretization of the stabilized formulation.We now discuss the numerical discretization of (5.1).The implementation of the discrete stabilized mixed formulation follows  the lines of Section 6.There are only two important changes: the appearance of the parameter α in (5.1) (compare with (3.1)) and the possibility of choosing different finite elements.In fact, the latter feature is the main motivation for introducing a stabilized formulation, since it allows us to work with finite element spaces which may not satisfy the inf-sup property, uniformly with respect to h.The stabilization formulation (5.1) coupled with the P 1 /P  Accordingly, we only discuss a P 1 /P 1 based approximation and we define (8.5) with W 1 h and W 2 h as follows We also define the finite dimensional space (8.6) 2).Tables 9 and 10 give the results for α = 0.5 and α = h 2 respectively.As expected, the stabilized α-terms allows to recover the convergence of the approximation with respect to h.
Compared with the richer P 1 /P 2 approximation discussed in the previous sections, we also check that the rate of convergence is reduced : precisely, we observe the following rates of convergence for R 1 (h) = v − q h L 2 (0,T ) and R 2 (h) = M (w h , q h ) L 2 (Q T ) using the non uniform mesh and with α = 0.5:

Conclusion -Perspective
We have introduced a space-time variational formulation which allows to characterize null controls for the wave equation posed in R n , restate as a first order system of n + 1 variables.A generalized observability inequality holds true for this system, leading to the wellposedness of the mixed formulation.Appropriate conformal finite dimensional discretization based on H 1 approximation both in time and space provides a convergent approximation of the solution and therefore allows to construct a sequence λ h (1, •) strongly convergent in L 2 (0, T ) toward the null control of minimal L 2 norm.An example of appropriate conformal approximation is obtained by considering piecewise P 2 finite element spaces for the primal variables and piecewise P 1 finite elements spaces for the dual (or lagrange multipliers) variables.Moreover, a stabilized procedure à la Barbosa-Hugues allows to consider any conformal approximation, for instance the one based on piecewise P 1 finite elements, both for the primal and dual variables.The numerical experiments based on a very stiff example for which the initial condition to be controlled is discontinuous support the analysis and show the robustness of the approach.
Much of the methods presented here could be applied to the numerical controllability of more general wave equations.For instance, we must mention that wave equations of the form u tt −div(A(x)∇u) = f , can be handled in an analogous way as long as its corresponding observability observability is available.In that case the natural choice for the variables (v, p) in order to rewrite the equation as a first order system would be v = u t and p = A∇u; in that case, we would have the following equivalent system v t − div p = f, p t − A∇v = 0.In this direction, we mention the elasticity system for which an observability holds true (see [20] chapter 3) if n controls are introduced and for which first order formulations are available (see for instance [4]).Another direction in which we could generalize this work is in considering a wave equation with a zero order term u tt − ∆u + qu = 0 leading to the equivalent first order system v t − div p + q t 0 v(•, τ )dτ = q u 0 , p t − ∇v = 0, involving a non local term.Eventually, we mention that this kind of methodology based on conformal H 1 approximation may be employed to address inverse problems as done in [22] in a parabolic context, following [10].

Appendix -Semigroup generation
Let A be the closed and densely defined operator A : D(A) ⊆ H → H with D(A) = V given by A(v, p) = (div p, ∇v).Lemma 9. A and −A generate strongly continuous semigroups of contractions.
i) Let λ > 0, then for any G ∈ H we must show that there exists a unique U ∈ D(A) such that λU − AU = G; this amounts to show that for any λ > 0 and (f, F) ∈ H there are (v, p) ∈ D(A) such that (10.1) λv − div p = f, λp − ∇v = F.
ii) In order to show that R λ ≤ λ −1 , we multiply the equations in (10.1) by v and p respectively, we integrate in Ω and sum the two expressions to get which implies the uniqueness of (v, p) and the desired estimate for R λ .
The structure of the operator A allows to repeat the same proof for −A.

Figure 4 .
Figure 4. Iterative refinement of the triangular mesh over Q T with respect to the variable λ h : 110, 2 880 and 8 636 triangles.

Table 3 .
δ h w.r.t.r and h, T = 2, for the W h -M h finite elements and non uniform mesh.

Table 4 .
δ h w.r.t.r and h, T = 2, for the W h -M h finite elements and uniform mesh.

Table 8 .
Non uniform mesh -Conjugate gradient method -Number of iterates for r = 1 (top), r = 10 −2 and r = h (bottom).approximation (7.1)-(7.2) leads exactly to the same results.This is expected since such P 1 /P 2 discretization passes the inf-sup test.