Hong Diep Nguyen

Nathalie Revol
email: nathalie.revol@ens-lyon.fr

Philippe Théveny
email: philippe.theveny@ens-lyon.fr

Tradeoffs between Accuracy and Efficiency for Optimized and Parallel Interval Matrix Multiplication

Keywords: interval arithmetic, verified computation, inf-sup representation, mid-rad representation, floating-point arithmetic, accuracy, overestimation, efficiency, BLAS, parallel implementation

Interval arithmetic is mathematically defined as set arithmetic. For implementation issues, it is necessary to detail the representation of intervals and to detail formulas for the arithmetic operations. Two main representations of intervals are considered here: inf-sup and midrad. Formulas for the arithmetic operations, using these representations, are studied along with formulas that trade off accuracy for efficiency. This tradeoff is particularly blatant on the example of interval matrix multiplication, implemented using floating-point arithmetic: according to the chosen formulas, the efficiency as well as the accuracy can vary greatly in practice, and not necessarily as predicted by the theory. Indeed, theoretical predictions are often based on exact operations, as opposed to floating-point operations, and on operations count, as opposed to measured execution time. These observations and the recommendations that ensue are further obfuscated by considerations on memory usage, multithreaded computations. . . when these algorithms are implemented on parallel architectures such as multicores.

Introduction

Interval arithmetic is a mean to implement operations on sets and not only on single real values or points. It is also considered as a tool of choice for verified computations. Indeed, the fundamental property of interval arithmetic is the inclusion property: every computed result encloses the exact result. However, it suffers from two main drawbacks: lack of accuracy and lack of efficiency. We show here that a panel of formulas exists for the implementation of interval arithmetic. This panel offers a choice ranging from maximal accuracy to acceptable runtime.

More precisely, the problems one has to face are a lack of:

accuracy: the computed interval result is an overestimation of the exact result, and this overestimation may be too crude; efficiency: each interval operation entails several operations on single values (e.g., 1 interval addition requires to perform 2 real additions) and this seems to be the necessary price to pay. Another source of slowdown lies in the fact that implementing interval arithmetic using floating-point arithmetic requires to change frequently the rounding modes, in order to compute an "outward" enclosure of the result, i.e. to preserve the inclusion of the exact result into the computed one. However, on most processors and programming languages, changing the rounding modes incurs flushing the pipeline and leads to severe slowdowns, up to 2 orders of magnitude, when compared to the (unverified) corresponding floating-point computation.

To remedy the problem of frequent changes of the rounding modes, to enable the use of optimized routines in linear algebra, several variants of the arithmetic operations have been proposed [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF][START_REF] Rump | Fast And Parallel Interval Arithmetic[END_REF][START_REF] Diep | Efficient implementation of interval matrix multiplication, Proceedings of PARA[END_REF][START_REF] Diep | Efficient algorithms for verified scientific computing: numerical linear algebra using interval arithmetic[END_REF][START_REF] Rump | Fast Interval Matrix Multiplication[END_REF]. The main goal of this paper is to show that there is a wide choice of formulas, each satisfying different requirements in terms of accuracy and efficiency, at least as far as the problem of multiplying two interval matrices is concerned. These variants are listed in Section 2, they are based on two representations for intervals, the inf-sup and the mid-rad representations which are recalled. The precise meaning of "the best representation in floating-point arithmetic" will be defined. The use of these variants for interval matrix multiplication is presented in Section 3; the questions are: how far it is possible to reduce the number of rounding modes' changes, and whether it is possible to use floatingpoint BLAS routines. A comparison of the theoretical complexity (the number of floating-point operations), the observed efficiency and the attained accuracy (the overestimation between the best result in floating-point arithmetic and the computed result) is performed experimentally. Finally, some remarks and (preliminary) results on the parallelization of these algorithms for interval matrix multiplication are presented in Section 4: optimizing the memory usage enters the scene.

2 Interval arithmetic: inf-sup and mid-rad representations The "Hull" means that the result is the smallest enclosing interval. Formulas have been devised that translate this abstract definition into implementable expressions. These formulas rely on the chosen representation for intervals. We will restrict ourselves to the inf-sup and mid-rad representations.

Another limit in this paper is the use of bounded intervals only. We refer the interested reader to a discussion on the mailing list of the IEEE 1788 working group for the standardization of interval arithmetic [5, January-February 2012] about the mid-rad representation of unbounded intervals.

In what follows, the notations below will be used:

-a in boldface corresponds to an interval quantity (whether it is an interval, a vector with interval components or a matrix with interval components will be made clear by the context); -[a, ā] with square brackets corresponds to an inf-sup representation with a as the infimum and ā as the supremum; -< a m , a r > with angular brackets corresponds to a mid-rad representation, with a m as the midpoint and a r as the radius.

Arithmetic operations: formulas for the inf-sup and mid-rad representations

Formulas for arithmetic operations +, -and × for both representations can be found in [7, pp. 7-8 and pp.22-23]. Each interval addition or subtraction requires 2 real additions/subtractions. Each interval multiplication requires 4 real multiplications with the inf-sup representation and 6 real multiplications with the mid-rad representation. However, these formulas yield the same result when exact real arithmetic is employed.

Other formulas have been proposed [START_REF] Diep | Efficient implementation of interval matrix multiplication, Proceedings of PARA[END_REF][START_REF] Diep | Efficient algorithms for verified scientific computing: numerical linear algebra using interval arithmetic[END_REF][START_REF] Rump | Fast And Parallel Interval Arithmetic[END_REF][START_REF] Rump | Fast Interval Matrix Multiplication[END_REF] for the multiplication: they reduce the number of real operations, at the cost of an increase of the radius of the result. The ratio between the radius of the resulting interval and the radius of the exact interval is called overestimation factor. For the formulas mentioned above, these overestimation factors have been bounded, the bounds being 4 -2 √ 2 1.18 for the inf-sup formulas in [START_REF] Diep | Efficient implementation of interval matrix multiplication, Proceedings of PARA[END_REF] and the mid-rad formulas in [START_REF] Diep | Efficient algorithms for verified scientific computing: numerical linear algebra using interval arithmetic[END_REF] and 1.5 for the midrad formulas in [START_REF] Rump | Fast And Parallel Interval Arithmetic[END_REF]. In [START_REF] Diep | Efficient implementation of interval matrix multiplication, Proceedings of PARA[END_REF] and [START_REF] Rump | Fast And Parallel Interval Arithmetic[END_REF] only 3 real products have to be performed, in [START_REF] Diep | Efficient algorithms for verified scientific computing: numerical linear algebra using interval arithmetic[END_REF] 5 real products are required.

Implementation using floating-point arithmetic

The overestimation ratios given above are established under the assumption that exact arithmetic is used. Let us focus now on what happens to the result of an interval operation when floating-point arithmetic is used for the implementation.

Floating-point arithmetic has been standardized as the IEEE-754 standard [START_REF]ANSI/IEEE 754-1985, Standard for Binary Floating-Point Arithmetic[END_REF] and IEEE-754 2008 revision of the standard [START_REF]Standard for Floating-Point Arithmetic[END_REF]. In particular, the standard specifies the result of an arithmetic operation between (one or) two floatingpoint operands: the operands are considered as exact values and the result is the exact result of the operation between these exact values, rounded to the closest floating-point number according to the prescribed rounding mode. Four rounding modes are defined by the IEEE-754 standard, three of them will be of use in this work: RN (rounding-to-nearest), RD (rounding downwards) and RU (rounding upwards).

A natural question when one implements interval arithmetic using floatingpoint arithmetic is "what is the best possible floating-point interval that represents a given interval a?". Here, "best possible" means the smallest with regards to enclosure. The second result below is new.

Proposition 1 (Best floating-point interval for the inf-sup representation). For the inf-sup representation, [RD(a), RU(ā)] is the best floating-point interval that represents [a, ā] in the sense that for each floating-point interval

[b, b] such that [a, ā] ⊂ [b, b] we have [RD(a), RU(ā)] ⊂ [b, b].
Proposition 2 (Best floating-point interval for the mid-rad representation). For the mid-rad representation, the best floating-point interval that represents the interval < a m , a r > is the interval < c m , c r > with

c m = RN(a m ), c r = RU(a r + |c m -a m |).
Every floating-point interval in mid-rad representation that contains < a m , a r > either contains < c m , c r > or its midpoint is not a m rounded to nearest.

Note: In the case where a m is the middle of two consecutive floating-point numbers, there are two best intervals in terms of smallest radius, but the tie rule determines the midpoint of one of them only.

These two theorems put in evidence the fact that, contrary to what happens in R where any representation corresponds to the same interval, in floating-point arithmetic the choice of the representation modifies the result. Indeed, when a mathematical a is rounded in floating-point arithmetic, the resulting interval depends on the chosen representation. In general, it holds that:

[RD(a), RU(ā)] =< RN(a m ), RU(a r + |RN(a m ) -a m |) >
and it is generally not straightforward to determine a priori which interval is best. One can exhibit examples where one representation is better than the other for both cases.

In [START_REF] Rump | Fast Interval Matrix Multiplication[END_REF], formulas for the multiplication of two floating-point intervals in midrad representation are given. These formulas require 3 multiplication of floatingpoint numbers. They explicitely incorporate bounds on the roundoff errors, see notably [START_REF] Nicholas | Accuracy and Stability of Numerical Algorithms[END_REF] for an accessible and comprehensive analysis of roundoff errors in numerical computations.

Interval Matrix Multiplication

Notations: matrices denoted with capital letters. Matrices with interval coefficients denoted with boldface capital letters.

Why is it interesting to use the variants presented above?

Let us cite again [START_REF] Rump | Fast Interval Matrix Multiplication[END_REF] to explain why it is especially interesting to try different variants of the formulas for the arithmetic operations, when it comes to multiply two interval matrices. "The classical approach (i.e. the three nested loops using the classical formulas for the addition and multiplication) is very slow on todays computers because i) it requires 2n 3 times switching the rounding mode, thus ii) jeopardizing compiler optimization and iii) putting the burden of well-known acceleration techniques such as blocked code, parallel implementations etc. on the user." In other words, gain in efficiency is expected when trying other formulas, by factoring the changes of rounding modes, in order to insert changes of rounding modes around large blocks of computations which can be optimized using classical approaches. These large blocks of computations typically correspond to the routines implemented in BLAS, and thus the use of optimized BLAS, such as GotoBLAS or Atlas, directly benefits to the performances.

In what follows, we consider only algorithms for the product of two interval matrices. They are good candidates to put into practice the panel of existing formulas for arithmetic operations, as they are simple algorithms. Indeed, they offer a wide range of opportunities for optimization and parallelization, and still have a polynomial complexity.

Let us also insist on the fact that only the product of matrices based on the classical approach is considered, as fast matrix products are not suited to interval arithmetic, because they suffer from overestimation due to variable dependency [START_REF] Ceberio | Fast Multiplication of Interval Matrices (Interval Version of Strassen's Algorithm)[END_REF]

Algorithms based on the inf-sup representation

The two algorithms considered here are on the one hand the so-called classical approach, where the three nested loops are explicitely written and a change of rounding mode is performed for each operation, at the very heart of these three loops. On the other hand, we also consider the algorithm presented in [START_REF] Diep | Efficient implementation of interval matrix multiplication, Proceedings of PARA[END_REF]. This algorithm adapts the formula for interval multiplication, with 3 real multiplications, to the case of interval matrix multiplication. The two notable features of this algorithm is that only 3 changes of the rounding modes are required, and that 9 products of floating-point matrices are used.

Algorithms for the mid-rad representation

The algorithms considered here are:

a counterpart of the classical approach, where each operation between two intervals is performed using Neumaier's formulas [7, p. 22-23], which incur no overestimation in exact arithmetic; the algorithm presented in [START_REF] Diep | Efficient algorithms for verified scientific computing: numerical linear algebra using interval arithmetic[END_REF], which uses 7 products of floating-point matrices and has a theoretical overestimation ratio (i.e. without considering roundoff errors) of 1.18; the algorithm presented in [START_REF] Rump | Fast And Parallel Interval Arithmetic[END_REF], which uses 4 products of floating-point matrices and has a theoretical overestimation ratio (i.e. without considering roundoff errors) of 1.5; the algorithm presented in [START_REF] Rump | Fast Interval Matrix Multiplication[END_REF], which uses 3 products of floating-point matrices and explicitely incorporates a bound on roundoff errors, based on the computation of an upper bound for the condition number.

The last three references first present a formula for interval multiplication, with 3 to 5 real multiplications required for one interval multiplication. Then it is shown how to adapt this formula to the case of interval matrix multiplication.

Experiments: efficiency and accuracy.

Experiments have been performed using MatLab, matrices are generated using the randn function of MatLab, which takes as argument n the dimension of the matrices (all dimensions are equal). This function returns a matrix of size n × n whose components are uniformly distributed random numbers between 0 and 1. Experimentally, one observes that the execution time of the more elaborate algorithms is roughly n times the execution time of a floating-point matrix multiplication, if the corresponding algorithm performs n punctual matrix multiplication. This is true when the dimension is large enough (typically above 500): the decomposition of the matrices thus becomes negligible (it costs O(n 2 ) operations).

The accuracy of each algorithm is related for one part to its theoretical overestimation ratio, but also for another part on the number of matrix operations it performs: a higher number of matrix operations entails a loss of accuracy. For the algorithm in [START_REF] Rump | Fast Interval Matrix Multiplication[END_REF], as stated by its author, the overestimation ratio is also related to the condition number and this can be experimentally observed.

Optimization and Parallelization of the Product of Interval Matrices

To improve further the execution time, to take benefit from the existing architectures present even in our laptops, we turn now to the parallelization of interval matrix multiplication.

The algorithms presented in Section 3 have some parts which have a high potential for parallelization:

conversion to and fro a given representation to the other one: each coefficient can be converted independently of the other ones; with the mid-rad representation, the computation of the midpoint of the result can be usually performed using any optimized, parallel, routine; then the radius can be computed, again using any optimized, parallel, routine; the same phenomenon appears for the variant based on the inf-sup representation.

This potential should easily be exploited on parallel architectures.

However, the hope that algorithms for interval matrix multiplication, that ultimately resort to floating-point matrix multiplication, would benefit for free of optimizations made for this last computation, is not entirely grounded when it comes to parallelization. Indeed, it does not suffice to use parallel BLAS to gain performances. What follows is based on a preliminary work and we do not have yet any firm conclusion nor any recommendation on this topic.

Let us first make clear that our target architectures in this work are multicores, supporting multithreaded computations.

Rounding modes and multithreading

A first issue is quite technical: during multithreaded computations, it may be difficult to access the device that controls the current rounding mode, to modify it for the duration of the whole computation of a thread and to restore it, or at least to release it, for other threads. Mutual exclusion is a possibility, but it seems contradictory with the concurrent execution of different computations. Maybe a solution as was adopted in fi lib [START_REF] Lerch | filib++, a Fast Interval Library Supporting Containment Computations[END_REF] (and even as soon as in 1998) where no reference was made to the rounding mode and where one ulp was systematically added or subtracted to the endpoints or radius, is a promising one.

Memory usage

Even if some blocks of computations are independent and thus parallelizable, they operate on common data and this raises the issue of memory usage (are the right data available in the cache).

Another issue is that the algorithms mentioned above use temporary matrices. Cache misses, reduction of data transfer is already an issue for developers of floating-point matrix multiplication. This problem is amplified for developers of interval matrix multiplication. Possible hints are to reduce the size of blocks to accomodate extra matrix variables, or maybe to re-compute -rather than store -temporary results to reduce memory usage.

Conclusion

Several formulas are at hand when it comes to implement interval arithmetic and in particular interval matrix multiplication. Choosing the formulas or algorithm that is most suited to a particular application should be guided by the requirements in terms of accuracy: are input intervals thin or large? to what extent is overestimation acceptable? efficiency: what is the acceptable slowdown?

It is a nice result that interval computations are no more doomed to prohibitive slowdowns, thanks mainly to the continued advances by S. M. Rump during the last decade. However, the user still has to trade off one criterion for the other. The goal of this paper was in part to show that several possibilities exist, and in the other part to illustrate the differences between these possibilities, both on sequential and on parallel architectures.

Further work include the parallelization of the product of interval matrices on a wider range of architectures, notably including GPU.

2. 1

 1 Inf-sup and mid-rad representationsThe mathematical definition of binary arithmetic operations, in interval arithmetic, considers intervals as sets: for any operation , for any sets A and B, A B = Hull{a b, a ∈ A, b ∈ B}.


Stéphane Gerbi
email: stephane.gerbi@univ-savoie.fr

Belkacem Said-Houari
email: saidhouarib@yahoo.fr

Existence and exponential stability of a damped wave equation with dynamic boundary conditions and a delay term

Keywords: Damped wave equations, boundary delay, global solutions, exponential stability, Kelvin-Voigt damping, dynamic boundary conditions

In this paper we consider a multi-dimensional wave equation with dynamic boundary conditions related to the Kelvin-Voigt damping and a delay term acting on the boundary. If the weight of the delay term in the feedback is less than the weight of the term without delay or if it is greater under an assumption between the damping factor, and the difference of the two weights, we prove the global existence of the solutions. Under the same assumptions, the exponential stability of the system is proved using an appropriate Lyapunov functional. More precisely, we show that even when the weight of the delay is greater than the weight of the damping in the boundary conditions, the strong damping term still provides exponential stability for the system.

Introduction

In this paper we consider the following linear damped wave equation with dynamic boundary conditions and a delay boundary term:

                       u tt -∆u -α∆u t = 0,
x ∈ Ω, t > 0 , u(x, t) = 0,

x ∈ Γ 0 , t > 0 , u tt (x, t) = -∂u ∂ν (x, t) + α∂u t ∂ν (x, t) + µ 1 u t (x, t) + µ 2 u t (x, t -τ )

x ∈ Γ 1 , t > 0 , u(x, 0) = u 0 (x) x ∈ Ω , u t (x, 0) = u 1 (x) x ∈ Ω , u t (x, t -τ ) = f 0 (x, t -τ ) x ∈ Γ 1 , t ∈ (0, τ ) , (1)
where u = u(x, t) , t ≥ 0 , x ∈ Ω , ∆ denotes the Laplacian operator with respect to the x variable, Ω is a regular and bounded domain of R N , (N ≥ 1), ∂Ω = Γ 0 ∪ Γ 1 , mes(Γ 0 ) > 0, Γ 0 ∩ Γ 1 = ∅ and ∂ ∂ν denotes the unit outer normal derivative, α, µ 1 and µ 2 are positive constants. Moreover, τ > 0 represents the time delay and u 0 , u 1 , f 0 are given functions belonging to suitable spaces that will be precised later.

This type of problems arise (for example) in modelling of longitudinal vibrations in a homogeneous bar in which there are viscous effects. The term ∆u t , indicates that the stress is proportional not only to the strain, 1 but also to the strain rate. See [START_REF] Caroll | Singular and Degenerate Cauchy Problems[END_REF]. From the mathematical point of view, these problems do not neglect acceleration terms on the boundary. Such type of boundary conditions are usually called dynamic boundary conditions. They are not only important from the theoretical point of view but also arise in several physical applications. For instance in one space dimension, problem (1) can modelize the dynamic evolution of a viscoelastic rod that is fixed at one end and has a tip mass attached to its free end. The dynamic boundary conditions represents the Newton's law for the attached mass, (see [START_REF] Budak | A collection of problems on mathematical physics[END_REF][START_REF] Andrews | Second order evolution equations with dynamic boundary conditions[END_REF][START_REF] Conrad | stabilization of a flexible beam with a tip mass[END_REF] for more details). In the two dimension space, as showed in [START_REF] Goldstein | Derivation and physical interpretation of general boundary conditions[END_REF] and in the references therein, these boundary conditions arise when we consider the transverse motion of a flexible membrane Ω whose boundary may be affected by the vibrations only in a region. Also some dynamic boundary conditions as in problem (1) appear when we assume that Ω is an exterior domain of R 3 in which homogeneous fluid is at rest except for sound waves. Each point of the boundary is subjected to small normal displacements into the obstacle (see [START_REF] Beale | Spectral properties of an acoustic boundary condition[END_REF] for more details). This type of dynamic boundary conditions are known as acoustic boundary conditions.

In the absence of the delay term (i.e. µ 2 = 0) problem ( 1) has been investigated by many authors in recent years (see, e.g., [START_REF] Gerbi | Asymptotic stability and blow up for a semilinear damped wave equation with dynamic boundary conditions[END_REF], [START_REF] Gerbi | Local existence and exponential growth for a semilinear damped wave equation with dynamic boundary conditions[END_REF], [START_REF] Grobbelaar-Van Dalsen | On fractional powers of a closed pair of operators and a damped wave equation with dynamic boundary conditions[END_REF], [START_REF] Grobbelaar-Van Dalsen | On the solvability of the boundary-value problem for the elastic beam with attached load[END_REF], [START_REF] Pellicer | Large time dynamics of a nonlinear spring-mass-damper model[END_REF], [START_REF] Pellicer | Analysis of a viscoelastic spring-mass model[END_REF]).

Among the early results dealing with the dynamic boundary conditions are those of Grobbelaar-Van Dalsen [START_REF] Grobbelaar-Van Dalsen | On fractional powers of a closed pair of operators and a damped wave equation with dynamic boundary conditions[END_REF][START_REF] Grobbelaar-Van Dalsen | On the initial-boundary-value problem for the extensible beam with attached load[END_REF] in which the author has made contributions to this field.

In [START_REF] Grobbelaar-Van Dalsen | On fractional powers of a closed pair of operators and a damped wave equation with dynamic boundary conditions[END_REF] the author introduced a model which describes the damped longitudinal vibrations of a homogeneous flexible horizontal rod of length L when the end x = 0 is rigidly fixed while the other end x = L is free to move with an attached load. This yields to a system of two second order equations of the form

             u tt -u xx -u txx = 0, x ∈ (0, L), t > 0, u(0, t) = u t (0, t) = 0, t > 0, u tt (L, t) = -[u x + u tx ] (L, t), t > 0, u (x, 0) = u 0 (x) , u t (x, 0) = v 0 (x) , x ∈ (0, L), u (L, 0) = η, u t (L, 0) = µ. (2)
By rewriting problem (2) within the framework of the abstract theories of the so-called B-evolution theory, an existence of a unique solution in the strong sense has been shown. An exponential decay result was also proved in [START_REF] Grobbelaar-Van Dalsen | On the initial-boundary-value problem for the extensible beam with attached load[END_REF] for a problem related to [START_REF] Beale | Spectral properties of an acoustic boundary condition[END_REF], which describe the weakly damped vibrations of an extensible beam. See [START_REF] Grobbelaar-Van Dalsen | On the initial-boundary-value problem for the extensible beam with attached load[END_REF] for more details. Subsequently, Zang and Hu [START_REF] Zhang | Energy decay for a nonlinear viscoelastic rod equations with dynamic boundary conditions[END_REF], considered the problem

u tt -p (u x ) xt -q (u x ) x = 0, x ∈ (0, 1) , t > 0 with u (0, t) = 0, p (u x ) t + q (u x ) (1, t) + ku tt (1, t) = 0, t ≥ 0.
By using the Nakao inequality, and under appropriate conditions on p and q, they established both exponential and polynomial decay rates for the energy depending on the form of the terms p and q.

It is clear that in the absence of the delay term and for µ 1 = 0, problem (2) is the one dimensional model of [START_REF] Andrews | Second order evolution equations with dynamic boundary conditions[END_REF]. Similarly, and always in the absence of the delay term, Pellicer and Solà-Morales [START_REF] Pellicer | Analysis of a viscoelastic spring-mass model[END_REF] considered the one dimensional problem of (1) as an alternative model for the classical spring-mass damper system, and by using the dominant eigenvalues method, they proved that their system has the classical second order differential equation

m 1 u ′′ (t) + d 1 u ′ (t) + k 1 u(t) = 0,
as a limit, where the parameter m 1 , d 1 and k 1 are determined from the values of the spring-mass damper system. Thus, the asymptotic stability of the model has been determined as a consequence of this limit. But they did not obtain any rate of convergence. See also [START_REF] Pellicer | Large time dynamics of a nonlinear spring-mass-damper model[END_REF][START_REF] Pellicer | Spectral analysis and limit behaviours in a spring-mass system[END_REF] for related results.

Recently, the present authors studied in [START_REF] Gerbi | Local existence and exponential growth for a semilinear damped wave equation with dynamic boundary conditions[END_REF] and [START_REF] Gerbi | Asymptotic stability and blow up for a semilinear damped wave equation with dynamic boundary conditions[END_REF] a more general situation of (1). They considered problem (1) with µ 2 = 0, a nonlinear damping of the form g (u t ) = |u t | m-2 u t instead of µ 1 u t and a nonlinear source term f (u) = |u| p-2 u t in the right hand side of the first equation of problem [START_REF] Andrews | Second order evolution equations with dynamic boundary conditions[END_REF]. A local existence result was obtained by combining the Faedo-Galerkin method with the contraction mapping theorem. Concerning the asymptotic behavior, the authors showed that the solution of such problem is unbounded and grows up exponentially when time goes to infinity if the initial data are large enough and the damping term is nonlinear. The blow up result was shown when the damping is linear. Also, we proved in [START_REF] Gerbi | Asymptotic stability and blow up for a semilinear damped wave equation with dynamic boundary conditions[END_REF] that under some restrictions on the exponents m and p, we can always find initial data for which the solution is global in time and decay exponentially to zero. The main difficulty of the problem considered is related to the non ordinary boundary conditions defined on Γ 1 . Very little attention has been paid to this type of boundary conditions. We mention only a few particular results in the one dimensional space [START_REF] Grobbelaar-Van Dalsen | Boundary stabilization for the extensible beam with attached load[END_REF][START_REF] Pellicer | Analysis of a viscoelastic spring-mass model[END_REF][START_REF] Doronin | Global solvability for the quasilinear damped wave equation with nonlinear second-order boundary conditions[END_REF][START_REF] Kirane | Blow-up for some equations with semilinear dynamical boundary conditions of parabolic and hyperbolic type[END_REF].

The purpose of this paper is to study problem [START_REF] Andrews | Second order evolution equations with dynamic boundary conditions[END_REF], in which a delay term acted in the dynamic boundary conditions. In recent years one very active area of mathematical control theory has been the investigation of the delay effect in the stabilization of hyperbolic systems and many authors have shown that delays can destabilize a system that is asymptotically stable in the absence of delays (see [START_REF] Datko | An example on the effect of time delays in boundary feedback stabilization of wave equations[END_REF] for more details).

In [START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF], Nicaise and Pignotti examined the wave equation with a linear boundary damping term with a delay. Namely, they looked to the following problem:

u tt -∆u = 0, x ∈ Ω, t > 0, ( 3
)
where Ω is a bounded domain with smooth boundary ∂Ω = Γ 0 ∪ Γ 1 . On Γ 0 , they considered the Dirichlet boundary conditions. While on Γ 1 they assumed the following boundary conditions:

∂u ∂ν (x, t) = µ 1 u t (x, t) + µ 2 u t (x, t -τ ), x ∈ Γ 1 , t > 0. ( 4
)
They proved under the assumption

µ 2 < µ 1 (5)
that the solution is exponentially stable. On the contrary, if [START_REF] Caroll | Singular and Degenerate Cauchy Problems[END_REF] does not hold, they found a sequence of delays for which the corresponding solution of (3) will be unstable. The main approach used in [START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF], is an observability inequality obtained with a Carleman estimate. The same results were showed if both the damping and the delay are acting in the domain. We also recall the result by Xu, Yung and Li [START_REF] Xu | Stabilization of the wave system with input delay in the boundary control[END_REF], where the authors proved the same result as in [START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF] for the one space dimension by adopting the spectral analysis approach. We point out that problem (1) has been already studied by Nicaise and Pignotti in [START_REF] Nicaise | Interior feedback stabilization of wave equations with time dependent delay[END_REF] for µ 1 = 0 and a time-varying delay. They find the same condition as the one used in this paper when µ 1 = 0 by a different way. However our result and our Lyapunov functional are slightly different here. (See Remark 3.2 for more details), and we want to point out that this paper may be viewed as a continuation of the work of Nicaise and Pignotti [START_REF] Nicaise | Interior feedback stabilization of wave equations with time dependent delay[END_REF] in which an additional damping term acts on the boundary and the study of the competition between these two damping terms is very interesting.

As it has been proved by Datko [9,Example 3.5], systems of the form

w tt -w xx -aw xxt = 0, x ∈ (0, 1), t > 0, w (0, t) = 0, w x (1, t) = -kw t (1, t -τ ) , t > 0, (6)
where a, k and τ are positive constants become unstable for an arbitrarily small values of τ and any values of a and k. In [START_REF] Conrad | stabilization of a flexible beam with a tip mass[END_REF] and even in the presence of the strong damping -aw xxt , without any other damping, the overall structure can be unstable. This was one of the main motivations for considering problem [START_REF] Andrews | Second order evolution equations with dynamic boundary conditions[END_REF]. (Of course the structure of problem ( 1) and ( 6) are different due to the nature of the boundary conditions in each problem). Subsequently, Datko et al [START_REF] Datko | An example on the effect of time delays in boundary feedback stabilization of wave equations[END_REF] treated the following one dimensional problem:

         u tt (x, t) -u xx (x, t) + 2au t (x, t) + a 2 u(x, t) = 0, 0 < x < 1, t > 0, u(0, t) = 0, t > 0, u x (1, t) = -ku t (1, t -τ ), t > 0, (7)
which models the vibrations of a string clamped at one end and free at the other end, where u(x, t) is the displacement of the string. Also, the string is controlled by a boundary control force (with a delay) at the free end. They showed that, if the positive constants a and k satisfy

k e 2a + 1 e 2a -1 < 1,
then the delayed feedback system [START_REF] Grobbelaar-Van Dalsen | On fractional powers of a closed pair of operators and a damped wave equation with dynamic boundary conditions[END_REF] is stable for all sufficiently small delays. On the other hand if

k e 2a + 1 e 2a -1 > 1,
then there exists a dense open set D in (0, ∞) such that for each τ ∈ D, system (7) admits exponentially unstable solutions.

As a consequence of what we have said before, two main questions naturally arise here:

• Is it possible for the damping term -∆u t to stabilize system (1) when the weight of the delay is greater than the weight of the boundary damping (i.e. when µ 2 ≥ µ 1 )?

• Does the particular structure of the problem prevents the instability result obtained in [START_REF] Datko | Two questions concerning the boundary control of certain elastic systems[END_REF] for problem (6)?

One of the main purpose of this paper is to give positive answers to the above two questions. More precisely, we study the asymptotic behavior (as t → ∞) and related decay rates for the corresponding solutions of system [START_REF] Andrews | Second order evolution equations with dynamic boundary conditions[END_REF] where the question to be addressed here is whether the delay term µ 2 u t (x, t -τ ) can destroy the stability of the system, which is exponentially stable in the absence of that delay [START_REF] Gerbi | Asymptotic stability and blow up for a semilinear damped wave equation with dynamic boundary conditions[END_REF]. As we shall see below, the presence of the strong damping term α∆u t in (1) plays a decisive role in the stability of the whole system if (5) does not hold. Thanks to the energy method, we built appropriate Lyapunov functionals lead to stability results.

The paper is organized as follows: in the next section, we prove the global existence of the solutions by using the Lumer-Phillips' theorem in the same way as in [START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF]. In section 3, we show that if the weight of the delay is less than the weight of the damping, then the energy defined by (38) decays exponentially to zero. We also prove that even if the weight of the delay is greater than the weight of the damping, the solution still decays to zero exponentially provided that the damping parameter α satisfies an appropriate condition. Let us mention that without the damping factor α, Nicaise and Pignotti [START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF] proved the instability of the null stationary solution in the case µ 2 ≥ µ 1 , whereas we will show that if µ 2 ≥ µ 1 , by adding a condition of the form α > (µ 2 -µ 1 )B 2 (with B a constant defined later), we are able to prove the stability of the null stationary state thanks to a suitable choice of a Lyapunov function.

2 Well-posedness of Problem [START_REF] Andrews | Second order evolution equations with dynamic boundary conditions[END_REF].

In this section we will first transform the delay boundary conditions by adding a new unknown. Then as in [START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF], we will use the Lumer-Phillips' theorem to prove the existence and uniqueness of the solution of problem (1).

Setup and notations

We denote

H 1 Γ0 (Ω) = u ∈ H 1 (Ω)/ u Γ0 = 0 . We set γ 1 the trace operator from H 1 Γ0 (Ω) on L 2 (Γ 1 ) and H 1/2 (Γ 1 ) = γ 1 H 1 Γ0 (Ω)
. We denote by B the norm of γ 1 namely:

∀u ∈ H 1 Γ0 (Ω) , u 2,Γ1 ≤ B ∇u 2 .
We recall that H 1/2 (Γ 1 ) is dense in L 2 (Γ 1 ) (see [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF]).

We denote E(∆, L 2 (Ω)) = {u ∈ H 1 (Ω) such that ∆u ∈ L 2 (Ω)} and recall that for a function u ∈ E(∆, L 2 (Ω)) , ∂u ∂ν ∈ H -1/2 (Γ 1 ) and the next Green's formula is valid (see [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF]):

Ω ∇u(x)∇v(x)dx = Ω -∆u(x)v(x)dx + ∂u ∂ν ; v Γ1 , ∀v ∈ H 1 Γ0 (Ω), (8)
where .; . Γ1 means the duality pairing between H -1/2 (Γ 1 ) and H 1/2 (Γ 1 ).

By (., .) we denote the scalar product in L 2 (Ω) i.e. (u, v) = Ω u(x)v(x)dx. Also we mean by . q the L q (Ω) norm for 1 ≤ q ≤ ∞, and by . q,Γ1 the L q (Γ 1 ) norm.

Throughout the paper, we use the standard notations as in the book [START_REF] Brezis | Analyse fonctionnelle[END_REF] for example.

In order to prove the local existence of the solution of problem (1), we consider the following two cases :

case 1: µ 2 < µ 1 .
We may define a positive real number ξ such that:

τ µ 2 ≤ ξ ≤ τ (2µ 1 -µ 2 ) . ( 9
)
case 2: µ 2 ≥ µ 1 . We will suppose that the damping parameter α verifies:

α > (µ 2 -µ 1 )B 2 . ( 10
)
In this case, we may define a positive real number ξ satisfying the two inequalities:

ξ ≥ τ µ 2 , (11)
α > µ 2 2 + ξ 2τ -µ 1 B 2 > 0 . (12)

Semigroup formulation of the problem

In this section, we prove the global existence and the uniqueness of the solution of problem (1). To overcome the problem of the boundary delay, we introduce, as in [START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF], the new variable:

z (x, ρ, t) = u t (x, t -τ ρ) , x ∈ Γ 1 , ρ ∈ (0, 1) , t > 0. (13)
Then, we have τ z t (x, ρ, t) + z ρ (x, ρ, t) = 0, in Γ 1 × (0, 1) × (0, +∞) .

Therefore, problem (1) is equivalent to:

                               u tt -∆u -α∆u t = 0, x ∈ Ω, t > 0 , τ z t (x, ρ, t) + z ρ (x, ρ, t) = 0, x ∈ Γ 1 , ρ ∈ (0, 1) , t > 0 , u(x, t) = 0, x ∈ Γ 0 , t > 0 , u tt (x, t) = - ∂u ∂ν (x, t) + α ∂u t ∂ν (x, t) + µ 1 u t (x, t) + µ 2 z(x, 1, t) , x ∈ Γ 1 , t > 0 , z(x, 0, t) = u t (x, t), x ∈ Γ 1 , t > 0 , u(x, 0) = u 0 (x), x ∈ Ω , u t (x, 0) = u 1 (x), x ∈ Ω , z(x, ρ, 0) = f 0 (x, -τ ρ), x ∈ Γ 1 , ρ ∈ (0, 1) . (15)
The first natural question is the existence of solutions of the problem [START_REF] Grobbelaar-Van Dalsen | On the solvability of the boundary-value problem for the elastic beam with attached load[END_REF]. In this section we will give a sufficient condition that guarantees the well-posedness of the problem.

For this purpose, as in [START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF], we will use a semigroup formulation of the initial-boundary value problem [START_REF] Grobbelaar-Van Dalsen | On the solvability of the boundary-value problem for the elastic beam with attached load[END_REF]. If we denote V := (u, u t , γ 1 (u t ), z)

T , we define the energy space:

H = H 1 Γ0 (Ω) × L 2 (Ω) × L 2 (Γ 1 ) × L 2 (Γ 1 × (0, 1)).
Clearly, H is a Hilbert space with respect to the inner product

V 1 , V 2 H = Ω ∇u 1 .∇u 2 dx + Ω v 1 v 2 dx + Γ1 w 1 w 2 dσ + ξ Γ1 1 0 z 1 z 2 dρdσ for V 1 = (u 1 , v 1 , w 1 , z 1 ) T , V 2 = (u 2 , v 2 , w 2 , z 2 )
T and ξ is defined by ( 9) or [START_REF] Doronin | Global solvability for the quasilinear damped wave equation with nonlinear second-order boundary conditions[END_REF]. Therefore, if V 0 ∈ H and V ∈ H , the problem ( 15) is formally equivalent to the following abstract evolution equation in the Hilbert space H :

V ′ (t) = A V (t), t > 0, V (0) = V 0 , (16)
where ′ denotes the derivative with respect to time t, V 0 := (u 0 , u 1 , γ 1 (u 1 ), f 0 (., -.τ ))

T and the operator A is defined by:

A        u v w z        =          v ∆u + α∆v - ∂u ∂ν -α ∂v ∂ν -µ 1 v -µ 2 z (., 1) - 1 τ z ρ          .
The domain of A is the set of V = (u, v, w, z) T such that:

(u, v, w, z) T ∈ H 1 Γ0 (Ω) × H 1 Γ0 (Ω) × L 2 (Γ 1 ) × L 2 Γ 1 ; H 1 (0, 1) , (17)
u + αv ∈ E(∆, L 2 (Ω)) , ∂(u + αv) ∂ν ∈ L 2 (Γ 1 ), (18)
w = γ 1 (v) = z(., 0) on Γ 1 . (19)
The well-posedness of problem ( 15) is ensured by: Theorem 2.1. Suppose that µ 2 ≥ µ 1 and α > (µ 2 -µ 1 )B 2 or µ 2 < µ 1 . Let V 0 ∈ H , then there exists a unique solution V ∈ C (R + ; H ) of problem [START_REF] Grobbelaar-Van Dalsen | Boundary stabilization for the extensible beam with attached load[END_REF]. Moreover, if V 0 ∈ D (A ), then

V ∈ C (R + ; D (A )) ∩ C 1 (R + ; H ) .
Proof. To prove Theorem 2.1, we use the Lumer-Phillips' theorem. For this purpose, we show firstly that the operator A is dissipative. Indeed, let V = (u, v, w, z) T ∈ D (A ). We have

A V, V H = Ω ∇u.∇vdx + Ω v (∆u + α∆v) dx + Γ1 w - ∂u ∂ν -α ∂v ∂ν -µ 1 v -µ 2 z (σ, 1) dσ - ξ τ Γ1 1 0 zz ρ dρdσ. But since u + αv ∈ E(∆, L 2 (Ω)) and ∂(u + αv) ∂ν ∈ L 2 (Γ 1 )
, we may apply Green's formula [START_REF] Grobbelaar-Van Dalsen | On the initial-boundary-value problem for the extensible beam with attached load[END_REF] where the duality pairing .; . Γ1 is simply the L 2 (Γ 1 ) inner product (because w = γ 1 (v) ∈ L 2 (Γ 1 )) and obtain:

A V, V H = -µ 1 Γ1 w 2 dσ -µ 2 Γ1 z (σ, 1) wdσ -α Ω |∇v| 2 dx - ξ τ Γ1 1 0 z ρ zdρdx. ( 20
)
At this point, we have to distinguish the following two cases: Case 1: We suppose that µ 2 < µ 1 . Let us choose then ξ that satisfies inequality [START_REF] Datko | Two questions concerning the boundary control of certain elastic systems[END_REF]. Using Young's inequality, (20) leads to

A V, V H + α Ω |∇v| 2 dx + µ 1 - ξ 2τ - µ 2 2 Γ1 w 2 dσ + ξ 2τ - µ 2 2 Γ1 z 2 (σ, 1) dσ ≤ 0.
Consequently, by using ( 9), we deduce that

A V, V H ≤ 0. ( 21
)
Case 2: We suppose that µ 2 ≥ µ 1 and α > (µ 2 -µ 1 )B 2 . Let us choose then ξ that satisfies the two inequalities [START_REF] Doronin | Global solvability for the quasilinear damped wave equation with nonlinear second-order boundary conditions[END_REF] and [START_REF] Gerbi | Asymptotic stability and blow up for a semilinear damped wave equation with dynamic boundary conditions[END_REF]. Using Young's inequality and the definition of the constant B, we can again prove that the inequality [START_REF] Pellicer | Large time dynamics of a nonlinear spring-mass-damper model[END_REF] holds. This means that in both cases A is dissipative. Now we show that λI -A is surjective for all λ > 0.

For F = (f 1 , f 2 , f 3 , f 4 ) T ∈ H , let V = (u, v, w, z) T ∈ D (A ) solution of (λI -A ) V = F, which is: λu -v = f 1 , (22) λv -∆(u + αv) = f 2 , ( 23
)
λw + ∂(u + αv) ∂ν + µ 1 v + µ 2 z(., 1) = f 3 , (24)
λz + 1 τ z ρ = f 4 . ( 25
)
To find V = (u, v, w, z) T ∈ D (A ) solution of the system ( 22), ( 23), ( 24) and ( 25), we proceed as in [START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF]. Suppose u is determined with the appropriate regularity. Then from [START_REF] Pellicer | Analysis of a viscoelastic spring-mass model[END_REF], we get:

v = λu -f 1 . (26)
Therefore, from the compatibility condition on Γ 1 , ( 19), we determine z(., 0) by:

z(x, 0) = v(x) = λu(x) -f 1 (x), for x ∈ Γ 1 . (27)
Thus, from [START_REF] Trench | Introduction to real analysis[END_REF], z is the solution of the linear Cauchy problem:

z ρ = τ f 4 (x) -λz(x, ρ) , for x ∈ Γ 1 , ρ ∈ (0, 1), z(x, 0) = λu(x) -f 1 (x). ( 28
)
The solution of the Cauchy problem (28) is given by:

z(x, ρ) = λu(x)e -λρτ -f 1 e -λρτ + τ e -λρτ ρ 0 f 4 (x, σ)e λστ dσ for x ∈ Γ 1 , ρ ∈ (0, 1). (29)
So, we have at the point ρ = 1,

z(x, 1) = λu(x)e -λτ + z 1 (x), for x ∈ Γ 1 (30)
with

z 1 (x) = -f 1 e -λτ + τ e -λτ 1 0 f 4 (x, σ)e λστ dσ, for x ∈ Γ 1 . Since f 1 ∈ H 1 Γ0 (Ω) and f 4 ∈ L 2 (Γ 1 ) × L 2 (0, 1), then z 1 ∈ L 2 (Γ 1
). Consequently, knowing u, we may deduce v by [START_REF] Xu | Stabilization of the wave system with input delay in the boundary control[END_REF], z by (29) and using (30), we deduce w = γ 1 (v) by [START_REF] Goldstein | Derivation and physical interpretation of general boundary conditions[END_REF].

In view of equations ( 23) and ( 24), we set, as in [START_REF] Nicaise | Interior feedback stabilization of wave equations with time dependent delay[END_REF], u = u + αv. Then, from (26), we have

v = λu -f 1 = λ(u -αv) -f 1 .
Since λ > 0 and α > 0 , 1 + λα = 0; thus we have:

v = λ 1 + λα u - f 1 1 + λα . (31)
But since u = u -αv, we have:

u = 1 1 + λα u + α 1 + λα f 1 . (32)
From equations ( 23) and ( 24), u must satisfy:

λ 2 1 + λα u -∆u = f 2 + λ 1 + λα f 1 , in Ω (33)
with the boundary conditions

u = 0, on Γ 0 (34) ∂u ∂ν = f 3 -λγ 1 (v) -µ 1 γ 1 (v) -µ 2 z(. , 1), on Γ 1 (35)
the last equation at least formally since we don't have yet found the regularity of u. Replacing u by its expression (32) and inserting it in equation (30), we get:

z(x, 1) = λ 1 + λα u(x)e -λτ + λα 1 + λα f 1 (x)e -λτ + z 1 (x), for x ∈ Γ 1 .
Using the preceding expression of z(., 1) and the expression of v given by (31), we have:

∂u ∂ν = - λ µ 2 e -λτ + (λ + µ 1 1 + λα u + f (x), for x ∈ Γ 1 (36)
with

f (x) = f 3 (x) + (λ + µ 1 ) -µ 2 λαe -λτ 1 + λα f 1 (x) -µ 2 z 1 (x), for x ∈ Γ 1 .
From the regularity of

f 1 , f 2 , z 1 , we get f ∈ L 2 (Γ 1 ).
The variational formulation of problem (33), (34),( 36) is to find u ∈ H 1 Γ0 (Ω) such that:

Ω λ 2 1 + λα uω + ∇u∇ωdx + Γ1 λ µ 2 e -λτ + (λ + µ 1 1 + λα u(σ)ω(σ)dσ, (37)
= Ω f 2 + λ 1 + λα f 1 ωdx + Γ1 f (σ)ω(σ)dσ,
for any ω ∈ H 1 Γ0 (Ω). Since λ > 0 , µ 1 > 0 , µ 2 > 0, the left hand side of (37) defines a coercive bilinear form on H 1 Γ0 (Ω). Thus by applying the Lax-Milgram theorem, there exists a unique u ∈ H 1 Γ0 (Ω) solution of (37). Now, choosing ω ∈ C ∞ c , u is a solution of (33) in the sense of distribution and therefore u ∈ E(∆, L 2 (Ω)). Thus using the Green's formula [START_REF] Grobbelaar-Van Dalsen | On the initial-boundary-value problem for the extensible beam with attached load[END_REF] in (37) and exploiting the equation (33) on Ω, we obtain finally:

Γ1 λ µ 2 e -λτ + (λ + µ 1 1 + λα u(σ)ω(σ)dσ + ∂u ∂ν ; ω Γ1 = Γ1 f (σ)ω(σ)dσ ∀ω ∈ H 1 Γ0 (Ω) .
So u ∈ E(∆, L 2 (Ω)) verifies (36) and by equation ( 32) and (31) we recover u and v and thus by (29), we obtain z and finally setting w = γ 1 (v), we have found V = (u, v, w, z) T ∈ D (A ) solution of (Id -A ) V = F . Thus, the proof of Theroem 2.1, follows from the Lumer-Phillips' theorem.

3 Asymptotic behavior 3.1 Exponential stability for µ 2 < µ 1

In this subsection, we show that under the assumption µ 2 < µ 1 , the solution of problem [START_REF] Grobbelaar-Van Dalsen | On the solvability of the boundary-value problem for the elastic beam with attached load[END_REF] decays to the null steady state with an exponential decay rate. For this goal, we use the energy method combined with the choice of a suitable Lyapunov functional. For a positive constant ξ satisfying the strict inequality ( 9), (i.e. < instead of ≤) we define the functional energy of the solution of problem [START_REF] Grobbelaar-Van Dalsen | On the solvability of the boundary-value problem for the elastic beam with attached load[END_REF] as

E(t) = E(t, z, u) = 1 2 ∇u(t) 2 2 + u t (t) 2 2 + u t (t) 2 2,Γ1 + ξ 2 Γ1 1 0 z 2 (σ, ρ, t) dρ dσ = 1 2 E 1 (t) + ξ 2 Γ1 1 0 z 2 (σ, ρ, t) dρ dσ, (38)
where E 1 (t) = ∇u(t) 2 2 + u t (t) 2 2 + u t (t) 2 2,Γ1 . Let us first remark that this energy is greater than the usual one of the solution of problem [START_REF] Andrews | Second order evolution equations with dynamic boundary conditions[END_REF], namely E 1 (t). Now, we prove that the above energy E (t) is a decreasing function along the trajectories. More precisely, we have the following result: Lemma 3.1. Assume that µ 1 > µ 2 , then the energy defined by ( 38) is a non-increasing positive function and there exists a positive constant C such that for (u, z) solution of [START_REF] Grobbelaar-Van Dalsen | On the solvability of the boundary-value problem for the elastic beam with attached load[END_REF], and for any t ≥ 0, we have:

dE (t) dt ≤ -C Γ1 u 2 t (σ, t) dσ + Γ1 z 2 (σ, 1, t) dσ -α Ω |∇u t (x, t)| 2 dx . (39)
Proof. We multiply the first equation in ( 15) by u t and perform integration by parts to get:

1 2 d dt ∇u(t) 2 2 + u t (t) 2 2 + u t (t) 2 2,Γ1 + α ∇u t (t) 2 2 +µ 1 u t (t) 2 2,Γ1 + µ 2 Γ1 u t (σ, t)u t (σ, t -τ )dσ = 0 . ( 40
)
We multiply the third equation in [START_REF] Grobbelaar-Van Dalsen | On the solvability of the boundary-value problem for the elastic beam with attached load[END_REF] by ξz, integrate the result over Γ 1 × (0, 1), we obtain:

ξ τ Γ1 1 0 z ρ z(σ, ρ, t) dρ dσ = ξ 2τ Γ1 1 0 ∂ ∂ρ z 2 (σ, ρ, t) dρ dσ = ξ 2τ Γ1 z 2 (σ, 1, t) -z 2 (σ, 0, t) dσ . (41)
Using the definition (13) of z in the equality (40) and using the same technique as in the first step of the proof of Theorem 2.1, where we proved that A is dissipative, inequality (39) holds.

The asymptotic stability result reads as follows:

Theorem 3.1. Assume that µ 2 < µ 1 . Then there exist two positive constants C and γ independent of t such that for (u, z) solution of problem ( 15), we have:

E(t) ≤ Ce -γt , ∀ t ≥ 0 . (42)
Proof. The proof of Theorem 3.1 relies on the construction of a Lyapunov functional. For a small positive constant ε to be chosen later, we define:

L(t) = E(t) + ε Ω u(x, t)u t (x, t) dx + ε Γ1 u(σ, t)u t (σ, t) dσ + εα 2 Ω |∇u(x, t)| 2 dx (43) + εξ Γ1 1 0 e -2τ ρ z 2 (σ, ρ, t) dρ dσ.
Let us say that the introduction of the last term in the Lyapunov functional L is inspired by the work of Nicaise and Pignotti [START_REF] Nicaise | Interior feedback stabilization of wave equations with time dependent delay[END_REF].

It is straightforward to see that for ε > 0, L(t) and E(t) are equivalent in the sense that there exist two positive constants β 1 and β 2 > 0 depending on ε such that for all t ≥ 0

β 1 E(t) ≤ L(t) ≤ β 2 E(t) . (44)
By taking the time derivative of the function L defined by (43), using the equations in problem [START_REF] Grobbelaar-Van Dalsen | On the solvability of the boundary-value problem for the elastic beam with attached load[END_REF], several integration by parts, and exploiting (39), we get:

dL(t) dt ≤ -C Γ1 u 2 t (σ, t) dσ + Γ1 z 2 (σ, 1, t) dσ -α ∇u t 2 2 -ε ∇u 2 2 + ε u t 2 2 + ε u t 2 2,Γ1 -εµ 1 Γ1 u t (σ, t)u(σ, t) dσ -εµ 2 Γ1 z(σ, 1, t)u(σ, t)dσ (45)
+ εξ d dt Γ1 1 0 e -2τ ρ z 2 (σ, ρ, t) dρ dσ .

By using the second equation in [START_REF] Grobbelaar-Van Dalsen | On the solvability of the boundary-value problem for the elastic beam with attached load[END_REF], the last term in (45) can be treated as follows:

εξ d dt Γ1 1 0 e -2τ ρ z 2 (σ, ρ, t) dρ dσ = - 2εξ τ Γ1 1 0 e -2τ ρ z(σ, ρ, t)z ρ (σ, ρ, t) dρ dσ = - εξ τ Γ1 1 0 e -2τ ρ ∂ ∂ρ z 2 (σ, ρ, t) dρ dσ.
Then, by using an integration by parts and the definition of z, the above formula can be rewritten as:

εξ d dt Γ1 1 0 e -2τ ρ z 2 (σ, ρ, t) dρ dσ = - εξ τ e -2τ Γ1 z 2 (1, ρ, t) dρ dσ + εξ τ Γ1 u 2 t (σ, t) dσ (46) -2εξ Γ1 1 0 e -2τ ρ z 2 (σ, ρ, t) dρ dσ .
Proof of Theorem 3.2. We use the same Lyapunov function as in the previous section, namely, for a small positive constant ε to be chosen later, we define: e -2τ ρ z 2 (σ, ρ, t) dρ dσ.

L(t) = E(t) + ε Ω u(x,
By taking the time derivative of the function L, using the equations in problem [START_REF] Grobbelaar-Van Dalsen | On the solvability of the boundary-value problem for the elastic beam with attached load[END_REF], several integration by parts, and exploiting (52), we get: e -2τ ρ z 2 (σ, ρ, t) dρ dσ .

dL(t) dt ≤ -κ Γ1 z 2 (
The remaining part of the proof is similar to the one of the proof of Theorem 3.1: by choosing firstly δ and then ε, we may find γ > 0 independent of t such that: dL(t) dt ≤ -γL(t) , ∀t ≥ 0 .

This inequality permits us to conclude the proof of Theorem 3.2.

Remark 3.2. After our work had been submitted, we noticed that a similar problem has been already studied by Nicaise and Pignotti [START_REF] Nicaise | Interior feedback stabilization of wave equations with time dependent delay[END_REF]. However our problem is slightly different from the one considered in [START_REF] Nicaise | Interior feedback stabilization of wave equations with time dependent delay[END_REF]:

• First, we have the extra term µ 1 u t on the boundary conditions, which makes the stability analysis independent of the strong damping α∆u t , for µ 1 > µ 2 . This is not the case in [START_REF] Nicaise | Interior feedback stabilization of wave equations with time dependent delay[END_REF]. See the assumption (2.56) in [START_REF] Nicaise | Interior feedback stabilization of wave equations with time dependent delay[END_REF]. Moreover, the paper [START_REF] Nicaise | Interior feedback stabilization of wave equations with time dependent delay[END_REF] does not cover the case α = 0 at all.

• Secondly, our Lyapunov functional is different from the one used in [START_REF] Nicaise | Interior feedback stabilization of wave equations with time dependent delay[END_REF]. We choose to use this one to have a strong control of the boundary delay term.

Finally, let us remark that our assumption (µ 2 -µ 1 )B 2 < α is exactly the same than the one obtained in [20, Condition (2.56)], when µ 1 = 0. So this work can be viewed as a continuation of the works of Nicaise and Pignotti [START_REF] Nicaise | Interior feedback stabilization of wave equations with time dependent delay[END_REF]

  t)u t (x, t) dx + ε Γ1 u(σ, t)u t (σ, t) dσ

	+	εα 2 Ω	|∇u(x, t)|	2 dx
			1
	+ εξ
		Γ1	0

  σ, 1, t) dσ + -2τ ρ z 2 (σ, ρ, t) dρ dσ .Let us remark that equality (46), and inequalities (47) and (48) used in the the proof of Theorem 3.1 are still valid. Inserting (46), (47) and (48) into (56) and using Poincaré's inequality, we have:

						|∇u t (x, t)|	2 dx
					Ω
	-α ∇u t	2 2 -ε ∇u 2 2 + ε u t	2 2 + ε u t	2 2,Γ1
	-εµ 1	u t (σ, t)u(σ, t) dσ -εµ 2	z(σ, 1, t)u(σ, t)dσ	(56)
	Γ1						Γ1
	+ εξ e dL(t) 1 d dt Γ1 0 dt ≤ -κ -ε ξ τ	e -2τ +	µ 2 4δ	Γ1	z 2 (σ, 1, t) dσ
	-κ -ε B 2 + 1 +	B 2 ξ τ	+	µ 1 4δ	∇u t	2 2
	-ε 1 -C(Ω) 2 δ (µ 1 + µ 2 ) ∇u 2 2	(57)
			1
	-2εξ
		Γ1	0

Acknowledgments. The second author was supported by MIRA 2007 project of the Région Rhône-Alpes. This author wishes to thank Univ. de Savoie of Chambéry for its kind hospitality. Moreover, the two authors wish to thank the referees for their remarks and the careful reading of the proofs presented in this paper.

Applying Young's inequality, and the trace inequality, we obtain, for any δ > 0:

We choose now δ small enough in (49) such that

.

Once δ is fixed, using once again Poincaré's inequality in (49), we may pick ε small enough to obtain the existence of η > 0, such that:

On the other hand, by virtue of (44), setting γ = -ηε/β 2 , the last inequality becomes:

Hence, integrating the previous differential inequality (51) between 0 and t, we get

for some positive constant C * . Consequently, by using (44) once again, we conclude that it exists C > 0 such that:

This completes the proof of Theorem 3.1 .

Exponential stability for

As, we have said in the Introduction, and it is clearly observed in Theorem 3.1, that the strong internal damping compensates the destabilizing effect of the delay in the boundary condition.

In this section, we assume that µ 2 > µ 1 and α > (µ 2 -µ 1 )B 2 . As we will see, we cannot directly perform the same proof as for the case where µ 2 ≤ µ 1 , since the boundary delay term µ 2 u t (x, t -τ ) is greater than the normal one µ 1 u t (x, t), i.e. (µ 2 ≥ µ 1 ). So we have to control this term by the damping term α∆u t in the equation.

Remark 3.1. In the case 2, namely µ 2 > µ 1 the condition α > (µ 2 -µ 1 )B 2 permits us to find ξ satisfying (11)-( 12). This choice of ξ is essential in the proofs of Lemma 3.2 and Theorem 3.2 below. Lemma 3.2. Assume that µ 2 > µ 1 and α > (µ 2 -µ 1 )B 2 . For any ξ satisfying ( 11)-( 12), the energy defined by ( 38) is a non-increasing positive function and there exists a positive constant κ such that for (u, z) solution of ( 15), and for any t ≥ 0, we have:

Proof. Let us first recall from ( 40) and (41) the following identity

Now, using Young's inequality, then (53) takes the form:

Since,

then, using the trace inequality, we obtain:

Using the two inequalities ( 11)-( 12) that ξ satisfy , we may find κ > 0 such that the inequality (52) holds.

We can now state that under the same assumption as in Lemma 3.2, the system (15) is also exponentially stable. The second stability result reads as follows: Theorem 3.2. Assume that µ 2 > µ 1 and α > (µ 2 -µ 1 )B 2 . For any ξ satisfying (11)- [START_REF] Gerbi | Asymptotic stability and blow up for a semilinear damped wave equation with dynamic boundary conditions[END_REF], there exist two positive constants C and γ independent of t such that for (u, z) solution of problem [START_REF] Grobbelaar-Van Dalsen | On the solvability of the boundary-value problem for the elastic beam with attached load[END_REF], we have:

(55)