Umberto Biccari
email: umberto.biccari@deusto.es

† Víctor Hernández-Santamaría

Controllability of a one-dimensional fractional heat equation: theoretical and numerical aspects

We analyze the controllability problem for a one-dimensional heat equation involving the fractional Laplacian (-d 2

x ) s on the interval (-1, 1). Using classical results and techniques, we show that, acting from an open subset ω ⊂ (-1, 1), the problem is null-controllable for s > 1/2 and that for s ≤ 1/2 we only have approximate controllability. Moreover, we deal with the numerical computation of the control employing the penalized Hilbert Uniqueness Method (HUM) and a finite element (FE) scheme for the approximation of the solution to the corresponding elliptic equation. We present several experiments confirming the expected controllability properties.

Introduction and main results

Let ω ⊂ (-1, 1) be an open and nonempty subset. In this work, we consider the following nonlocal one-dimensional heat equation defined on the domain (-1, 1) × (0, T )

             z t + (-d 2
x ) s z = g1 ω , (x, t) ∈ (-1, 1) × (0, T ) z = 0, (x, t) ∈ [ R \ (-1, 1) ] × (0, T ) z(x, 0) = z 0 (x),

x ∈ (-1, 1),

where z 0 ∈ L 2 (-1, 1) is a given initial datum. In (1.1), for all s ∈ (0, 1), (-d 2 x ) s denotes the one-dimensional fractional Laplace operator, which is defined as the following singular integral (-d 2

x ) s u(x) = c 1,s P.V. R u(x)u(y) |x -y| 1+2s dy.

(1.2)

Here, c 1,s is a normalization constant given by

c 1,s = s2 2s Γ 1+2s 2 √ πΓ(1 -s) ,
where Γ is the usual Gamma function. Moreover, we have to mention that, for having a completely rigorous definition of the fractional Laplace operator, it is necessary to introduce also the class of functions u for which computing (-d 2 x ) s u makes sense. We postpone this discussion to the next section.

The analysis of non-local operators and non-local PDEs is a topic in continuous development. A motivation for this growing interest relies in the large number of possible applications in the modeling of several complex phenomena for which a local approach turns up to be inappropriate or limiting. Indeed, there is an ample spectrum of situations in which a non-local equation gives a significantly better description than a PDE of the problem one wants to analyze. Among others, we mention applications in turbulence ( [START_REF] Bakunin | Turbulence and diffusion: scaling versus equations[END_REF]), anomalous transport and diffusion ( [START_REF] Bologna | Anomalous diffusion associated with nonlinear fractional derivative Fokker-Planck-like equation: exact time-dependent solutions[END_REF][START_REF] Meerschaert | Fractional calculus, anomalous diffusion, and probability[END_REF]), elasticity ( [START_REF] Dipierro | Dislocation dynamics in crystals: a macroscopic theory in a fractional Laplace setting[END_REF]), image processing ( [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF]), porous media flow ( [START_REF] Vázquez | Nonlinear diffusion with fractional Laplacian operators[END_REF]), wave propagation in heterogeneous high contrast media ( [START_REF] Zhu | Modeling acoustic wave propagation in heterogeneous attenuating media using decoupled fractional Laplacians[END_REF]). Also, it is well known that the fractional Laplacian is the generator of s-stable processes, and it is often used in stochastic models with applications, for instance, in mathematical finance ( [START_REF] Levendorskii | Pricing of the American put under Lévy processes[END_REF][START_REF] Pham | Optimal stopping, free boundary, and American option in a jump-diffusion model[END_REF]).

One of the main differences between these non-local models and classical Partial Differential Equations is that the fulfillment of a non-local equation at a point involves the values of the function far away from that point.

In the present paper, we are interested in studying control properties for equation (1.1). In more detail, we aim to give an answer to the following question:

given any T > 0 and any initial datum z 0 ∈ L 2 (-1, 1), can we find a control function g ∈ L 2 (ω×(0, T )) such that the corresponding solution to (1.1) satisfies z(x, T ) = 0? It is well-known that the classical local heat equation (as well as many other general variants) is null-controllable for any time T > 0 (see, e.g., [START_REF] Fattorini | Exact controllability theorems for linear parabolic equations in one space dimension[END_REF][START_REF] Fursikov | Controllability of evolution equations[END_REF][START_REF] Lebeau | Contrôle exact de léquation de la chaleur[END_REF]). Nevertheless, to the best of our knowledge, there are few results in the literature on the null-controllability of the fractional heat equation, and none of them is for a problem involving the fractional Laplacian in its integral form (1.2). The existing results ( [START_REF] Micu | On the controllability of a fractional order parabolic equation[END_REF][START_REF] Miller | On the controllability of anomalous diffusions generated by the fractional Laplacian[END_REF]), instead, deal with the so-called spectral fractional Laplace operator, whose definition will be given later.

In this paper, we deal with the controllability of (1.1), both from the theoretical and the numerical point of view. Employing spectral analysis techniques based on the works [START_REF] Kulczycki | Spectral properties of the Cauchy process on half-line and interval[END_REF][START_REF] Kwaśnicki | Eigenvalues of the fractional Laplace operator in the interval[END_REF], the first main result that we obtain is the following Theorem 1.1. Given any z 0 ∈ L 2 (-1, 1), the parabolic problem (1.1) is null-controllable at time T > 0 with a control function g ∈ L 2 (ω × (0, T )) if and only if s > 1/2. Furthermore, even if for s ≤ 1/2 null controllability for (1.1) fails, we still have the following result of approximate controllability, as a consequence of unique continuation properties for the fractional Laplace operator ( [START_REF] Fall | Unique continuation property and local asymptotics of solutions to fractional elliptic equations[END_REF]).

Theorem 1.2. Let s ∈ (0, 1). Given any z 0 ∈ L 2 (-1, 1), there exists a control function g ∈ L 2 (ω × (0, T )) such that the unique solution z to the parabolic problem (1.1) is approximately controllable at time T > 0.

Theorems 1.1 and 1.2 will then find a confirmation in the study of the corresponding numerical control problem. For this purpose, we will employ the penalized Hilbert Uniqueness Method, which relies in the classical works of Glowinski and Lions ([28,[START_REF] Glowinski | Exact and approximate controllability for distributed parameter systems: a numerical approach[END_REF]). This method is very general and it can be applied to a broad class of PDE control problems, see, for instance, [START_REF] Boyer | On the penalised HUM approach and its applications to the numerical approximation of null-controls for parabolic problems[END_REF][START_REF] Boyer | Insensitizing controls for a semilinear parabolic equation: a numerical approach[END_REF][START_REF] Boyer | Approximate controllability conditions for some linear 1D parabolic systems with space-dependent coefficients[END_REF][START_REF] Khodja | Partial null controllability of parabolic linear systems[END_REF]. However, all of the previous works have one thing in common: they are devoted to the study of equations of local nature. Here, we will see that, when dealing with a nonlocal equation as (1.1), new issues arise during the numerical implementation.

In this context, for the resolution of the numerical control problem, we will need a finite element (FE) approximation of the solution to the following non-local Poisson equation

       (-d 2
x ) s u = f, x ∈ (-1, 1) u ≡ 0,

x ∈ R \ (-1, 1). (1.3) In the recent past, the fractional Laplacian has been widely analyzed also from the point of view of numerical analysis. We refer, for instance, to the works [START_REF] Acosta | A short FE implementation for a 2d homogeneous Dirichlet problem of a Fractional Laplacian[END_REF][START_REF] Acosta | A fractional Laplace equation: regularity of solutions and Finite Element approximations[END_REF][START_REF] Borthagaray | Laplaciano fraccionario: regularidad de soluciones y aproximaciones por elementos finitos[END_REF]. There, the authors present a FE scheme for implementing the solution of (1.3) in a bounded domain Ω ⊂ R 2 . In particular, they provide appropriate quadrature rules in order to solve numerically the variational formulation associated to the problem. Moreover, in [START_REF] Acosta | A fractional Laplace equation: regularity of solutions and Finite Element approximations[END_REF][START_REF] Borthagaray | Laplaciano fraccionario: regularidad de soluciones y aproximaciones por elementos finitos[END_REF] it is also developed an accurate analysis of the efficiency of the FE method, employing several existing results. The techniques of the aforementioned works have then been applied in [START_REF] Acosta | Finite element approximations for fractional evolution problems[END_REF], combined with a convolution quadrature approach, for solving evolution equations involving the fractional Laplacian. For the sake of completeness, we also mention [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF], where it is presented a discretization of the spectral fractional Laplacian and its application to the evolutionary case [START_REF] Bonito | The approximation of parabolic equations involving fractional powers of elliptic operators[END_REF], and [START_REF] Nochetto | A PDE approach to fractional diffusion in general domains: a priori error analysis[END_REF], where the same problem is treated applying the well known extension of Caffarelli and Silvestre ([16]).

In the present paper, we propose a FE approximation for the fractional Poisson equation (1.3) which does not require any quadrature rule. Indeed, exploiting the one-dimensional nature of the problem, each entry of the stiffness matrix can be computed explicitly in terms of its position, the parameter s and the mesh size. This, in particular, allows a quick and simple implementation of the control problem.

This paper is organized as follows. In Section 2, we briefly present the functional setting and some existing theory related to the problems that we are going to analyze. In particular, we give a more accurate definition of the fractional Laplace operator, we discuss the known results on the controllability of fractional parabolic equations and we introduce the variational formulation associated to (1.3) (needed for the development of the FE scheme). Section 3 is devoted to the proof of Theorems 1.1 and 1.2. In Section 4, we describe the penalized HUM with its application to our control problem. Moreover, we present our FE method for the elliptic equation (1.3), that then will be employed for the numerical control of (1.1). In Section 5 we present and comment the results of our numerical simulations. Finally, in Appendix A we include the complete details for computing the stiffness matrix associated to our FE scheme.

Preliminary results on the fractional Laplace operator

In this Section, we introduce some preliminary result that will be useful in the remainder of the paper.

We start by giving a more rigorous definition of the fractional Laplace operator, as we have anticipated in Section 1. Let

L 1 s (R) := u : R -→ R : u measurable , R |u(x)| (1 + |x|) 1+2s dx < ∞ .
For any u ∈ L 1 s and ε > 0 we set

(-d 2 x ) s ε u(x) = c 1,s |x-y|>ε u(x) -u(y) |x -y| 1+2s dy, x ∈ R.
The fractional Laplacian is then defined by the following singular integral

(-d 2 x ) s u(x) = c 1,s P.V. R u(x) -u(y) |x -y| 1+2s dy = lim ε→0 + (-d 2 x ) s ε u(x), x ∈ R, (2.1)
provided that the limit exists. We notice that if 0 < s < 1/2 and u is a smooth function, for example bounded and Lipschitz continuous on R, then the integral in (2.1) is in fact not really singular near x (see e.g. [START_REF] Di Nezza | Hitchhiker's guide to the fractional Sobolev spaces[END_REF]Remark 3.1]). Moreover, L 1 s (R) is the right space for which v := (-d 2

x ) s ε u exists for every ε > 0, v being also continuous at the continuity points of u. It is by now well-known (see, e.g., [START_REF] Di Nezza | Hitchhiker's guide to the fractional Sobolev spaces[END_REF]) that the natural functional setting for problems involving the Fractional Laplacian is the one of the fractional Sobolev spaces. Since these spaces are not as familiar as the classical integral order ones, for the sake of completeness, we recall here their definition.

Given s ∈ (0, 1), the fractional Sobolev space H s (-1, 1) is defined as

H s (-1, 1) :=        u ∈ L 2 (-1, 1) : |u(x) -u(y)| |x -y| 1 2 +s ∈ L 2 (-1, 1) × (-1, 1)        .
It is classical that this is a Hilbert space, endowed with the norm (derived from the scalar product)

u H s (-1,1) := u 2 L 2 (-1,1) + |u| 2 H s (-1,1) 1 2 ,
where the term

|u| H s (-1,1) := 1 -1 1 -1 |u(x) -u(y)| 2 |x -y| 1+2s dxdy 1 2
is the so-called Gagliardo seminorm of u. We set

H s 0 (-1, 1) := C ∞ 0 (-1, 1) H s (-1,1)
the closure of the continuous infinitely differentiable functions with compact support in (-1, 1) with respect to the H s (-1, 1)-norm. The following facts are well-known.

• For 0 < s ≤ 1 2 , the identity H s 0 (-1, 1) = H s (-1, 1) holds. This is because, in this case, the C ∞ 0 (-1, 1) functions are dense in H s (-1, 1) (see, e.g., [START_REF] Lions | Non-homogeneous boundary value problems and applications I[END_REF]Theorem 11.1]). [START_REF] Fiscella | Density properties for fractional Sobolev spaces[END_REF]). Finally, in what follows we will indicate with H -s (-1, 1) = (H s (-1, 1)) the dual space of H s (-1, 1) with respect to the pivot space L 2 (-1, 1).

• For 1 2 < s < 1, we have H s 0 (-1, 1) = {u ∈ H s (R) : u = 0 in R \ (-1, 1)} ([
A more exhaustive description of fractional Sobolev spaces and of their properties can be found in several classical references (see, e.g., [START_REF] Adams | Sobolev spaces[END_REF][START_REF] Di Nezza | Hitchhiker's guide to the fractional Sobolev spaces[END_REF][START_REF] Lions | Non-homogeneous boundary value problems and applications I[END_REF]).

Let us now discuss the parabolic equation (1.1). First of all, we mention that the issues of existence, uniqueness and regularity of the solutions have been studied by several authors. Among others, we mention the works [START_REF] Biccari | Local regularity for fractional heat equations[END_REF][START_REF] Fernández-Real | Boundary regularity for the fractional heat equation[END_REF][START_REF] Leonori | Basic estimates for solutions of a class of nonlocal elliptic and parabolic equations[END_REF]. In particular, in [START_REF] Leonori | Basic estimates for solutions of a class of nonlocal elliptic and parabolic equations[END_REF]Theorem 26] it is showed that, assuming z 0 ∈ L 2 (Ω) and g ∈ L 2 (0, T ; H -s (Ω)), then equation (1.1) admits a unique weak solution z ∈ L 2 (0, T ;

H s 0 (Ω)) ∩ C([0, T ]; L 2 (Ω)) with z t ∈ L 2 (0, T ; H -s (Ω)).
Notice that taking as in our case g ∈ L 2 (ω × (0, T )), the same result holds due to the continuous injection of L 2 into H -s .

In this paper we are mainly interested in the study of control properties for the parabolic system (1.1). For the sake of completeness, we include below the definitions of null and approximate controllability. Definition 2.1. System (1.1) is said to be null-controllable at time T > 0 if, for any z 0 ∈ L 2 (-1, 1), there exists g ∈ L 2 (ω × (0, T )) such that the corresponding solution z satisfies z(x, T ) = 0. Definition 2.2. System (1.1) is said to be approximately controllable at time T > 0 if, for any z 0 , z T ∈ L 2 (-1, 1) and any δ > 0, there exists g ∈ L 2 (ω × (0, T )) such that the corresponding solution z satisfies z(x, T )z T L 2 (-1,1) ≤ δ.

We already mentioned that, to the best of our knowledge, there are no results in the literature concerning the controllability of the fractional heat equation involving the integral operator (2.1). The existing ones deal with the spectral definition of the fractional Laplace operator, which is given as follows.

Let {ψ k , λ k } k∈N ⊂ H 1 0 (-1, 1) × R + be the set of normalized eigenfunctions and eigenvalues of the Laplace operator in (-1, 1) with homogeneous Dirichlet boundary conditions, so that {ψ k } k∈N is an orthonormal basis of L 2 (-1, 1) and

       -d 2 x ψ k = λ k ψ k , x ∈ (-1, 1), ψ k (-1) = ψ k (1) = 0.
Then, the spectral fractional Laplacian (-d 2

x ) s S is defined by

(-d 2 x ) s S u(x) = k≥1 u, ψ k λ s k ψ k (x), (2.2)
firstly for u ∈ C ∞ 0 (-1, 1) and then for u ∈ H s 0 (-1, 1) employing a density argument. It is important to notice that the spectral fractional Laplacian and the fractional Laplacian defined as in (2.1) are two different operators. For instance, definition (2.2) depends on the choice of the domain, while the integral definition does not. For a complete discussion on the differences of these two operators, we refer to [START_REF] Servadei | On the spectrum of two different fractional operators[END_REF].

The control problem for the fractional heat equation involving the operator (-d 2 x ) s S has been analyzed in [START_REF] Micu | On the controllability of a fractional order parabolic equation[END_REF], where the authors proved null controllability provided that s > 1/2. For s ≤ 1/2, instead, null controllability does not hold, not even for T large. This negative result is based on the equivalence (consequence of Müntz Theorem, see, e.g., [START_REF] Schwartz | Etude des sommes d'exponentielles[END_REF]Page 24]) between the controllability property (more specifically, the possibility of proving an observability inequality), and the following condition for the eigenvalues of the considered operator

k≥1 1 λ k < ∞, (2.3)
which is clearly not satisfied for the spectral fractional Laplacian when s ≤ 1/2, since in that case the eigenvalues are λ k = (kπ) 2s . Finally, in [START_REF] Miller | On the controllability of anomalous diffusions generated by the fractional Laplacian[END_REF], the same result as in [START_REF] Micu | On the controllability of a fractional order parabolic equation[END_REF] is obtained in a multi-dimensional setting, by means of a spectral observability condition for a negative self-adjoint operator, which allows to prove the null-controllability of the semi-group that it generates. As we anticipated in Section 1, the same null controllability result holds for the parabolic equation (1.1). This will be obtained by means of classical tools ( [START_REF] Fattorini | Exact controllability theorems for linear parabolic equations in one space dimension[END_REF]) and by an explicit approximations of the eigenvalues and the eigenfunctions the fractional Laplacian with homogeneous Dirichlet boundary conditions.

We stress that (1.1) is a different model with respect to the ones analyzed in [START_REF] Micu | An introduction to the controllability of partial differential equations[END_REF][START_REF] Miller | On the controllability of anomalous diffusions generated by the fractional Laplacian[END_REF], since the operators (2.1) and (2.2) are not equivalent.

We conclude this section by introducing the variational formulation associated to equation (1.3), which is at the heart of our numerical method. That is, find u ∈ H s 0 (-1, 1) such that

a(u, v) = 1 -1 f v dx,
for all v ∈ H s 0 (-1, 1), where the bilinear form a(•, •) :

H s 0 (-1, 1) × H s 0 (-1, 1) → R is given by a(u, v) = c 1,s 2 R R (u(x) -u(y))(v(x) -v(y))
|x -y| 1+2s dxdy.

(2.4)

Since the bilinear form a is continuous and coercive, Lax-Milgram Theorem immediately implies existence and uniqueness of solutions to the Dirichlet problem (1.3). In more detail, if f ∈ H -s (-1, 1), then (1.3) admits a unique weak solution u ∈ H s 0 (-1, 1) (see, e.g., [6, Proposition 2.1]). Furthermore, in the literature it is possible to find improved regularity results for the solution to (1.3), both in Hölder and Sobolev spaces. The interested reader may refer, for instance, to [START_REF] Acosta | A fractional Laplace equation: regularity of solutions and Finite Element approximations[END_REF][START_REF] Biccari | Local elliptic regularity for the Dirichlet fractional Laplacian[END_REF][START_REF] Leonori | Basic estimates for solutions of a class of nonlocal elliptic and parabolic equations[END_REF][START_REF] Ros-Oton | The Dirichlet problem for the fractional Laplacian: regularity up to the boundary[END_REF][START_REF] Ros-Oton | The extremal solution for the fractional Laplacian[END_REF].

Proof of the controllability properties

This section is devoted to study the control properties for the parabolic system (1.1). We begin by proving the null controllability result.

Proof of Theorem 1.1. First of all, for all ϕ T ∈ L 2 (-1, 1), let ϕ(x, t) be the unique weak solution to the adjoint system

             -ϕ t + (-d 2 x ) s ϕ = 0, (x, t) ∈ (-1, 1) × (0, T ) ϕ = 0, (x, t) ∈ [ R \ (-1, 1) ] × (0, T ) ϕ(x, T ) = ϕ T (x),
x ∈ (-1, 1).

(3.1)

Multiplying (1.1) by ϕ and integrating over (-1, 1) × (0, T ), it is straightforward to check that z(x, T ) = 0 if and only if

T 0 1 -1 ϕ(x, t)g(x, t)1 ω (x) dxdt = - 1 -1 z 0 (x)ϕ(x, 0) dx, (3.2)
In turn, it is classical that (3.2) is equivalent to the existence of a constant C > 0 such that the following observability inequality holds

ϕ(x, 0) 2 L 2 (-1,1) ≤ C T 0 1 -1 ϕ(x, t)g(x, t)1 ω (x) dx 2 dt, (3.3)
Notice that ϕ can be expressed in the basis of the eigenfunctions of the fractional Laplacian on (-1, 1) with zero Dirichlet boundary conditions. Namely,

ϕ(x, t) = k≥1 ϕ k e -λ k (T -t) k (x), (3.4)
where ϕ k = ϕ T , k and, for k ≥ 1, k (x) are the solutions to the following eigenvalue problem

       (-d 2 x ) s k = λ k k , x ∈ (-1, 1), k ∈ N k = 0, x ∈ R \ (-1, 1).
Now, plugging (3.4) into (3.3), using the orthonormality of the eigenfunctions k as a basis of L 2 (-1, 1) and employing the change of variables Tt → t, the observability inequality becomes

k≥1 |ϕ k | 2 e -2λ k T ≤ C T 0 k≥1 ϕ k g k (t)e -λ k t 2 dt, (3.5)
where

g k = g1 ω , k .
By means of the classical moment method ( [START_REF] Fattorini | Exact controllability theorems for linear parabolic equations in one space dimension[END_REF]), inequalities of the form (3.5) are well known to be true if and only if (2.3) holds and the eigenfunctions k satisfy the following lower bound:

k L 2 (ω) ≥ C > 0, ∀ k ≥ 1, (3.6)
where the constant C is independent of k. The proof of (3.6) is an easy adaptation of the one of [32, Lemma 2]. Moreover, according to [START_REF] Kulczycki | Spectral properties of the Cauchy process on half-line and interval[END_REF][START_REF] Kwaśnicki | Eigenvalues of the fractional Laplace operator in the interval[END_REF] we have

λ k = kπ 2 - (1 -s)π 4 2s + O 1 k .
Therefore, we easily see that the condition (2.3) is satisfied if and only if s > 1/2. If s ≤ 1/2, instead, the series diverges, since it behaves like an harmonic series. In conclusion, the observability inequality (3.3) holds true when s > 1/2, but it is false when s ≤ 1/2. This concludes the proof. Even if for s ≤ 1/2 null controllability for (1.1) fails, Theorem 1.2 ensures that, for all s ∈ (0, 1), we still have approximate controllability. This is consequence of a unique continuation property for the fractional Laplacian, which has been obtained in [START_REF] Fall | Unique continuation property and local asymptotics of solutions to fractional elliptic equations[END_REF].

Proof of Theorem 1.2. It is classical (see, e.g., [START_REF] Micu | An introduction to the controllability of partial differential equations[END_REF]Theorem 5.2]) that the result is true as soon as one has the following unique continuation property for the solution to the adjoint equation (3.1).

Given s ∈ (0, 1) and ϕ T ∈ L 2 (-1, 1), let ϕ be the unique solution to the system (3.1). Let ω ⊂ (-1, 1) be an arbitrary open set. If ϕ = 0 on ω × (0, T ), then ϕ = 0 on (-1, 1) × (0, T ).

(P)

Therefore, we are reduced to the proof of the property (P). To this end, let us recall that ϕ can be expressed in the form (3.4) and let us assume that ϕ = 0 in ω × (0, T ).

(3.7)

Let {ψ k j } 1≤k≤m k be an orthonormal basis of ker(λ k -(-d 2 x ) s ). Then, (3.4) can be rewritten as

ϕ(x, t) = k≥1         m k j=1 ϕ k j ψ k j (x)         e -λ k (T -t) , (x, t) ∈ (-1, 1) × (-∞, T ). Let z ∈ C with η := (z) > 0 and let N ∈ N. Since the functions ψ k j , 1 ≤ j ≤ m k , 1 ≤ k ≤ N are orthonormal, we have that N k=1         m k j=1 ϕ k j ψ k j (x)         e z(t-T ) e -λ k (T -t) 2 L 2 (-1,1) ≤ N k=1         m k j=1 |ϕ k j | 2         e 2η(t-T ) e -2λ k (T -t) ≤ k≥1         m k j=1 |ϕ k j | 2         e 2η(t-T ) e -2λ k (T -t) ≤ Ce 2η(t-T ) ϕ T 2 L 2 (-1,1) .
Hence, letting

w N (x, t) := N k=1         m k j=1 ϕ k j ψ k j (x)         e z(t-T ) e -λ k (T -t) ,
we have shown that

w N (x, t) L 2 (-1,1) ≤ Ce η(t-T ) ϕ T L 2 (-1,1) . Moreover, we have T -∞ e η(t-T ) ϕ T L 2 (-1,1) dt = 1 η ϕ T L 2 (-1,1) +∞ 0 e -τ dτ = 1 η ϕ T L 2 (-1,1)
.

Therefore, we can apply the Dominated Convergence Theorem, obtaining

lim N→+∞ T -∞ w N (x, t) dt = T -∞ lim N→+∞ w N (x, t) dt = T -∞ e z(t-T ) k≥1         m k j=1 ϕ k j ψ k j (x)         e -λ k (T -t) dt = +∞ k≥1 m k j=1 ϕ k j ψ k j (x) T -∞ e z(t-T ) e -λ k (T -t) dt = +∞ k≥1 m k j=1 ϕ k j ψ k j (x) +∞ 0 e -(z+λ k )τ dτ = +∞ k≥1 m k j=1 ϕ k j z + λ k ψ k j (x), x ∈ (-1, 1), (z) > 0. (3.8)
It follows from (3.7) and (3.8) that

+∞ k≥1 m k j=1 ϕ k j z + λ k ψ k j (x) = 0, x ∈ ω, (z) > 0.
This holds for every z ∈ C \ {-λ k } k∈N , using the analytic continuation in z. Hence, taking a suitable small circle around -λ not including {-λ k } k and integrating on that circle we get that

w := m j=1 ϕ j ψ j (x) = 0, x ∈ ω.
According to [START_REF] Fall | Unique continuation property and local asymptotics of solutions to fractional elliptic equations[END_REF]Theorem 1.4] , (-d 2 x ) s has the unique continuation property in the sense that if λ k is an eigenvalue of (-d 2

x ) s on (-1,1) with Dirichlet boundary conditions, and ((-

d 2 x ) s -λ k ) k = 0 in (-1, 1) with k = 0 in ω, then k = 0 in (-1, 1
). This can applied to w , in order to conclude w = 0 in (-1, 1) for every . Since {ψ j } 1≤ j≤m are linearly independent in L 2 (-1, 1), we get ϕ j = 0, 1 ≤ j ≤ m k , ∈ N. It follows that ϕ T = 0 and hence, ϕ = 0 in (-1, 1) × (0, T ), meaning that ϕ enjoys the property (P). As an immediate consequence, we have that our original equation (1.1) is approximately controllable. Our proof is then concluded. Remark 3.1. According to [START_REF] Fall | Unique continuation property and local asymptotics of solutions to fractional elliptic equations[END_REF], the elliptic unique continuation property for the fractional Laplacian holds in any space dimension. In view of that, Theorem 1.2 may be extended also to the case N > 1. On the other hand, the same does not applies to Theorem 1.1. Indeed, the proof of this result use arguments that are designed specifically for one-dimensional problems ( [START_REF] Fattorini | Exact controllability theorems for linear parabolic equations in one space dimension[END_REF]). If one would like to analyze the null-controllability in a general multi-dimensional setting, other tools (for instance Carleman estimates) are needed. As far as we know, these techniques have not been fully developed yet for problems involving the fractional Laplacian on a domain. Remark 3.2. For the sake of simplicity, in the results presented above we focused on the interval (-1, 1). Nevertheless, everything that we did in this section actually holds in the more general case x ∈ (-L, L) and the extension is immediate.

The penalized HUM and its implementation

We devote this Section to the description of the numerical scheme that we are going to employ for solving the control problem. Let us start with a brief description of the so called penalized Hilbert Uniqueness Method (HUM in what follows) that we shall employ for computing the controls for equation (1.1). Here, we will mostly refer to the work of Boyer [START_REF] Boyer | On the penalised HUM approach and its applications to the numerical approximation of null-controls for parabolic problems[END_REF].

Let (E, •, • ) be a Hilbert space whose norm is denoted by • . Let (A, D(A)) be an unbounded operator in E such that -A generates an analytic semi-group in E that we indicate by t → e -tA . Also, we denote (A * , D(A * )) the adjoint of this operator and by t → e -tA * the corresponding semi-group.

Let (U, [•, •]) be another Hilbert space whose norm is denoted by • . Let B be an unbounded operator from U to D(A * ) and let B * : D(A * ) → U be its adjoint. Let T > 0 be given and, for any y 0 ∈ E and v ∈ L 2 (0, T ; U), let us consider the non-homogeneous evolution problem

       y t + Ay = Bv, t ∈ [0, T ] y(0) = y 0 . (4.1)
The well posedness of (4.1) is guaranteed by [START_REF] Coron | Control and nonlinearity[END_REF]Theorem 2.37]. From now on, we will denote the solution at time T corresponding to the initial datum y 0 and the control v by

y v,y 0 (T ) = L T (v, y 0 ). (4.2)
The linear operator

L T (•, •) is then continuous from L 2 (0, T ; U) × E into E.
In the framework of both controllability notions that we introduced in Section 2, if one control exists it is certainly not unique. In the penalized version of the HUM, we look for a control that is solution to a suitable optimization problem. In particular, for any ε > 0, we shall find

v ε = min v∈L 2 (0,T ;U) F ε (v) (4.3)
where

F ε (v) := 1 2 T 0 v(t) 2 dt + 1 2ε L T (v, y 0 ) 2 , ∀v ∈ L 2 (0, T ; U).
Observe that, for any ε > 0, the functional F ε has a unique minimizer in L 2 (0, T ; U) that we denote by v ε . This is due to the fact that F ε is strictly convex, continuous and coercive.

However, the space L 2 (0, T ; U) in which one has to minimize F ε is a quite big one and it depends on the time variable. This makes the optimization problem computationally expensive. This issue can be circumvented by considering the following problem, defined on the smaller space E. Namely, we consider the minimization problem

q T ε = min q∈E J ε (q T ) (4.4)
where

J ε (q T ) := 1 2 T 0 B * e -(T -t)A * q T 2 dt + ε 2 q T 2 + L T (y 0 , 0), q T , ∀q T ∈ E. (4.5)
In fact, it is classical to prove that (4.3) and (4.4) are equivalent since, for any ε > 0, the minimizers v ε and q T ε of the functionals F ε and J ε , respectively, are related through the formula

v ε = B * e -(T -t)A * q T ε , for a.e. t ∈ (0, T ).
Notice also that we can express the approximate and null controllability properties of the system, for a given initial datum y 0 , in terms of the behavior of the penalized HUM approach described above. In particular, we have Theorem 4.1 (Theorem 1.7 of [START_REF] Boyer | On the penalised HUM approach and its applications to the numerical approximation of null-controls for parabolic problems[END_REF]). Problem (4.1) is approximately controllable from the initial datum y 0 if and only if we have

L T (v ε , y 0 ) = y v ε ,y 0 (T ) → 0, as ε → 0.
Problem (4.1) is null-controllable from the initial datum y 0 if and only if we have

M 2 y 0 := 2 sup ε>0 inf L 2 (0,T ;U) F ε < +∞.
In this case, we have

v ε L 2 (0,T ;U) ≤ M y 0 , L T (v ε , y 0 ) ≤ M y 0 √ ε.
Observe that the fractional Laplacian (-d 2 x ) s satisfies the properties required for the operator A (see, e.g., [25, Theorem 2.14]). Therefore, the penalized HUM approach that we just described can be applied to the control problem (1.1). Following the discussion in [START_REF] Boyer | On the penalised HUM approach and its applications to the numerical approximation of null-controls for parabolic problems[END_REF], we expect that, under discretization, the solution z to (1.1) retains the properties of controllability stated in Theorem 4.1. This will be in accordance with the theoretical results obtained in Theorems 1.1 and 1.2.

To this end, let us study the fully-discrete version of (1.1). For any given mesh M and any integer M > 0, we set δt = T/M and we consider an implicit Euler method, with respect to the time variable. More precisely, we consider

         M h z n+1 -z n δt + A h z n+1 = v n+1 h 1 ω , ∀n ∈ {1, . . . , M -1} z 0 = z 0 , (4.6)
where z 0 ∈ R M , M h is the classical mass matrix and A h is a suitable stiffness matrix approximating the fractional Laplacian. We are going to present more details on the construction of (4.6) in the next section.

In system (4.6), v h,δt = (v n h ) 1≤n≤M is a fully-discrete control function whose cost, that is the discrete L 2 δt (0, T ; R M )norm, is defined by

v δt L 2 δt (0,T ;R M ) :=        M i=1 δt|v n | 2 L 2 (R M )        1/2
, and where

| • | L 2 (R M ) stands for the norm associated to the L 2 -inner product on R M (u, v) L 2 (R M ) = h N i=1 u i v i .
With the above notation and according to the penalized HUM strategy, we introduce, for some penalization parameter ε > 0, the following primal fully-discrete functional

F ε,h,δt (v δt ) = 1 2 v δt 2 L 2 δt (0,T ;R M ) + 1 2ε |L h T (v δt , y 0 ) | 2 L 2 (R M ) , ∀ v δt ∈ L 2 δt (0, T ; R M ),
that we wish to minimize onto the whole fully-discrete control space L 2 δt (0, T ; R M ) and where z M is the final value of the controlled problem (4.6). Here L h T (•, •) stands for the discrete version of the operator (4.2). We can apply Fenchel-Rockafellar theory results to obtain the corresponding dual functional, which reads as follows

J ε,h,δt (ϕ T ) = 1 2 L h T (•, 0) ϕ T 2 L 2 δt (0,T ;R M ) + ε 2 |ϕ T | 2 L 2 (R M ) + ϕ T , L h T (0, y 0 ) L 2 (R M ) , ∀ϕ T ∈ L 2 h,δt (-1, 1) (4.7)
where

L h T (•, •) is the adjoint of L h T (•, •) and L h T (•, 0) ϕ T = (1 ω ϕ n ) 1≤n≤M+1
with ϕ = (ϕ n ) 1≤n≤M+1 solution to the adjoint system

         M h ϕ n -ϕ n+1 δt + A h ϕ n = 0, ∀n ∈ {1, . . . , M} ϕ M+1 = ϕ T . (4.8)
Notice that (4.7) is the fully-discrete approximation of (4.5). Moreover, it can be readily verified that this functional has a unique minimizer without any additional assumption on the problem. Therefore, by minimizing (4.7), and from duality theory, we obtain a control function

v ε,h,δt = 1 ω ϕ n ε,h,δt 1≤n≤M
, where ϕ is the solution to (4.8) evaluated in the optimal datum ϕ T ε . Thus, the optimal penalized control always exists and is unique. Deducing controllability properties amounts to study the behavior of this control with respect to the penalization parameter ε, in connection with the discretization parameters.

It is well known that, in general, we cannot expect for a given bounded family of initial data that the fully-discrete controls are uniformly bounded when the discretization parameters h, δt and the penalization term ε tend to zero independently. Instead, we expect to obtain uniform bounds by taking the penalization parameter ε = φ(h) that tends to zero in connection with the mesh size not too fast (see [START_REF] Boyer | On the penalised HUM approach and its applications to the numerical approximation of null-controls for parabolic problems[END_REF]) and a time step δt verifying some weak condition of the kind δt ≤ ζ(h) where ζ tends to zero logaritmically when h → 0 (see [START_REF] Boyer | Uniform controllability properties for space/time-discretized parabolic equations[END_REF]). This fact will be confirmed by the numerical simulations that we are going to present in Section 5.1 below, by observing the behavior of the norm of the control, the optimal energy inf F ε , and the norm of the solution at time T . In this way, with the help of Theorem 4.1, we obtain numerical evidences of the properties of null and approximate controllability for equation (1.1), which are in accordance with the theoretical results in Section 2.

Finite element approximation

According to the formulation of systems (4.6) and (4.8), in order to solve numerically our control problem, a proper approximation A h of the operator (-d 2

x ) s is needed. Furthermore, observe that, for every every time step n, (4.6) and (4.8) are actually discrete equations of the form

I + δtM -1 h A h u h = f h .
The matrices appearing in the above equation will be computed employing a finite element scheme on a uniform mesh.

Computation of the stiffness matrix

A h .
First of all, we recall the variational formulation associated to the elliptic equation (1.3): find u ∈ H s 0 (-1, 1) such that

a(u, v) = 1 -1 f v dx, (4.9)
for all v ∈ H s 0 (-1, 1), where a(u, v) is the bilinear form introduced in (2.4). Let us take a partition of the interval (-1, 1) as follows:

-1 = x 0 < x 1 < . . . < x i < x i+1 < . . . < x N+1 = 1 , with x i+1 = x i + h, i = 0, . . . N. We call M the mesh composed by the points {x i : i = 1, . . . , N}, while the set of the boundary points is denoted ∂M := {x 0 , x N+1 }. Now, define K i := [x i , x i+1 ] and consider the discrete space

V h := v ∈ H s 0 (-1, 1) v | K i ∈ P 1 , (4.10)
where P 1 is the space of the continuous and piece-wise linear functions. Hence, we approximate (4.9) with the following discrete problem: find u h ∈ V h such that

c 1,s 2 R R (u h (x) -u h (y))(v h (x) -v h (y)) |x -y| 1+2s dxdy = 1 -1 f v h dx,
for all v h ∈ V h . If now we indicate with φ i N i=1 a basis of V h , it will be sufficient that the above equality is satisfied for all the functions of the basis, since any element of V h is a linear combination of them. Therefore the problem takes the following form

c 1,s 2 R R (u h (x) -u h (y))(φ i (x) -φ i (y)) |x -y| 1+2s dxdy = 1 -1 f v h dx, i = 1, . . . , N. (4.11)
Clearly, since u h ∈ V h , we have

u h (x) = N j=1 u j φ j (x),
where the coefficients u j are, a priori, unknown. In this way, (4.11) is reduced to solve the linear system A h u = F, where the stiffness matrix A h ∈ R N×N has components

a i, j = c 1,s 2 R R (φ i (x) -φ i (y))(φ j (x) -φ j (y)) |x -y| 1+2s dxdy, (4.12)
while the vector F ∈ R N is given by F = (F 1 , . . . , F N ) with

F i = f, φ i = 1 -1 f φ i dx, i = 1, . . . , N.
Moreover, the basis φ i N i=1 that we will employ is the classical one in which each φ i is the tent function with supp(φ i ) = (x i-1 , x i+1 ) and verifying φ i (x j ) = δ i, j . In particular, for x ∈ {x i-1 , x i , x i+1 } the i th function of the basis is explicitly defined as (see Figure 1)

φ i (x) = 1 - |x -x i | h . (4.13) • • • (x i-1 , 0) (x i , 1) (x i+1 , 0) (x i , 0) y x φ i (x)
Figure 1: Basis function φ i (x) on its support (x i-1 , x i+1 ).

Let us now describe our algorithm for computing the coefficients a i, j . Before that, we shall make the following preliminary comments. 1. The constant c 1,s in the definition (2.1) of the fractional Laplacian is common for all the entries of the matrix.

For the sake of simplicity, we will drop this constant in the following computations.

2. It is evident from the definition (4.12) that A h is symmetric. Therefore, in our algorithm we will only need to compute the values a i, j with j ≥ i.

3. Due to the non-local nature of the problem, the matrix A h is full. However, while computing its components, we will encounter many simplifications, due to the fact that supp(φ i ) ∩ supp(φ j ) = ∅ for j ≥ i + 2.

4. While computing the values a i, j , we will only work on the mesh M, not considering the points of the set ∂M. In this way, we will ensure that the basis functions φ i satisfy the zero Dirichlet boundary conditions. In other words, in our FE approximation we are considering only the functions from φ 1 to φ N . Instead, if we considered the points x 0 and x N+1 , then we would need to introduce in our discretization also the basis functions φ 0 and φ N+1 , which take value one at the boundary, and this would not be consistent with the continuous problem. Figure 2 provides a graphical explanation of this last discussion.

x 0 = -1

φ 1 x 2 x 1 φ 2 x 3 φ 3 x 4 x N-1 φ N x N+1 = 1 x N y x 1 . . . . . . . . .

• •

Figure 2: Basis functions φ i (x) on the whole interval (-L, L).

We now start building the stiffness matrix A h . This will be done it in three steps, since the values of the matrix can be computed differentiating among three well defined regions: the upper triangle, corresponding to j ≥ i + 2, the upper diagonal corresponding to j = i + 1 and the diagonal, corresponding to j = i (see Figure 3). In fact, as it will be clear during our computations, in each of these regions the intersections among the support of the basis functions are different, thus generating different values of the bilinear form. In what follows, we will briefly present which will be the contributions to the matrix in each of these three steps, including the complete computations as an appendix at the end of the paper.

Step 1: j ≥ i + 2 As we mentioned in Remark 4.1 [START_REF] Acosta | A fractional Laplace equation: regularity of solutions and Finite Element approximations[END_REF], in this case we have supp(φ i ) ∩ supp(φ j ) = ∅ (see also Figure 4). Hence, (4.12) is reduced to computing only the integral

a i, j = -2 x j+1 x j-1 x i+1 x i-1 φ i (x)φ j (y) |x -y| 1+2s dxdy. (4.14)
Taking into account the definition of the basis function (4.13), from (4.14) we obtain a i, j = -2

x j+1

x j-1

x i+1

x i-1

1 -|x-x i | h 1 - |y-x j | h |x -y| 1+2s dxdy.
Finally, this last integral can be computed explicitly employing the following change of variables: (x i-1 , 0)

x -x i h = x, y -x i h = ŷ. ( 4
(x i , 1) (x i+1 , 0) (x i , 0) (x j-1 , 0) (x j , 1) (x j+1 , 0) (x j , 0) y x φ i (x) φ j (x)
Figure 4: Basis functions φ i (x) and φ j (x) for j ≥ i + 1. In this case, the supports are disjoint.

In this way, for the elements a i, j , j ≥ i + 2, we get the following values:

a i, j =                                          -h 1-2s 4(k + 1) 3-2s + 4(k -1) 3-2s -6k 3-2s -(k + 2) 3-2s -(k -2) 3-2s 2s(1 -2s)(1 -s)(3 -2s) , k = j -i, s 1 2 -4( j -i + 1) 2 log( j -i + 1) -4( j -i -1) 2 log( j -i -1) s = 1 2 , j > i + 2 +6( j -i) 2 log( j -i) + ( j -i + 2) 2 log( j -i + 2) + ( j -i -2) 2 log( j -i -2), 56 ln(2) -36 ln(3), s = 1 2 , j = i + 2.
Step 2:

j = i + 1
This is the most cumbersome case, since it is the one with the most interactions between the basis functions (see Figure 5). According to (4.12), and using the symmetry of the integral with respect to the bisector y = x, we have

a i,i+1 = R R (φ i (x) -φ i (y))(φ i+1 (x) -φ i+1 (y)) |x -y| 1+2s dxdy = +∞ x i+1 +∞ x i+1 . . . dxdy + 2 +∞ x i+1 x i+1 x i . . . dxdy + 2 +∞ x i+1 x i -∞ . . . dxdy + x i+1 x i x i+1 x i . . . dxdy + 2 x i+1 x i x i -∞ . . . dxdy + x i -∞ x i -∞ . . . dxdy :=Q 1 + Q 2 + Q 3 + Q 4 + Q 5 + Q 6 .
These contributions will be calculated separately, employing changes of variables analogous to (4.15). After several

(x i-1 , 0) (x i , 1) (x i+1 , 0) (x i , 0) (x j-1 , 0) (x j , 1) (x j+1 , 0) (x j , 0) y x φ i (x) φ j (x)
Figure 5: Basis functions φ i (x) and φ i+1 (x). In this case, the intersection of the supports is the interval [x i , x i+1 ]. computations, we obtain

a i,i+1 =              h 1-2s 3 3-2s -2 5-2s + 7 2s(1 -2s)(1 -s)(3 -2s)
, s 1 2 9 ln 3 -16 ln 2, s = 1 2 .

Step 3: j = i

As a last step, we fill the diagonal of the matrix A h , which collects the values corresponding to the case φ i (x) = φ j (x) (see Figure 6). We have

a i,i = R R (φ i (x) -φ i (y)) 2 |x -y| 1+2s dxdy = +∞ x i+1 +∞ x i+1 . . . dxdy + 2 +∞ x i+1
x i+1

x i-1

. . . dxdy + +∞ x i+1

x i-1 -∞

. . . dxdy

+ x i+1
x i-1

x i+1

x i-1

. . . dxdy + 2

x i-1 -∞ x i+1
x i-1

. . . dxdy + +

x i-1 -∞ +∞ xi+1
. . . dxdy

+ x i-1 -∞ x i-1 -∞ . . . dxdy := R 1 + R 2 + R 3 + R 4 + R 5 + R 6 + R 7 . (x i-1 , 0) (x i , 1) = (x j , 1) (x i+1 , 0) (x i , 0) (x j-1 , 0) (x j+1 , 0) (x j , 0) y x φ i (x) = φ j (x) • • • • • • • • • • • Figure 6:
Basis functions φ i (x) and φ j (x). In this case, the two functions coincide.

Once again, the terms R i , i = 1, . . . , 7 will be computed separately, and summing them we will obtain

a i,i =                    h 1-2s 2 3-2s -4 s(1 -2s)(1 -s)(3 -2s) , s 1 2 8 ln 2, s = 1 2 .

Conclusion

Summarizing, we have the following values for the elements of the stiffness matrix A h : for s 1/2

a i, j = -h 1-2s                                      4(k + 1) 3-2s + 4(k -1) 3-2s -6k 3-2s -(k + 2) 3-2s -(k -2) 3-2s 2s(1 -2s)(1 -s)(3 -2s) , k = j -i, k ≥ 2 3 3-2s -2 5-2s + 7 2s(1 -2s)(1 -s)(3 -2s) , j = i + 1 2 3-2s -4 s(1 -2s)(1 -s)(3 -2s) , j = i.
For s = 1/2, instead, we have

a i, j =                                            -4( j -i + 1) 2 log( j -i + 1) -4( j -i -1) 2 log( j -i -1) +6( j -i) 2 log( j -i) + ( j -i + 2) 2 log( j -i + 2) + ( j -i -2) 2 log( j -i -2), j > i + 2 56 ln(2) -36 ln(3), j = i + 2.
9 ln 3 -16 ln 2,

j = i + 1 8 ln 2, j = i.
Remark 4.2. We point out the following facts:

1. The matrix A h has the structure of a N-diagonal matrix, meaning that value of its elements remain constant along its diagonals. This is in analogy with the tridiagonal matrix approximating the classical Laplace operator. Notice, however, that in our case we obtain a full matrix. This is consistent with the nonlocal nature of the operator that we are discretizing.

2. The value of each element a i, j is given explicitly, and it only depends on i, j, s and h. In other words, when approximating the left hand side of (4.11), no numerical integration is needed.

3. For s = 1/2, the elements a i, j do not depend on the value of h which, in turn, is a function of N. This implies that, in this particular case, no matter how many points we consider in our mesh, the matrix A h will always have the same entries.

Computation of the mass matrix M h .

The matrix M h appears when deriving the FE scheme for the parabolic problem (1.1). Adapting the same construction that we presented before for the elliptic equation (1.3), the approximate solution z h shall be written in the form

z h (x, t) := N j=1 z j (t)φ j (x),
where the basis functions {φ j } N j=1 are the same as before while the coefficients {z j (t)} N j=1 have to be computed. Then, the corresponding variational formulation reads as

N j=1                  z j (t) 1 -1 φ j φ i dx m i, j +z j (t) a(φ j , φ i ) a i, j                  = ω gφ i dx g i (t) , i = 1, . . . , N.
Therefore, denoting g(t) := (g 1 (t), . . . , g N (t)) T , the unknown vector z(t) := (z 1 (t), . . . , z N (t)) T is obtained by solving the following linear system

M h z (t) + A h z(t) = g(t), (4.16)
where M h is the mass matrix with entries m i, j and A h is the stiffness matrix composed by the entries a i, j previously computed. Solving (4.16) by means of an implicit Euler method in time, we obtain exactly the formulation (4.6). Moreover, (4.8) is obtained with the same procedure, taking into account that this time, the original problem being homogeneous, the vector g is actually zero.

Numerical results

In this Section, we present the numerical simulations corresponding to the algorithm previously described, and we provide a complete discussion of the results obtained.

First of all, we test numerically the accuracy of our method for the resolution of the elliptic equation (1.3), by applying it to the following problem

(-d 2 x ) s u = 1, x ∈ (-1, 1) u ≡ 0, x ∈ R \ (-1, 1). (5.1)
In this particular case, the unique solution to (5.1) can be computed exactly and it is given in [START_REF] Getoor | First passage times for symmetric stable processes in space[END_REF]. It reads as follows,

u(x) = 2 -2s √ π Γ 1+2s 2 Γ(1 + s) 1 -x 2 s • 1 (-1,1) . (5.2)
In Figure 7, we show a comparison for different values of s between the exact solution (5.2) and the computed numerical approximation. Here we consider N = 50. One can notice that when s = 0.1 (and also for other small values of s), the computed solution is to a certain extent different from the exact solution. However, one should be careful with such result and a more precise analysis of the error should be carried.

In the same spirit as in [START_REF] Acosta | A short FE implementation for a 2d homogeneous Dirichlet problem of a Fractional Laplacian[END_REF], the computation of the error in the space H s 0 (-1, 1) can be readily done by using the definition of the bilinear form, namely where have used the orthogonality condition a(v h , uu h ) = 0 ∀v h ∈ V h . For this particular test, since f ≡ 1 in (-1, 1), the problem is therefore reduced to

u -u h 2 H s 0 (-1,1) = a(u -u h , u -u h ) = a(u, u -u h ) = 1 -1 f (x) (u(x) -u h (x)) dx,
u -u h H s 0 (-1,1) = 1 -1 (u(x) -u h (x)) dx 1/2
where the right-hand side can be easily computed, since we have the closed formula

1 -1 u dx = π 2 2s Γ(s + 1 2 )Γ(s + 3 2 )
and the term corresponding to 1 -1 u h can be carried out numerically. In Figure 8, we present the computational errors evaluated for different values of s and h. The rates of convergence shown are of order (in h) of 1/2. This is in accordance with the following result: Theorem 5.1 ([2, Theorem 4.6]). For the solution u of (4.9) and its FE approximation u h given by (4.11), if h is sufficiently small, the following estimates hold where C is a positive constant not depending on h.

u -u h H s 0 (-1,1) ≤ Ch 1/2 |ln h| f C 1/2-s (-1,1) , if s < 1/2, u -u h H s 0 (-1,1) ≤ Ch 1/2 |ln h| f L ∞ (-1,1) , if s = 1/2 u -u h H s 0 (-1,1) ≤ C 2s-1 h 1/2 |ln h| f C β (-1,1) , if s > 1/2, β > 0 10 -3 10 -2
Moreover, Figure 8 shows that the convergence rate is maintained also for small values of s. This confirms that the behavior shown in Figure 7a is not in contrast with the known theoretical results. Indeed, since it is well-known that the notion of trace is not defined for the spaces H s (-1, 1) with s ≤ 1/2 (see [START_REF] Lions | Non-homogeneous boundary value problems and applications I[END_REF][START_REF] Tartar | An introduction to Sobolev spaces and interpolation spaces[END_REF]), it is somehow natural that we cannot expect a point-wise convergence in this case.

Control experiments

To address the actual computation of fully-discrete controls for a given problem, we use the methodology described, for instance, in [START_REF] Glowinski | Exact and approximate controllability for distributed parameter systems: a numerical approach[END_REF]. We apply an optimization algorithm to the dual functional (4.7). Since these functionals are quadratic and coercive, the conjugate gradient is a natural and quite simple choice.

First of all, starting from the expression (4.7), the gradient of the functional J ε,h,δt (ϕ T ) can be easily computed and it reads as

∇J ε,h,δt (ϕ T ) = L h T L h T (•, 0) ϕ T , 0 + εϕ T + L h T (0, y 0 )
. Hence, the computation of this gradient at each iteration amounts to solve first the homogeneous equation (4.8). Then, set v n = ϕ n 1 ω and finally solve (4.6) with zero initial datum. In this way, the procedure to compute the gradient of (4.7) basically requires to solve two parabolic equations: a homogeneous backward one associated with the final data ϕ T , and a non-homogeneous forward problem with zero initial data.

We present now some results obtained with the described methodology. In accordance with the discussion in Section 4, we use the finite-element approximation of (-d 2

x ) s for the space discretization and the implicit Euler scheme in the time variable. We denote by N the number of points in the mesh and by M the number of time intervals. As discussed in [START_REF] Boyer | Uniform controllability properties for space/time-discretized parabolic equations[END_REF], the results in this kind of problems does not depend too much in the time step, as soon as it is chosen to ensure at least the same accuracy as the space discretization. The same remains true here, and therefore we always take M = 2000 in order to concentrate the discussion on the dependency of the results with respect to the mesh size h and the parameter s.

As we mentioned, we choose the penalization term ε as a function of h. A reasonable practical rule ( [START_REF] Boyer | On the penalised HUM approach and its applications to the numerical approximation of null-controls for parabolic problems[END_REF]) is to systematically choose ε = φ(h) ∼ h 2p where p is the order of accuracy in space of the numerical method employed for the discretization of the spatial operator involved (in this case the fractional Laplacian (2.1)).

We recall that the solution z to (1.1) belongs to the space L 2 (0, T ; H s 0 (-1, 1)) ∩ C([0, T ]; L 2 (-1, 1)). In view of that, we immediately have that z(•, T ) ∈ L 2 (-1, 1). Therefore, we shall choose the value of p as the convergence rate in the L 2 -norm for the discretization of the elliptic problem (1.3). This convergence rate is given by the following result. Theorem 5.2 ([11, Proposition 3.3.2]). Let s ∈ (0, 1), f ∈ L 2 (-1, 1) and u be the solution to (1.3). Given a uniform mesh M with mesh size h, and the space V h defined as in (4.10), let u h be the finite element solution to the corresponding discrete problem. Then, it holds that

u -u h L 2 (-1,1) ≤ C(s, α)h 2α f L 2 (-1,1) ,
where α := min{s, 1/2 -ε}, for all ε > 0.

By virtue of Theorem 5.2, the appropriate value of p that we shall employ is

p = 2α =        2s, for s < 1 2 1 -2ε, for s ≥ 1 2
. We present below the numerical experiments obtained applying our method. We begin by plotting on Figure 9 the time evolution of the uncontrolled solution as well as the controlled solution. Here, we set s = 0.8, ω = (-0.3, 0.8) and T = 0.3, and as an initial condition we take z 0 (x) = sin(πx). The control domain is represented as highlighted zone on the plane (t, x). As expected, we observe that the uncontrolled solution is damped with time, but does not reach zero at time T , while the controlled solution does. In Figure 10, we present the computed values of various quantities of interest when the mesh size goes to zero. More precisely, we observe that the control cost v δt L 2 δt (0,T ;R M ) and the optimal energy remain bounded as h → 0. On the other hand, we see that

|y M | L 2 (R M ) ∼ C φ(h) = Ch 1/2 .
(5.3)

We know that, for s = 0.8, system (1.1) is null controllable. This is now confirmed by (5.3), according to Theorem 4.1. In fact, the same experiment can be repeated for different values of s > 1/2, obtaining the same conclusions.

According to the discussion in Section 2, one can prove that null controllability does not hold for system (1.1) in the case s ≤ 1/2. However approximate controllability can be proved by means of the unique continuation property of the operator (-d 2

x ) s . We would like to illustrate this property in Figure 11. We observe that the results are different from what we obtained in Figure 10. In fact, the cost of the control and the optimal energy increase in both cases, while the target y M tends to zero with a slower rate than h 1/2 . This seems to confirm that a uniform observability estimate for (1.1) does not hold and that we can only expect to have approximate controllability (see Theorem 4.1).

A Explicit computations of the elements of the matrix A h

We present here the explicit computations for each element a i, j of the stiffness matrix, completing the discussion that we started in Section 4. Step 1:

j ≥ i + 2
We recall that, in this case, the value of a i, j is given by the integral

a i, j = -2 x j+1 x j-1 x i+1 x i-1 φ i (x)φ j (y)
|x -y| 1+2s dxdy.

(A.1)

In Figure 12, we give a scheme of the region of interaction (marked in gray) between the basis functions in this case. These are the only regions in which (A.1) will be different than zero. Now, taking into account the definition of the basis function (4.13), the integral (A.1) becomes a i, j = -2

x j+1

x j-1

x i+1

x i-1

1 -|x-x i | h 1 - |y-x j | h |x -y| 1+2s dxdy.
Let us introduce the following change of variables:

x -x i h = x, y -x i h = ŷ.
Then, rewriting (with some abuse of notations since there is no possibility of confusion) x = x and ŷ = y, we get

a i, j = -2h 1-2s 1 -1 1 -1 (1 -|x | )(1 -|y | ) |x -y + i -j | 1+2s dxdy. (A.2)
The integral (A.2) can be computed explicitly in the following way. First of all, for simplifying the notation, let us define k = ji. We have Figure 11: Convergence properties of the method for s < 1/2. Same legend as in Figure 10 = -2h 1-2s (B 1 + B 2 + B 3 + B 4 ).

a i, j = -2h 1-2s 1 -1 1 -1 (1 -|x | )(1 -|y | ) |x -y + i -j | 1+2s dxdy = -2h 1-2s 1 -1 1 -1 (1 -|x | )(1 -|y | ) |x -y -k | 1+2s dxdy = -2h 1-2s 1 0 1 0 (1 -x)(1 -y) (y -x + k) 1+2s dxdy -2h 1-2s 1 0 0 -1 (1 + x)(1 -y) (y -x + k) 1+2s dxdy -2h 1-2s 0 -1 1 0 (1 -x)(1 + y) (y -x + k) 1+2s dxdy -2h 1-2s 0 -1 0 -1 (1 + x)(1 + y) (y -x + k) 1+2s dxdy
These terms B i , i = 1, 2, 3, 4, can be computed integrating by parts several times. In more detail, we have

B 1 = 1 4s(1 -2s) 2k 1-2s - (k + 1) 2-2s -(k -1) 2-2s 1 -s - 2k 3-2s -(k + 1) 3-2s -(k -1) 3-2s (1 -s)(3 -2s) B 2 = 1 4s(1 -2s) -2k 1-2s + 2(k + 1) 2-2s -2k 2-2s 1 -s + 2(k + 1) 3-2s -k 3-2s -(k + 2) 3-2s (1 -s)(3 -2s) B 3 = 1 4s(1 -2s) -2k 1-2s + 2k 2-2s -2(k -1) 2-2s 1 -s + 2(k -1) 3-2s -k 3-2s -(k -2) 3-2s (1 -s)(3 -2s) B 4 = 1 4s(1 -2s) 2k 1-2s - (k + 1) 2-2s -(k -1) 2-2s 1 -s - 2k 3-2s -(k + 1) 3-2s -(k -1) 3-2s (1 -s)(3 -2s) .
Therefore, we obtain

a i, j = -h 1-2s 4(k + 1) 3-2s + 4(k -1) 3-2s -6k 3-2s -(k + 2) 3-2s -(k -2) 3-2s 2s(1 -2s)(1 -s)(3 -2s) . (A.3)
We notice that, when s = 1/2, both the numerator and the denominator of the expression above are zero. Hence, in this particular case, it would not be possible to introduce the value that we just encountered in our code. Nevertheless, this difficulty can be overcome at least in two ways:

1. by setting s = 1/2 in B i , i = 1, 2, 3, 4, before computing these integrals; 2. by computing the limit s → 1/2 in (A.3).

However, since we already have the expressions of B i , i = 1, 2, 3, 4, for a general s, the second approach is actually straightforward and quicker. Indeed, we can easily compute lim

s→ 1 2 -h 1-2s 4(k + 1) 3-2s + 4(k -1) 3-2s -6k 3-2s -(k + 2) 3-2s -(k -2) 3-2s 2s(1 -2s)(1 -s)(3 -2s) = -4(k + 1) 2 log(k + 1) -4(k -1) 2 log(k -1) + 6k 2 log(k) + (k + 2) 2 log(k + 2) + (k -2) 2 log(k -2), O x i-1 x i-1 x i x i x i+1
x i+1

x j-1

x j-1

x j

x j

x j+1

x j+1 y x Figure 12: Interactions between the basis function φ i and φ j when

j ≥ i + 2. if k 2. When k = 2, instead, since lim k→2 (k -2) 2 log(k -2) = 0,
the corresponding value a i, j = a i,i+2 is given by a i,i+2 = 56 ln(2) -36 ln(3).

Step 2: j = i + 1 This is the most cumbersome case, since it is the one with the most interactions between the basis functions (see Figure 5). According to (4.12), and using the symmetry of the integral with respect to the bisector y = x, we have

a i,i+1 = R R (φ i (x) -φ i (y))(φ i+1 (x) -φ i+1 (y)) |x -y| 1+2s dxdy = +∞ x i+1 +∞ x i+1 . . . dxdy + 2 +∞ x i+1 x i+1 x i . . . dxdy + 2 +∞ x i+1 x i -∞ . . . dxdy + x i+1 x i x i+1 x i . . . dxdy + 2 x i+1 x i x i -∞ . . . dxdy + x i -∞ x i -∞ . . . dxdy :=Q 1 + Q 2 + Q 3 + Q 4 + Q 5 + Q 6 .
In Figure 13, we give a scheme of the regions of interactions between the basis functions φ i and φ i+1 enlightening the domain of integration of the Q i . The regions in grey are the ones that produce a contribution to a i,i+1 , while on the regions in white the integrals will be zero.

Le us now compute the terms Q i , i = 1, . . . , 6, separately.

Computation of Q 1

Since φ i = 0 on the domain of integration we have

Q 1 = +∞ x i+1 +∞ x i+1 φ i+1 (x) -φ i+1 (y) |x -y| 1+2s dxdy O x i-1 x i-1 x i x i x i+1 x i+1 x i+2 x i+2 y x Q 1 Q 2 Q 2 Q 3 Q 3 Q 4 Q 5 Q 5 Q 6
Figure 13: Interactions between the basis function φ i and φ i+1 .

= +∞ x i+1 +∞ x i+1 φ i+1 (x) |x -y| 1+2s dxdy - +∞ x i+1 +∞ x i+1 φ i+1 (y) |x -y| 1+2s dxdy = 0.
The fact that this integral is zero is, actually, not surprising since, according to Figure 13, the region of integration for Q 1 is outside of the region of interaction of the basis functions.

Computation of Q 2

We have

Q 2 = 2 +∞ x i+1 x i+1 x i φ i (x)(φ i+1 (x) -φ i+1 (y))
|x -y| 1+2s dxdy.

Now, using Fubini's theorem we can exchange the order of the integrals, obtaining

Q 2 = 2 x i+1 x i φ i (x)φ i+1 (x) +∞ x i+1 dy |x -y| 1+2s dx -2 x i+2 x i+1 x i+1 x i φ i (x)φ i+1 (y) |x -y| 1+2s dxdy = 1 s x i+1 x i φ i (x)φ i+1 (x) (x i+1 -x) 2s dx -2 x i+2 x i+1 x i+1 x i φ i (x)φ i+1 (y) |x -y| 1+2s dxdy = 1 s x i+1 x i 1 -|x-x i | h 1 -|x-x i+1 | h (x i+1 -x) 2s dx -2 x i+2 x i+1 x i+1 x i 1 -|x-x i | h 1 -|y-x i+1 | h |x -y| 1+2s dxdy := Q 1 2 + Q 2 2 .
The two integrals above can be computed explicitly. Indeed, employing the change of variables

x i+1 -x h = x,
and then renaming x = x, Q 1 2 becomes

Q 1 2 = h 1-2s s 1 0 x 1-2s (1 -x) dx = h 1-2s s(2 -2s)(3 -2s) .
For computing Q 2 2 , instead, we introduce the change of variables

x i -x h = x, y -x i+1 h = ŷ, (A.4)
and we obtain

Q 2 2 = -2h 1-2s 1 0 1 0 (1 -x)(1 -y) (y -x + 1) 1+2s dxdy = h 1-2s 2 1-2s + s -2 s(1 -s)(3 -2s) .
Adding the two contributions, we get the following expression for the term Q 2

Q 2 = h 1-2s 2 2-2s + 2s -3 s(2 -2s)(3 -2s) .
Computation of Q 3

In this case, we simply take into account the intervals in which the basis functions are supported, so that we obtain

Q 3 = -2 x i+2 x i+1 x i x i-1 φ i (x)φ i+1 (y) |x -y| 1+2s dxdy = -2 x i+2 x i+1 x i x i-1 1 -|x-x i | h 1 -|y-x i+1 | h |x -y| 1+2s dxdy.
This integral can be computed applying again (A.4), and we get

Q 3 = -2h 1-2s 1 0 0 -1 (1 + x)(1 -y) (y -x + 1) 1+2s dxdy = h 1-2s 13 -5 • 2 3-2s + 3 3-2s + s(2 4-2s -14) + 4s 2 2s(1 -2s)(1 -s)(3 -2s) , (A.5) if s 1/2. If s = 1/2, instead, we have Q 3 = -2 1 0 0 -1 (1 + x)(1 -y) (y -x + 1) 2 dxdy = 1 + 9 ln 3 -16 ln(2).
Notice that this last value could have been computed directly from (A.5), by taking the limit as s → 1/2 in that expression, being this limit exactly 1 + 9 ln 3 -16 ln(2).

Computation of Q 4

In this case, we are in the intersection of the supports of φ i and φ i+1 . Therefore, we have

Q 4 = x i+1 x i x i+1 x i (φ i (x) -φ i (y))(φ i+1 (x) -φ i+1 (y))
|x -y| 1+2s dxdy.

Moreover, we notice that, this time, it is possible that x = y, meaning that Q 4 could be a singular integral. To deal with this difficulty, we will exploit the explicit definition of the basis function. We have (see also Figure 14)

φ i (x) = 1 - x -x i h , φ i+1 (x) = x i+1 -x h , x ∈ (x i , x i+1 ).
Therefore,

(φ i (x) -φ i (y))(φ i+1 (x) -φ i+1 (y)) = y -x h x -y h = - |x -y| 2 h 2 ,
and the integral becomes

Q 4 = - x i+1 x i x i+1 x i |x -y| 1-2s dxdy = - h 1-2s (1 -s)(3 -2s) . (x i , 0) (x i+1 , 1) (x i , 1) (x i+1 , 0) y x φ i (x) φ i+1 (x)
Figure 14: Functions φ i (x) and φ i+1 (x) on the interval (x i , x i+1 ).

Computation of Q 5

Here the procedure is analogous to the one for Q 2 before. Using again Fubini's theorem we have

Q 5 = 2 x i+1
x i φ i (y)φ i+1 (y) Applying again (A.4), it is now easy to check that

Q 5 = Q 2 .
Computation of Q 6

In analogy with what we did for Q 1 , we can show that also Q 6 = 0.

Conclusion

The elements a i,i+1 are now given by the sum 2Q 2 + Q 3 + Q 4 , according to the corresponding values that we computed.

In particular, we have

a i,i+1 =              h 1-2s 3 3-2s -2 5-2s + 7 2s(1 -2s)(1 -s)(3 -2s)
, s 1 2 9 ln 3 -16 ln 2, s = 1 2 .

Step 3: j = i

As a last step, we fill the diagonal of the matrix A h . In this case we have x i-1

x i+1

x i-1

. . . dxdy + 2

x i-1 -∞

x i+1

x i-1

. . . dxdy + +

x i-1 -∞ +∞ xi+1
. . . dxdy

+ x i-1 -∞ x i-1 -∞ . . . dxdy := R 1 + R 2 + R 3 + R 4 + R 5 + R 6 + R 7 . O x i-1 x i-1 x i x i x i+1 x i+1 y R 1 R 2 R 2 R 3 R 1 4 R 2 4 R 3 4 R 4 4 R 5 R 5 R 6 R 7
Figure 15: Interactions between the basis function φ i (x) and φ i (y).

In Figure 15, we give a scheme of the regions of interactions between the basis functions φ i (x) and φ i (y) enlightening the domain of integration of the R i . The regions in grey are the ones that produce a contribution to a i,i , while on the regions in white the integrals will be zero.

Le us now compute the terms R i , i = 1, . . . , 7, separately. First of all, according to Figure 15 we have that R 1 = R 3 = R 6 = R 7 = 0. This is due to the fact that the corresponding regions are all away from the support of the basis functions.

Computation of R 2

Since φ i (y) = 0 on the domain of integrations we have (1 -|x| ) 2 1x dx = 2 ln 16 -4.

Computation of R 4

In this case, we are in the intersection of the supports of φ i (x) and φ i (y). Therefore, we have

R 4 = x i+1
x i-1

x i+1

x i-1

(φ i (x) -φ i (y)) 2 |x -y| 1+2s dxdy.

In order to compute this integral, we divide it in four components as follows:

R 4 =

x i x i-1

x i

x i-1 Moreover, we notice that, due to symmetry reason, we have R 2 4 = R 3 4 . Therefore, we can compute only one of this two terms and add its value twice when building the matrix A h . Also, notice that in these two region it cannot happen that x = y. On the other hand, R 1 4 and R 4 4 may be singular integrals, and we shall deal with them as we did before.

Computation of R 1

Using again the explicit expression of the basis functions we find

(φ i (x) -φ i (y)) 2 = |x -y| 2 h 2 ,
and the integral becomes

R 1 4 = x i x i-1
x i

x i-1

|x -y| 1-2s dxdy = h 1-2s (1s) (3 -2s) .

Computation of R 2 4 In this case, we simply have

R 2 4 =
x i

x i-1

x i+1

x i (φ i (x) -φ i (y)) 2 |x -y| 1+2s dxdy.

Employing (4.15) we obtain Computation of R 4 4 Also in this case we can use the explicit expression of the basis functions and the integral becomes

R 4 4 = x i+1 x i x i+1 x i |x -y| 1-2s dxdy = h 1-2s (1 -s)(3 -2s) = R 1 4 .
Adding the values that we just computed, we therefore obtain

R 4 = 2(R 1 4 + R 2 4 ) =              h 1-2s 8 -8s -2 3-2s 2s(1 -2s)(1 -s)(3 -2s)
, s 1 2 8 ln 3 -8 ln 2, s = 1 2 .

Computation of R 5

Since, once again, φ i (y) = 0 on the domain of integration we have

R 5 = 2 x i-1 -∞ x i+1
x i-1 φ 2 i (x) |x -y| 1+2s dxdy = 2 (1 -|x| ) 2 1 + x dx = 8 ln 2 -4.

Conclusion

The elements a i,i are now given by the sum 2R 2 + R 4 , according to the corresponding values that we computed. In particular, we have

a i,i =                    h 1-2s 2 3-2s -4 s(1 -2s)(1 -s)(3 -2s)
, s 1 2 8 ln 2, s = 1 2 .

Remark 4 . 1 .

 41 The following fact are worth noticing.

Figure 3 :

 3 Figure 3: Structure of the stiffness matrix A h .

8 Figure 7 :

 87 Figure 7: Plot for different values of s.

9 Figure 8 :

 98 Figure 8: Convergence of the error.

Figure 9 :

 9 Figure 9: Time evolution of system (4.6).

Figure 10 :

 10 Figure 10: Convergence properties of the method for the controllability of the fractional heat equation for s = 0.8.

φ 1 φ

 1 i (x)φ i+1 (y) |x -y| 1+2s dxdy i (y)φ i+1 (y) (yx i ) 2s dy -2 i (x)φ i+1 (y)|x -y| 1+2s dxdy.

R 2 2 = h 1 -2s s 1 - 1 ( 1 R 2 = 2 1 - 1

 221111211 x) 2s dxdy.This integral is computed employing (4.15), and we obtainR -|x| ) 2 (1x) 2s dx = h 1-2s 4s -6 + 2 3-2s s(1 -2s)(1s)(3 -2s) , if s 1/2. If s = 1/2, instead, we have

  y) 2 (xy) 1+2s dxdy = h 1-2s 2s 2 -5s + 4 -2 2-2s s(1 -2s)(1s)(3 -2s) , if s 1/2. If s = 1/2, instead, we get y) 2 (xy) 2 dxdy = 3 -4 ln 2.

5 = h 1 -2s s 1 - 1 ( 1 -|x| ) 2 ( 1 2 , if s 1 / 2 .R 5 = 2 1 - 1

 5111121212511 x i-1 ) 2s dxdy.Employing one last time (4.15), we getR + x) 2s dx = h 1-2s 4s -6 + 2 3-2s s(1 -2s)(3 -2s)(1s) = R If s = 1/2, instead, we have

A quick computation shows that taking the limit s → 1 -in A h we recover the tridiagonal matrix of the classical FE approximation of the local Poisson equation. This is in accordance with the same behavior of the continuous operator, as shown, for instance, in[START_REF] Di Nezza | Hitchhiker's guide to the fractional Sobolev spaces[END_REF] Proposition 4.4].

Acknowledgements

The authors wish to acknowledge F. Boyer (Institut de Mathématiques de Toulouse) for fruitful discussion on the numerical implementation of the control problem. Moreover, a special thanks goes to J. Lohéac (Laboratoire de Sciences Numériques de Nantes), for helping with some computations in Section 4.

Funding

The work of Umberto Biccari was partially supported by the European Research Council [694126 DYCON (Dynamic control and numerics of partial differential equations)], by the Ministerio de Economía, Industria y Competividad (MINECO, Spain) [MTM2014-52347] and by the Air Force Office of Scientific Research [FA9550-15-1-0027]. The work of Víctor Hernández-Santamaría was partially supported by the European Research Council [694126 DYCON (Dynamic control and numerics of partial differential equations)].


Langlois David
email: langlois@loria.fr

Raybaud Sylvain
email: raybauds@loria.fr

Smaïli Kamel

LORIA System for the WMT12 Quality Estimation Shared Task

In this paper we present the system we submitted to the WMT12 shared task on Quality Estimation. Each translated sentence is given a score between 1 and 5. The score is obtained using several numerical or boolean features calculated according to the source and target sentences. We perform a linear regression of the feature space against scores in the range [1:5]. To this end, we use a Support Vector Machine. We experiment with two kernels: linear and radial basis function. In our submission we use the features from the shared task baseline system and our own features. This leads to 66 features. To deal with this large number of features, we propose an in-house feature selection algorithm. Our results show that a lot of information is already present in baseline features, and that our feature selection algorithm discards features which are linearly correlated.

Introduction

Machine translation systems are not reliable enough to be used directly. They can only be used to grasp the general meaning of texts or help human translators. Confidence measures detect erroneous words or sentences. Such information could be useful for users to decide whether or not to post-edit translated sentences [START_REF] Specia | Exploiting objective annotations for measuring translation post-editing effort[END_REF][START_REF] Specia | Predicting machine translation adequacy[END_REF] or select documents mostly correctly translated [START_REF] Soricut | Trustrank: Inducing trust in automatic translations via ranking[END_REF]. Moreover, it is possible to use confidence measures to compare outputs from different systems and to recommend the best one [START_REF] He | Bridging SMT and TM with translation recommendation[END_REF]. One can also imagine that confidence measures at word-level could be also useful for a machine translation system to automatically correct parts of output: for example, a translation system translates the source sentence, then, this output is translated with another translation system [START_REF] Simard | Rule-based translation with statistical phrase-based post-editing[END_REF]. This last step could be driven by confidence measures.

In previous works [START_REF] Raybaud | Detecting errors in machinetranslated sentences[END_REF]Raybaud et al., 2009a;Raybaud et al., 2009b) we used stateof-the-art features to predict the quality of a translation at sentence-and word-level. Moreover, we proposed our own features based on previous works on cross-lingual triggers [START_REF] Lavecchia | Discovering phrases in machine translation by simulated annealing[END_REF][START_REF] Latiri | Phrase-based machine translation based on text mining and statistical language modeling techniques[END_REF]. We evaluated our work in terms of Discrimination Error Trade-off, Equal Error Rate and Normalised Mutual Information.

In this article, we compare the features used in the shared task baseline system and our own features. This leads to 66 features which will be detailed in sections 3 and 4. We therefore deal with many features. We used a machine learning approach to perform regression of the feature space against scores given by humans. Machine learning algorithms may not efficiently deal with high dimensional spaces. Moreover, some features may be less discriminant descriptors and then in some cases could add more noise than information. That is why, in this article we propose an in-house feature selection algorithm to remove useless features.

The article is structured as follows. In Section 2, we give an overview of our quality estimation system. Then, in Sections 3 and 4, we describe the features we experimented with. In section 6, we describe the algorithm we propose for feature selec-tion. Then we give the results of several configurations in Section 7.

Overview of our quality estimation submission

Each translated sentence is assigned a score between 1 and 5. 5 means that the machine translation output is perfectly clear and intelligible and 1 means that it is incomprehensible. The score is calculated using several numerical or boolean features extracted according to the source and target sentences. We perform a regression of the feature space against [1 : 5].

The baseline features

The quality estimation shared task organizers provided a baseline system including several interesting features. Among them, several are yet used in [START_REF] Raybaud | Detecting errors in machinetranslated sentences[END_REF] but we give below a brief review of the whole baseline features set1 :

• Source and target sentences lengths: there is a correlation between the sizes of source and target sentences.

• Average source token length: this is the average number of letters of the words in the sentence. We guess that this feature can be useful because short words have more chance to be tool words.

• Language model likelihood of source and target sentences: a source sentence with low likelihood is certainly far from training corpus statistics. There is a risk it is badly translated. A target sentence with low likelihood is not suitable in terms of target language.

• Average number of occurrences of the words within the target sentence: too many occurrences of the same word in the target sentence may indicate a bad translation.

• Average number of translations per source word in the sentence: for each word in the source sentence, the feature indicates how many words of the target sentence are indeed translations of this word in the IBM1 table (with probability higher than 0.2).

• Weighted average number of translations per source word in the sentence: this feature is similar to the previous one, but a frequent word is given a low weight in the averaging.

• n-gram frequency based features: the baseline system proposes to group the n-gram frequencies into 4 quartiles. The features indicate how many n-gram (unigram to trigram) in source sentence are in quartiles 1 and 4. These features indicate if the source sentence contains n-grams relevant to the training corpus.

• Punctuation based features: there may exist a correlation between punctuation of source and target sentences. The count of punctuation marks in both sentences may then be useful.

Overall, the baseline system proposes 17 features.

The LORIA features

In a previous work [START_REF] Raybaud | Detecting errors in machinetranslated sentences[END_REF], we tested several confidence measures. The Quality Measure Task campaign constitutes a good opportunity for us to compare our approach to others. We give below a brief review of our features (we cite again features which are yet presented in baseline features because sometimes, we use a variant of them):

• lengths: three features are generated, lengths of source and target sentences (already presented in baseline features), and ratio of target over source length

• n-gram based features [START_REF] Duchateau | Confidence scoring based on backward language models[END_REF]: each word in the source and target sentences is given its 5-gram probability. Then, the sentence-level score is the average of the scores across all words in the sentence. There are 4 features: one for each language (source and target) and one for each direction (left-to-right and right-to-left 5-gram).

• backoff n-gram based features: in the same way, a score is assigned to a word according to how many times the language model had to back off in order to assign a probability to the sequence [START_REF] Uhrik | Confidence metrics based on n-gram language model backoff behaviors[END_REF]. Here too, word scores are averaged and we get 4 scores.

• averaged features: a common property of all ngram based and backoff based features is that a word can get a low score if it is actually correct but its neighbours are wrong. To compensate for this phenomenon we took into account the average score of the neighbours of the word being considered. More precisely, for every relevant feature x . defined at word level we also computed:

x lef t . (w i ) = x . (w i-2 ) * x . (w i-1 ) * x . (w i ) x centred . (w i ) = x . (w i-1 ) * x . (w i ) * x . (w i+1 ) x right . (w i ) = x . (w i) * x . (w i+1 ) * x . (w i+2 )
A sentence level feature is then calculated according to the average of each new "averaged feature".

• intra-lingual features: the intra-lingual score of a word in a sentence is the average of the mutual information between that word and the other words in that sentence. Mutual information is defined by: I(w 1 , w 2 ) = P (w 1 , w 2 )×log P (w 1 , w 2 ) P (w 1 )P (w 2 )

(1) The intra-lingual score of a sentence is the average of the intra-lingual scores of the words in this sentence. There are two features, one for each language.

• cross-lingual features: the cross-lingual score of a word in a target sentence is the average of the mutual information between that word and the words in the source sentence. The crosslingual score of a target sentence is the average of its constituents.

• IBM1 features: the score of the target sentence is the average translation probability provided by the IBM1 model.

• basic parser: this produces two scores, a binary flag indicating whether any bracketing inside the target sentence is correct, and one indicating if the sentence ends with an end of sentence symbol (period, colon, semi-colon, question/exclamation/quotation mark, comma, apostrophe, close parenthese)

• out-of-vocabulary: this generates two scores, the number of out-of-vocabulary words in the sentence, and the same one but normalized by the length of the sentence. These scores are used for both sides.

This leads to 49 features. A few ones are equivalent to or are strongly correlated to baseline ones.

As we want to be able to integrate several sets of features without prior knowledge, our system is able to discard redundant features (see Section 6).

Regression

Our system predicts a score between 1 and 5 for each test sentence. For that, we used the training corpus to perform the linear regression of the input features against scores given by humans. We used SVM algorithm to perform this regression (LibSVM toolkit [START_REF] Chang | LIBSVM: A library for support vector machines[END_REF]). We experimented two kernels: linear, and radial basis function. For the radial basis function, we used grid search to optimise parameters.

Feature Selection

We experimented with many features. Some of them may be very poor predictors. Then, these features may disturb the convergence of the training algorithm of SVM. To prevent this drawback, we applied an in-house feature selection algorithm. A feature selection algorithm selects the most relevant features by maximizing a criterion. Feature selection algorithms can be divided into two classes: backward and forward [START_REF] Guyon | An introduction to variable and feature selection[END_REF]. Backward algorithms remove useless features from a set. Forward algorithms start with an empty feature set and insert useful features. We implemented a greedy backward elimination algorithm for feature selection. It discards features until a quality criterion stops to decrease. The criterion used is the Mean Average Error (MAE) calculated on the development corpus:

M AE(s, r) = n i=1 |s i -r i | n (2)
where s is the list of scores predicted by the system, r is the list of scores given by experts, n is the size of these lists.

The algorithm is described below: For calculating the MAE for a feature set, several steps are necessary: performing the regression between the features and the expert scores on the training corpus, using this regression to predict the scores on the development corpus, calculate the MAE between the predicted scores and the expert scores on this development corpus.

Results

We used the data provided by the shared task on Quality Estimation 2 , without additional corpus. This data is composed of a parallel English-Spanish training corpus. This corpus is made of the concatenation of europarl-v5 and news-commentary10 corpora (from WMT-2010), followed by tokenization, cleaning (sentences with more than 80 tokens removed) and truecasing. It has been used for baseline models provided in the baseline package by the shared task organizers. We used the same training corpus to train additional language models (forward and backward 5-gram with kneyser-ney discounting, obtained with the SRILM toolkit) and triggers required for our features. For feature extrac-2 http://dl.dropbox.com/u/6447503/resources.tbz tion, we used the files provided by the organizers: 1832 source english sentences, their translations by the baseline translation system, and the score given by humans to these translations. We split these files into a training part (1000 sentences) and a development part (832 sentences). We used the train part to perform the regression between the features and the scores. We used the development corpus to optimise the parameters of the regression and for feature selection. We did not use additional provided information such as phrase alignment, word alignment, word graph, etc.

Table 1 presents our results in terms of MAE and Root Mean Squared Error (RMSE). MAE is described in Formula 2, and RMSE is defined by:

RM SE(s, r) = n i=1 (s i -r i ) 2 n (3)
Each line of Table 1 gives the performance for a set of features. BASELINE+LORIA constitutes the union of both features BASELINE (Section 3) and LORIA (Section 4). the 'feature selection' column indicates if feature selection algorithm is applied. We experimented the SVM with two kernels: linear (LIN in Table 1) and radial basis function (RBF in Table 1). As the radial basis function uses parameters, we proposed results with default values (DEF) and with values optimised by grid search on the development corpus (OPT). MAE and RMSE are given for development corpus and for the test corpus. This test corpus (and its reference scores given by humans) is the one released for the shared 2012 task3 . MAE and RMSE has been computed against the scores given by humans to the translations in this test corpus 4 .

The results show that the performance on development corpus are always confirmed by those of the test corpus. The BASELINE features alone achieve already good performance, better than ours. Although the differences are well inside the confidence interval, the fusion of both sets outperforms slightly the BASELINE. The feature selection algorithm allows to gain 0.01 point. The gain is the same for the optimisation of the radial basis function parameters. Surprisingly, the linear kernel, simpler than other kernels, yields the same performance as radial basis function.

In addition to MAE and RMSE results, we studied the linear correlations between features: our objective is to check if BASELINE and LORIA complement each other. We computed the linear correlation between all features (BASELINE+LORIA). This leads to 2145 values. Table 2 shows in line +/the number of features pairs which correlate with an absolute score higher than thresholds 0.9, 0.8 or 0.7. Among these pairs we give in line + the number of pairs with positive correlation, and in line -the number of pairs with negative correlation. For lines + and -, we give 4 numbers: number of pairs, number of LORIA-LORIA (e.g. the number of correlations between a LORIA feature and another LORIA feature) pairs, number of BASELINE-BASELINE pairs, number of LORIA-BASELINE pairs. We remark that only 6% of the pairs correlates (column 0.7, line +/-) and that the correlations are mostly between LORIA features. This last point is not surprising because there are more LORIA features than BASELINE ones. There are very few correlations between LORIA and BASELINE features. We studied precisely the correlated pairs. There is a strong (more than 0.9) positive correlation between n-gram and backoff based features and their averaged feature versions. Sometimes, there is also a strong correlation between 'forward' and 'backward' features. Source and target sentences lengths linearly correlate (0.98). This is the same case for source and target language model likelihoods. There is also a high correlation between forward and backward 5-gram scores (0.89). There are very few negative correlations between features. As they are not numerous, one can list these pairs with correlation between -1 and -0.7: target sentence length and target language model probability; source sentence length and source language model probability; ratio of OOV words over sentence length in source sentence and percentage of unigrams in the source sentence seen in the SMT training corpus; and number of OOV words in source sentence and percentage of unigrams in the source sentence seen in the SMT training corpus. These correlations are not surprising. First, language model probability is not normalized ≥ 0.9 ≥ 0.8 ≥ 0. by the number of tokens: the more tokens, the lower probability. Second, the more OOV in the sentence, the fewer known unigrams. Last, we present the set of features discarded by our feature selection algorithm. We give only this description for the LORIA+BASELINE set, with linear kernel. The algorithm discards 18 LORIA features out of 49 (37%) and 3 BASELINE out of 17 (18%). The features discarded from LORIA are mostly averaged features based on n-gram and backoff. This is consistent with the fact that these features are strongly correlated with n-gram and backoff features. We remark that very few BASELINE features are discarded: lengths of source and target language because these features are yet included in LORIA features, and "average number of translations per source word in the sentence" maybe because the LORIA feature giving the average IBM1 probabilities is more precise. Last, we remark that the target length feature is discarded, and only ratio between target and source length is kept.

Conclusion

In this paper, we present our system to evaluate the quality of machine translated sentences. A sentence is given a score between 1 and 5. This score is predicted using a machine learning approach. We use the training data provided by the organizers to perform the regression between numerical features calculated from source and target sentences and scores given by human experts. The features are the baseline ones provided by the organizers and our own features. We proposed a feature selection algorithm to discard useless features. Our results show that baseline features contain already the main part of information for prediction. Concerning our own features, a study of the linear correlations shows that averaged features do not provide new information compared to n-gram and backoff features. This last remark is confirmed by our feature selection algorithm. Our feature selection algorithm seems to discard features linearly correlated with others while keeping relevant features for prediction. Last, we remark that the choice of kernel, optimisation of parameters and feature selection have not a strong effect on performance. The main effort may have to be concentrated on features in the future.

Table 2 :

 2 Statistics on the linear correlations between LO-RIA+BASELINE features

	7

Table 1 :

 1 Results of the various sets of features in terms of MAE and RMSE

	Dev	Test

Indeed, our system takes into input a set of features, and is able to discard redundant features (see Section 6).

https://github.com/lspecia/QualityEstimation/blob/master/ test set.tar.gz

available at https://github.com/lspecia/QualityEstimation/ blob/master/test set.likert