Michel Duprez
email: michel.duprez@inria.fr

Vanessa Lleras

Alexei Lozinski

A new φ-FEM approach for problems with natural boundary conditions

We present a new finite element method, called φ-FEM, to solve numerically elliptic partial differential equations with natural (Neumann or Robin) boundary conditions using simple computational grids, not fitted to the boundary of the physical domain. The boundary data are taken into account using a level-set function, which is a popular tool to deal with complicated or evolving domains. Our approach belongs to the family of fictitious domain methods (or immersed boundary methods) and is close to recent methods of cutFEM/XFEM type. Contrary to the latter, φ-FEM does not need any non-standard numerical integration on cut mesh elements or on the actual boundary, while assuring the optimal convergence orders with finite elements of any degree and providing reasonably well conditioned discrete problems. In the first version of φ-FEM, only essential (Dirichlet) boundary conditions was considered. Here, to deal with natural boundary conditions, we introduce the gradient of the primary solution as an auxiliary variable. This is done only on the mesh cells cut by the boundary, so that the size of the numerical system is only slightly increased . We prove theoretically the optimal convergence of our scheme and a bound on the discrete problem conditioning, independent of the mesh cuts. The numerical experiments confirm these results.

Introduction

We consider a second order elliptic partial differential equation with Neumann boundary conditions

-∆u + u = f in Ω, ∂u ∂n = 0 on Γ (1)
in a bounded domain Ω ⊂ R d (d = 2, 3) with smooth boundary Γ assuming that Ω and Γ are given by a level-set function φ: Ω := {φ < 0} and Γ := {φ = 0}.

(

) 2
Such a representation is a popular and useful tool to deal with problems with evolving surfaces or interfaces [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF]. In the present article, the level-set function is supposed known on R d , smooth, and to behave near Γ similar to the signed distance to Γ. Our goal is to develop a finite element method for (1) using a mesh which is not fitted to Γ, i.e. we allow the boundary Γ to cut the mesh cells in an arbitrary manner. The existing finite elements methods on non-matching meshes, such as the fictitious domain/penalty method [START_REF] Glowinski | Error estimates for fictitious domain/penalty/finite element methods[END_REF], XFEM [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Moës | Imposing dirichlet boundary conditions in the extended finite element method[END_REF][START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite-element method[END_REF][START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF], CutFEM [START_REF] Burman | Fictitious domain finite element methods using cut elements: Ii. A stabilized Nitsche method[END_REF][START_REF] Burman | Cutfem: discretizing geometry and partial differential equations[END_REF], CutIGA [START_REF] Elfverson | A new least squares stabilized nitsche method for cut isogeometric analysis[END_REF] (see also [START_REF] Mittal | Immersed boundary methods[END_REF] for a review on immersed boundary methods) contain the integrals over the physical domain Ω and thus necessitate non-standard numerical integration on the parts of mesh cells cut by Γ. In this article, we propose a finite element method, based on an alternative variational formulation on an extended domain matching the computational mesh, thus avoiding any non-standard quadrature while maintaining the optimal accuracy and controlling the conditioning uniformly with respect to the position of Ω over the mesh. 1In the recent article [START_REF] Duprez | φ-FEM: a finite element method on domains defined by level-sets[END_REF], we have proposed such a method for the Poisson problem with homogeneous Dirichlet boundary conditions u = 0 on Γ. The idea behind this method, baptised φ-FEM, is to put u = φw so that u = 0 on Γ for whatever w since φ = 0 there. We then replace φ and w by the finite element approximations φ h and w h , substitute u ≈ φ h w h into an appropriate variational formulation and get an easily implementable discretization in terms of the new unknown w h . Such a simple idea cannot be used directly to discretize the Neumann boundary conditions in [START_REF] Balay | PETSc users manual[END_REF]. Indeed, multiplication by φ works well to strongly impose the essential Dirichlet boundary conditions whereas Neumann conditions are natural, i.e. they come out of the usual variational formulation without imposing them into the functional spaces. We want thus to reformulate Problem (1) so that Neumann conditions become essential. The way to go is the dualization of this problem, in the terminology of [START_REF] Boffi | Mixed finite element methods and applications[END_REF], consisting in introducing an auxiliary (vectorvalued) variable for the gradient ∇u. In the present article, we want to use the usual conforming scalar finite elements as much as possible. Accordingly, we do not pursue the classical route of mixed methods, as in Chapter 7 of [START_REF] Boffi | Mixed finite element methods and applications[END_REF]. We shall rather introduce the additional unknowns only where they are needed, i.e. in the vicinity of boundary Γ.

More specifically, let us assume that Ω lies inside a simply shaped domain O (typically a box in R d ) and introduce a quasi-uniform simplicial mesh T O h on O (the background mesh). Let T h be a submesh of T O h obtained by getting rid of mesh elements lying entirely outside Ω (the definition of T h will be slightly changed afterwords). Denote by Ω h the domain covered by mesh T h (Ω h only slightly larger than Ω) and by Ω Γ h the domain covered by mesh elements of T h cut by Γ (a narrow strip of width ∼ h around Γ). Assume that the right-hand side f is actually well defined on Ω h and imagine for the moment that the solution u of eq. ( 1) can be extended to a function on Ω h , still denoted by u, which solves the same equation, now on Ω h :

-∆u

+ u = f, in Ω h . (3)
As announced above, we now introduce an auxiliary vector-valued unknown y on Ω Γ h , setting y = -∇u there, so that u, y satisfy the dual form of the original equation

y + ∇u = 0 , div y + u = f, in Ω Γ h . (4)
This allows us to rewrite the natural boundary condition ∂u ∂n = 0 on Γ as the essential condition on y: y • n = 0 on Γ. The latter can now be imposed using the idea of multiplication by the level-set φ. To this end, we note that the outward-looking unit normal n is given on Γ by n = 1 |∇φ| ∇φ . Hence, we have

y • n = 0 on Γ if we put y • ∇φ + pφ = 0, in Ω Γ h , (5)
where p is yet another (scalar-valued) auxiliary unknown on Ω Γ h . The variable p is artificial and has no physical meaning (roughly speaking, it represents the second derivative in the normal direction on Γ, extended somehow to Ω Γ h ). The only purpose of its introduction is to avoid the terms with integral on Γ in our variational formulation.

Our finite element method, cf. ( 6) below, will be based on a variational formulation of system (3)-(5) treating eqs. ( 4)-( 5) in a least squares manner and adding a stabilization in the vein of the Ghost penalty [START_REF] Burman | Ghost penalty[END_REF][START_REF] Elfverson | A new least squares stabilized nitsche method for cut isogeometric analysis[END_REF]. As in [START_REF] Duprez | φ-FEM: a finite element method on domains defined by level-sets[END_REF], we coin our method φ-FEM in accordance with the tradition of denoting the level-sets by φ. Contrary to [START_REF] Duprez | φ-FEM: a finite element method on domains defined by level-sets[END_REF], we need here additional finite element unknowns discretizing y and p on Ω Γ h . Since, the latter represents only a small portion of the whole computational domain Ω h , the extra cost induced by these unknowns is negligible as h → 0. We want to emphasize that the reformulation (3)-( 5) is very formal and will serve only as a motivation for our discrete scheme [START_REF] Burman | Ghost penalty[END_REF]. The system (3)-( 5) itself is clearly over-determined and may well be ill-posed (the "boundary" conditions hidden in [START_REF] Boiveau | Fictitious domain method with boundary value correction using penalty-free Nitsche method[END_REF] are actually not on the boundary of domain Ω h where the problem is now posed). We shall assume neither the existence of a continuous solution to (3)-( 5), nor any properties of such a solution in the theoretical analysis of our scheme, cf. Theorem 2.1.

The article is organized as follows: our φ-FEM method is presented in the next section. We also give there the assumptions on the level-set φ and on the mesh, and announce our main result: the a priori error estimate for φ-FEM in the Neumann case. We work with standard continuous P k finite elements (k ≥ 1) on a simplicial mesh and prove the optimal order h k for the error in the H 1 norm and the (slightly) suboptimal order h k+1/2 for the error in the L 2 norm. We note in passing that employing finite elements of any order is quite straightforward in our approach contrary to more traditional schemes of CutFEM type, cf. [START_REF] Boiveau | Fictitious domain method with boundary value correction using penalty-free Nitsche method[END_REF][START_REF] Lehrenfeld | High order unfitted finite element methods on level set domains using isoparametric mappings[END_REF] for a special treatment of the case k > 1. The proofs of the error estimates are the subject of Section 3. Moreover, we show in Section 4 that the associated finite element matrix has the condition number of order 1/h 2 , i.e. of the same order as that of a standard finite element method on a matching grid of comparable size. In particular, the conditioning of our method does not suffer from arbitrarily bad intersections of Γ with the mesh. Numerical illustrations are given in Section 5.

Definitions, assumptions, description of φ-FEM, and the main result

Assume Ω ⊂ O and let T O h be a quasi-uniform simplicial mesh on O with h = max T ∈T h diam T and ρ(T ) ≥ βh for all T ∈ T O h with the mesh regularity parameter β > 0 fixed once for all (here ρ(T ) is the radius of the largest ball inscribed in T ). Fix integers k, l ≥ 1 and let φ h be the FE interpolation of φ on T O h by the usual continuous finite elements of degree l. 2 Let Γ h := {φ h = 0} and introduce the computational mesh T h (approximately) covering Ω and the auxiliary mesh T Γ h covering Γ h :

T h = {T ∈ T O h : T ∩ {φ h < 0} = ∅} and Ω h = (∪ T ∈T h T ) • , T Γ h = {T ∈ T h : T ∩ Γ h = ∅} and Ω Γ h = (∪ T ∈T Γ h T ) • .
We shall also denote by Ω i h = Ω h \ Ω Γ h the domain of mesh elements completely inside Ω and set Γ i h = ∂Ω i h . We now introduce the finite element spaces

V (k) h = {v h ∈ H 1 (Ω h ) : v h | T ∈ P k (T ) ∀T ∈ T h }, Z (k) h = {z h ∈ H 1 (Ω Γ h ) d : z h | T ∈ P k (T ) d ∀T ∈ T Γ h }, Q (k) h = {q h ∈ L 2 (Ω Γ h ) : q h | T ∈ P k-1 (T ) ∀T ∈ T Γ h }, W (k) h = V (k) h × Z (k) h × Q (k) h and the finite element problem: Find (u h , y h , p h ) ∈ W (k) h such that a h (u h , y h , p h ; v h , z h , q h ) = Ω h f v h + γ div Ω Γ h f (div z h + v h ), (6)
for all (v h , z h , q h ) ∈ W (k)

h , where

a h (u, y, p; v, z, q) = Ω h ∇u • ∇v + Ω h uv + ∂Ω h y • nv + γ div Ω Γ h (div y + u)(div z + v) + γ u Ω Γ h (y + ∇u) • (z + ∇v) + γ p h 2 Ω Γ h (y • ∇φ h + 1 h pφ h )(z • ∇φ h + 1 h qφ h ) + σh Γ i h ∂u ∂n ∂v ∂n
with some positive numbers γ div , γ u , γ p , and σ properly chosen in a manner independent of h. We have assumed here that f is well defined on Ω h , rather than on Ω only. The finite element problem ( 6) is inspired by (3)-( 5). The first line in the definition of a h comes from multiplying (3) by a test function v, integrating by parts [START_REF] Boffi | Mixed finite element methods and applications[END_REF]. Equations ( 4)-( 5) are then added in least squares manner, introducing the test functions z and q corresponding to y and p respectively (a similar idea is also used in the construction of CutIGA method [START_REF] Elfverson | A new least squares stabilized nitsche method for cut isogeometric analysis[END_REF]). Note that we replace p by 1 h p in the term stemming from (5). This rescaling does not affect the discretization of u (which is the only quantity that interests us) and will be crucial to control the conditioning of the method. Finally, the terms multiplied by σh is the Ghost penalty from [START_REF] Burman | Ghost penalty[END_REF] (we need to penalize the jumps only on Γ i h because some continuity of ∇u h on the facets inside Ω Γ h is already enforced by assimilating ∇u h to y h which is continuous). We now recall some technical assumptions on the domain and the mesh, the same as in [START_REF] Lozinski | CutFEM without cutting the mesh cells: a new way to impose Dirichlet and Neumann boundary conditions on unfitted meshes[END_REF][START_REF] Duprez | φ-FEM: a finite element method on domains defined by level-sets[END_REF]. These assumptions hold true for smooth domains and sufficiently refined meshes. Assumption 1. There exists a neighborhood of Γ, a domain Ω Γ , which can be covered by open sets O i , i = 1, . . . , I and one can introduce on every O i local coordinates ξ 1 , . . . , ξ d with ξ d = φ such that all the partial derivatives ∂ α ξ/∂x α and ∂ α x/∂ξ α up to order k + 1 are bounded by some C 0 > 0. Thus, φ is of class C k+2 on Ω Γ . Moreover, |∇φ| ≥ m on Ω Γ with some m > 0.

Ω h ∇u • ∇v + Ω h uv - ∂Ω h ∇u • nv = Ω h f v and noting that -∇u • n = y • n on ∂Ω h by

Assumption 2. Ω Γ

h ⊂ Ω Γ and |∇φ h | ≥ m 2 on all the mesh elements of Ω Γ h .

Assumption 3. The approximate boundary Γ h can be covered by element patches {Π k } k=1,...,NΠ having the following properties:

• Each Π k is composed of a mesh element T k lying inside Ω and some elements cut by Γ, more precisely andΠ Γ k contains at most M mesh elements; • Each mesh element in a patch Π k shares a facet with at least one other mesh element in the same patch (so that Π k is a connected set). In particular, T k shares a facet F k with an element in Π Γ k ;

Π k = T k ∪ Π Γ k where T k ∈ T h , T k ⊂ Ω, Π Γ k ⊂ T Γ h ,
• T Γ h = ∪ NΠ k=1 Π Γ k and Γ i h = ∪ NΠ k=1 F k ; • Π k and Π l are disjoint if k = l.
Assumption 3 prevents strong oscillations of Γ on the length scale h. It can be reformulated by saying that each cut element T ∈ T Γ h can be connected to an uncut element T ⊂ Ω i h by a path consisting of a small number of mesh elements adjacent to one another; see [START_REF] Lozinski | CutFEM without cutting the mesh cells: a new way to impose Dirichlet and Neumann boundary conditions on unfitted meshes[END_REF] for a more detailed discussion and an illustration (Fig. 2). Theorem 2.1. Suppose that Assumptions 1-3 hold true, l ≥ k + 1, Ω ⊂ Ω h and f ∈ H k (Ω h ). Let u ∈ H k+2 (Ω) be the solution to (1) and (u h , y h , p h ) ∈ W (k) h be the solution to [START_REF] Burman | Ghost penalty[END_REF]. Provided γ div , γ u , γ p , σ are sufficiently big, it holds

|u -u h | 1,Ω ≤ Ch k f k,Ω h and u -u h 0,Ω ≤ Ch k+1/2 f k,Ω h (7)
with C > 0 depending on the constants in Assumptions 1, 3 (and thus on the norm of φ in C k+2 ), on the mesh regularity, on the polynomial degrees k and l, and on Ω, but independent of h, f , and u.

Remark 1 ((Condition Ω ⊂ Ω h )). The assumptions of Theorem 2.1 include Ω ⊂ Ω h . Note that one would automatically have Ω ⊂ Ω h , were Ω h defined as the set of mesh cells having a non empty intersection with Ω = {φ < 0}. However, Ω h is based on the intersections with {φ h < 0} which can result in some rare situation where tiny portions of Ω lie outside Ω h . In such a case, the a priori estimates [START_REF] Burman | Cutfem: discretizing geometry and partial differential equations[END_REF]will control the error only on Ω ∩ Ω h .

Remark 2 ((non-homogeneous Neumann and Robin conditions))

. We can also treat the case of more general boundary conditions:

(i) non-homogeneous Neumann boundary conditions ∂u ∂n = g on Γ by adding the term

- γ p h 2 Ω Γ h g|∇φ h |(z h • ∇φ h + 1 h q h φ h )
in the right-hand side of (6) where g ∈ H k+1 (Ω Γ h ) is lifting of g from Γ to a vicinity of Γ.

(ii) Robin boundary condition ∂u ∂n + αu = g on Γ (α ∈ R) by replacing the penultimate term in a h by

γ p h 2 Ω Γ h (y • ∇φ h -|∇φ h |αu + 1 h pφ h )(z • ∇φ h -|∇φ h |αv + 1 h qφ h )
and by adding the term

- γ p h 2 Ω Γ h g|∇φ h |(z h • ∇φ h -|∇φ h |αv + 1 h q h φ h )
in the right-hand side of (6) where g ∈ H k+1 (Ω Γ h ) is defined as before.

Theorem 2.1 remains valid, adding g k+1,Ω Γ h to f k,Ω h in [START_REF] Burman | Cutfem: discretizing geometry and partial differential equations[END_REF]. This framework will be used in first test case of the numerical simulations performed in Section 5: Fig. 2345678for (i) and Fig. 9 for (ii).

Proof of the a priori error estimates

From now on, we shall use the letter C for positive constants (which can vary from one line to another) that depend only on the regularity of the mesh and on the constants in Assumptions 1-3.

We shall begin with some technical results, mostly adapted from [START_REF] Lozinski | CutFEM without cutting the mesh cells: a new way to impose Dirichlet and Neumann boundary conditions on unfitted meshes[END_REF] and [START_REF] Duprez | φ-FEM: a finite element method on domains defined by level-sets[END_REF] to be used later in the proofs of the coercivity of a h (Section 3.2) and the a priori error estimates (Sections 3.3 and 3.4).

Technical lemmas

We recall first a lemma from [START_REF] Duprez | φ-FEM: a finite element method on domains defined by level-sets[END_REF]: Lemma 3.1. Let T be a triangle/tetrahedron, E one of its sides and p a polynomial on T such that p = a on E for some a ∈ R, ∂p ∂n = 0 on E, and ∆p = 0 on T . Then p = a on T .

We now adapt a lemma from [START_REF] Lozinski | CutFEM without cutting the mesh cells: a new way to impose Dirichlet and Neumann boundary conditions on unfitted meshes[END_REF]: Lemma 3.2. Let B h be the strip between ∂Ω h and Γ h . For any β > 0, there exist 0 < α < 1 and δ > 0 depending only on the mesh regularity and geometrical assumptions such that, for all

v h ∈ V (k) h , z h ∈ Z (k) h B h z h • ∇v h ≤ α|v h | 2 1,Ω h +δ z h +∇v h 2 0,Ω Γ h +βh ∂v h ∂n 2 0,Γ i h +βh 2 div z h +v h 2 0,Ω Γ h +βh 2 v h 2 0,Ω Γ h . (8)
Proof. The boundary Γ can be covered by element patches {Π k } k=1,...,NΠ as in Assumption 3. Choose any β > 0 and consider α := max

Π k ,(z h ,v h ) =(0,0) F (Π k , z h , v h ) (9)
with

F (Π k , z h , v h ) = z h 0,Π Γ k |v h | 1,Π Γ k -β z h + ∇v h 2 0,Π Γ k -βh ∂v h ∂n 2 0,F k -β 2 h 2 div z h 2 0,Π Γ k 1 2 z h 2 0,Π Γ k + 1 2 |v h | 2 1,Π k ,
where the maximum is taken over all the possible configurations of a patch Π k allowed by the mesh regularity and over all

v h ∈ V (k) h and z h ∈ Z (k) h restricted to Π k . Note that F (Π k , z h , v h ) is invariant under the scaling transformation x → 1 h x, v h → 1 h v h , z h → z h .
We can thus assume h = 1 when computing the maximum in [START_REF] Duprez | φ-FEM: a finite element method on domains defined by level-sets[END_REF]. Moreover,

F (Π k , z h , v h ) is homogeneous with respect to v h , z h , i.e. F (Π k , z h , v h ) = F (Π k , µz h ,
µv h ) for any µ = 0. Thus, the maximum in ( 9) is indeed attained since it can be taken over a closed bounded set in a finite dimensional space (all the admissible patches on a mesh with

h = 1 and all v h , z h such that |v h | 2 1,Π k + z h 2 0,Π Γ k = 1
). Clearly, α ≤ 1. Supposing α = 1 leads to a contradiction. Indeed, if α = 1, we can then take Π k , v h , z h yielding this maximum (in particular,

|v h | 2 1,Π k + z h 2 0,Π Γ k > 0). We observe then 1 2 |v h | 2 1,Π k -z h 0,Π Γ k |v h | 1,Π Γ k + 1 2 z h 2 0,Π Γ k + β z h + ∇v h 2 0,Π Γ k + βh ∂v h ∂n 2 0,F k + β 2 h 2 div z h 2 0,Π Γ k = 0
and consequently (recall

|v h | 2 1,Π k = |v h | 2 1,T k + |v h | 2 1,Π Γ k ) 1 2 |v h | 2 1,T k + β z h + ∇v h 2 0,Π Γ k + βh ∂v h ∂n 2 0,F k + β 2 h 2 div z h 2 0,Π Γ k = 0. ( 10
)
This implies |v h | 1,T k = 0 so that v h = const on T k . Moreover, z h + ∇v h 0,Π Γ k = 0 so that ∇v h = -z h on Π Γ k , hence ∇v h is continuous on Π Γ k and ∆v h = 0 on Π Γ k since div z h = 0 there. The jump ∂v h ∂n vanishes also on the facet F k separating T k from Π Γ
k , as implied directly by [START_REF] Elfverson | A new least squares stabilized nitsche method for cut isogeometric analysis[END_REF]. Combining these observations with Lemma 3.1, starting from T k and its neighbor in Π Γ k and then propagating to other elements of Π Γ k , we see that v h = const on the whole Π k . We have thus

∇v h = 0 on Π k and z h = 0 on Π Γ k , which is in contradiction with |v h | 2 1,Π k + z h 2 0,Π Γ k > 0. Thus α < 1 and z h 0,Π Γ k |v h | 1,Π Γ k ≤ α 2 z h 2 0,Π Γ k + α 2 |v h | 2 1,Π k + β z h + ∇v h 2 0,Π Γ k + βh ∂v h ∂n 2 0,∂T k ∩∂Π Γ k + β 2 h 2 div z h 2 0,Π Γ k
for all v h , z h and all admissible patches Π k . We now observe

B h z h • ∇v h ≤ k B h ∩Π Γ k z h • ∇v h ≤ k z h 0,Π Γ k |v h | 1,Π Γ k ≤ α 2 z h 2 0,Ω Γ h + α 2 |v h | 2 1,Ω h + β z h + ∇v h 2 0,Ω Γ h + βh ∂v h ∂n 2 0,Γ i h + β 2 h 2 div z h 2 0,Ω Γ h .
We now use the Young inequality with any ε > 0 to obtain

z h 2 0,Ω Γ h = z h +∇v h 2 0,Ω Γ h + ∇v h 2 0,Ω Γ h -2(z h +∇v h , ∇v h ) 0,Ω Γ h ≤ 1 + 1 ε z h +∇v h 2 0,Ω Γ h +(1+ε)|v h | 2 1,Ω h , which leads to B h z h • ∇v h ≤ α 1 + ε 2 |v h | 2 1,Ω h + β + α 2 + α 2ε z h + ∇v h 2 0,Ω Γ h + βh ∂v h ∂n 2 0,Γ i h + βh 2 div z h 2 0,Ω Γ h .
Taking ε sufficiently small, redefining α as α 1 + ε 2 and putting δ = β + α 2 + α 2ε we obtain

B h z h • ∇v h ≤ α|v h | 2 1,Ω h + δ z h + ∇v h 2 0,Ω Γ h + βh ∂v h ∂n 2 0,Γ i h + βh 2 div z h 2 0,Ω Γ h .
This leads to (8) by the triangle inequality div

z h 0,Ω Γ h ≤ div z h + v h 0,Ω Γ h + v h 0,Ω Γ h . Lemma 3.3. For all v ∈ H 1 (Ω Γ h ), v 0,Ω Γ h ≤ C √ h v 0,Γ i h + h|v| 1,Ω Γ h and for all v ∈ H 1 (Ω h \Ω), v 0,Ω h \Ω ≤ C √ h v 0,Γ + h|v| 1,Ω h \Ω .
We refer to [START_REF] Lozinski | CutFEM without cutting the mesh cells: a new way to impose Dirichlet and Neumann boundary conditions on unfitted meshes[END_REF] for the first inequality. The second one can be treated similarly.

The following lemma is borrowed from [START_REF] Duprez | φ-FEM: a finite element method on domains defined by level-sets[END_REF]. It's a partial generalization of Lemma 3.3 to derivatives of higher order. Lemma 3.4. Under Assumption 1, it holds for all v ∈ H s (Ω h ) with integer

1 ≤ s ≤ k + 1, v vanishing on Ω, v 0,Ω h \Ω ≤ Ch s v s,Ω h \Ω .
Lemma 3.5. For all piecewise polynomial (possibly discontinuous) functions

v h on T Γ h , v h 0,Γ h ≤ C √ h v h 0,Ω Γ
h with a constant C > 0 depending on the maximal degree of polynomials in v h and on the constants in Assumptions 1-3.

Proof. A scaling argument on all T ∈ T Γ

h .

Finally, we recall a Hardy-type lemma, cf. [START_REF] Duprez | φ-FEM: a finite element method on domains defined by level-sets[END_REF].

Lemma 3.6. Assume that the domain Ω Γ is a neighborhood of Γ, given by ( 2), and satisfies Assumption 1. Then, for any u ∈ H s+1 (Ω Γ ) vanishing on Γ and an integer s ∈ [0, k], it holds u φ s,Ω Γ ≤ C u s+1,Ω Γ with C > 0 depending only on the constants in Assumption 1 and on s.

Coercivity of the bilinear form a

It will be convenient to rewrite the bilinear form a h in a manner avoiding the integral on ∂Ω h . To this end, we recall that B h is the strip between ∂Ω h and Γ h and observe for any

y ∈ H 1 (B h ) d , v ∈ H 1 (B h ), q ∈ L 2 (Γ h ): ∂Ω h y • nv = ∂Ω h y • nv - Γ h 1 |∇φ h | (y • ∇φ h )v + Γ h 1 |∇φ h | (y • ∇φ h + 1 h qφ h )v = B h (v div y + y • ∇v) + Γ h 1 |∇φ h | (y • ∇φ h + 1 h qφ h )v.
Indeed, φ h = 0 on Γ h and the unit normal to Γ h , looking outward from B h , is equal to -∇φ h /|∇φ h |. Thus,

a h (u, y, p; v, z, q) = Ω h ∇u • ∇v + Ω h uv + B h (v div y + y • ∇v) + Γ h 1 |∇φ h | (y • ∇φ h + 1 h qφ h )v + γ div Ω Γ h (div y + u)(div z + v) + γ u Ω Γ h (y + ∇u) • (z + ∇v) + σh Γ i h ∂u ∂n ∂v ∂n + γ p h 2 Ω Γ h (y • ∇φ h + 1 h pφ h )(z • ∇φ h + 1 h qφ h ). ( 11
)
Proposition 1. Provided γ div , γ u , γ p , σ are sufficiently big, there exists an h-independent constant c > 0 such that

a h (v h , z h , q h ; v h , z h , q h ) ≥ c|||v h , z h , q h ||| 2 h , ∀(v h , z h , q h ) ∈ W (k) h with |||v, z, q||| 2 h = v 2 1,Ω h + div z + v 2 0,Ω Γ h + z + ∇v 2 0,Ω Γ h + h ∂v ∂n 2 0,Γ i h + 1 h 2 z • ∇φ h + 1 h qφ h 2 0,Ω Γ h .
Proof. Using the reformulation of the bilinear form a h given by ( 11), we have for all

(v h , z h , q h ) ∈ W (k) h , a h (v h , z h , q h ; v h , z h , q h ) = |v h | 2 1,Ω h + v h 2 0,Ω h + B h (v h div z h + z h • ∇v h ) + Γ h 1 |∇φ h | (z h • ∇φ h + 1 h q h φ h )v h + γ div div z h + v h 2 0,Ω Γ h + γ u z h + ∇v h 2 0,Ω Γ h + σh ∂v h ∂n 2 0,Γ i h + γ p h 2 z h • ∇φ h + 1 h q h φ h 2 0,Ω Γ h .
Since B h ⊂ Ω Γ h , we remark that the integral of v h div z h can be combined with that of v h on Ω Γ h to give

v h 2 0,Ω Γ h + B h v h div z h ≥ B h v h (div z h + v h ) ≥ -v h 0,Ω Γ h div z h + v h 0,Ω Γ h .
We also use an inverse inequality from Lemma 3.5 and the fact that 1/|∇φ h | is uniformly bounded by Assumption 2, to estimate

Γ h 1 |∇φ h | (z h • ∇φ h + 1 h q h φ h )v h ≤ C h z h • ∇φ h + 1 h q h φ h 0,Ω Γ h v h 0,Ω Γ h .
Applying the Young inequality (for any ε > 0) to the last two bounds and combining this with (8) yields

a h (v h , z h , q h ; v h , z h , q h ) ≥ (1-α)|v h | 2 1,Ω h + v h 2 0,Ω i h -(ε+βh 2 ) v h 2 0,Ω Γ h + γ div - 1 2ε -βh 2 div z h +v h 2 0,Ω Γ h + (γ u -δ) z h + ∇v h 2 0,Ω Γ h + (σ -β)h ∂v h ∂n 2 0,Γ i h + γ p h 2 - C 2 2εh 2 z h • ∇φ h + 1 h q h φ h 2 0,Ω Γ h .
To bound further from below the first 3 terms we note, using Lemma 3.3 and the trace inverse inequality,

v h 2 0,Ω Γ h ≤ C(h v h 2 0,Γ i h + h 2 |v h | 2 1,Ω Γ h ) ≤ C( v h 2 0,Ω i h + h 2 |v h | 2 1,Ω h )
so that, introducing any κ ≥ 0 and observing h ≤ h 0 := diam(Ω),

(1 -α)|v h | 2 1,Ω h + v h 2 0,Ω i h -(ε + βh 2 ) v h 2 0,Ω Γ h ≥ (1 -α)|v h | 2 1,Ω h + v h 2 0,Ω i h + κ v h 2 0,Ω Γ h -(ε + βh 2 0 + κ) v h 2 0,Ω Γ h ≥ (1 -α -C(ε + βh 2 0 + κ)h 2 0 )|v h | 2 1,Ω h + (1 -C(ε + βh 2 0 + κ)) v h 2 0,Ω i h + κ v h 2 0,Ω Γ h .
Taking ε, κ, β sufficiently small and γ u , γ p , γ div sufficiently big, gives the announced lower bound for a h (v h , z h , q h ; v h , z h , q h ).

3.3 Proof of the H 1 error estimate in Theorem 2.1

Under the Theorem's assumptions, the solution to (1) is indeed in H k+2 (Ω) and it can be extended to a function ũ ∈ H k+2 (Ω h ) such that ũ = u on Ω and

ũ k+2,Ω h ≤ C( f k,Ω + g k+1/2,Γ ) ≤ f k,Ω . (12)
Introduce y = -∇ũ and p = -h φ y • ∇φ on Ω Γ h . Then, y ∈ H k+1 (Ω Γ h ) and p ∈ H k (Ω Γ h ) by Lemma 3.6. Moreover,

y k+1,Ω Γ h ≤ C ũ k+2,Ω h ≤ C f k,Ω and p k,Ω Γ h ≤ Ch y k+1,Ω Γ h ≤ Ch f k,Ω . (13)
Clearly, ũ, y, p satisfy

a h (ũ, y, p; v h , z h , q h ) = Ω h f v h + γ div Ω Γ h f (div z h + v h ) + γ p h 2 Ω Γ h (y • ∇φ h + 1 h pφ h )(z h • ∇φ h + 1 h q h φ h ), ∀(v h , z h , q h ) ∈ W (k) h
with f := -∆ũ + ũ. It entails a Galerkin orthogonality relation

a h (ũ -u h , y -y h , p -p h ; v h , z h , q h ) = Ω h ( f -f )v h + γ div Ω Γ h ( f -f )(div z h + v h ) + γ p h 2 Ω Γ h (y • ∇φ h + 1 h pφ h )(z h • ∇φ h + 1 h q h φ h ), ∀(v h , z h , q h ) ∈ W (k) h . ( 14
)
Introducing the standard nodal interpolation I h or, if necessary, a Clément interpolation (recall that p is only in

H 1 (Ω Γ h ) if k = 1
), we then have by Proposition 1,

c u h -I h ũ, y h -I h y, p h -I h p h ≤ sup (v h ,z h ,q h )∈W (k) h a h (u h -I h ũ, y h -I h y, p h -I h p; v h , z h , q h ) v h , z h , q h h ≤ sup (v h ,z h ,q h )∈W (k) h I -II -III v h , z h , q h h ,
where

I = a h (e u , e y , e p ; v h , z h , q h ), II = Ω h ( f -f )v h + γ div Ω Γ h ( f -f )(div z h + v h ), III = γ p h 2 Ω Γ h (y • ∇φ h + 1 h pφ h )(z h • ∇φ h + 1 h q h φ h ),
with e u = ũ -I h ũ, e y = y -I h ỹ and e p = p -I h p. We now estimate each term separately. Recalling [START_REF] Glowinski | Error estimates for fictitious domain/penalty/finite element methods[END_REF], we have

I ≤ e u 1,Ω h v h 1,Ω h + dive y 0,B h v h 0,B h + e y 0,B h |v h | 1,B h + 1 |∇φ h | (e y • ∇φ h + 1 h e p φ h ) 0,Γ h v h 0,Γ h + γ div div e y + e u 0,Ω Γ h div z h + u h 0,Ω Γ h + γ u e y + ∇e u 0,Ω Γ h z h + ∇v h 0,Ω Γ h + σh ∂e u ∂n 0,Γ i h ∂v h ∂n 0,Γ i h + γ p h 2 e y • ∇φ h + 1 h e p φ h 0,Ω Γ h z h • ∇φ h + 1 h q h φ h 0,Ω Γ h .
Applying Lemma 3.5 to the L 2 norms on Γ h , recalling that 1/|∇φ h | is uniformly bounded on Ω Γ h (cf. Assumption 2), and recombining the terms, we get

I ≤ C e u 2 1,Ω h + e y 2 1,Ω Γ h + h ∂e u ∂n 2 0,Γ i h + 1 h 2 e y • ∇φ h + 1 h e p φ h 2 0,Ω Γ h 1/2 v h , z h , q h h .
The usual interpolation estimates give

e u 2 1,Ω h + e y 2 1,Ω Γ h + h ∂eu ∂n 2 0,Γ i h ≤ Ch 2k ( ũ 2 k+1,Ω h + y 2 k+1,Ω Γ h ) .
Moreover, recalling that |∇φ h | and

1 h |φ h | are uniformly bounded on Ω Γ h , we get 1 h 2 e y • ∇φ h + 1 h e p φ h 2 0,Ω Γ h ≤ C h 2 e y 2 0,Ω Γ h + e p 2 0,Ω Γ h ≤ Ch 2k (|y| 2 k+1,Ω Γ h + 1 h 2 |p| 2 k,Ω Γ h ).
Thus, by regularity estimates [START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF], I ≤ Ch k f k,Ω v h , z h , q h h . We now estimate the second term

|II| ≤ C( f -f 0,Ω h v h 0,Ω h + f -f 0,Ω Γ h div z h + v h 0,Ω Γ h ) ≤ C f -f 0,Ω h v h , z h , q h h ≤ Ch k f k,Ω∪Ω h v h , z h , q h h .
Indeed, thanks to Lemma 3.4 and

f = f on Ω, f -f 0,Ω h = f -f 0,Ω h \Ω ≤ Ch k f -f k,Ω h \Ω ≤ Ch k f k,Ω∪Ω h . (15)
Finally,

|III| ≤ C h y • ∇φ h + 1 h pφ h 0,Ω Γ h v h , z h , q h h and, recalling y • ∇φ + 1 h pφ = 0 on Ω Γ h , 1 h y • ∇φ h + 1 h pφ h 0,Ω Γ h = 1 h y • ∇(φ h -φ) + 1 h p(φ h -φ) 0,Ω Γ h ≤ 1 h y 0,Ω Γ h ∇(φ h -φ) ∞ + 1 h 2 p 0,Ω Γ h φ h -φ ∞ ≤ Ch k ( y 0,Ω Γ h + p 0,Ω Γ h ) ≤ Ch k f k,Ω by
regularity estimates [START_REF] Hecht | New development in FreeFem++[END_REF]. Note that the optimal order is achieved here since φ is assumed of regularity C k+2 and it is approximated by finite elements of degree at least k + 1.

Combining the estimate for the terms I-III leads to

u h -I h ũ, y h -I h y, p h -I h p h ≤ Ch k f k,Ω∪Ω h ,
so that, by the triangle inequality together with interpolation estimate, we get

u h -ũ, y h -y, p h -p h ≤ Ch k f k,Ω∪Ω h . (16)
This implies the announced H 1 error estimate for u -u h .

3.4 Proof of the L 2 error estimate in Theorem 2.1

Since Ω ⊂ Ω h , we can introduce w : Ω → R such that

-∆w + w = u -u h in Ω, ∂w ∂n = 0 on Γ.
By elliptic regularity, w 2,Ω ≤ C u -u h 0,Ω . Let w be an extension of w from Ω to Ω h preserving the H 2 norm estimate and set w h = I h w. We observe

u -u h 2 0,Ω = Ω ∇(u -u h ) • ∇(w -w h ) + Ω (u -u h )(w -w h ) + Ω ∇(u -u h ) • ∇w h + Ω (u -u h )w h ≤ Ch k+1 f k,Ω h | w| 2,Ω h + Ω ∇(u -u h ) • ∇w h + Ω (u -u h )w h
by the already proven H 1 error estimate and interpolation estimates for I h w (recall also Ω ⊂ Ω h ). Taking v h = w h , z h = 0 and q h = 0 in the Galerkin orthogonality relation ( 14), we obtain, thanks to [START_REF] Glowinski | Error estimates for fictitious domain/penalty/finite element methods[END_REF],

Ω h ∇(ũ -u h ) • ∇w h + Ω h (ũ -u h )w h + B h (w h div(y -y h ) + (y -y h ) • ∇w h ) + Γ h 1 |∇φ h | ((y -y h ) • ∇φ h + 1 h (p -p h )φ h )w h + γ div Ω Γ h (div(y -y h ) + ũ -u h )w h + γ u Ω Γ h ((y -y h ) + ∇(ũ -u h )) • ∇w h + σh Γ i h ∂(ũ -u h ) ∂n ∂w h ∂n = (1 + γ div ) Ω h ( f -f )w h .
Using the last relation in the bound for u -u h 2 0,Ω , we can further bound it as

u -u h 2 0,Ω Ch k+1 f k,Ω h | w| 2,Ω h + Ω h \Ω ∇(ũ -u h ) • ∇w h + Ω h \Ω (ũ -u h )w h + B h (w h div(y -y h ) + (y -y h ) • ∇w h ) + Γ h 1 |∇φ h | ((y -y h ) • ∇φ h + 1 h (p -p h )φ h )w h + γ div Ω Γ h (div(y -y h ) + ũ -u h )w h + γ u Ω Γ h ((y -y h ) + ∇(ũ -u h )) • ∇w h + σh Γ i h ∂(ũ -u h ) ∂n ∂w h ∂n + (1 + γ div ) Ω h ( f -f )w h Ch k+1 f k,Ω h | w| 2,Ω h + C ũ -u h , y -y h , p -p h h × w h 1,Ω h \Ω + w h 1,Ω Γ h + h w h 0,Γ h + √ h [∇w h ] 0,Γ i h + C f -f 0,Ω h \Ω w h 1,Ω h \Ω .
It remains to bound different norms of w h featuring in the estimate above. By Lemma 3.3 and interpolation estimates

w h 0,Ω h \Ω ≤ w-I h w 0,Ω h \Ω + w 0,Ω h \Ω ≤ Ch 2 | w| 2,Ω h \Ω +C √ h w 0,Γ + h| w| 1,Ω h \Ω ≤ C √ h w 2,Ω h .
Similarly,

∇w h 0,Ω h \Ω ≤ ∇( w -I h w) 0,Ω h \Ω + ∇ w 0,Ω h \Ω ≤ Ch| w| 2,Ω h \Ω + C √ h ∇ w 0,Γ + h|∇ w| 1,Ω h \Ω ≤ C √ h w 2,Ω h .
Analogous estimates also hold for w h 1,Ω Γ h . Moreover, by interpolation estimates,

[∇w h ] 0,Γ i h = [∇( w -I h w)] 0,Γ i h C √ h| w| 2,Ω h
and, by Lemma 3.5,

h w h 0,Γ h ≤ C √ h w h 0,Ω Γ h ≤ C √ h w 2,Ω h .
Hence,

u -u h 2 0,Ω ≤ Ch k+1 f k,Ω h | w| 2,Ω h + C √ h( ũ -u h , y -y h , p -p h h + f -f 0,Ω h \Ω ) w 2,Ω h .
This implies, by ( 15) and ( 16), u -

u h 2 0,Ω ≤ Ch k+ 1 2 f k,Ω h w 2,Ω h , which entails the announced error estimate in L 2 (Ω) since w 2,Ω h ≤ C u -u h 0,Ω .

Conditioning

We are now going to prove that the condition number of the finite element matrix associated to the bilinear form a h is of order 1/h 2 . Theorem 4.1. Under Assumptions 1-3 and recalling that the mesh T h is quasi-uniform, the condition number defined by κ(A) := A 2 A -1 2 of the matrix A associated to the bilinear form a h on W

(k) h satisfies κ(A) ≤ Ch -2 .
Here, • 2 stands for the matrix norm associated to the vector 2

-norm | • | 2 .
Proof. The proof is divided into 4 steps:

Step 1. We shall prove for all q h ∈ Q

(k) h q h φ h 0,Ω Γ h ≥ Ch q h 0,Ω Γ h . ( 17
)
We have min

T,q h =0,φ h =0 q h φ h 0,T h T q h 0,T ∇φ h ∞,T C, (18)
where the minimum is taken over all simplexes T with h T = diam (T ) satisfying the regularity assumptions and all polynomials q h of degree k and φ h of degree l, with φ h vanishing at at least one point on T . Note that this excludes ∇φ h ∞,T = 0 because φ h would then vanish identically on T . The minimum in ( 18) is indeed attained since, by homogeneity, it can be taken over the compact set q h 0,T = ∇φ h ∞,T = 1 and simplexes with h T = 1. Hence, ( 18) is valid with some C > 0. Applying [START_REF] Moës | Imposing dirichlet boundary conditions in the extended finite element method[END_REF] on any mesh element

T ∈ T Γ h to any q h ∈ Q (k)
h and φ h approximation to φ satisfying Assumption 2 leads to q h φ h 0,T Ch T m 2 q h 0,T . Taking the square on both sides and summing over all T ∈ T Γ h yields (17).

Step 2. We shall prove for all (v h , z h , q h ) ∈ W

(k) h a h (v h , z h , q h ; v h , z h , q h ) ≥ c v h , z h , q h 2 0 (19) with v h , z h , q h 2 0 = v h 2 0,Ω h + z h 2 0,Ω Γ h + q h 2 0,Ω Γ h . Indeed, by Lemma 1, a h (v h , z h , q h ; v h , z h , q h ) ≥ c|||v h , z h , q h ||| 2 h ≥ c( v h 2 1,Ω h + z h + ∇v h 2 0,Ω Γ h + z h • ∇φ h + 1 h q h φ h 2 0,Ω Γ h ).
We have assumed here (without loss of generality) h ≤ 1. By Young's inequality with any 1 ∈ (0, 1),

z h + ∇v h 2 0,Ω Γ h = z h 2 0,Ω Γ h + ∇v h 2 0,Ω Γ h + 2(z h , ∇v h ) 0,Ω Γ h ≥ (1 -1 ) z h 2 0,Ω Γ h - 1 -1 1 ∇v h 2 0,Ω Γ h . (20)
Similarly, for any 2 ∈ (0, 1), using that ∇φ h is uniformly bounded,

z h • ∇φ h + 1 h q h φ h 2 0,Ω Γ h ≥ 1 -2 h 2 φ h q h 2 0,Ω Γ h -C 1 -2 2 z h 2 0,Ω Γ h . (21)
Thus, combining ( 20), ( 21) and ( 17),

a h (v h , z h , q h ; v h , z h , q h ) ≥ c 1 - 1 -1 1 v h 2 1,Ω h + 1 -1 -C 1 -2 2 z h 2 0,Ω Γ h + C(1 -2 ) q h 2 0,Ω Γ h .
Taking 1 , 2 close to 1, we get [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF].

Step 3. We shall prove for all (v h , z h , q h ) ∈ W

(k) h a h (v h , z h , q h ; v h , z h , q h ) ≤ C h 2 v h , z h , q h 2 0 . ( 22
)
By definition of a h [START_REF] Glowinski | Error estimates for fictitious domain/penalty/finite element methods[END_REF] and Cauchy-Schwarz inequality,

a h (v h , z h , q h ; v h , z h , q h ) ≤ C v h 2 1,Ω h + z h 2 1,Ω Γ h + h ∂v h ∂n 2 0,Γ i h + 1 h 2 z h • ∇φ h 2 0,Ω Γ h + 1 h 4 q h φ h 2 0,Ω Γ h .
This leads to (22) thanks to inverse inequalities and to the fact that both ∇φ h and 1 h φ h are uniformly bounded on Ω Γ h .

Step 4. We combine ( 19) and ( 22), and observe that the norm v h , z h , q h 0 is equivalent to the 2-norm of the vector representing (v h , z h , q h ). This leads to the desired result as at the end of the proof of Theorem 4.1 from [START_REF] Duprez | φ-FEM: a finite element method on domains defined by level-sets[END_REF].

Numerical simulations

In this section, we illustrate φ-FEM on three different test cases, cf. Fig. 1, exploring the errors with respect to exact "manufactured" solutions. The numerical results for the 1st test case (in 2D) confirm the predicted theoretical estimates (in fact, better than theoretically predicted convergence rate is observed for the L 2 error). In the 2nd test case (also in 2D), we show that the optimal convergence is recovered even when the level-set function φ is less regular than assumed by the theory. Our method is also compared with CutFEM [START_REF] Burman | Fictitious domain finite element methods using cut elements: Ii. A stabilized Nitsche method[END_REF] in the last case. Finally, a 3D example is given in the 3rd test case.

In the first test case, we will treat some examples with non-homogeneous Neumann condition (Fig. 2345678) and Robin condition (Fig. 9) thanks to the modification of the scheme given in Remark 2. In the two last test cases, we will consider homogeneous Neumann conditions.

The surrounding domains O are always chosen as boxes aligned with the Cartesian coordinates and the background meshes T O h are obtained from uniform Cartesian grids, dividing the cells into the simplexes (semi-cross meshes in 2D) . We always use the numerical quadrature of a high enough order so that all the integrals in ( 6) are computed exactly.

We have implemented φ-FEM both in FreeFEM [START_REF] Hecht | New development in FreeFem++[END_REF] and in multiphenics [START_REF] Ballarin | multiphenics[END_REF]. Both implementations give the same results in our test cases and we present here only those obtained with FreeFEM. The implementation scripts can be consulted on GitHub. 3Remark 3. Our method (6) features "mixed" terms, such as

γ 1 Ω Γ h (y h +∇u h )•(z h +∇v h ) involving u h , v h
defined on mesh T h and y h , z h defined on mesh T Γ h , a submesh of T h . Such integrals cannot be implemented in the current version of FEniCS since it requires all the finite elements involved in a problem to be defined on the same mesh. This is why we have turned to the multiphenics library, a spin-off of FEniCS, that does not have such a restriction. On the other hand, FreeFEM features interpolations between meshes in a user-friendly manner. However, we have discovered that a straightforward implementation of ( 6) in FreeFEM involving an implicit interpolation from T h to T Γ h can lead to some spurious oscillations in the error curves. Much better results (reported below) are obtained if we introduce explicitly the interpolation matrix from V (k) h to its restriction on T Γ h , using the FreeFEM function interpolate, or, equivalently but more efficiently, the renumbering between the degrees of freedom of these two finite element spaces, , using the FreeFEM function restrict.

1st test case

Domain Ω (see Fig. 1 left) is defined by the level-set function φ given in the polar coordinates (r, θ) by φ(r, θ) = r 4 (5 + 3 sin(7(θ -θ 0 ) + 7π/36))/2 -R 4 , where R = 0.47 and θ 0 ∈ [0, 2π). The surrounding domain O is fixed to (-0.5, 0.5) 2 . Varying the angle θ 0 results in a rotation of Ω, so that the boundary Γ cuts the triangles of the background mesh in a different manner, creating sometimes the "dangerous" situations when certain mesh triangles of T h have only a tiny portion inside the physical domain Ω.

L 2 relative error u -u h 0,Ω i h / u 0,Ω i h ; Right: H 1 relative error |u -u h | 1,Ω i h /|u| 1,Ω i h .
We use φ-FEM to solve numerically Poisson-Neumann problem [START_REF] Balay | PETSc users manual[END_REF] with non-homogeneous boundary conditions ∂u ∂n = g adjusting f and g so that the exact solution is given by u(x, y) = sin(x) exp(y). The Neumann boundary condition is extrapolated to a vicinity of Γ by g = ∇u•∇φ |∇φ| + uφ, cf. Remark 2. The addition of uφ here does not perturb g on Γ. Its purpose is to mimick the real life situation where g is known on Γ only and g is some extension of g, not necessarily the natural one ∇u • ∇φ/|∇φ|.

L 2 relative error u -u h 0,Ω i h / u 0,Ω i h ; Right: H 1 relative error |u -u h | 1,Ω i h /|u| 1,Ω i h .
We report at Figs. 2 and 3 the evolution of the relative error under the mesh refinement for a fixed position of Ω (θ 0 = 0), using finite element spaces W (k) h with k = 1 (P 1 FE for u h ) and k = 2 (P 2 FE for u h ). We also try there different values of l, the degree of finite element used to approximate the level-set φ, recalling that it should be chosen as k+1 or greater. The experiments reported in these figures confirms the optimal convergence order of the method in both H 1 and L 2 norms (orders k and k + 1 respectively). The convergence order in the L 2 norm is thus better than in theory. An interesting experimental observation comes from exploring the degree l: while the lowest possible value l = k + 1 ensures indeed the optimal convergence orders, it seems advantageous to increase the degree to l = k + 2, leading to more accurate results, especially in the L 2 norm. Another series of experiments is reported at Figs. 4 and5. We explore there the errors with respect to the rotation of Ω over the background mesh (varying θ 0 ). We restrict ourselves here with finite elements degree k = 1 but compare two different values of l: l = k + 1 = 2 at Fig. 4 vs. l = k + 2 = 3 at Fig. 5. We observe again an advantage of the choice l = k + 2: the oscillations on any given background mesh become less important when increasing l and fade away under the mesh refinement in the case l = k +2 (this concerns mostly the L 2 errors; the H 1 errors are pretty much the same in both cases). The influence of the parameters σ, γ div , γ u , γ p on the accuracy of the method is explored by the numerical experiments reported at Figs. 6 and7. Although a full assessment of the role of all the 4 parameters is difficult (we have chosen, somewhat arbitrarily, two scenarios of parameter variations out of endless other possibilities), the conclusion of our numerical experiments seems clear: the method is not sensible to variation of the parameters in the wide range from 10 -6 to 10, and there is no need to take these parameters greater than 10. Finally, we report at Fig. 8 evolution of the condition number of the φ-FEM matrix under the mesh refinement and also its sensitivity with respect to the rotations of Ω. The theoretically predicted behavior of ∼ 1/h 2 is confirmed. The conditioning of the method is also found to be rather insensitive to the position of Ω over the mesh. We end this section by given an example with Robin boundary condition with α = 1 thanks to modification of main scheme presented in Remark 2. We consider the same domain Ω, level-set function φ and solution u as before. The Robin condition is extrapolated by g = ∇u•∇φ |∇φ| + αu + uφ. In Fig. 9, we report the L 2 errors and the H 1 error (left) and the condition number (right) for k = 1 and l = 3. We observe that optimal convergence order and standard condition number remain valid for our Robin formulation. 10 -2

L 2 relative error u -u h 0,Ω i h / u 0,Ω i h ; Right: H 1 relative error |u -u h | 1,Ω i h /|u|
-u h 0,Ω i h / u 0,Ω i h ; Right: H 1 relative error |u -u h | 1,Ω i h /|u|
10 -1 1 1 1 2 h u -u h 0,Ω i h / u 0,Ω i h |u -u h | 1,Ω i h /|u| 1,Ω i h 2 × 10 -1 10

2nd test case

In this test case, the domain Ω is the rectangle (-1, 1) × (-2, 2) rotated by an angle θ 0 counter-clockwise around the origin. It is defined by the level-set function φ given by φ(x, y) = Φ • Π(x, y), with Φ(x, y) = max(|x|, |y|/2) -1 and Π

x y = cos(θ 0 ) -sin(θ 0 ) sin(θ 0 ) cos(θ 0 )

x y . The surrounding domain is taken as The results are presented at Figs. 10 (left) and 11, first choosing a fixed inclination angle θ 0 = π/8, and then varying θ 0 from 0 to 2π/7. The numerical tests show again the optimal convergence of φ-FEM with P 1 finite elements in the L 2 and H 1 norms, notwithstanding the fact that the level-set function φ is less regular than assumed in our theoretical results. Note that we have used here the FE of degree l = 3 to represent the level-set, which is higher than the minimal degree k + 1 = 2 suggested by the theory. The situation is here similar to that of the tests case 1-2: the implementation using the lower degree l = 2 elements (not reported here) is also optimally convergent but turns out to be less robust than l = 3 wih respect to the placement of Ω over the mesh (higher oscillations, especially in the L 2 error, when varying θ 0 ).

O = (-R, R) 2 , with R = 1.1 √ 5, cf.
We have also compared our method with CutFEM [START_REF] Burman | Fictitious domain finite element methods using cut elements: Ii. A stabilized Nitsche method[END_REF]:

Find u h ∈ V (k) h s.t. Ω ∇u h • ∇v h + Ω u h v h + σh E∈F Γ E ∂u h ∂n ∂v h ∂n = Ω f v h + Γ gv h ∀ v h ∈ V (k) h ,
where F Γ = {E(internal facet of T h ) such that ∃T ∈ T h : T ∩ Γ = ∅ and E ∈ ∂T }.

The results are reported at Figs. 10 (right, the simulation at fixed inclination angle θ 0 ) and 12 (simulations with the rotating domain Ω). Comparing two parts of Fig. 10, we conclude that φ-FEM and CutFEM are both optimally convergent and produce very similar results. However, looking closer at Figs. 11 and12, we can point out an advantage of the φ-FEM over the CutFEM: the former seems more robust with respect to the position of Ω over the background mesh, the oscillations of the L 2 errors with rotating the domain are more pronounced for the latter method (the H 1 errors are almost the same in both cases).

3rd test case

We here take Ω ⊂ R 3 as the ball of radius R = 0.75 centered at the origin encapsulated into the box O = (-1, 1) 3 . Ω is defined by the level-set function φ(x, y, z) = x 2 + y 2 + z 2 -R 2 . Fig. 1 (right) gives an example of mesh T h for this test case. We choose the exact solution as u(x, y, z) = cos

x 2 + y 2 + z 2 .
The normal derivative ∂u ∂n turns out to be constant on Γ in this case, equal to -sin(R). Accordingly, we choose to extrapolate the Neumann boundary condition to a vicinity of Γ by g = -sin(R), cf. Remark 2. Again, we observe in Fig. 13 (left) the optimal orders of convergence for the L 2 and H 1 errors. We have done the experiments on moderate meshes only, using a direct solver (MUMPS) for the linear systems. The condition numbers are reported at Fig. 13 (right). Estimating numerically the condition numbers is tricky even on these relatively coarse meshes in 3D since their direct evaluation, by eg. the standard function of NumPy used in our 2D tests, is too expensive on some of our meshes in 3D. We have computed the condition numbers reported here by tracking the largest and smallest singular values on iterations of the GMRES method without any preconditioner, as provided by PETSc [START_REF] Balay | PETSc users manual[END_REF]. The condition numbers obtained like this can be thus under-estimated. This can explain why they seem to behave more like ∼ 1/h rather than the theoretically expected ∼ 1/h 2 . Another possible explanation is that the meshes explored here are yet too coarse to see the asymptotic regime.

Figure 1 :

 1 Figure 1: Domains and meshes considered in φ-FEM for the test case 1 (left), test case 2 (center) and test case 3 (right).

  l = 2 k = 1 and l = 3 k = 1 and l = 4

Figure 2 :

 2 Figure 2: φ-FEM for the test case 1, θ 0 = 0, σ = 0.01 and γ u = γ p = γ div = 10, k = 1 and different values of l. Left:L 2 relative error u -u h 0,Ω i h / u 0,Ω i h ; Right: H 1 relative error |u -u h | 1,Ω i h /|u| 1,Ω i h .

  and l = 3 k = 2 and l = 4

Figure 3 :

 3 Figure 3: φ-FEM for the test case 1, θ 0 = 0, σ = 0.01 γ u = γ p = γ div = 10, k = 2 and different values of l. Left:L 2 relative error u -u h 0,Ω i h / u 0,Ω i h ; Right: H 1 relative error |u -u h | 1,Ω i h /|u| 1,Ω i h .

Figure 4 :

 4 Figure 4: Sensitivity of the relative error with respect to θ 0 in φ-FEM for the test case 1, σ = 0.01 andγ u = γ p = γ div = 20, k = 1 and l = 2. Left: L 2 relative error u -u h 0,Ω i h / u 0,Ω i h ; Right: H 1 relative error |u -u h | 1,Ω i h /|u| 1,Ω i h .

Figure 5 :

 5 Figure 5: Sensitivity of the relative error with respect to θ 0 in φ-FEM for the test case 1, σ = 0.01 and γ u = γ p = γ div = 10, k = 1 and l = 3. Left:L 2 relative error u -u h 0,Ω i h / u 0,Ω i h ; Right: H 1 relative error |u -u h | 1,Ω i h /|u| 1,Ω i h .

Figure 6 :Figure 7 :

 67 Figure 6: Sensitivity of the relative error in φ-FEM with respect to σ with γ u = γ p = γ div = 10 being fixed for the test case 1, θ 0 = 0, k = 1 and l = 3. Left: L 2 relative error; Right: H 1 relative error.

Figure 8 :

 8 Figure 8: Condition number in φ-FEM for the test case 1, σ = 0.01, γ u = γ p = γ div = 10, θ 0 = 0, k = 1 and l = k + 2. Left: θ 0 = 0; Right: different values of θ 0 .

Figure 9 :

 9 Figure 9: φ-FEM for the test case 1 and Robin boundary conditions, σ = 0.01, γ u = γ p = γ div = 10, α = 1, k = 1 and l = 3. Left:L 2 relative error u -u h 0,Ω i h / u 0,Ω i h and H 1 relative error |u -u h | 1,Ω i h /|u|1,Ω i h . ; Right: Condition number.

Fig. 1

 1 middle. We use φ-FEM to solve numerically Poisson-Neumann problem[START_REF] Balay | PETSc users manual[END_REF] with the exact solution given by u(x, y) = U • Π(x, y), where U (x, y) = cos(πx) cos(πy/2).

uFigure 10 :Figure 11 :

 1011 Figure 10: L 2 and H 1 relative error for the test case 2. Left: φ-FEM with σ = 0.01, γ u = γ p = γ div = 10, k = 1 and l = 3; Right: CutFEM , θ 0 = π/8, σ = 0.01.

Figure 12 :

 12 Figure12: Sensitivity of the relative error with respect to θ 0 in CutFEM for the test case 2, σ = 0.01. Left: L 2 relative error u -u h 0,Ω i h / u 0,Ω i h ; Right:H 1 relative error |u -u h | 1,Ω i h /|u| 1,Ω i h .

Figure 13 :

 13 Figure 13: φ-FEM for the test case 3, σ = 0.01, γ u = γ p = γ div = 10, k = 1 and l = 3. Left: L 2 and H 1 relative error; Right: condition number.

  i h / u 0,Ω i h ; Right: H 1 relative error |u -u h | 1,Ω i h /|u| 1,Ω i h .

	10 -1				h 0.217				h 0.217
					0.109 0.0543	10 -1			0.109 0.0544
	10 -2				0.0271				0.0272
					0.0136				0.0136
	10 -3
	10 -4					10 -2
	0	0.2	0.4	0.6	0.8	0	0.4	0.6	0.8
			θ0				θ0

Another approach that can bypass non-standard quadrature is the Fat Boundary Method[START_REF] Maury | A fat boundary method for the poisson problem in a domain with holes[END_REF][START_REF] Bertoluzza | Analysis of the fully discrete fat boundary method[END_REF] where a non-matching global mesh is combined with a local matching mesh around the boundary. Our method, on the other hand, avoids completely the construction of boundary fitted meshes.

The integer k is the degree of finite elements which will be used to approximate the principal unknown u while φ is approximated by finite elements of degree l. We shall require l ≥ k + 1 in our convergence Theorem 2.1. Note, that we cannot set l = k unlike the Dirichlet case in[START_REF] Duprez | φ-FEM: a finite element method on domains defined by level-sets[END_REF]. This is essentially due to the fact that φ h is used here to approximate the normal on Γ in addition to approximating Γ itself.

https://github.com/michelduprez/PhiFEM-Neumann


Xinyong Wang

Christophe Fiter

Ying Tang

Laurentiu Hetel

Sampled-data Control for a Class of Linear Hyperbolic Systems via the Lyapunov-Razumikhin Technique*

This work investigates the stability for a class of linear hyperbolic systems with distributed sampled-data controllers. First, we convert the original system into an equivalent system in which the sampling induced error is modeled as a reset integrator. Then by means of an appropriate Lyapunov function coupled with the Razumikhin technique, sufficient conditions are given for the Rε -stability of the system. Finally, our results are validated by a numerical example.

I. INTRODUCTION

The application of digital computer in control system has become a general trend, which makes sampled-data control an active field of research in the past decades [START_REF] Åström | Computer-controlled systems: theory and design[END_REF], [START_REF] Chen | Optimal sampled-data control systems[END_REF], [START_REF] Kuo | Discrete-data control systems[END_REF]. Stability and control design for finite-dimentional systems have been considered in many research works: see e.g. the survey [START_REF] Hetel | Recent developments on the stability of systems with aperiodic sampling: An overview[END_REF], [START_REF] Laila | Sampled-data control of nonlinear systems[END_REF], [START_REF] Monaco | Issues on nonlinear digital control[END_REF]. Compared with the research method of finite dimensional system, the analysis and control of infinite dimensional systems is more challenging. Few results exist for sampled-data infinite dimensional system [START_REF] Logemann | Stability of infinitedimensional sampled-data systems[END_REF], [START_REF] Logemann | Generalized sampled-data stabilization of well-posed linear infinite-dimensional systems[END_REF].

In general, sampled-data systems can be analyzed using discrete-time, time-delay and Input-Output methods (see [START_REF] Hetel | Recent developments on the stability of systems with aperiodic sampling: An overview[END_REF] and references therein). For the class of partial differential equations (PDEs), using discrete-time finite-dimensional approximate models, [START_REF] Tan | Dynamic practical stabilization of sampled-data linear distributed parameter systems[END_REF] proposed a methodology for the design of sampled-data controller with practical stability guarantees. In references [START_REF] Fridman | Robust sampled-data control of a class of semilinear parabolic systems[END_REF], [START_REF] Kang | Distributed sampled-data control of Kuramoto-Sivashinsky equation[END_REF], [START_REF] Selivanov | Sampled-data relay control of diffusion pdes[END_REF], the time-delay approach has been used for the analysis of parabolic PDEs. Hold boundary feedback control in one-dimensional linear hyperbolic systems were considered in [START_REF] Karafyllis | Sampled-data boundary feedback control of 1-D linear transport pdes with non-local terms[END_REF]. In [START_REF] Espitia | Event-based control of linear hyperbolic systems of conservation laws[END_REF], [START_REF] Espitia | Stabilization of boundary controlled hyperbolic pdes via Lyapunov-based event triggered sampling and quantization[END_REF], event-triggered sampled-data control with controller on the boundaries was developed. The boundary feedback control of a 2 × 2 hyperbolic system was implemented by backstepping method in [START_REF] Davó | Stability analysis of a 2 × 2 linear hyperbolic system with a sampled-data controller via backstepping method and looped-functionals[END_REF], [START_REF] Espitia | Event-based boundary control of a linear 2 × 2 hyperbolic system via backstepping approach[END_REF].

It can be seen from the literature review that the analysis of sampled-data controller for hyperbolic PDEs is a wideopen area of research, and there are still many topics worth studying. The present paper aims at studying the distributed sampled-control for a class of hyperbolic PDEs. The idea is to generalize the Input-Output approach [START_REF] Fujioka | Stability analysis of systems with aperiodic sample-andhold devices[END_REF], [START_REF] Kao | Stability analysis of systems with uncertain time-varying delays[END_REF], [START_REF] Omran | Stability analysis of bilinear systems under aperiodic sampled-data control[END_REF] for finite dimensional systems, to the case of hyperbolic PDEs. An interconnected equivalent system consisting of a continuous-time PDE and a reset-integral operator is derived from the original system. In our previous work [START_REF] Wang | Stability Analysis for A Class of Linear Hyperbolic System of Balance Laws with Sampled-data Control[END_REF], the stability of linear hyperbolic systems with sampled-data controller has been ensured for a sufficiently small sampling period. In the present paper, new stability conditions are 1 Univ. Lille, CNRS, Centrale Lille, UMR 9189 CRIStAL, F-59000 Lille, France xinyong.wang, christophe.fiter, ying.tang, laurentiu.hetel@univ-lille.fr proposed by using Lyapunov-Razumikhin stability criteria (e.g. [START_REF] Fridman | Lyapunov-based stability analysis[END_REF], [START_REF] Gu | Survey on recent results in the stability and control of time-delay systems[END_REF]), the estimation of the maximum sampling interval for the system stability is improved.

The paper is structured as follows: Section II presents the systems and the problem under study. In Section III, we propose the equivalent remodelling of system, followed by the concrete stability analysis process. A numerical example is given to illustrate the feasibility of our method in Section IV. The paper is ended with conclusions and perspectives.

Notation: N is the set of nonnegative integers from 0 to infinity, R + is the set of nonnegative reals, R n is used to denote the set of n-dimensional Euclidean space with the norm | • |. L 2 (0, L) stands for the Hilbert space of square integrable scalar functions on (0, L) with the norm

• L 2 (0,L) , defined by τ L 2 (0,L) = L 0 |τ (x)| 2 dx.
The associated norm to Sobolev space H 1 (0, L) is defined

as τ H 1 (0,L) = L 0 |τ (x)| 2 + |τ x (x)| 2 dx. Given a functional V : H 1 ([0, L]; R n ) → R + such that L V ≤C = y ∈ H 1 ([0, L]; R n ) : V (y) ≤ C .
The notation W ≤ 0 denotes that W is symmetric and negative semidefinite. The symmetric elements are denoted by * in the symmetric matrix. The identity matrix is denoted by I and λ min (Θ) and λ max (Θ) are the minimum and maximum eigenvalues of the matrix Θ. C 0 is the space of continuous functions, whereas C 1 is the space of continuously differentiable functions. ⌈•⌉ is the ceiling function.

II. SYSTEM DESCRIPTION AND PROBLEM FORMULATION

A. System Description

We consider the following sampled-data controlled hyperbolic system (1)

         ∂ t z (t, x) + Λ∂ x z (t, x) + Γz (t, x) + u (t, x) = 0, u (t, x) = ̥z (t k , x) , ∀t ∈ [t k , t k+1 ), k ∈ N, z(t, 0) = 0, ∀t ≥ 0, z(0, x) =z 0 (x), ∀x ∈ [0, L] , (1a) (1b) (1c) (1d)
where z : [0, +∞) × [0, L] → R n , Λ = diag {λ 1 , λ 2 , ..., λ n } with λ 1 , λ 2 , ..., λ n > 0, Γ and ̥ are real n × n constants matrices. The sampling instants are defined as a sequence {t k } k∈N where

t 0 = 0, t k+1 -t k ∈ [h, h], (2)
and h, h are the given bounds of the sampling intervals satisfying h ≥ h > 0.

To address the issue under consideration, we need the compatibility condition given below: Condition 1. The initial condition z 0 (x), satisfies

z 0 (0) = 0, ∀x ∈ [0, L] . (3)
Remark 1. We explain the concept of the solution and rewrite the system (1)-( 2) as a first order system

dz(t) dt = Υz (t) + f (z (t k )) , t ∈ [t k , t k+1 ) , k ∈ N, z (0) = z 0 ,
where f (z(t k )) = -̥z(t k ), and Υ is an operator defined by Υz = -Λ∂ x z (t, x) -Γz (t, x), with domain

D(Υ) = z ∈ H 1 (0, L; R n ) z(0) = 0. (4)
A stable C 0 semigroup is produced by the operator Υ (see the proof of theorem A.1. in [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF]). Moreover, the fact is that

f k : H 1 (0, L) → H 1 (0, L) is continuously differentiable for t ∈ [t k , t k+1 ). If z 0 ∈ D(Υ)
, then in the light of Theorem 6.1.5 of [START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF], there is a classical solution for each t ∈ [t k , t k+1 ), k ∈ N. Consequently, a solution can be constructed by selecting the last value of the previous sampling interval as the initial condition for the next sampling interval so that it is continuous at each sampling instant.

B. Problem Formulation

In the present work, we adopt the Rε-stability for the system (1)-( 2), which is defined as follows.

Definition 1. Rε-stability [START_REF] Polyakov | Practical stabilization via relay delayed control[END_REF] Consider positive scalars R and ε, such that 0 < ε < R, and a Lyapunov function V :

H 1 ([0, L]; R n ) → R + .
If for all solutions of system (1) with z 0 (x) ∈ L V <R , the trajectory of the state z(t, x) converges to L V ≤ε as t goes to infinity, then, system (1) is called Rεstable from L V <R to L V ≤ε .

Our main goal is to ensure that the closed-loop system (1)-( 2) is Rε-stable due to Input-Output method.

III. MAIN RESULT

This section consists of two parts. Firstly, the sampled-data system is equivalently expressed as a continuous hyperbolic PDE with sampling induced error as disturbances in the input. Secondly, we provide constructive Rε-stability criteria based on the provided model.

A. System Remodelling

System (1) can be rewritten equivalently as

         ∂ t z (t, x) + Λ∂ x z (t, x) + (Γ + ̥)z (t, x) + ̥̟ (t, x) = 0, ∀t ∈ [t k , t k+1 ), k ∈ N, z(t, 0) = 0, ∀t ≥ 0, z(0, x) = z 0 (x), ∀x ∈ [0, L] , (5a)
(5b) (5c)
with the sampling error

̟ (t, x) = z (t k , x) -z (t, x) . (6)
ϕ P ̟ J Define the function ϕ

(t, x) = ∂z(t,x) ∂t , ∀t ≥ 0, x ∈ [0, L] . Note that for all t ∈ [t k , t k+1 ) , k ∈ N, x ∈ [0, L], we have ̟ (t, x) = - t t k ∂z (θ, x) ∂θ dθ = - t t k ϕ (θ, x) dθ. (7)
Therefore, the closed-loop system can be seen as the interconnection of two systems P and J shown in Fig. 1, where the operator P : L 2 (0, L) → L 2 (0, L) is defined by

P :                ∂ t z (t, x) = -Λ∂ x z (t, x) -(Γ + ̥) z (t, x) -̥̟ (t, x) , z(t, 0) = 0,∀t ≥ 0, z(0, x) =z 0 (x), ∀x ∈ [0, L] , ϕ (t, x) = -Λ∂ x z (t, x) -(̥ + Γ) z (t, x) -̥̟ (t, x) = ∂ t z (t, x) , (8)
and the operator J :

L 2 (0, L) → L 2 (0, L) is defined by J : ̟ (t, x) = (J z)(t, x) = - t t k ϕ (θ, x) dθ, ∀t ∈ [t k , t k+1 ) , k ∈ N, x ∈ [0, L] . (9)
Remark 2. The operator P is a nominal continuous-time control-loop since we use continuous-time sampling error instead of sampled-data controller in the system (8), J is an integral operator representing the sampling error. For the simplicity of the closed-loop system structure, we choose this form of J so that we need one output ϕ instead of two z(t, x) and z(t k , x), which is why we wrote (6) as [START_REF] Hetel | Discrete and intersample analysis of systems with aperiodic sampling[END_REF].

B. Stability Analysis

In the following, we present our primary results. Proposition 1. Consider systems ( 8)-( 9) with (2) and a function V : H 1 ([0, L] ; R n ) → R + which is differentiable w.r.t. its argument and such that there exists

0 < ι 1 < ι 2 satisfying ι 1 ̺ 2 H 1 ([0,L]; R n ) ≤ V (̺) ≤ ι 2 ̺ 2 H 1 ([0,L]; R n ) .
Suppose that along the trajectories of the system (8)-( 9), the corresponding solution z(t, •) satisfies V (z) + 2δV (z) ≤ 0, for some δ > 0, whenever

1) R > V (z(t, •)) ≥ max{ε, V (z(t k , •))/α}, with some α > 1, 2) z(t k , •) ∈ L V <R .
Then the system is Rε-stable from L V <R to L V ≤ε . The proof of Proposition 1 can be found in the appendix. Remark 3. Proposition 1 is based on the generalization of the Razumikhin technique to get the Rε-stability for hyperbolic systems. In the following theorem, we will show how it can be used in constructed manner.

W (x) =      -e -2µx (̥ + Γ) T Θ 1 + Θ 1 (̥ + Γ) -e -2µx Θ 1 ̥ 0 0 * -γI 0 0 * * -e -2µx Γ T Θ 2 + Θ 2 Γ + βΘ 2 -e -2µx Θ 2 ̥ * * * -γI      (11) N (x) = e -2µx     (α -1)Θ 1 -Θ 1 0 0 * -Θ 1 0 0 * * αΘ 2 0 * * * -Θ 2     (12)
Theorem 1. Consider systems ( 8)-( 9) with ( 2) and an initial condition satisfying (3):

1) Let λ = min i∈{1,...,n} λ i . Assume that there exist µ, γ, κ > 0, α > 1 and symmetric positive matrices Θ 1 ∈ R n×n , Θ 2 ∈ R n×n satisfying the commutativity conditions:

ΛΘ 1 = Θ 1 Λ, ΛΘ 2 = Θ 2 Λ and W (0) + κN (0) 0, W (L) + κN (L) 0, (10)
with W (x) and N (x) defined for all x ∈ [0, L] as ( 11)-( 12).

2) For given decay rate δ > 0,

∃ε ∈ R + , R ∈ R + s.t. 0 < ε < R and it holds γ3 h |Λ| 2 Ω 1 + |Γ| 2 + |̥| 2 Ω 2 + γΩ 1 ≤ (2σ -β)ε -2δR, (13)
for some 0

< β < 2σ withΩ 1 = R λmin(Θ2)e -2µL , Ω 2 = R λmin(Θ1)e -2µL , σ = µλ.
Then the considered system (1) is Rε-stable from L V <R to L V ≤ε for any sampling sequence satisfying (2), with the Lyapunov function defined by

V (z) = V 1 (z) + V 2 (z), (14)
where

V 1 (z) = L 0 z T e -2µx Θ 1 zdx, V 2 (z) = L 0 z T x e -2µx Θ 2 z x dx.

Proof. Consider the Lyapunov function (14). It can be bounded as

Φ z (t, •) 2 H 1 ( [0,L] ; R n ) ≤ V (z (t, •)) ≤ Ψ z (t, •) 2 H 1 ([0,L]; R n ) , where Φ = min{λ min (Θ 1 ) , λ min (Θ 2 )}e -2µL , Ψ = max{λ max (Θ 1 ) , λ max (Θ 2 )}.
Step 1: In this step, we clarify that the function V defined in ( 14) is continuous by using the construction method [START_REF] Wang | Stability Analysis for A Class of Linear Hyperbolic System of Balance Laws with Sampled-data Control[END_REF].

Remark 4. V 1 is used to bound z, and V 2 is used to deal with the term z x that appears in the derivative of V 1 .

Step 2: In this step we study the time derivative of V (z) defined in [START_REF] Fridman | Robust sampled-data control of a class of semilinear parabolic systems[END_REF]. Thanks to commutativity condition:

ΛΘ 1 = Θ 1 Λ, we first compute the time derivative of V 1 (z) along the solutions to (8)-(9),∀t ∈ [t k , t k+1 ) , k ∈ N, V1 (z) = L 0 ∂ t z T e -2µx Θ 1 z + z T e -2µx Θ 1 ∂ t z dx = L 0 (-Λ∂ x z -(̥ + Γ)z -̥̟) T e -2µx Θ 1 z + z T e -2µx Θ 1 (-Λ∂ x z -(̥ + Γ)z -̥̟) dx = -z T Λe -2µx Θ 1 z L 0 + L 0 -z T (̥ + Γ) T e -2µx Θ 1 +e -2µx Θ 1 (̥ + Γ) z -̟ T ̥ T e -2µx Θ 1 z -z T e -2µx Θ 1 ̥̟ dx -2µ L 0 z T Λe -2µx Θ 1 zdx. ( 15
)
In order to get the time derivative of z x in V 2 , we refer to the original system (1). Since z

: [0, +∞) × [0, L] → R n has consecutive partial derivatives in [0, +∞) × [0, L],
according to Schwartz's theorem [START_REF] James | Advanced calculus belmont[END_REF] we can obtain ∀t ∈ (t k , t k+1 )

∂ xt z (t, x) = ∂ tx z (t, x) = -Λ∂ xx z (t, x) -Γ∂ x z (t, x) -̥∂ x z (t k , x) . ( 16
)
For the next calculation of the time derivative of V 2 , we use Lemma 1 in the appendix. According to [START_REF] Hetel | Recent developments on the stability of systems with aperiodic sampling: An overview[END_REF] and Lemma 1, we have

∂ x z(t, 0) = 0, ∀t ∈ [t k , t k+1 ), k ∈ N, z 0 (0) = 0, ∂ x z 0 (0) = 0. (17a) (17b)
Similarly to the computation of V1 , by using the commutativity condition: ΛΘ 2 = Θ 2 Λ, the time derivative of V 2 (z) along the solutions to ( 16)-( 17), ∀t ∈ [t k , t k+1 ) , k ∈ N is shown as follows

V2 (z) = -∂ x z T Λe -2µx Θ 2 ∂ x z L 0 + L 0 -∂ x z T Γ T e -2µx Θ 2 + e -2µx Θ 2 Γ ∂ x z -∂ x z T (t k , •) ̥ T e -2µx Θ 2 ∂ x z -∂ x z T e -2µx Θ 2 ̥∂ x z (t k , •) dx -2µ L 0 ∂ x z T Λe -2µx Θ 2 ∂ x zdx. ( 18
)
Adding γ ̟ (s, •) [START_REF] Kang | Distributed sampled-data control of Kuramoto-Sivashinsky equation[END_REF] for some γ > 0, β > 0, and using boundary condition (5b) and (17a) we have

2 L 2 ( [0,L]; R n ) -γ ̟ (s, •) 2 L 2 ( [0,L]; R n ) to (15) and γ ∂ x z (t k , •) 2 L 2 ( [0,L]; R n ) - γ ∂ x z (t k , •) 2 L 2 ( [0,L]; R n ) , β L 0 z T x e -2µx Θ 2 z x dx - β L 0 z T x e -2µx Θ 2 z x dx to
V (z) = V1 (z) + V2 (z) ≤ -2σV 1 (z) -(2σ -β)V 2 (z) + L 0 η T W (x)ηdx + γ ̟ (s, •) 2 L 2 ( [0,L]; R n ) + γ ∂ x z (t k , •) 2 L 2 ( [0,L]; R n ) . (19)
with

σ = µλ, η = [z T ̟ T (∂ x z) T (∂ x z(t k , •)) T ]
T , and W (x) defined in [START_REF] Wang | Stability Analysis for A Class of Linear Hyperbolic System of Balance Laws with Sampled-data Control[END_REF].

Step 3: In this step, we show that V (z) + 2δV (z) ≤ 0, whenever

R > V (z(t, •)) ≥ max{ε, V (z(t k , •))/α}, z(θ, •) ∈ L V <R , ∀θ ∈ [t k , t), k ∈ N. (20a) (20b)
Let us assume that conditions (20) hold. Since condition [START_REF] Espitia | Event-based control of linear hyperbolic systems of conservation laws[END_REF] is linear in e -2µx and 0 ≤ x ≤ L, by convexity, we have

W (x) + κN (x) ≤ 0, for x ∈ [0, L]. Therefore, we get L 0 η T (W (x) + κN (x))ηdx ≤ 0, (21)
with W (x) and N (x) given in ( 11) and ( 12). Now, consider t ∈ [t k , t k+1 ) and a trajectory z satisfying [START_REF] Karafyllis | Sampled-data boundary feedback control of 1-D linear transport pdes with non-local terms[END_REF]. Since condition (20a) is satisfied, we have V (z(t, •)) ≥ V (z(t k , •))/α with some α > 1, which can be rewritten as

L 0 η T N (x)ηdx ≥ 0. (22)
In view of ( 21), ( 22) and κ > 0, by S-procedure, it implies

L 0 η T W (x)ηdx ≤ 0. (23)
According to condition (20b), the following inequalities are further derived for all θ ∈ [t k , t]:

z (θ, •) 2 L 2 ([0,L];R n ) < R λ min (Θ 1 ) e -2µL = A 1 , ∂ x z (θ, •) 2 L 2 ([0,L];R n ) < R λ min (Θ 2 ) e -2µL = A 2 . (24)
Using ( 7) and [START_REF] Gu | Survey on recent results in the stability and control of time-delay systems[END_REF]

for t ∈ [t k , t k+1 ), k ∈ N, the upper bound of ̟ (s, •) 2 L 2 ([0,L]; R n ) can be calculated ̟ (s, •) 2 L 2 ([0,L]; R n ) = L 0 |̟ (s, x)| 2 dx = L 0 t t k ∂z (θ, x) ∂θ dθ 2 dx ≤ 3 L 0 t t k |Λ| 2 |∂ x z (θ, x)| 2 + |Γ| 2 |z (θ, x) | 2 + |̥| 2 |z (t k , x)| 2 dθdx ≤ 3 h |Λ| 2 A 2 + |Γ| 2 + |̥| 2 A 1 = ω. (25)
In addition, since condition (20a) is satisfied, we have

-2σV 1 (z) -(2σ -β)V 2 (z) ≤ -(2σ -β)(V 1 (z) + V 2 (z)) < -(2σ -β)ε. (26)
Therefore, substituting [START_REF] Fridman | Lyapunov-based stability analysis[END_REF], [START_REF] Logemann | Stability of infinitedimensional sampled-data systems[END_REF] and ( 26) into [START_REF] Kao | Stability analysis of systems with uncertain time-varying delays[END_REF], we have for all t ∈ [t k , t k+1 ), k ∈ N,

V (z) < -(2σ -β)ε + γω + γ R λ min (Θ 2 ) e -2µL . (27)
Since ( 13) holds, we deduce from [START_REF] Mirkin | Some remarks on the use of time-varying delay to model sample-and-hold circuits[END_REF],

V (z) < -2δR ≤ -2δV (z). (28)
Therefore, we have shown that V (z)+2δV (z) ≤ 0, whenever conditions [START_REF] Karafyllis | Sampled-data boundary feedback control of 1-D linear transport pdes with non-local terms[END_REF] are satisfied.

Step 4: In this step, we show that if [START_REF] Karafyllis | Sampled-data boundary feedback control of 1-D linear transport pdes with non-local terms[END_REF] are going to be satisfied for any t ∈ [t k , T • ). From step 3, we know that V is going to decrease during that time interval, either continuously, or until V reaches below max{ε, V (z(t k , •))/α} and when it reaches that region, it never gets back out. Therefore, we have V (z(T • , •)) < V (z(t k , •)) < R, which contradicts the assumption that there exists t • ∈ (t k , t k+1 ) such that V (z(t • , •)) ≥ R.

z(t k , •) ∈ L V <R , then z(t, •) ∈ L V <R , ∀t ∈ [t k , t k+1 ). Consider z such that z(t k , •) ∈ L V <R , assume that ∃ t • ∈ (t k , t k+1 ) s.t. V (z(t • , •)) ≥ R. Let us then call T • the minimum of such t • , then ∀t ∈ [t k , T • ), V (z(t, •)) < R. Therefore conditions
Summary: From step 3 and step 4, it is clear that

V (z) + 2δV (z) ≤ 0 wherever R > V (z(t, •)) ≥ max{ε, V (z(t k , •))/α}, z(t k , •) ∈ L V <R , (29a) (29b)
and therefore, the conditions of Theorem 1 are satisfied, which concludes the proof of Rε-stability.

Remark 5. It is worth pointing out that we use several parameters and now we summarize each parameter in detail. For Rε-stability, R is the domain of attraction for a given Lyapunov function, ε specifies the positive invariant level set of V . They satisfy 0 < ε < R. In this paper, we can fix R then compute ε or vice versa. α is a parameter introduced in the Lyapunov-Razumikhin method to define level set in which the time derivative of V (z(t, •)) should be negative between two sampling interval, we choose it greater than 1. The closer α is to 1, the greater the values of V (z(t k , •))/α are, and the less conservative the conditions of V-convergence are. µ is related to the decay rate of V 1 , V 2 , and δ is related to the decay rate of V . γ and κ are found by line search to realize the conditions given in Theorem 1. First, the algorithm of Theorem 1 is implemented in Matlab using Yalmip [START_REF] Lofberg | YALMIP: A toolbox for modeling and optimization in MATLAB[END_REF] to solve the condition 1). Then we use the same parameters to test the condition 2). Due to (13), we adjust γ, β to be the smallest possible and µ to be the largest possible. In the numerical section, we sort out their relationship: µ > 0, δ > 0, λ = min i∈{1,...,n}

λ i > 0, σ = µλ, 0 < β < 2σ.

IV. NUMERICAL SIMULATION

In this section, we present a numerical example to illustrate the method we proposed in Section III. According to Remark 5, the parameters in condition (10) are selected as: µ = 0.09, κ = 1.8, γ = 0.001, α = 1.001, then we choose β = 0.01, δ = 0.001 satisfying condition [START_REF] Espitia | Stabilization of boundary controlled hyperbolic pdes via Lyapunov-based event triggered sampling and quantization[END_REF]. We fix R = 20, and choose appropriate Θ 1 , Θ 2 to observe the evolution of states. The simulation results are introduced in Figs. 234. Figs. 2-3 present that both state trajectories converge to near the origin with the controller and the initial conditions satisfying the compatibility condition (3). As can be seen from Fig. 4, the time-evolution of Lyapunov function V (z(t, •)) decreases when R > V (z(t, •)) ≥ max{ε, V (z(t k , •))/α}, α > 1. Please note that using the method proposed in [START_REF] Wang | Stability Analysis for A Class of Linear Hyperbolic System of Balance Laws with Sampled-data Control[END_REF] is not feasible in the case presented here.

V. CONCLUSIONS

The main work of this paper is to use a sampling controller for distributed control of linear hyperbolic balance laws. The closed-loop system is reformulated based on Input-Output approach. New stability condition has been obtained by means of the Lyapunov-Razumikhin method. In the future, we will consider the global stability with controller discretized both in time and in space.

Fig. 1 .

 1 Fig. 1. Alternative representation of the closed-loop system.

1 -

 1 cos 2πx) sin 4πx 0.15(1 -cos 4πx) sin 2πx .

Fig. 2 .Fig. 3 .Fig. 4 .

 234 Fig. 2. Response of state z 1 .

APPENDIX Lemma 1. Consider the system (1)-( 2) with initial condition z 0 satisfying Condition 1. Then ∀t ∈ [t k , t k+1 ), k ∈ N, ∂ x z(t, 0) = 0.

Proof: We recall system (1), the time derivative of the boundary condition leads to

The proof of Proposition 1: During a sampling interval [t k , t k+1 ) with initial state z(t k , •):

(1

(

Then we will discuss two possibilities in case (b):

))/α} otherwise, according to the same principle, it would contradict the proposition in Theorem 1. So, over the whole sampling interval t ∈ [t k , t k+1 ), we can get inequality [START_REF] Omran | Stability analysis of bilinear systems under aperiodic sampled-data control[END_REF].

(b2) When V (z(t, •)) > max{ε, V (z(t k , •))/α}, ∀t ∈ [t k , t k+1 ), since V (z) + 2δV (z) ≤ 0, we have V (z(t, •)) ≤ e -2δ(t-t k ) V (z(t k , •)). Then it is not hard to get [START_REF] Omran | Stability analysis of bilinear systems under aperiodic sampled-data control[END_REF].

Consider

with ξ = max{1/α, e -2δ(t-t k ) } ≤ 1, then we can derive

, where h is the lower bound of the sampling interval. By recursion, the following inequality holds if z(t 0 , •) ∈ L V <R , ∀k ∈ N, we have [START_REF] Polyakov | Practical stabilization via relay delayed control[END_REF] Then combining [START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF] and [START_REF] Polyakov | Practical stabilization via relay delayed control[END_REF], we get that

when k is large enough. Therefore, there ∃ k = ⌈log ζ ε / V (z(t 0 , •)) ⌉, such that z(t, •) ∈ L V ≤ε , ∀ t ≥ tk , which leads the proof of Rε-stability.


Xinyong Wang
email: xinyong.wang@univ-lille.fr

Ying Tang
email: ying.tang@univ-lille.fr

Christophe Fiter
email: christophe.fiter@univ-lille.fr

Laurentiu Hetel
email: laurentiu.hetel@univ-lille.fr

Hetel A Laurentiu

Hybrid

A Hybrid System Approach to Exponential Stability with Sampled-data Control for a Class of Linear Hyperbolic Systems

Keywords: Linear hyperbolic systems, sampled-data control, hybrid approach, exponential stability, Lyapunov function

come

A Hybrid System Approach to Exponential Stability

with Sampled-data Control for a Class of Linear Hyperbolic Systems ⋆ Xinyong Wang * Ying Tang * Christophe Fiter * Laurentiu Hetel *

INTRODUCTION

Many physical and chemical phenomena can be described by partial differential equations (PDEs). These dynamics take place in an infinite dimensional space, which makes it harder to analyze than finite dimensional systems. This prompted many scholars to study the control of PDEs (See [START_REF] Krstic | Boundary control of PDEs: A course on backstepping designs[END_REF]; [START_REF] Necas | Weak and measure-valued solutions to evolutionary PDEs[END_REF]; [START_REF] Majda | Introduction to PDEs and Waves for the Atmosphere and Ocean[END_REF]). In practice, controllers are implemented numerically with algorithms on computers. At present, it is a challenging research topic to analyze infinite dimensional systems with sampled-data controllers [START_REF] Logemann | Stabilization of well-posed infinitedimensional systems by dynamic sampled-data feedback[END_REF]; [START_REF] Ke | Adaptive sampled-data integral control of stable infinite-dimensional linear systems[END_REF]). Here we study the sampled-data controller for hyperbolic systems.

Typically, sampled-data control can be handled by discrete time [START_REF] Kao | On stability of systems with aperiodic sampling devices[END_REF]), time-delay [START_REF] Fridman | Robust sampled-data stabilization of linear systems: an input delay approach[END_REF]), Input/Output [START_REF] Fujioka | Stability analysis of systems with aperiodic sample-and-hold devices[END_REF]), and hybrid system [START_REF] Postoyan | A framework for the observer design for networked control systems[END_REF]) methods (see e.g., [START_REF] Hetel | Recent developments on the stability of systems with aperiodic sampling: An overview[END_REF] for a survey). For infinite dimensional systems, fewer results are available. In [START_REF] Logemann | Discrete-time and sampled-data low-gain control of infinite-dimensional linear systems in the presence of input hysteresis[END_REF], the sampled-data low-gain control is studied for systems with input hysteresis.

(see also [START_REF] Logemann | Stability of infinite-dimensional sampled-data systems[END_REF][START_REF] Logemann | Generalized sampled-data stabilization of well-posed linear infinitedimensional systems[END_REF] for other sampled-data control laws). For parabolic PDEs, the time-delay method has been used in Selivanov andFridman (2016, 2017); [START_REF] Fridman | Robust sampled-data control of a class of semilinear parabolic systems[END_REF]; [START_REF] Kang | Distributed sampled-data control of Kuramoto-Sivashinsky equation[END_REF] for systems with distributed sampled-data control. Event-triggered control of hyperbolic PDEs was developed in [START_REF] Espitia | Event-based control of linear hyperbolic systems of conservation laws[END_REF][START_REF] Espitia | Stabilization of boundary controlled hyperbolic PDEs via Lyapunov-based event triggered sampling and quantization[END_REF]. In Espitia et al. (2017a); [START_REF] Davó | Stability analysis of a 2 × 2 linear hyperbolic system with a sampled-data controller via backstepping method and looped-functionals[END_REF], the authors introduced the backstepping approach to stabilize a class of event-triggered hyperbolic systems. Compared with the existing paper addressing the boundary control case [START_REF] Diagne | Lyapunov exponential stability of 1-D linear hyperbolic systems of balance laws[END_REF]; [START_REF] Safi | Tractable sufficient stability conditions for a system coupling linear transport and differential equations[END_REF]; [START_REF] Tang | Stability analysis of coupled linear ode-hyperbolic pde systems with two time scales[END_REF]), few results addressing the case of distributed sampled-data control for hyperbolic systems.

Recently, in our previous work [START_REF] Wang | Stability analysis for a class of linear hyperbolic system of balance laws with sampled-data control[END_REF] and [START_REF] Wang | Sampleddata distributed control for homo-directional linearhyperbolic system with spatially sampled state measurements[END_REF]), based on the Lyapunov-Razumikhin method, new stability conditions have been provided. The method there-⋆ This work was partially supported by the French "Règion Hauts-de-France" grant.

in allows to check local practical stability for hyperbolic PDEs with distributed sampled-data control.

In summary, it can be seen for the literature survey that there is a wide open research space for the analysis of sampled hyperbolic PDEs. In this paper, we proposed a new hybrid system approach for the analysis of hyperbolic PDEs with distributed sampled-data control. The idea of this paper is to use an augmented state model with an impulsive form (see [START_REF] Haddad | Impulsive and hybrid dynamical systems: stability, dissipativity, and control[END_REF]; [START_REF] Naghshtabrizi | Exponential stability of impulsive systems with application to uncertain sampled-data systems[END_REF]) for the finite dimensional case in order to derive Integral Quadratic Constraints (IQCs) [START_REF] Megretski | System analysis via integral quadratic constraints[END_REF]) which characterise the effect of the sampling. Using the obtained IQCs, we derive numerical criteria for analyzing stability. Compared with our previous local practical stability conditions [START_REF] Wang | Stability analysis for a class of linear hyperbolic system of balance laws with sampled-data control[END_REF]), this work provides global exponential stability conditions in the form of linear matrix inequalities (LMIs).

The structure of this paper is as follows. Section 2 re-models the system to be studied into an augmented hybrid system and states the purpose of our research. In Section 3, a useful preliminary result is proposed in order to obtain an IQC characteristic of the sampling effect. Next the main stability result is proposed. A numerical example is shown in Section 4. Finally, we summarize the paper with a conclusion.

Notations: N is a nonnegative integer from 0 to infinity, the set of real numbers is denoted by R, R + is the set of non-negative reals, R n is used to denote the set of n-dimensional Euclidean space with the norm | • |. L 2 (0, L) stands for the Hilbert space of square integrable scalar functions on (0, L) with the correspond-

ing norm ∥ • ∥, defined by ∥ ρ ∥ L 2 (0,L) = √ ( ∫ L 0 |ρ (x)| 2 dx ) . The set of functions ψ : [0, L] → R n such that ∫ L 0 |ψ (x)| 2 dx < ∞ is denoted by L 2 ([0, L]; R n ).
The notation M ≼ 0 denotes that M is a symmetric negative semidefinite matrix, and the symmetric elements are denoted by ⋆. The identity matrix is denoted by I.

PROBLEM STATEMENT

We consider the following hyperbolic systems:

         ∂ t y (t, x) + Λ∂ x y (t, x) + ϒy (t, x) + u (t, x) = 0, u (t, x) = Ky (t k , x) , ∀t ∈ [t k ,t k+1 ), k ∈ N, y(t, 0) = y(t, L), ∀t ≥ 0, y(0, x) =y 0 (x), ∀x ∈ [0, L] , (1a) (1b) (1c) (1d)
where y :

[0, +∞) × [0, L] → R n , Λ = diag {Λ + , Λ -}, Λ + = diag {λ 1 , • • • , λ m }, Λ -= diag {λ m+1 , • • • , λ n } with λ 1 > • • • λ m > 0 > λ m+1 > • • • > λ n , K and ϒ are real n × n constants matrices.
The sampling sequence is defined as {t k } k∈N where

t 0 = 0,t k+1 -t k ∈ (0, h]. ( 2
               ∂ t η (t, x) + A∂ x η (t, x) + Bη (t, x) = 0, ∀t ∈ (t k ,t k+1 ), k ∈ N, η (t k , x) = Cη ( t - k , x ) ,t = t k , k ∈ N, η(t, 0) = η(t, L), ∀t ≥ 0, η(0, x) = η 0 (x) = [y 0 (x) y 0 (x)] T , ∀x ∈ [0, L] . (3a) (3b) (3c) (3d) with A = [ Λ 0 n×n 0 n×n 0 n×n ] , B = [ ϒ K 0 n×n 0 n×n ] , C = [ I n×n 0 n×n I n×n 0 n×n ] , η ( t - k , • ) = lim t↑t k η (t, •).
Note that the initial condition (3d) is chosen the same as the condition (1d), so that all the solution of (1) are characterized by the first component of the augment state. As a result, the closed-loop system can be regarded as an augmented impulsive system of the y (t, x)-variable and the y (t k , x)-variable. In this article, we intend to find the exponential stability criteria of the original system with the help of the augmented one.

STABILITY ANALYSIS

Before the statement of the main result, a technical lemma is first given based on model (3) to characterize the link between the system state y of system (1) and its sampled vector ŷ by an IQC. The idea of Lemma 1 is to use a norm Φ depending on the augmented system state and study its growth along the solution of hybrid system (3) during one sampling interval [t k ,t k+1 ).

Then we overbound the growth of exponential function with maximum growth rate α. Next, stability conditions for system (1) are derived using the obtained IQC.

IQC Condition

Lemma 1: Consider system (1) and augmented system (3). Let α ∈ R and Θ ∈ R 2n×2n a diagonal matrix satisfying 2αΘ -B T Θ -ΘB ≼ 0.

(4) Then, the inequality

∫ L 0 η T (t, x)N(t -t k )η (t, x) dx ≥ 0,t ∈ [t k ,t k+1 ), k ∈ N (5) holds along the solutions η ∈ L 2 ( [0, L] ; R 2n )
which is the solution of (3), with N(h) defined for all h ∈ [0, h] as

N(h) = e -2αh [ 0 n×n I n×n 0 n×n I n×n ] T Θ [ 0 n×n I n×n 0 n×n I n×n ] -Θ. (6)
Proof: Let us consider the following functional Φ defined for all η ∈ L 2 ([0, L]; R 2n ):

Φ(η) = ∫ L 0 η T Θηdx, (7)
with η ∈ L 2 ( [0, L] ; R 2n ) solution of (3).
Computing the time derivative of Φ along the solution to (3)

Φ(η) = ∫ L 0 ( η t T Θη + η T Θη t ) dx = ∫ L 0 ( (-A∂ x η -Bη) T Θη + η T Θ (-A∂ x η -Bη) ) dx = ∫ L 0 -∂ x [ η T AΘη ] dx + ∫ L 0 ( -η T ( B T Θ + ΘB ) η) dx,
(8) using the boundary condition (3c), we get

Φ(η) = ∫ L 0 ( η T ( -B T Θ -ΘB ) η) dx. ( 9
)
Since the condition (4) holds, then we have

Φ(η) + 2αΦ(η) = ∫ L 0 ( η T ( 2αΘ -B T Θ -ΘB ) η) dx ≤ 0. (10)
Then, according to the comparison lemma we have

Φ(η(t, •)) ≤ e -2α(t-t k ) Φ(η(t k , •)), ∀t ∈ [t k ,t k+1 ), k ∈ N. (11)
The inequality (11) can be rewritten as:

e -2α(t-t k ) Φ(η(t k , •)) -Φ(η(t, •)) = ∫ L 0 η T (t, x)N(t -t k )η (t, x) dx ≥ 0, (12)
with N(tt k ) defined in (6).

Remark 1. The parameter α in the above lemma represents an upper bound on the growth rate of the function Φ (norm of the state η) between two sampling points. This upper bound is shown in (11). The following example provides a more intuitive explanation. The function Φ captures the growth of the norm of η.

Example 1. Consider system (3) with

L = 1, A =    1.2 0 0 0 ⋆ -0.8 0 0 ⋆ ⋆ 0 0 ⋆ ⋆ ⋆ 0    , B =    1.3 1.5 1.1 0.5 -0.5 2.5 1 0.9 0 0 0 0 0 0 0 0    y 0 (x) = [ 2.5 -1.5 cos 4πx -0.8 sin 6πx -1 ] ,
and a maximum sampling interval h = 0.1. Condition (4) in Lemma 1 is an LMI that can be used in order to characterize the link between the system state and sampled version based on the IQC (5). Using Lemma 1 with α = 3, we can derive an IQC of the form ( 5) and ( 6) with ( 13) given below:

Θ =    -0.25 0 0
Fig. 1. Time-evolution of the norm Φ along the trajectory of system (3) for the system describing Example 1 (green) and its upper bound in each sampling interval (red).

The bounded growth performance of function Φ(η(t, •)) and the upper bound e -2α(t-t k ) Φ(η(t k , •)) as stated in ( 11) for all sampling intervals [t k ,t k+1 ), k ∈ N are shown in Fig. 1.

Main Stability Result

Next, the following theorem gives two LMIs allowing to check the exponential stability of system (1).

Theorem 1: Consider the system (1)-( 2). Assume that there exist ε > 0, λ > 0, α ∈ R, a diagonal positive definite matrix Q ∈ R n×n and a diagonal matrix Θ ∈ R 2n×2n satisfying (4) in Lemma 1

2αΘ -B T Θ -ΘB ≼ 0, and

M + εN(0) ≼ 0, M + εN( h) ≼ 0, ( 14
)
with M defined as

M = [ 2λ Q -ϒ T Q -Qϒ -QK ⋆ 0 n×n ] , ( 15
)
and N(h) defined for all h ∈ [0, h] as in (6):

N(h) = e -2αh [ 0 n×n I n×n 0 n×n I n×n ] T Θ [ 0 n×n I n×n 0 n×n I n×n ] -Θ.
Then system (1) is exponentially stable in L 2 -norm for any sampling sequence satisfying (2), with a decay-rate larger than λ .

Proof: We consider the original system (1) and choose a Lyapunov function as

V (y) = ∫ L 0 y T Qydx. ( 16
)
For simplicity, we use the notation y instead of y(t, x). For t ∈ [t k ,t k+1 ), k ∈ N, the time derivative of V (y) along the trajectories of (1) is

V (y) = ∫ L 0 ( y t T Qy + y T Qy t ) dx = ∫ L 0 ( (-Λ∂ x y -ϒy -Ky (t k , •)) T Qy + y T Q (-Λ∂ x y -ϒy -Ky (t k , •)) ) dx = - [ y T ΛQy ] L 0 + ∫ L 0 ( -y T ( ϒ T Q + Qϒ ) y -y T (t k , •) K T Qy -y T QKy (t k , •) ) dx. ( 17
)
Using the boundary condition (1c) and adding the term 2λV (y) to both sides of ( 17), then, we have

V (y) + 2λV (y) = ∫ L 0 ( y T ( 2λ Q -ϒ T Q -Qϒ ) y (18) -y T (t k , •) K T Qy -y T QKy (t k , •) ) dx = ∫ L 0 η T Mηdx (19) with M defined in (15) and η = [ y T (t, x) ŷT (t, x) ] T , ∀t ∈ [t k ,t k+1 ), k ∈ N.
Since condition ( 14) is linear in e -2αh and tt k < h, by using a convexity argument, we have

M + εN(t -t k ) ≼ 0, for t ∈ [t k ,t k+1 ), k ∈ N and t k+1 -t k ∈ (0, h]. Therefore, we get ∫ L 0 η T (M + εN(t -t k ))ηdx ≤ 0, (20)
recalling the condition ( 5)

∫ L 0 η T N(t -t k )ηdx ≥ 0.
Since condition (14) holds and the Integral S-procedure in the appendix implies that if (20) and ( 5) are satisfied then we obtain

∫ L 0 η T Mηdx ≤ 0. ( 21
)
In view of ( 18) and ( 21), we have V (y) + 2λV (y) ≤ 0 (22) for ∀t ∈ [t k ,t k+1 ), k ∈ N, which means that the following inequality is obtained:

V (y(t, •)) ≤ e -2λ (t-t k ) V (y(t k , •)), ∀t ∈ [t k ,t k+1 ), k ∈ N, by recursion, we get V (y(t, •)) ≤ e -2λt V (y 0 ) , i.e. ∀t ∈ [t k ,t k+1 ) , k ∈ N, we have ∥ y(t, •) ∥ 2 L 2 ([0,L];R n ) ≤ λ max (Q) λ min (Q) e -2λt ∥y 0 ∥ 2 L 2 ([0,L];R n ) .
Hence, according to Proposition 1, we can conclude that the system is exponentially stable.

Remark 2. Theorem 1 provides a numerical method for checking the stability of system (1). In order to ease the applicability of the results, we now summarize each parameter in detail. α is a scalar which is an upper bound of the growth rate of functional Φ used in Lemma 1. λ is a positive scalar which is an lower bound of the decay rate of V . ε is found by line search to realize ( 14). With a fixed α, Θ, Q can be found by solving the LMIs in ( 4) and ( 14) in Matlab using Yalmip [START_REF] Lofberg | Yalmip: A toolbox for modeling and optimization in matlab[END_REF]. An example illustrating the tuning of these different parameters is given in section 4.

NUMERICAL EXAMPLE

Consider system (1) and ( 2) with

L = 1, Λ = [ 1.2 0 0 -0.8 ] , ϒ = [ 1.3 1.5 -0.5 2.5 ] , y 0 (x) = [ 2.5 -1.5 cos 4πx -0.8 sin 6πx -1 ] .
Then, we consider the closed-loop system under the sampleddata controller with

K = [ 1.1 0.5 1 0.9 ] .
Following Remark 2, Theorem 1 was tested for several values of α and h. The results are illustrated in Fig. 2, where the growth rate α of Φ can be infinite, however, there is a minimum value and as the growth rate α decreases, h increases. This is consistent with Lemma 1. Fig. 3 illustrates that the sampling interval h decreases as the decay rate λ increases for some constant parameters: α = 3, ε = 1.

Such figure can be used in order to find a tradeoff between system performance (in terms of decay rate) and robustness which satisfy the conditions (4) and ( 14). The results for the closed-loop system are presented in Figs. [4][5][6] show that the states converge to the origin with the controller. The time-evolution of Lyapunov function V is shown in Fig. 6.

CONCLUSION

This paper focused on the distributed sampled-data control for a class of hyperbolic systems using hybrid system approach. The closed-loop system is first represented as an augmented impulsive hybrid system. In addition, by means of the IQCs, we prove the exponential stability of the system. In the future, we will pay more attention to sampled-data control using space and time discretization for different types of hyperbolic systems.

Fig. 2 .

 2 Fig. 2. Feasible region (blue area) and nonfeasible region (red area) guaranteed by Theorem 1 with λ = 10 -5 , ε = 1.

Fig. 4 .Fig. 6 .

 46 Fig. 4. Response of state y 1 for the closed-loop system.

  ) and h > 0. Let ŷ indicate a piecewise constant signal representing the latest state measurement of the plant available at the controller, ŷ(t, x) = y(t k , x), for all t ∈ [t k ,t k+1 ), k ∈ N. Using the augmented system state η = [y

T (t, x), ŷT (t, x)] T ∈ R 2n , we recast system (1) into an augmented hybrid model with the following structure:

⋆ -0.39 0 0 ⋆ ⋆ -0.19 0 ⋆ ⋆ ⋆ -0.19    . (13)

Appendix A. A LEMMA AND PROOF

Lemma 2 (Integral S-procedure): The property Integral Sprocedure is given as follows: Let F and G be symmetric matrices. Assume that the strict integral inequality ∫ L 0 y T (x)Gy(x)dx ≥ 0, holds. Then the implication ∫ L 0 y T (x)Fy(x)dx ≤ 0 holds if there exists some nonnegative number ϖ such that F + ϖG ≼ 0 Proof: It is immediately proved by rewriting the inequality above as F ≼ -ϖG and multiplying on the left by y T (x) and on the right by y(x), then taking the integral operation on both sides of the inequality. We have

Since ∫ L 0 y T (x)Gy(x)dx ≥ 0 and ϖ ≥ 0 hold, one can obtain

Hence, this lemma is proved.

Appendix B. A PROPOSITION AND PROOF

Proposition 1 : Consider the systems (1)-( 2) and a candidate Lyapunov function

Assume that along the trajectories of the system (1)-( 2), the corresponding solution y(t,

2) for some λ > 0. Then the system is exponentially stable in L 2 -norm with a decay-rate larger than λ , that is for any initial condition

Proof: First, we declare that V defined in (B.1) is continuous (One can consult [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF]). Since y(t, x) is continuous with respect to t for all t ∈ [t k ,t k+1 ), k ∈ N, and continuous at sampling instants by choosing the last value of the previous sampling interval as the initial condition of the following sampling interval, then V is continuous for all t ≥ 0.

Then, we consider the differentiable Lyapunov function V :

Applying the comparison lemma, we have

) by recursion, the following inequality holds

Then instituting (B.7) into (B.5), we obtain

=e -2λ (t-t 0 ) V (y(t 0 , •))

=e -2λt V (y 0 ). (B.8) Combining (B.1) and (B.8), we get

This concludes the proof of Proposition 1.