Internal feedback stabilization of multi-dimensional wave equations with a boundary delay: a numerical study

In this paper, we consider two internal stabilization problems for the multi-dimensional wave equation with a boundary time-delay. We prove that the first problem is well-posed in an appropriate functional space. Subsequently, we numerically study the exponential stability in a two-dimensional case under Geometric Control Condition (GCC) derived by Lebeau. In addition, we provide a numerical investigation of the second wave system, which corresponds to the two-dimensional variant of the system studied by Datko et al.

Delay effects arise in many applications and practical problems. In particular, time delays are ubiquitous in control systems since sensors and actuators are commonly used. Based on numerous articles, it has been noticed that even an arbitrarily small delay may destabilize a system which was uniformly asymptotically stable in the absence of a time-delay (see, e.g., [20,21]). Therefore, it is crucial to investigate the effect of a delay either on the stability of the system or on the performance of the controller. Several articles appeared in this research avenue aiming to eliminate the delay or compensate it, or at least reduce its negative impact. In contrast to this predominant view, a number of unstable systems can be stabilized via a purposeful introduction of a time-delay in the systems [11,24], while others are not affected by the presence of a delay [23].
Now, let us go back to the system under consideration (1.1)-(1.5). First, we consider the boundary control U(u t (x, tτ ) = ∂Gu t ∂ν (x, tτ ), (x, t) ∈ 1 × R + , where ν stands for the unit normal vector of pointing towards the exterior of , and ∂ ∂ν is the normal derivative and in which G = (-) -1 : H -1 ( ) − → H 1 0 ( ). This gives rise to our first problem: The exponential stability of (1.6)-(1.10) with k = 0 has been studied in [28], where it has been shown that the system is exponentially stable if the support of a satisfies some geometric control condition (GCC). In turn, if τ = 0 and k < 0, that is, in the absence of any time-delay, then the above system (1.1)-(1.5) is exponentially stable when the support of a or 1 satisfies a geometrical control condition (GCC) (see [9,10,27,28]) even if k = 0 (see, e.g., [28] and [36]). Nevertheless, if a delay term occurs in the system, then instability phenomena may arise as shown in [31] for Neumann feedback boundary control.
In order to contrast the destabilizing effect of the boundary time delay in a wave equation, a "good" (undelayed) damping term is introduced in (1.6) (see [31]). More precisely, the problem considered in [31] is with a 0 , a > 0 are real constants, and the initial data belong to suitable spaces. If a 0 > a, it was shown in [31] that system (1.11)-(1.15) is uniformly exponentially stable (see also [1,3,4,6,7,12,22,32,33] for related results and similar situations).
On the other hand, we note that studying the system (1.16)-(1.20) is mainly motivated by the work [21], where one can find an interesting analysis conducted for the onedimensional case of the system (1.16)-(1.20). Indeed, a positive answer to the above problem (concern) is provided in [21]. Specifically, the authors in [21] considered the following one-dimensional system where b and c are positive real numbers. Moreover, the initial data and the function f are in appropriate functional spaces. Through careful spectral analysis, the authors have shown in [21] that, for any a > 0, and if c satisfies then the spectrum of the system (1.21)-(1.23) satisfies the following property: Re(ω) ≤ -β for each eigenvalue ω of the system and some positive constant β depending on the delay τ . Notwithstanding, as far as we know, there are no results in the literature for the multi-dimensional case described by (1.16)- (1.20). This has motivated us to numerically solve this problem in the case when the space variable x belongs to R 2 .
Finally, we would like to point out that one variant of the one-dimensional version of the system (1.6)-(1.10), where a = f = 0 and for τ = 2, has the following form: where k is a real number. Arguing as in [5] (see the appendix), one can obtain the exponential stability of the above system for (u 0 , u 1 ) ∈ L 2 (0, 1) × H -1 (0, 1) provided that k ∈ (-1, 0). The paper is organized as follows: The second section is devoted to the well-posedness of the problem (1.6)-(1.10). In the third section, we numerically study in R 2 the exponential stability of both delayed systems (1. Next, we define a bounded linear operator B ∈ L(H -1 ( )) as follows: Then, let C ∈ L(L 2 ( 1 ); H -1 2 ) such that 2 )) (the duality is in the sense of H), while A -1 is the extension of A to H, namely, for all h ∈ H and ϕ ∈ D(A), A -1 h is the unique element in H -1 = (D(A)) (the duality is in the sense of H) such that (see, for instance, [35]) Here and in the sequel, D ∈ L(L 2 ( 1 ); L 2 ( )) (called the Dirichlet mapping), ∀v ∈ L 2 ( 1 ), Dv is the unique solution (transposition solution) of To study the well-posedness of the systems (1.6)-(1.10), we write it as an abstract Cauchy problem in a product space and use the semigroup approach. For this purpose, take the Hilbert space H := H 1 2 × H and the unbounded linear operators It is well known that the operator (A, D(A)) defined by (2.1), generates a strongly continuous semigroup of contractions on H denoted by (T (t)) t≥0 . We also denote (T -1 (t)) t≥0 the extension of (T (t)) t≥0 to (D(A)) : 2 )) , and the duality is in the sense of H.
We have the following result: Theorem 2. 1 The system (1.6)-(1.10) is well-posed. More precisely, for every (u 0 , u 1 ) ∈ H, and f ∈ L 2 (-τ , 0; L 2 ( 1 )), the solution of (1.6)-(1.10) is given by The proof requires several steps (see [6,8] for other types of problems). First, consider the evolution problems and A natural question is the regularity of y j when v j ∈ L 2 (τ j, (j + 1)τ ; L 2 ( 1 )), j ∈ N * . By applying standard energy estimates, we can easily check that ). However, if C satisfies a certain admissibility condition, then y j is more regular. More precisely, the following result, which is a version of the general transposition method (see, for instance, Lions and Magenes [29]) holds true: The system (2.4)-(2.5) admits a unique solution φ having the regularity Moreover, according to [2], C * φ(·) ∈ H 1 (0, τ ; L 2 ( 1 )) and for all t ∈ (0, τ ), there exists a constant C > 0 such that Next, we have the following lemma: admits a unique solution having the regularity and 3) can be written as After simple calculations, one can check that the operator C * : D(A) → L 2 ( 1 ) is given by This implies that in which φ satisfies (2.4)-(2.5). Amalgamating the last claim with (2.6), we deduce that there exists a constant C > 0 such that for all T ∈ (0, τ ) In view of [13, Theorem 3.1, p.187], the last estimate leads to the interior regularity (2.7) (see also [35]).

Numerical study
Numerical solutions for the multi-dimensional wave systems (1.6)-(1.10) as well as (1.16)-(1.20) with and without the presence of a boundary time-delay are simulated using COM-SOL Multiphysics software. The domain is taken to be the square (0, 1) × (0, 1) (see Fig. 1), while different values of k and various functions a(x), x = (x 1 , x 2 ) T are considered. It should also be noted that for the numerical experiments we assume that 0 ∩ 1 = ∅, which is different from the assumption in the theoretical discussion.
The finite element method with the quadratic Lagrange shape functions as basis functions is used with a fixed fine element mesh size, where the minimum element mesh size equals to 1.25 × 10 -4 , and the maximum element mesh size is 0.037, and the generalized-α method as the time-stepping scheme is used with dt = 0.001 (see Fig. 1). The generalizedα method was introduced by Chung and Hulbert [19] to solve the following second-order ODEs: where A, B, and C are the mass, damping, and stiffness matrices, respectively, H(t) is the vector applied loads, and Y (t) is the displacement vector, where the superimposed dots indicate differentiation with respect to time; d and v are given vectors of initial displacements and velocities, respectively. By approximating Y (t i ),Ẏ (t i ), andŸ (t i ) by d i , v i , and a i , respectively, the generalized-α algorithm is given by:  (3.11) in which i ∈ [0, . . . , N -1], N is the number of time steps, and t = t i+1t i is the time step size. The algorithm parameters α f , α m , β, and γ are used to control the degree of damping high frequencies while minimizing unwanted low-frequency dissipation. The values of the parameters α f , α m , β, and γ depend on the user-specified high-frequency dissipation ρ ∞ .
It should be noted that ρ ∞ was taken to be 0.75 in all numerical simulations conducted herein. Therefore, the values of the parameters α f , α m , β, and γ are 0.428, 0.285, 0.326, and 0.642, respectively.

The 2-d wave system without a boundary time-delay
In this subsection, we consider the 2-d wave problem (1.6)-(1.10) without delay. Letting x = (x 1 , x 2 ) T and taking the initial data u 0 (x) = 0 and u 1 (x) = sin(πx 1 ), the system is: First, we choose a(x) = 1 and vary the values of k. It has been observed that the dynamics of the system (3.13)-(3.16) is exponentially stable if k ≤ 0 and is unstable for k > 0.   Figure 3b is a plot of the energy E 0 (t) for different negative values of k, indicating that the energy exponentially converges faster as the value of k decreases. On the other hand, Fig. 3c shows that the dynamics of the system is unstable when k = 0.1.

The 2-d wave equation (1.6)-(1.10) with a boundary time-delay
In this subsection, we go back to the 2-d wave equation (1.6)-(1.10) with a boundary time-delay τ . The inial data are taken as follows: u 0 (x) = 0, u 1 (x) = sin(πx 1 ), and f (x, t) = sin(πx 1 ). The system is simulated for different parameters k, functions a(x), and time-delays τ . First, we choose the same value of k = -1 and function a(x) = 1 as in the case without delay presented in Sect. 3.1.1, but with time-delay τ = 1. Figure 4 presents the dynamics of u(x, t) at different time, and Fig. 5 shows the corresponding energy, E 0 (t), versus time.   Figure 7b shows that the dynamics stabilizes as the value of k increases from -0.2 to -0.1. On the other hand, if the values of k increase from 0.1 to 0.2, the dynamics becomes more unstable (see Fig. 7c). These observations indicate that if a(x) = 1 and |k| < 0.2, then the 2-d wave equation with a boundary time-delay τ = 1 is stable. However, this is not true if we change the function a(x) = 3, where it is shown that the dynamics becomes exponentially stable (see Fig. 7d). Figure 7e indicates similar observations to Fig. 7d when the time-delay is  Figure 8a shows that as we increase the value of the time-delay τ while fixing k = -0.2 and a(x) = 1, the system destabilizes, and the rate of divergence of the energy increases as the time-delay τ increases. Similar observations are deduced for the case k = 0.2 and a(x) = 1 (see Fig. 8b). However, increasing the value of a(x) = 3 exponentially stabilizes the dynamics for the previous two cases: k = -0.2, and k = 0.2 (see Figs. 8c and 8d).
Then, we let k = 1 and modify a(x). Figure 11d shows the convergence of the energy for various functions of a(x) (i.e., a(x) = 1, a(x) = x 1 , and a(x) = e x 1 ). The figure shows that faster convergence is achieved when a(x) = e x 1 , and this is due to the fact that e x 1 > 1 > x 1 for all x 1 ∈ [0, 1].  Figure 15b shows that the dynamics becomes more unstable as the value of k increases from 0.1 to 0.4. On the other hand, if the values of k increase from -0.4 to -0.1, the dynamics becomes more stable (see Fig. 15c). These observations indicate that if a(x) = 1 and |k| < 0.3, then the 2-d wave equation with a boundary delay τ = 1 is stable. However, this is not true if we change the function a(x) = 3, where it is shown that the dynamics becomes more exponentially stable (see Fig. 15d). Figure 15e indicates similar observations to Fig. 7e when the delay is increased to τ = 2. Now, to investigate the effect of the choice of the delay τ on the stability of the system (1.16)-(1.20), we fix k and a(x) and vary τ . Figure 15f shows that as we increase the value of the delay τ while fixing k = -0.4 and a(x) = 1, the system is still unstable. However, the rate of divergence of the energy decreases as the delay τ increases. Similar observations are deduced for the cases k = 0.4 and a(x) = 1 (see Fig. 15g) and k = 1 and a(x) = 1 (see Fig. 15h).