Benjamin Boutin

Thi Hoai

Thuong Nguyen

Nicolas Seguin

Thị Hoài

AND Thương Nguyễn

A

Keywords: AMS classification: 35F46, 35L50, 65M06, 65M12 hyperbolic relaxation system, damped wave equation, stiff stability, discrete transparent boundary condition, central schemes, implicit scheme, energy estimates, Z-transform

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1. Introduction 1.1. Context and motivation. Hyperbolic systems of partial differential equations with relaxation terms [START_REF] Benzoni-Gavage | Multi-dimensional hyperbolic partial differential equations, First-order systems and applications[END_REF] are important in many physical situations, such as kinetics theories [START_REF] Cercignani | The Boltzmann Equation[END_REF], gases not in thermodynamic equilibrium [START_REF] Vincenti | Introduction to Physical Gas Dynamics[END_REF], phase transitions with small transition time [START_REF] Landau | Statistical Physics[END_REF], water waves [START_REF] Stoker | Water waves. The mathematical theory with applications[END_REF][START_REF] Whitham | Linear and nonlinear waves[END_REF], reactive flows [START_REF] Colella | Theoretical and Numerical Structure for Reacting Shock Waves[END_REF], river flows traffics flows and more general waves [START_REF] Whitham | Linear and nonlinear waves[END_REF]. The study of the zero relaxation limit for such systems has caught much interest, both from a theoretical and numerical point of view, after the works of Liu [START_REF] Liu | Hyperbolic conservation laws with relaxation[END_REF], Chen, Levermore and Liu [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF], Hanouzet and Natalini [START_REF] Hanouzet | Global Existence of Smooth Solutions for Partially Dissipative Hyperbolic Systems with a Convex Entropy[END_REF], Yong [START_REF] Yong | Boundary conditions for hyperbolic systems with stiff source terms[END_REF][START_REF] Yong | Entropy and Global Existence for Hyperbolic Balance Laws[END_REF]. The major issue in the theory of the relaxation approximations to equilibrium system of conservation laws is the appearance of stiff boundary layers in the presence of physical or numerical boundary conditions due to the additional characteristic speeds introduced in the relaxation systems. In this article, we are concerned the development of grid algorithms for solving initial boundary value problem, which involves the question of formulating boundary conditions to get a stable scheme. Due to the presence of boundary layers and to the possible interaction of the boundary and initial layers, numerical schemes have to be properly designed so as to provide accurate approximations and consistent behaviors. One of the simplest linear hyperbolic systems with relaxation is linear damped wave equation in one space dimension (1.1)

∂ t u ε (x, t) + ∂ x v ε (x, t) = 0, ∂ t v ε (x, t) + a∂ x u ε (x, t) = -ε -1 v ε (x, t),
where a > 0 and the relaxation time ε > 0 may be introduced to characterize the stiffness of the relaxation. When ε goes to zero, the model may be simplified. We expect indeed that for any position x and time t, the solution (u ε , v ε )(x, t) tends to some (u(x), 0), which is the solution of the corresponding equilibrium system [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Xin | Stiff Well-Posedness and Asymptotic Convergence for a Class of Linear Relaxation Systems in a Quarter Plane[END_REF].

In order to determine a unique solution to the problem (1.1) in the quarter plane x > 0, t > 0, the values of the solution at initial time are prescribed

u ε (x, 0) = u 0 (x), v ε (x, 0) = v 0 (x). (1.2)
In some cases, the suitable boundary condition comes from physical considerations. At a solid wall that bounds the flow of a fluid, for example, one sets the normal component of the fluid velocity equal to zero (if the effects of viscosity are considered, the tangential component must also vanish). In other situations, Date: May 7, 2020. The authors B. Boutin and T.H.T. Nguyễn were partially supported by ANR project NABUCO, ANR-17-CE40-0025.

1 the choice of boundary conditions is not so obvious. This is the case when considering artificial boundaries, which do not correspond to a well-identifies physical phenomenon. In general, not any boundary condition is suitable for a given hyperbolic problem. In the case of the problem (1.1), which is a particular case of the Jin-Xin relaxation model in one space dimension [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF], the solution at the boundary is imposed B u u ε (0, t) + B v v ε (0, t) = b(t), (1.3) where B u and B v are constants. For simplicity, we also assume the initial data f (x) = u 0 , v 0 (x) and the boundary data b(t) to be compatible at the space-time corner x = 0, t = 0, i.e.

f (0) = f (0) = 0, b(0) = b (0) = 0.
It is easy to see that the hyperbolic structure is related to the Riemann invariants √ au ε ± v ε and to the characteristic velocities ± √ a. Therefore, the boundary condition (1.3) has to satisfy the Uniform Kreiss Condition (UKC)

B u + √ aB v = 0. (1.4)
Only under this assumption, the incoming flow √ au + v at the boundary x = 0 can be deduced from the outgoing flow √ au ε -v ε and the data b(t). Therefore, the initial boundary value problem (IBVP) (1.1)-(1.3) is well-posed for each fixed ε (see [START_REF] Benzoni-Gavage | Multi-dimensional hyperbolic partial differential equations, First-order systems and applications[END_REF][START_REF] Xin | Stiff Well-Posedness and Asymptotic Convergence for a Class of Linear Relaxation Systems in a Quarter Plane[END_REF][START_REF] Yong | Boundary conditions for hyperbolic systems with stiff source terms[END_REF]).

In [START_REF] Xin | Stiff Well-Posedness and Asymptotic Convergence for a Class of Linear Relaxation Systems in a Quarter Plane[END_REF], Xin and Xu study the asymptotic equivalence of a general linear system of one-dimensional conservation laws and the corresponding relaxation model proposed by Jin and Xin [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF] in the limit of a small relaxation rate ε. The main issue is to extend and precise this asymptotic equivalence in the presence of physical boundaries. Within the same problematic, Yong in [START_REF] Yong | Boundary conditions for hyperbolic systems with stiff source terms[END_REF] proposed a Generalized Kreiss Condition (GKC) for general multi-dimensional linear constant coefficient relaxation systems, or one-dimensional nonlinear systems, with non-characteristic boundaries. This condition enables uniform stability estimates and a reduced boundary condition for the corresponding equilibrium system. For the special Jin-Xin system (1.1) with boundary condition (1.3) but with stiff source terms of the form ε -1 (λu ε -v ε ) for some λ ∈ R, Xin and Xu identify and rigorously justify a necessary and sufficient condition (which they call the Stiff Kreiss Condition, or SKC in short) on the boundary condition to guarantee the uniform well-posedness of the IBVP, independently of the relaxation parameter. In addition to the work in [START_REF] Yong | Entropy and Global Existence for Hyperbolic Balance Laws[END_REF], their study also covers the characteristic case and provides optimal asymptotic expansions for the limit process, handling with boundary and/or initial layers. In the case of our system (1.1), the parameter λ = 0 so that the boundary is characteristic for limit equation, and the SKC simply reduces to

B v = 0, or B u B v / ∈ - √ a, 0 .
It is classical to notice that, by linearity, the IBVP (1.1)-(1.3) can be broken up into two simpler problems, one with homogeneous initial condition

(1.5)      ∂ t U (x, t) + A∂ x U (x, t) = ε -1 SU (x, t), x > 0, t > 0, U (x, 0) = 0, x > 0, BU (0, t) = b(t), t > 0,
and the other with homogeneous boundary condition

(1.6)      ∂ t U (x, t) + A∂ x U (x, t) = ε -1 SU (x, t), x > 0, t > 0, U (x, 0) = f (x), x > 0, BU (0, t) = 0, t > 0,
where

U (x, t) = u(x, t) v(x, t) , A = 0 1 a 0 , S = 0 0 0 -1 , B = B u B v .
The main difficult part of the proof of stability under the SKC appears with the study of (1.6). This is due to the complex interactions between the initial data, the boundary condition and the stiff relaxation term. By means of the energy method, Xin and Xu [START_REF] Xin | Stiff Well-Posedness and Asymptotic Convergence for a Class of Linear Relaxation Systems in a Quarter Plane[END_REF] show that the time evolution of the energy

U (t) L 2 (R + ,R 2 ) decreases if the boundary condition (1.3) satisfies B u B v > 0,
which is only the subclass of the SKC, so that this analysis is not satisfactory. Actually, under the SKC, they also find explicitly the solution U (x, t) of the IBVP (1.6) by using the Laplace transform. The solution is decomposed into two parts, by assembling the solution for the case of the Cauchy problem without boundary and another one for the case of the IBVP but with homogeneous initial condition [START_REF] Xin | Stiff Well-Posedness and Asymptotic Convergence for a Class of Linear Relaxation Systems in a Quarter Plane[END_REF]Section 5]. It means that the issue of stiff well-posedness of the IBVP (1.1)-(1.3) is solved if we can prove the Cauchy problem and the IBVP with homogeneous initial condition, that is to say (1.5), themselves are stiffly well-posed.

The motivation of the present study is to analyze the counterpart of the above results for the fully discrete difference approximation of the IBVP (1.1)- (1.3). Because of the propagation of waves of the problem (1.1) with the characteristic velocities ± √ a, one of the waves at the boundary x = 0 is coming into the computation domain while the other one is going out of it. Thus, the way of formulating boundary conditions for the relaxation systems so as to guarantee the uniform stability and to minimize the artificial boundary layer is a crucial issue to the success of the schemes. In [START_REF] Boutin | A stiffly stable semi-discrete scheme for the characteristic linear hyperbolic relaxation with boundary[END_REF], the boundary is approximated using a summation-by-parts method. Using energy estimates and Laplace transforms, the semi-discrete approximation of the IBVP (1.1)-(1.3) is proved to be stiffly stable if the boundary condition (1.3) satisfies B u B v > 0, which is only the subclass of the SKC. In this article, we consider the discrete transparent boundary technique to construct a stiffly stable boundary condition. The technique and its analysis has been proposed by Arnold and Ehrhardt in [START_REF] Arnold | Discrete transparent boundary conditions for the Schrödinger equation: fast calculation, approximation, and stability[END_REF]. Besse, Noble and co-authors apply the tools to dispersive problems [START_REF] Besse | Discrete artificial boundary conditions for the linearized Korteweg-de Vries equation. Numerical Methods for Partial Differential Equations[END_REF][START_REF] Besse | Discrete transparent boundary conditions for the mixed KDV-BBM equation[END_REF][START_REF] Kazakova | Discrete Transparent Boundary Conditions for the Linearized Green-Naghdi System of Equations[END_REF]. We also refer the reader to [START_REF] Grote | Exact Nonreflecting Boundary Conditions for the Time Dependent Wave Equation[END_REF] for non-reflecting methods in the context of wave problems. The recent work of Coulombel [START_REF] Coulombel | Transparent numerical boundary conditions for evolution equations: derivation and stability analysis[END_REF] proposes a systematic study of transparent boundary conditions for evolution problems. Our aim here is to prove that the SKC derived in [START_REF] Xin | Stiff Well-Posedness and Asymptotic Convergence for a Class of Linear Relaxation Systems in a Quarter Plane[END_REF] is a sufficient condition for the stiff stability of the proposed fully discrete of the IBVP (1.1)-(1.3).

1.2. Description of the numerical scheme. Let ∆t > 0 being the time step. The space step ∆x > 0 will always be chosen so that the parameter λ xt = ∆x∆t -1 is kept fixed. Letting now U n j = (u n j , v n j ) T denotes the approximation of the exact solution to (1.1)-(1.3) at the grid point (x j , t n ) = (j∆x, n∆t), for any (j, n) ∈ N × N (where we omit the explicit dependence on ε). We focus in this paper on the fully discrete approximation of the IBVP (1.1)-(1.3) obtained by the central differencing scheme in space and the implicit scheme in time.

A first step towards the fully discrete approximation of the IBVP (1.1)-(1.3) is the following system

(1.7)          U n+1 j -U n j ∆t + 1 2∆x U n+1 j+1 -U n+1 j-1 = 1 ε SU n+1 j , j ≥ 1, n ≥ 0, U 0 j = f j , j ≥ 0, BU n 0 = b n , n ≥ 0,
where the approximations of the initial condition f j and of the boundary data b n are defined for example by setting

f j = f (j∆x) for j ≥ 0 and b n = b(n∆t) for n ≥ 0.
Let us emphasize that the numerical scheme (1.7) still needs one more scalar equation at the boundary point j = 0 so as to be fully defined, due to the fact that the matrix B has rank one only. This is actually a discrete feature only, since in the continuous case this single equation is exactly complemented by the only incoming characteristic (at least under UKC). An additional relation to define U n+1 0 is thus needed. We want to use the central scheme at the boundary point j = 0, so that the modification of the ghost value U n+1 -1 can also be interpreted as the use of an extra boundary condition. From a mathematical point of view, the problem is set, in both cases, as follows: given an initial data compactly supported, one can construct boundary condition at j = 0 with the objective to approximate the exact solution of the whole space problem {j ∈ Z}, restricted to {j ∈ N}. If the approximate solution on {j ∈ N} coincides with the exact solution, one refers to these boundary conditions as transparent boundary conditions. Of course, these boundary condition should lead to a well-posed initial boundary value problem. It means that we use the discrete transparent boundary condition at j = 0 that determines a ghost value U n+1 -1 through the identity

U n+1 -1 = n+1 k=0 C n+1-k U k 0 ,
where the coefficients C k will be precised explicitly in the forthcoming Definition 2.4. The extra boundary condition determines U n+1 -1 as a linear function of U k 0 for past step times only: 0 ≤ k ≤ n + 1. We propose the following numerical approximation at the boundary:

1 ∆t Γ U n+1 0 -U n 0 + 1 2∆x ΓA U n+1 1 - n+1 k=0 C n+1-k U k 0 = 1 ε ΓSU n+1 0
with the matrix Γ = (-aB v B u ). Under the SKC, this choice for the matrix Γ will be useful to construct the numerical solution (U n j ) j∈N in the Propositions 2.3 and 3.1.

To summarize, we study all along the paper the following fully discrete approximation of the IBVP (1.1)-(1.3): 3), the UKC (1.4) is not enough and a more stringent restriction has to be imposed. Our aim is to prove that the SKC derived in [START_REF] Xin | Stiff Well-Posedness and Asymptotic Convergence for a Class of Linear Relaxation Systems in a Quarter Plane[END_REF] is then a sufficient condition for the stiff stability of the fully discrete IBVP (1.8), in other words the uniform stability with respect to the stiffness of the relaxation term.

(1.8)                      U n+1 j -U n j ∆t + 1 2∆x A U n+1 j+1 -U n+1 j-1 = 1 ε SU n+1 j , j ≥ 1, n ≥ 0, U 0 j = f j , j ≥ 0, BU n 0 = b n , n ≥ 0, 1 ∆t Γ U n+1 0 -U n 0 + 1 2∆x ΓA U n+1 1 - n+1 k=0 C n+1-k U k 0 = 1 ε ΓSU n+1 0 , n ≥ 0.
Theorem 1.1 (Main result). Assume that (B u , B v ) ∈ R 2 satisfies the SKC

B v = 0 or B u B v / ∈ - √ a, 0 . (1.9) Let λ xt ≤ 3
√ a/8 be a positive number. For any T > 0, there exists a constant C T > 0 such that for all ∆t > 0 and ∆x = λ xt ∆t, any (f j ) j∈N ∈ 2 (N, R 2 ) and (b n ) n∈N ∈ 2 (N, R), the solution (U n j ) j∈N to the scheme (1.8) satisfies N n=0 j≥0

∆x∆t|U n j | 2 + N n=0 ∆t|U n 0 | 2 ≤ C T   j≥0 ∆x|f j | 2 + N n=0 ∆t|b n | 2   , (1.10)
where N := T /∆t and C T is independent of ε ∈ (0, +∞).

In [START_REF] Xin | Stiff Well-Posedness and Asymptotic Convergence for a Class of Linear Relaxation Systems in a Quarter Plane[END_REF], Xin and Xu considered the IBVP for the Jin-Xin relaxation model [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF] and derived the SKC (1.9) to characterize its stiff well-posedness. They show in particular that the IBVP (1.1)-(1.3) is well-posed if and only if (1.9) holds. In the discrete IBVP (1.8), it seems that the SKC is also sufficient to derive uniform stability estimates. Besides, by linearity, the numerical scheme of the IBVP (1.8) can be broken up into two simpler problems, one with homogeneous initial condition (f j ) j∈N ≡ 0 and the other with homogeneous boundary b n ≡ 0, for any n ∈ N. The proof of Theorem 1.1 is based on two main ingredients, by assembling a result for the case of the following Cauchy problem

(1.11)        U I j n+1 -U I j n ∆t + 1 2∆x A U I j+1 n+1 -U I j-1 n+1 = 1 ε S U I j n+1 , j ∈ Z, n ≥ 0, U I j 0 = f j , j ∈ Z.
and another one for the problem (1.8) with homogeneous initial data. We state hereafter these two statements.

Proposition 1.2 (Cauchy problem). For any T > 0, there exists C T > 0 such that for all ∆t > 0, any

(f j ) j∈N ∈ 2 (N, R 2 ), the solution (U I j ) n j∈Z to (1.11) satisfies j∈Z ∆x|(U I j ) n | 2 ≤ C T j∈Z ∆x|f j | 2 , n ∈ N, (1.12)
where C T is independent of ε ∈ (0, +∞) and ∆x = λ xt ∆t. Proposition 1.3 (Homogeneous initial condition). Assume that the SKC (1.9) is satisfied. Then, there exists a constant C > 0 such that for any γ > 0 and any positive constant λ xt ≤ 3 √ a/8, the following property holds. For any ∆t > 0 together with ∆x = λ xt ∆t, and any boundary data

(b n ) n∈N ∈ 2 (N, R), the solution (U n j ) j∈N to (1.8) with (f j ) j∈N ≡ 0 satisfies γ γ∆t + 1 n≥0 j≥0 e -2γn∆t ∆t∆x|U n j | 2 + n≥0 e -2γn∆t ∆t|U n 0 | 2 ≤ C n≥0 e -2γn∆t ∆t|b n | 2 , (1.13)
where C is independent of ε ∈ (0, +∞).

To isolate the effects of a possible boundary layer and avoid the complicated interaction of boundary and initial layers, in Section 2, we consider the IBVP (1.8) with homogeneous initial data and nonzero boundary data b n , for any n ≥ 0. The numerical solution (U n j ) j∈N is constructed in Section 2.2 thanks to the Z-transform [START_REF] Jury | Theory and Application of the Z-Transform Method[END_REF][START_REF] Oppenheim | Signals and Systems[END_REF]. Furthermore, we follow the discrete transparent boundary condition at j = 0 as proposed in [START_REF] Arnold | Discrete transparent boundary conditions for the Schrödinger equation: fast calculation, approximation, and stability[END_REF][START_REF] Besse | Discrete transparent boundary conditions for the mixed KDV-BBM equation[END_REF][START_REF] Kazakova | Discrete Transparent Boundary Conditions for the Linearized Green-Naghdi System of Equations[END_REF] to find the explicit formula of the sequence (C m ) m≥0 . By using the Plancherel's theorem, under the SKC, the Proposition 1.3 is proved in Section 2.3. In order to illustrate the relevance of the SKC (1.9), we present in Section 2.4 some numerical results, for various values of the parameters (B u , B v ) and show that the numerical solution at the boundary x = 0 increases quickly if the SKC (1.9) does not hold. Besides, by the decrease of the error U (., t n ) -U n 2 2 (N,R 2 ) , we can observe the convergence of the discrete solution U n j to the exact one U (x j , t n ). After that, we observe the behavior of the energy terms U 2 (N×[0,T ),R 2 ) and U 2 ({0}×[0,T ),R 2 ) corresponding to whether or not the SKC (1.9) is valid. The nonzero initial data case is much more difficult with the sufficiency proof. This is due to the complicated interactions between the initial data, the boundary condition and the stiff relaxation term. Under the SKC, the numerical solution is again described by means of the Z-transform in Section 3.1. It is decomposed into three parts, by assembling the solution for the case of Cauchy problem (1.11), a numerical error term (U II j ) n j∈N and the solution for the case IBVP (1.8) with homogeneous initial data. Since the coefficients for computing the boundary value U n+1 -1 are defined for homogeneous initial data, this numerical error (U II j ) n j∈N is due to the interaction between the Cauchy problem and the IBVP with zero initial data. For the Cauchy problem, the Proposition 1.2 is studied in Section 3.2 by means of the discrete energy method. By an application of the Plancherel's theorem for Z-transform, the numerical error term (U II j ) n j∈N will be estimated in Section 3.3. In Section 3.4, we get the expected result of the Theorem 1.1 in the case IBVP with nonzero initial condition. In Section 3.5, we also look at the behavior of the numerical solution (U n j ) j∈Z and the energy terms U 2 (N×[0,T ),R 2 ) and U 2 ({0}×[0,T ),R 2 ) corresponding to whether or not the SKC (1.9) is valid. It seems that the SKC (1.9) is also necessary condition to guarantee the uniform stability of the IBVP (1.8) independent of the effect of the relaxation source term and the boundary dissipation.

Stiff stability of the IBVP with homogeneous initial condition

In this section, we consider the discrete IBVP (1.8) with nonzero boundary condition (b n ) n∈N ∈ 2 (N, R) and homogeneous Cauchy data (f j ) j∈N ≡ 0. Assuming that the SKC is satisfied, the numerical solution (U n j ) j∈N is obtained by using the Z-transform [START_REF] Jury | Theory and Application of the Z-Transform Method[END_REF][START_REF] Oppenheim | Signals and Systems[END_REF]. Thanks to the Plancherel's theorem, we then are able to get the expected result of the Proposition 1.3.

2.1. Notations and preliminary results. Before we enter the important proofs, let us introduce some notations and preliminary results. All along the paper, the complex values z and ξ are related through the formula

ξ = 1 -z -1 ∆t -1 , z = Re iθ , with R > 1, θ ∈ (-π, π].
Then, ξ obeys the inequalities

1 -R -1 ∆t -1 ≤ Re ξ ≤ 2∆t -1 . (2.1)
Besides, one also introduces the following matrix, already concerned with the continuous case [START_REF] Xin | Stiff Well-Posedness and Asymptotic Convergence for a Class of Linear Relaxation Systems in a Quarter Plane[END_REF]:

M (εξ) = A -1 (S -εξI) = 1 a 0 -(1 + εξ) -aεξ 0 ,
We recall that the eigenvalues and eigenvectors of M (εξ) can be easily found to be respectively

µ ± (εξ) = ± εξ(1 + εξ) a , r ± (εξ) =   1 aµ ∓ (εξ) 1 + εξ   .
In the above formula and all along the paper, the complex square root is defined with the branch cut along the negative real axis. Applying Lemma A.1 with the property Re ξ > 0 and ε > 0, we can prove

Re (µ -(εξ)) ≤ - εRe ξ √ a < 0, (2.2)
while, as a consequence,

Re (µ + (εξ)) ≥ εRe ξ √ a > 0.
Let us introduce We further define the following spectral projections

κ ± (εξ) = µ ± (εξ)λ xε + (µ + (εξ)λ xε ) 2 + 1, (2.3
Φ + (εξ) = 1 2g(εξ) 1 -g(εξ) g(εξ) -1 , Φ -(εξ) = 1 2g(εξ) 1 g(εξ) g(εξ) 1 , (2.4)
where we set

g(εξ) = aµ + (εξ) 1 + εξ . (2.5)
We also set Φ(εξ) the 2 × 2 matrix whose columns are composed by the components of the eigenvectors of the matrix M (εξ). We recall these matrices thus satisfy the following usefull identities

Φ + (εξ) = Φ(εξ) 0 0 0 1 Φ -1 (εξ) and Φ -(εξ) = Φ(εξ) 1 0 0 0 Φ -1 (εξ) (2.6) and Φ 2 + (εξ) = Φ + (εξ), Φ 2 -(εξ) = Φ -(εξ), Φ + (εξ)Φ -(εξ) = Φ -(εξ)Φ + (εξ) = 0. (2.7)
In order later on to construct and estimate the numerical solution (U n j ) j∈N by the Z-transform, the following lemmas are usefull: Lemma 2.1. [From [START_REF] Xin | Stiff Well-Posedness and Asymptotic Convergence for a Class of Linear Relaxation Systems in a Quarter Plane[END_REF]] Consider C + = {ζ ∈ C, Re ζ ≥ 0} the closed complex right half-plane. Under the SKC (1.9), the quantity g(ζ) is uniformly bounded in C + and the quantity B u + g(ζ)B v is uniformly bounded away from 0 in C + .

We omit the proof that the reader can find in the work by Xin and Xu [START_REF] Xin | Stiff Well-Posedness and Asymptotic Convergence for a Class of Linear Relaxation Systems in a Quarter Plane[END_REF].

Lemma 2.2. Let us consider the 4 × 4 matrix

M 1 (εξ) = 2λ xε M (εξ) I I 0 . (2.8)
Then, the k-th power of M 1 (εξ) reads also

M k 1 (εξ) = - 1 κ + (εξ) + κ -(εξ) κ k+1 (εξ) Ψ k (εξ) κ k (εξ) Ψ k+1 (εξ) κ k (εξ) Ψ k+1 (εξ) κ k-1 (εξ) Ψ k (εξ) , (2.9)
where

κ k (εξ) = (-1) k κ k + (εξ) -κ k -(εξ), Ψ k (εξ) = Φ -(εξ) + (-1) k Φ + (εξ).
(2.10) Proof. In this algebraic proof, we skip the dependence on εξ. Since the columns of the matrix Φ are composed by the components of the eigenvectors of the matrix M , the considered matrix M k 1 can be reformulated simply as

M k 1 = ΦM k 2 Φ -1 , (2.11)
where

Φ = Φ 0 0 Φ , M 2 = D 1 I I 0 , D 1 = 2λ xε diag (µ -, µ + ) .
Let Ψ is the 4×4 matrix whose columns are composed by the components of the eigenvectors of the matrix

M 2 Ψ =     -κ + κ - 0 0 0 0 -κ -κ + 1 1 0 0 0 0 1 1     so that M k 2 = ΨD k 2 Ψ -1 with D 2 = diag (-κ + , κ -, -κ -, κ + ). Therefore, the formula M k 1 in (2.11) reads M k 1 = ΦΨD k 2 Ψ -1 Φ -1 . (2.12)
By using the properties Φ ± in (2.6) and (2.12), one obtains

M k 1 = - 1 κ + + κ - κ k+1 Ψ k κ k Ψ k+1 κ k Ψ k+1 κ k-1 Ψ k
with κ k and Ψ k are the same as in (2.10).

2.2.

Solution by Z-transform. Firstly, we apply the Z-transform with respect to time index n ∈ N, which is discrete analogue of the Laplace transform in time t ∈ R + . This method enables the representation and estimations of the numerical solution (U n j ) j∈N . The definition reads as follows (see [START_REF] Jury | Theory and Application of the Z-Transform Method[END_REF][START_REF] Oppenheim | Signals and Systems[END_REF] for more details)

U j (z) = Z{U n j }(z) = n≥0 U n j z -n , |z| > 1.
Since we assume (U 0 j ) j∈N ≡ 0, observe that the Z-transform of the time-shifted numerical solution reads

n≥0 U n+1 j z -n = z U j (z) -zU 0 j = z U j (z).
Therefore, the IBVP (1.8) with zero initial data becomes

       U j+1 (z) -U j-1 (z) = 2λ xε M (εξ) U j (z), j ≥ 1, (2.13a) B U 0 (z) = b(z), (2.13b) ΓA U 1 (z) -z -1 Υ( U 0 (z)) -2λ xε M (εξ) U 0 (z) = 0, (2.13c)
where Υ( U 0 (z)) is the Z-transform of the sequence

(z) = Z{b n }(z) = n≥0 b n z -n .
Secondly, we look at the solution ( U j ) j∈N (z) to (2.13a)-(2.13c). This is the object of the following proposition:

Proposition 2.3. Assume that the SKC (1.9) is satisfied. Assume that Γ and Υ in the boundary condition (2.13c) are defined by

Γ = -aB v B u , Υ( U 0 (z)) = κ + (εξ)z U 0 (z). (2.14) Then the solution ( U j ) j∈N (z) ∈ 2 (N, C 2 ) to (2.13a)-(2.13c) takes the form U j (z) = b(z) B u + g(εξ)B v κ j -(εξ)r -(εξ). (2.15)
Proof. Before we prove the above result, let us notice that we omit the explicit dependence in εξ. Firstly, we look at the solution ( U j ) j∈N (z) to (2.13a) and consider the two-dimensional problem (2.13a) under the following one-step recurrence form

W j+1 (z) = M 1 W j (z), (2.16)
where M 1 is given by (2.8) and

W j (z) = U j (z) U j-1 (z) . (2.17)
The solution (W j ) j∈N (z) to (2.16) is simply W j (z) = M j 1 W 0 (z). Together with the the explicit formula of M j 1 in Lemma 2.2, the solution ( U j ) j∈N (z) to (2.13a) is therefore given by

U j (z) = - 1 κ + + κ - × κ j+1 Ψ j U 0 (z) + κ j Ψ j+1 U -1 (z) .
By using the definition of κ k and Ψ k in (2.10), the above formula is now equivalent to

U j (z) = - (-1) j κ j + κ + + κ - × Φ --κ + U 0 (z) + U -1 (z) + (-1) j+1 Φ + κ + U 0 (z) + U -1 (z) + κ j - κ + + κ - × κ -Φ -+ (-1) j Φ + U 0 (z) + Φ -+ (-1) j+1 Φ + U -1 (z) .
(2.18)

Since we expect ( U j ) j∈N (z) ∈ 2 (N, C 2 ), we need a natural boundary condition at x = +∞. Besides, one gets |κ + | > 1 and |κ -| < 1. Thus, the natural boundary condition takes the form

(2.19)    Φ --κ + U 0 (z) + U -1 (z) = 0, Φ + κ + U 0 (z) + U -1 (z) = 0.
By the definition of Φ ± in (2.4), the system (2. [START_REF] Oppenheim | Signals and Systems[END_REF]) is equivalent to

   (g, 1) -κ + U 0 (z) + U -1 (z) = 0, (g, -1) κ + U 0 (z) + U -1 (z) = 0.
Then, we have

U -1 (z) = κ + g × 0 1 g 2 0 U 0 (z).
Furthermore, we can see that

Φ --Φ + = 1 g 0 1 g 2 0 .
Thus,

U -1 (z) = κ + (Φ --Φ + ) U 0 (z). (2.20)
Plugging (2.20) into (2.18), we have

U j (z) = κ j - κ + + κ - κ -Φ -+ (-1) j Φ + + κ + Φ -+ (-1) j+1 Φ + Φ --Φ + U 0 (z).
Under the properties of Φ ± in (2.7), the above formula becomes

U j (z) = κ j -Φ -+ (-1) j Φ + U 0 (z). (2.21)
Secondly, we look at the boundary condition (2.13b) and (2.13c). Under the choice Υ( U 0 (z)) in (2.14), the boundary condition (2.13c) becomes

ΓA U 1 (z) -(κ + I + 2λ xε M ) U 0 (z) = 0. (2.22)
Indeed, we can compute separately

κ + I + 2λ xε M = κ -Φ -+ (κ + + 2λ xε µ + ) Φ + , U 1 (z) = κ -(Φ --Φ + ) U 0 (z). (2.
ΓAΦ + U 0 (z) = 0. (2.24)
Under the choice Γ in (2.14), we have

ΓAΦ + = a 2g × (B u + gB v ) × (g, -1).
Thus,

(B u + gB v ) × (g, -1) U 0 (z) = 0. (2.25)
From the Lemma 2.1, the equation (2.25) is equivalent under the SKC to

(g, -1) U 0 (z) = 0. (2.26)
Together with the boundary condition (2.13b), the value of U 0 (z) has to satisfy

B u B v g -1 U 0 (z) = b(z) 1 0 .
Then, again under the SKC, we have

U 0 (z) = b(z) B u + g(εξ)B v r -(εξ). (2.27)
Plugging the value of U 0 (z) in (2.27) into (2.21), the solution ( U j ) j∈N (z) to (2.13a)-(2.13c) is given by

U j (z) = b(z) B u + gB v κ j -Φ -+ (-1) j Φ + r -(εξ).
Since Φ + r -= 0 and Φ -r -= r -, the solution ( U j ) j∈N (z) finally is

U j (z) = b(z) B u + gB v κ j -r -.
This ends the proof of Proposition 2.3.

With ( U j ) j∈N (z) found in (2.15), the numerical solution (U n j ) j∈N to the IBVP (1.8) with nonzero boundary condition (b n ) n∈N ∈ 2 (N, R) and homogeneous Cauchy data (f j ) j∈N ≡ 0 can be obtained by inverting the Z-transform [START_REF] Jury | Theory and Application of the Z-Transform Method[END_REF][START_REF] Oppenheim | Signals and Systems[END_REF]

U n j = 1 2π π -π U j Re iθ R n e inθ dθ, R > 1.
Let us remark that an important assumption of Proposition 2.3 is Υ( U 0 (z)) = κ + (εξ)z U 0 (z). We now follow the discrete transparent boundary as proposed in [START_REF] Arnold | Discrete transparent boundary conditions for the Schrödinger equation: fast calculation, approximation, and stability[END_REF][START_REF] Besse | Discrete transparent boundary conditions for the mixed KDV-BBM equation[END_REF][START_REF] Kazakova | Discrete Transparent Boundary Conditions for the Linearized Green-Naghdi System of Equations[END_REF] to find the explicit formula for the sequence (C m ) m≥0 .

Definition 2.4. Let ε > 0, R > 1, θ ∈ (-π, π] and then κ + (εξ) be given by (2.3). The value of (C m ) m≥0 is defined as follows:

C m = 1 π π 0 Re κ + (εξ)R m e imθ dθ. (2.28)
Let us mention at this step that the values C m above are designed in the case of homogeneous initial data and are kept unchanged in the case of nonzero initial data in forthcoming Section 3. Thanks to the convolution property and inverting Z-transform, we now show that the definition of (C m ) m≥0 in (2.28) is the suitable choice to get the required identity Υ( U 0 (z)) = κ + (εξ)z U 0 (z) in Proposition 2.3. This is the object of the next lemma: Lemma 2.5. Let (C m ) m≥0 be defined from Definition 2.4, then

Υ( U 0 (z)) = Z n+1 k=0 C n+1-k U k 0 (z) = κ + (εξ)z U 0 (z). Proof. Since µ + (εξ) = µ + (εξ), one obtains κ + (εξ) = κ + (εξ). Then, Re κ + (εξ)R m e imθ = 1 2 κ + (εξ)R m e imθ + κ + (εξ)R m e -imθ = 1 2 κ + (εξ)R m e imθ + κ + (εξ)R m e -imθ .
Thus, the value of (C m ) m≥0 in (2.28) can be reformulated as

C m = 1 2π π 0 κ + (εξ)R m e imθ dθ + π 0 κ + (εξ)R m e -imθ dθ = 1 2π π -π κ + (εξ)R m e imθ dθ = Z -1 (κ + (εξ)) (m).
By the convolution property and inverting Z-transform, we can conclude that

Υ( U 0 (z)) = Z n+1 k=0 C n+1-k U k 0 (z) = κ + (εξ)z U 0 (z),
This ends the proof of Lemma 2.5. π, π] and then κ -(εξ) be given by (2.3). Then the following property holds

Let ε > 0, R > 1, θ ∈ (-
j≥0 |κ -(εξ)| 2j ≤ ∆t √ a ∆x(1 -R -1 ) . (2.30) Proof. Since the property of Re (µ -(εξ)) in (2.2), we can prove Re (µ -(εξ)) λ xε + (Re (µ -(εξ)) λ xε ) 2 + 1 2 ≤ η∆x + η 2 ∆x 2 + 1 2 , (2.31)
where η = -a -1/2 Re ξ. According to Lemma A.2 and the inequality (2.31), we have

|κ -(εξ)| 2 = µ -(εξ)λ xε + (µ -(εξ)λ xε ) 2 + 1 2 ≤ η∆x + η 2 ∆x 2 + 1 2 .
Then, we obtain the following estimate

j≥0 |κ -(εξ)| 2j = 1 -|κ -(εξ)| 2 -1 ≤ 1 -η∆x + η 2 ∆x 2 + 1 2 -1 .
Since Re ξ satisfies the property (2.1), we get

∆t √ a 2 ≤ - 1 η ≤ ∆t √ a 1 -R -1 .
If we assume now ∆x ≤ 3 √ a 8 ∆t ≤ -3 4η then we have

1 -η∆x + η 2 ∆x 2 + 1 2 -1 ≤ -η -1 ∆x -1 .
Thus, we conclude that

j≥0 |κ -(εξ)| 2j ≤ -η -1 ∆x -1 ≤ ∆t √ a ∆x(1 -R -1 )
.

This ends the proof of Lemma 2.6.

Secondly, by an application of the following Plancherel's theorem for

Z-transform n≥0 R -2n |U n j | 2 = 1 2π π -π | U j (Re iθ )| 2 dθ, R > 1, we have n≥0 R -2n |U n 0 | 2 = 1 2π π -π | U 0 (Re iθ )| 2 dθ = 1 2π π -π b(z) B u + g(εξ)B v 2 1 + |g(εξ)| 2 dθ.
From the Lemma 2.1, under the SKC, we then obtain n≥0

R -2n |U n 0 | 2 1 2π π -π | b(Re iθ )| 2 dθ n≥0 R -2n |b n | 2 . (2.32)
Similarly, by an application of the Plancherel's theorem for Z-transform, we have

n≥0 j≥0 R -2n |U n j | 2 = 1 2π j≥0 π -π | U j (Re iθ )| 2 dθ = 1 2π j≥0 π -π b(z) B u + g(εξ)B v 2 |κ -(εξ)| 2j 1 + |g(εξ)| 2 dθ.
Again, under the SKC, we get from Lemma 2.1

n≥0 j≥0 R -2n |U n j | 2 1 2π j≥0 π -π | b(z)| 2 |κ -(εξ)| 2j dθ. (2.33)
Following Lemma 2.6, if we assume (2.29) holds, then the inequality (2.33) becomes

R -1 R ∆x n≥0 j≥0 R -2n |U n j | 2 ∆t 2π π -π | b(Re iθ )| 2 dθ ∆t n≥0 R -2n |b n | 2 . (2.34)
According to (2.32) and (2.34), there exists a constant C > 0 such that

R -1 R n≥0 j≥0 R -2n ∆x|U n j | 2 + n≥0 R -2n ∆t|U n 0 | 2 ≤ C n≥0 R -2n ∆t|b n | 2 .
By setting in the above formula R = e γ∆t for γ > 0 and ∆t > 0, and using the classical lower bound e γ∆t ≥ 1 + γ∆t, we obtain that there exists a constant c > 0 such that γ γ∆t + 1

n≥0 j≥0 e -2γn∆t ∆t∆x|U n j | 2 + n≥0 e -2γn∆t ∆t|U n 0 | 2 ≤ c n≥0 e -2γn∆t ∆t|b n | 2 .
This ends the proof of the Proposition 1.3.

Let us observe that the scheme (1.8) together also with its boundary condition is closed to be forwardin-time, except it is one-step implicit. By this property, changing the data b to zero after some time T and unchanged before that time T , the discrete solution U n j is the same for n∆t < T . Therefore, there exists a constant C T > 0 such that

N n=0 j≥0 ∆x∆t U n j 2 + N n=0 ∆t |U n 0 | 2 ≤ C T N n=0 ∆t |b n | 2 , (2.35)
with N := T /∆t. This will be useful to prove the Theorem 1.1.

Numerical experiments.

In this paragraph, we first provide the behavior of the numerical solution (U n j ) j∈N according to whether or not the SKC (1.9) is valid. We also look at the degenerate case when the UKC (1.4) does not hold (and thus, none of the other stability conditions). Following the continuous case studied by Xin and Xu in [START_REF] Xin | Stiff Well-Posedness and Asymptotic Convergence for a Class of Linear Relaxation Systems in a Quarter Plane[END_REF], the solution of the IBVP (1.1)-(1.3) with homogeneous initial condition can be constructed by the method of Laplace transform. By inverting the Laplace transform, the solution U (x, t) has form

U (x, t) = L -1 U = 1 2π +∞ -∞ e ζt b(ζ) B u + g(εζ)B v e µ -(εζ)x/ε r -(εζ)dβ,
where ζ = α + iβ, α > 0. Then, we observe the error between the exact solution U (x j , t n ) and the numerical solution U n j of the numerical scheme (1.8) with homogeneous initial data at the grid point (x j , t n ) = (j∆x, n∆t). After that, we present some numerical experiments and observe the effective behavior of the energy terms U 2 (N×[0,T ),R 2 ) and U 2 ({0}×[0,T ),R 2 ) corresponding to whether or not the SKC (1.9) is valid.

As main parameters for the experiments, we choose a = 1, B v = 1, λ xt = 1/3 and let the relaxation rate ε and the boundary data B u vary. The test case we consider concerns the following data. The initial data is the homogeneous one (f j ) j∈N ≡ 0. The boundary data is

b(t) = t 2 sin(t).
Let us observe that these data are compatible in the corner (x, t) = (0, 0) in the sense that Bf (0) = b(0). Moreover, the Laplace transform of b(t) is

b(ζ) = ζ (ζ 2 + 1) 2 .
2.4.1. The behavior of the numerical solution. Let the space step ∆x = 10 -2 and the time step ∆t = λ -1 xt ∆x. Firstly, we choose the value of B u such that the SKC (1.9) is satisfied with ε = 10 -2 and also with ε = 10 In the first case, ε = 10 -2 , the incoming solution at the boundary x = 0 go slowly. This is due to the initial relaxation of solution to the equilibrium system. In the case ε = 10 2 , its solution seems to be faster. It is not so much influenced by relaxation source term but more by the boundary dissipation.

Secondly, we choose the value of B u such that the SKC (1.9) is not satisfied. Besides, we also present the numerical solution when the Uniform Kreiss Condition (1.4) is wrong. The Figures 2.3 and 2.4 show the numerical solution (U n j ) j∈N over the time interval t ∈ [0, 0.5). When the SKC (1.9) fails, we observe that the numerical solution at the boundary rise gradually. This is the case for example for ε = 10 -2 together with the parameters (B u , B v ) = (-1/2, 1). The behavior is even worse when the UKC (1.4) is not satisfied (see Figure 2.4).

2.4.2.

The error between the exact solution and the numerical solution. Let us begin with the notation

E(t n ) :=   ∆x j≥0 |U (x j , t n ) -U n j | 2   1/2 . (2.36)
We choose a set of values B u such that the SKC (1.9) is satisfied with the space step ∆x and the relation rate ε vary. The error, as measured in (2.36), are reported in the Tables 1 and2 According to the experiments in Tables 1 and2, for some ε ∈ (0, +∞) and (B u , B v ) satisfying the SKC (1.9), the observed convergence rate is 1 since going down ∆x by a factor 2 decreases the error of the same factor 2. It means that the behavior of the numerical solution U n j is the same as the evolution of the exact solution U (x j , t n ). This is the case for example for ε = 10 -2 together with the parameters (B u , B v ) = (-4, 1).

∆x ε = 10 -2 ε = 10 -1 ε = 1 ε = 10 ε = 10 2 5 × 10 -2 6.8 × 10 -3 1.2 × 10 -2 2.6 × 10 -2 3.1 × 10 -2 3.2 × 10 -2 25 × 10 -3 3 × 10 -3 5.9 × 10 -3 1.3 × 10 -2 1.6 × 10 -2 1.6 × 10 -2 125 × 10 -4 1.5 × 10 -3 2.9 × 10 -3 6.8 × 10 -3 8.2 × 10 -3 8.3 × 10 -3 625 × 10 -5 7.2 × 10 -4 1.5 × 10 -3 3.4 × 10 -3 4.1 × 10 -3 4.2 × 10 -3
Table 1. The error E(1.2) for B u = -4. Table 2. The error E(1.2) for B u = 3.

second to the 2 in time energy of the numerical trace at the boundary:

E 1 := U 2 2 (N×[0,T ),R 2 ) = N n=0 j≥0 ∆x∆t|U n j | 2 , E 2 := U 2 2 ({0}×[0,T ),R 2 ) = N n=0 ∆t|U n 0 | 2 ,
(2.37) which are shown in the Table 3 and Figures 2.5, 2.6.

B u ε = 10 -2 ε = 1 ε = 10 2 -4 35 × 10 -5 77 × 10 -5
47 × 10 -4 -2 15 × 10 -4 39 × 10 -4 43 × 10 -3 -1 2.6 × 10 17 1.69 × 10 31 3.17 × 10 56 -0.5 44047.9 418525.12 2837033.2 1 34 × 10 -4 5 × 10 -3 10 -2 3 48 × 10 -5 87 × 10 -5 2 × 10 -2

B u ε = 10 -2 ε = 1 ε = 10 2 -4 66 × 10 -4 8 × 10 -3 2 × 10 -2 -2 29 × 10 -3 4 × 10 -2

Energy

The SKC is no sa isfied wi h (B u , B v ) = (-1/2, 1) and a = 1 ε = 10 -2 ε = 1 ε = 10 2 • For some ε ∈ (0, +∞), the values of E 1 and E 2 rise gradually when the SKC (1.9) is not satisfied. This is the case for example for ε = 10 2 together with the parameters (B u , B v ) = (-1/2, 1). The behavior is even worse when the UKC (1.4) is not hold. • In the case ε = 10 -2 , the energy term E 1 and E 2 increase slowly. This is due to the effect of incoming solution at the boundary when the initial relaxation of solution tends to the equilibrium system. In the case ε = 10 2 , those values increase fairly rapidly. It is not so much influenced by relaxation source term but more by the boundary dissipation.

Clearly, the numerical results show that the numerical solution at the boundary x = 0 increase quickly as soon as the SKC (1.9) does not hold. If the UKC (1.4) is not satisfied, the behavior of numerical solution is even worse. Besides, by the decrease of the error E(t n ), we can see that the values of U n j tend to the exact solution U (x j , t n ). Indeed, it seems that the SKC (1.9) is also necessary to ensure the non-increase rapidly of the energy terms E 1 and E 2 under the effect of the relaxation source term and the boundary dissipation.

Stiff stability of the IBVP with homogeneous boundary condition

For convenience in the forthcoming discussions, we recall that the IBVP (1.8) with homogeneous boundary condition writes

(3.1)                      U n+1 j -U n j ∆t + 1 2∆x A U n+1 j+1 -U n+1 j-1 = 1 ε SU n+1 j , j ≥ 1, n ≥ 0, U 0 j = f j , j ≥ 0, BU n 0 = 0, n ≥ 0, 1 ∆t Γ U n+1 0 -U n 0 + 1 2∆x ΓA U n+1 1 - n+1 k=0 C n+1-k U k 0 = 1 ε ΓSU n+1 0 , n ≥ 0.
In [23, Section 5], under the SKC, Xin and Xu find explicitly the solution U (x, t) of the IBVP (1.6) by the method of Laplace transform. The solution is decomposed into two ingredients, by assembling a solution for the case of the Cauchy problem and another one for the case of the IBVP with homogeneous initial condition. In our case, assuming the SKC to hold, the numerical solution (U n j ) j∈N can be constructed by the method of Z-transform. Since the coefficients (C m ) m≥0 are defined for homogeneous initial data, the numerical solution (U n j ) j∈N of the IBVP (3.1) consists of not only the solutions for case of the Cauchy problem and for the IBVP with zero initial data but also another numerical error term (U II j ) n j∈N . To complete the proof of the Theorem 1.1 with homogeneous boundary condition, we first use the means of discrete energy method in order to prove the Proposition 1.2. By an application of the Plancherel's theorem for Z-transform [START_REF] Jury | Theory and Application of the Z-Transform Method[END_REF][START_REF] Oppenheim | Signals and Systems[END_REF], the numerical error term of (U II j ) n j∈N will be estimated in Section 3.3. After that, we get the expected result of the case IBVP with homogeneous initial condition.

3.1. Solution by Z-transform. Again, we follow the explicit solving of the IBVP (3.1) by using the Z-transform. With

U j (z) = Z{U n j }(z) = n≥0 U n j z -n , |z| > 1.
Importantly, we now have (f j ) j∈N = 0, and thus we get

n≥0 U n+1 j z -n = z U j (z) -zU 0 j = z U j (z) -zf j .
Therefore, (3.1) becomes

         U j+1 (z) -U j-1 (z) = 2λ xε M (εξ) U j (z) + f j+1 -2λ xε M (ε∆t -1 )f j -f j-1 , j ≥ 1, (3.2a) B U 0 (z) = 0, (3.2b) ΓA U 1 (z) -κ + (εξ)I + 2λ xε M (εξ) U 0 (z) -f 1 = 0. (3.2c) Let us recall that the Z-transform of n+1 k=0 C n+1-k U k 0 is given by κ + (εξ)z U 0 (z).
Firstly, we look at the solution ( U j ) j∈N (z) to (3.2a)-(3.2c). This is the object of the following proposition: Proposition 3.1. Assume that the SKC (1.9) is satisfied. Let (f j ) j∈N ∈ 2 (N, R 2 ) and denote V I k , w I (εξ) and w II (εξ) as follows:

V I k = f k+1 -2λ xε M (ε∆t -1 )f k -f k-1 , w I (εξ) = +∞ k=0 (-1) -k κ -k + (εξ)V I k , w II (εξ) = - +∞ k=0 κ -k + (εξ)V I k . (3.3)
Then, the solution

( U j ) j∈N (z) ∈ 2 (N, C 2 ) to (3.2a)-(3.2c) takes the form U j (z) = κ j+1 -(εξ) 4g(εξ) × B u -g(εξ)B v B u + g(εξ)B v × (g(εξ), 1) w I (εξ) -(g(εξ), -1) w II (εξ) × B u -g(εξ)B v B u + g(εξ)B v × r -(εξ) + (-1) j+1 r + (εξ) - κ j -(εξ) κ + (εξ) + κ -(εξ) × Φ -(εξ)w I (εξ) + (-1) j Φ + (εξ)w II (εξ) + 1 κ + (εξ) + κ -(εξ) × j-1 k=0 κ j-k -(εξ) Φ -(εξ) + (-1) j-k-1 Φ + (εξ) V I k + +∞ k=j (-1) j-k κ j-k + (εξ) Φ -(εξ) + (-1) j-k-1 Φ + (εξ) V I k . (3.4)
Proof. Before we prove the above result, let us notice that we omit the explicit dependence in εξ. Firstly, we look at the solution ( U j ) j∈N (z) to (3.2a) and consider the two-dimensional problem (3.2a) under the following nonhomogeneous one-step recurrence form

W j+1 (z) = M 1 W j (z) + V j , (3.5)
where W j is the same as in (2.17) and

V j = V I j 0 .
The solution (W j ) j∈N (z) to (3.5) is thus given by

W j (z) = M j 1 W 0 (z) + j-1 k=0 M j-1-k 1 V k . (3.6)
Together with the the explicit formula of M j 1 in Lemma 2.2, the solution ( U j ) j∈N (z) to (3.2a) is given by

U j (z) = - 1 κ + + κ - × κ j+1 Ψ j U 0 (z) + κ j Ψ j+1 U -1 + j-1 k=0 κ j-k Ψ j-k-1 V I k .
By the definition of κ k and Ψ k in (2.10), the above formula is equivalent to

U j (z) = - (-1) j κ j + κ + + κ - × -κ + Φ -+ (-1) j Φ + U 0 (z) + Φ -+ (-1) j+1 Φ + U -1 (z) + j-1 k=0 (-κ + ) -k Φ -+ (-1) j-k-1 Φ + V I k + κ j - κ + + κ - × κ -Φ -+ (-1) j Φ + U 0 (z) + Φ -+ (-1) j+1 Φ + U -1 (z) + j-1 k=0 κ -k - Φ -+ (-1) j-k-1 Φ + V I k . (3.7)
Thanks to the definition of w I and w II in (3.3), one has

j-1 k=0 (-κ + ) -k Φ -+ (-1) j-k-1 Φ + V I k = Φ -w I + (-1) j Φ + w II - +∞ k=j (-κ + ) -k Φ -+ (-1) j-k-1 Φ + V I k . (3.8)
Substituting (3.8) into (3.7), we have

U j (z) = - (-1) j κ j + κ + + κ - × Φ --κ + U 0 (z) + U -1 (z) + w I + (-1) j+1 Φ + κ + U 0 (z) + U -1 (z) -w II + 1 κ + + κ - × +∞ k=j (-1) j-k κ j-k + Φ -+ (-1) j-k-1 Φ + V I k + κ j - κ + + κ - × κ -Φ -+ (-1) j Φ + U 0 (z) + Φ -+ (-1) j+1 Φ + U -1 (z) + j-1 k=0 κ -k - Φ -+ (-1) j-k-1 Φ + V I k .
(3.9)

Since we expect By the definition of Φ ± in (2.4), the system (3.10) is equivalent to

( U j ) j∈N (z) ∈ 2 (N, C 2 ),
   g, 1 -κ + U 0 (z) + U -1 (z) + w I = 0, g, -1 κ + U 0 (z) + U -1 (z) -w II = 0.
Then, we have

U -1 (z) = κ + (Φ --Φ + ) U 0 (z) -Φ -w I + Φ + w II . (3.11)
Plugging (3.11) into (3.9), we get

U j (z) = 1 κ + + κ - × +∞ k=j (-1) j-k κ j-k + Φ -+ (-1) j-k-1 Φ + V I k + κ j - κ + + κ - × κ -Φ -+ (-1) j Φ + + κ + Φ -+ (-1) j+1 Φ + Φ --Φ + U 0 (z) + κ j - κ + + κ - × Φ -+ (-1) j+1 Φ + -Φ -w I + Φ + w II + 1 κ + + κ - × j-1 k=0 κ j-k - Φ -+ (-1) j-k-1 Φ + V I k .
Under the properties of Φ ± in (2.7), the above formula becomes

U j (z) = κ j -Φ -+ (-1) j Φ + U 0 (z) - κ j - κ + + κ - × Φ -w I + (-1) j Φ + w II + 1 κ + + κ - × j-1 k=0 κ j-k - Φ -+ (-1) j-k-1 Φ + V I k + +∞ k=j (-1) j-k κ j-k + Φ -+ (-1) j-k-1 Φ + V I k .
(3.12)

Secondly, we look at the boundary data U 0 (z) and extend the initial data (f j ) j∈N to the whole line by setting f j = 0 for j ≤ 0. Since f 0 = f -1 = 0 and Φ + + Φ -= I, we can see that

U 1 (z) = κ -Φ --Φ + U 0 (z) - κ - κ + + κ - × Φ -w I -Φ + w II + 1 κ + + κ - × κ -f 1 -κ + +∞ k=1 (-1) -k κ -k + Φ -+ (-1) -k Φ + V I k . (3.13)
On the other hand, by the definition of w I and w II in (3.3), we get the following property

+∞ k=1 (-1) -k κ -k + Φ -+ (-1) -k Φ + V I k = Φ -w I -Φ + w II -f 1 . (3.14)
Substituting (3.14) into (3.13), the value of U 1 (z) can be reformulated as

U 1 (z) = κ -(Φ --Φ + ) U 0 (z) -Φ -w I + Φ + w II + f 1 .
Thus, the equation (3.2c) becomes

ΓA κ -(Φ --Φ + ) U 0 (z) -Φ -w I + Φ + w II -(κ + I + 2λ xε M ) U 0 (z) = 0. (3.15)
Indeed, we observe that

κ + I + 2λ xε M = κ -Φ -+ (κ + + 2λ xε µ + )Φ + .
Then, the equation (3.15) can be represented as

ΓAΦ + U 0 (z) = - 1 2κ + ΓA Φ -w I -Φ + w II . (3.16)
Besides, one has

ΓAΦ + = a 2g × (B u + gB v ) × (g, -1).
From the Lemma 2.1, the equation (3.16) is equivalent under the SKC to

(g, -1) U 0 (z) = - g aκ + (B u + gB v ) ΓA Φ -w I -Φ + w II .
Together with the boundary condition (3.2b), one gets

B u B v g -1 U 0 (z) = - g aκ + (B u + gB v )   0 ΓA Φ -w I -Φ + w II   .
Then, under the SKC, we have

U 0 (z) = - g aκ + (B u + gB v ) 2 × ΓA Φ -w I -Φ + w II B v -B u . (3.17)
Let us extend the initial data (f j ) j∈Z to the whole line by setting f j = 0 for j ≤ 0. It is easy to verify that U I j (z) corresponds to the Z-transform of the solution (U I j ) n of the following extended Cauchy problem (3.22).

(3.22)        U I j n+1 -U I j n ∆t + 1 2∆x A U I j+1 n+1 -U I j-1 n+1 = 1 ε S U I j n+1 , j ∈ Z, n ≥ 0, U I j 0 = f j , j ∈ Z.
With ( U II j ) j∈N (z) found in (3.20), the value of (U II j ) n j∈N can be obtained by inverting the Z-transform

U II j n = 1 2π π -π U II j (Re iθ )R n e inθ dθ, R > 1. (3.23)
Indeed, the value of U III j (z) can be reformulated as

U III j (z) = -B U I 0 (z) + U II 0 (z) B u + g(εξ)B v × κ -(εξ)r -(εξ).
Following Section 2.2, U III j (z) corresponds to the Z-transform of the solution (U III j ) n of the IBVP with the homogeneous initial data (3.24)

                             U III j n+1 -U III j n ∆t + 1 2∆x A U III j+1 n+1 -U III j-1 n+1 = 1 ε S U III j n+1 , j ≥ 1, n ≥ 0, U III j 0 = 0, j ≥ 0, B U III 0 n = -B U I 0 n + U II 0 n , n ≥ 0, 1 ∆t Γ U III 0 n+1 -U III 0 n + 1 2∆x ΓA U III 1 n+1 - n+1 k=0 C n+1-k U III 0 k = 1 ε ΓS U III 0 n+1
, n ≥ 0.

3.2.

The energy method for the Cauchy problem. In this paragraph, we prove the Proposition 1.2 by means of the discrete energy method. The energy estimate in the continuous case are obtained using the integration by parts rule. Therefore, we need the corresponding summation by parts rules for the discrete approximations of ∂/∂ x [START_REF] Gustafsson | Time dependent problems and difference methods[END_REF]. The idea is to find a symmetric positive definite matrix H, such that HA is symmetric and HS is negative definite. Therefore, we choose

H = a 0 0 1 . Now,
(U I j ) n+1 -(U I j ) n , H(U I j ) n+1 + ∆t 2∆x j∈Z A (U I j+1 ) n+1 -(U I j-1 ) n+1 , H(U I j ) n+1 = ∆t ε j∈Z S(U I j ) n+1 , H(U I j ) n+1 , (3.25)
where ., . denotes the usual Euclidean inner product. Since H is a symmetric positive definite matrix, we have j∈Z

(U I j ) n+1 -(U I j ) n , H(U I j ) n+1 ≥ 1 2 j∈Z (U I j ) n+1 , H(U I j ) n+1 -(U I j ) n , H(U I j ) n .
Together with the symmetric matrix HA, the second flux term in (3.25) becomes j∈Z A (U I j+1 ) n+1 -(U I j-1 ) n+1 , H(U I j ) n+1 = 0.

Thus, we directly get the inequality j∈Z

(U I j ) n+1 , H(U I j ) n+1 -(U I j ) n , H(U I j ) n ≤ 2∆t ε j∈Z S(U I j ) n+1 , H(U I j ) n+1 . (3.26)
Let us remind that HS is negative definite. Then, from the inequality (3.26), for any n ∈ N, the following inequality holds j∈Z

(U I j ) n , H(U I j ) n ≤ j∈Z f j , Hf j . (3.27)
Furthermore, since H is a symmetric positive definitive matrix, the following inequality holds for some constants m, k > 0 To complete the proof of the Theorem 1.1 for the numerical scheme of the IBVP (3.1), observe that from (3.29) and setting f j = 0 for j < 0, for any T > 0, there exists C T > 0 such that

m (U I j ) n , H(U I j ) n ≤ (U I j ) n , (U I j ) n ≤ k (U I j ) n , H(U I j ) n . ( 3
(f j ) j∈N ∈ 2 (N, R 2 ), the values of (U II j ) j∈N defined in (3.23) satisfy N n=0 j≥0 ∆t∆x (U II j ) n 2 + N n=0 ∆t (U II 0 ) n 2 ≤ C T j≥0 ∆x|f j | 2 (3.32)
where N := T /∆t and C T is independent of ε ∈ (0, +∞).

Proof. By an application of the following Plancherel's theorem for Z-transform, we have

n≥0 R -2n (U II 0 ) n 2 = 1 2π π -π | U II 0 (Re iθ )| 2 dθ, R > 1 = 1 2π π -π |κ -(εξ)| 2 16|g(εξ)| 2 × B u -g(εξ)B v B u + g(εξ)B v 2 × (g(εξ), 1)w I (εξ) 2 × 1 + |g(εξ)| 2 dθ.
From the Lemma 2.

|κ -(εξ)| 2j ≤ ∆t √ a ∆x(1 -R -1 )
.

Together with λ xt = ∆x/∆t, the inequality (3.34) becomes

R -1 R n≥0 j≥0 R -2n ∆x (U II j ) n 2 k≥0 ∆x|f k | 2 . (3.35)
Assembling the estimates (3.33) and (3.35), there exists C > 0 such that

R -1 R n≥0 j≥0 R -2n ∆x (U II j ) n 2 + n≥0 j≥0 R -2n ∆t (U II 0 ) n 2 ≤ C k≥0 ∆x|f k | 2 .
By setting in the above formula R = e γ∆t for γ > 0 and ∆t > 0, and using the classical lower bound e γ∆t ≥ 1 + γ∆t, we obtain that there exists a constant c > 0 such that γ γ∆t + 1

n≥0 j≥0 e -2γn∆t ∆t∆x (U II j ) n 2 + n≥0 e -2γn∆t ∆t (U II 0 ) n 2 ≤ C k≥0 ∆x|f k | 2 .
Then, for all T > 0, there exists a constant C T > 0 such that

N n=0 j≥0 ∆t∆x (U II j ) n 2 + N n=0 ∆t (U II 0 ) n 2 ≤ C T k≥0 ∆x|f k | 2 ,
with N = T /∆t.

3.4.

Stiff stability analysis. Following Section 2.3, for any T > 0, there exists C T > 0 such that the solution (U III j ) n j∈N to (3.24) satisfies

N n=0 j≥0 ∆x∆t (U III j ) n 2 + N n=0 ∆t (U III 0 ) n 2 ≤ C T N n=0 ∆t B (U I 0 ) n + (U II 0 ) n 2 .
Furthermore, from the inequalities (3.31) and (3.32), one obtains

N n=0 ∆t B (U I 0 ) n + (U II 0 ) n 2 ≤ C T k≥0 ∆x|f k | 2 .
Therefore, we show the uniform estimate on (U

III j ) n j∈N N n=0 j≥0 ∆x∆t (U III j ) n 2 + N n=0 ∆t (U III 0 ) n 2 ≤ C T j≥0 ∆x|f j | 2 , (3.36)
with the positive constant C T independent of ε, ∆x and ∆t.

To complete the proof of the Theorem 1.1 for the numerical scheme of the IBVP (3.1), observe that from the inequalities (3.30) and (3.36), for any T > 0, there exists C T > 0 such that for any (f j ) j∈N ∈ 2 (N, R 2 ) and N := T /∆t, the solution (U n j ) j∈N to (3.1) satisfies

N n=0 j≥0 ∆x∆t|U n j | 2 + N n=0 ∆t|U n 0 | 2 ≤ C T j≥0 ∆x|f j | 2 , (3.37)
where the constant C T independent of ε, ∆x and ∆t. This is the last step to prove the Theorem 1.1.

Numerical experiments.

In this paragraph, we present some numerical experiments for the behavior of the numerical solution (U n j ) j∈N corresponding to whether or not the SKC (1.9) holds. We also look at the numerical solution when the UKC (1.4) is wrong. After that, we observe the effective behavior of the energy terms E 1 inside the domaine and E 2 along the boundary, which are defined in (2.37).

In our numerical experiments, we choose a = 1, B v = 1, λ xt = 1/3, fix the space step ∆x = 5 × 10 -3 , the time step ∆t = λ -1

xt ∆x, and let the relaxation rate ε and the boundary data B u vary. The boundary data is the homogeneous one b n ≡ 0, for any n ∈ N. The initial data is Let us first observe that these data are compatible in the corner (x, t) = (0, 0) in the sense that Bf 0 = 0. Moreover, the choice of an initial data with support in [1/4, 13/30] is motivated by the property of finite speed of propagation available at the continuous side (1.1). More precisely, the exact solution we approximate has characteristic velocities ±1 and therefore vanishes outside some space interval [0, 0.63] for small times in [0, 0.2]. Thus, we choose for our experiments the space interval [0, 1] and the time interval [0, T ) with T = 0.2. Let us mention that the numerical experiments are performed we another discrete right boundary condition at x = 1. This is chosen to be the classical homogeneous first order Neumann extrapolation boundary condition U n J+1 = U n J , for any n ∈ N, at the rightmost cell J. That boundary condition indeed exhibits convenient stability features for both the inflowing and the outflowing transport equation [START_REF] Goldberg | On a boundary extrapolation theorem by Kreiss[END_REF].

f j =      100 × 13 30 -x j x j - 1 4 × 1 -1 T , if x j ∈ 1 4 , 13
3.5.1. The behavior of the numerical solution. Firstly, we choose a set of values B u such that the SKC (1.9) is satisfied with ε = 10 -2 and also with ε = 10 2 . The Figures 3.1 and 3.2 show numerical solution (U n j ) j∈N over the time interval t ∈ [0, 0.2).

In the first case, ε = 10 -2 , due to the initial relaxation of solution to the equilibrium system, the numerical solution descends over time (see Figure 3.1). In the case ε = 10 2 , at time t < 0.2, its solution seems to translate to the left and the ghost solutions do not go backward in space for the implicit scheme. After that, the initial condition re-enters the domain from the left boundary (see Figure 3.2). It is not so much influenced by relaxation source term but more by the boundary dissipation.

Secondly, we choose the value of B u such that the SKC (1.9) is not satisfied. Besides, we also present the numerical solution when the Uniform Kreiss Condition (1.4) is wrong. The Figures 3.3 and 3.4 show numerical solution (U n j ) j∈N over the time interval t ∈ [0, 0.2). We can observe that the numerical solution at the boundary rise gradually when the SKC (1.9) fails. This is the case for example for ε = 10 -2 together with the parameters (B u , B v ) = (-1/2, 1). When the UKC (1.4) does not hold, the behavior is even worse (see Figure 3.4).

3.5.2. The effective behavior of the energy terms. We present hereafter the effective behavior of the energy terms E 1 and E 2 for ε ∈ (0, +∞), T = 0.2 and N = T /∆t.

According to Table 4 and Figures 3.5, 3.6, we can see that • For any ε ∈ (0, +∞), the values of E 1 and E 2 rise gradually when the SKC (1.9) is not satisfied. This is the case for example for ε = 10 2 together with the parameters (B u , B v ) = (-1/2, 1). The behavior is even worse when the UKC (1.4) is not hold. but since x < 0, one obtains by simple considerations the inequality x + √ x 2 + 1 < 1. Case 2: In the general case y = 0, let us begin with some notations: ζ 2 + 1 = p 1 + q 1 i, with p 1 = x 2 -y 2 + 1 and q 1 = 2xy,

ζ 2 + 1 = a 1 + b 1 i, with a 1 = p 1 + p 2 1 + q 2 1 2
and b 1 = sgn(q 1 ) -p 1 + p 2 1 + q 2 1 2 .

Together with these notations, some algebraic identities are available:

1. 3 .

 3 Main result. Dealing with the continuous IBVP (1.1)-(1.

  ) with the notation λ xε = ∆x/ε. According to Lemma A.2 together with the properties Re (µ -(εξ)) < 0 for ε > 0 and Re ξ > 0, we can prove |κ -(εξ)| < 1. Besides, since µ -(εξ) = -µ + (εξ), we get κ + (εξ)κ -(εξ) = 1. As a consequence, for any ε > 0 and Re ξ > 0, one has the separation property |κ + (εξ)| > 1.

  n+1 k=0 C n+1-k U k 0 n≥0 and b stands for the Ztransform of the scalar boundary data: b

Figure 2 . 1 .

 21 Figure 2.1. The numerical solution u(x, t) (left) and v(x, t) (right) for ε = 10 -2 . The SKC (1.9) holds with B u = -4.

Figure 2 . 3 .

 23 Figure 2.3. The numerical solution u(x, t) (left) and v(x, t) (right) for ε = 10 -2 . The SKC (1.9) does not hold with B u = -1/2.

2. 4 . 3 .Figure 2 . 4 .

 4324 Figure 2.4. The numerical solution u(x, t) (left) and v(x, t) (right) for ε = 10 2 . The UKC (1.4) is wrong with B u = -1.

5 ×

 5 10 -2 8.5 × 10 -3 1.3 × 10 -2 2.1 × 10 -2 2.3 × 10 -2 2.4 × 10 -2 25 × 10 -3 3.8 × 10 -3 6.4 × 10 -3 1.1 × 10 -2 1.2 × 10 -2 1.2 × 10 -2 125 × 10 -4 1.8 × 10 -3 3.2 × 10 -3 5.4 × 10 -3 6.1 × 10 -3 6.2 × 10 -3 625 × 10 -5 9.1 × 10 -4 1.6 × 10 -3 2.7 × 10 -3 3.1 × 10 -3 3.1 × 10 -3

Figure 2 . 5 .

 25 Figure 2.5. Energy evolution E 1 for B u = -4 (left) and B u = -0.5 (right).

Figure 2 . 6 .

 26 Figure 2.6. Energy evolution E 2 for B u = -4 (left) and B u = -0.5 (right).

3 Φ

 3 we need a natural boundary condition at x = +∞. Besides, one gets |κ + | > 1 and |κ -| < 1. Thus, the natural boundary condition takes the form (--κ + U 0 (z) + U -1 (z) + w I = 0, Φ + κ + U 0 (z) + U -1 (z) -w II = 0.

2 ≤

 2 .28) According to (3.27) and (3.28), there exists a constant C > 0 such that j∈Z ∆x U I j n C j∈Z ∆x|f j | 2 , for any n ∈ N, (3.29) with the constant C independent of ε and ∆x. This ends the proof of the Proposition 1.2.

Lemma 3 . 2 .

 32 with N = T /∆t. Furthermore, from the inequality (3.29) and ∆x = ∆tλ xt , one obtainsN n=0 ∆t (U I 0 ) n 2 ≤ C T j≥0 ∆x|f j | 2 . (3.31) 3.3. The uniform estimate on (U II j ) n .The following lemma concerns the estimate on (U II j ) n : Assume that the SKC (1.9) is satisfied and let λ xt ≤ 3 √ a/8 be a positive number. Then, for any T > 0, there exists a constant C T > 0 such that for any ∆t > 0 together with ∆x = λ xt ∆t, for any

Figure 3 . 3 .Figure 3 . 4 .

 3334 Figure 3.3. The numerical solution u(x, t) (left) and v(x, t) (right) for ε = 10 -2 . The SKC (1.9) is not valid with B u = -0.5.

(A. 1 )

 1 x 2 + b 2 1 + 1 = a 2 1 + y 2 and y = a 1 b 1 x ,Firstly, we prove the next inequality(A.2) a 1 x 2 + b 2 1 x + a 1 b 2 1 ≥ 0.

  2.3. Stiff stability analysis. Under the SKC, we now consider the Proposition 1.3 with nonzero boundary condition (b n ) n∈N ∈ 2 (N, R) and homogeneous Cauchy data (f j ) j∈N ≡ 0. In order to get the uniform estimate on (U n j ) j∈N , firstly, we prove the following lemma: Lemma 2.6. Assume that the parameters a, ∆x, ∆t > 0 satisfy

	(2.29)	∆x ≤	3	√ 8	a	∆t.

  2 . The Figures 2.1 and 2.2 show the numerical solution (U n j ) j∈N over the time interval t ∈ [0, 1.2).

			The numerical solution u(x,t)			The numerical solution v(x,t)
		0.00							0.0000
		-0.02							-0.0025
		-0.04							-0.0050
	u (x,t)	-0.08 -0.06						v (x,t)	-0.0100 -0.0075
		-0.10						t = 0 t = 0.24	-0.0125						t = 0 t = 0.24
								t = 0.48							t = 0.48
		-0.12						t = 0.72	-0.0150						t = 0.72
								t = 0.96							t = 0.96
		-0.14						t = 1.19	-0.0175						t = 1.19
		0.0	0.2	0.4	x	0.6	0.8	1.0	0.0	0.2	0.4	x	0.6	0.8	1.0

  .

				The numerical solution u(x,t)				The numerical solution v(x,t)
		0.000							0.000
		-0.025							-0.025
		-0.050							-0.050
	u (x,t)	-0.100 -0.075						v (x,t)	-0.100 -0.075
		-0.125					t = 0		-0.125				t = 0
								t = 0.24						t = 0.24
		-0.150					t = 0.48 t = 0.72	-0.150				t = 0.48 t = 0.72
								t = 0.96						t = 0.96
		-0.175					t = 1.19	-0.175				t = 1.19
			0.0	0.2	0.4	x	0.6	0.8	1.0		0.0	0.2	0.4	x	0.6	0.8	1.0
				The numerical solution u(x,t)				The numerical solution v(x,t)
								t = 0							t = 0
								t = 0.48						t = 0.48
			800								400
			600								300
		u (x,t)	400							v (x,t)	200
			200								100
			0								0
			0.0	0.2	0.4	x	0.6	0.8	1.0		0.0	0.2	0.4	x	0.6	0.8	1.0

Figure 2.2. The numerical solution u(x, t) (left) and v(x, t) (right) for ε = 10 2 . The SKC (1.9) holds with B u = -4.

  18 × 10 -2 -1 1.24 × 10 19 7.17 × 10 32 1.89 × 10 58

	-0.5 191420.5 9910910.98 106714237.85
	1	46 × 10 -3	55 × 10 -3	67 × 10 -3
	3	93 × 10 -4	94 × 10 -4	12 × 10 -3

Table 3 .

 3 The energy terms E 1 (left) and E 2 (right).

		The SKC is sa isfied wi h (B u , B v ) = (-4, 1) and a = 1
		0.004	ε = 10 -2 ε = 1					1400000
			ε = 10 2					1200000
		0.003						1000000
	Energy	0.002						600000 800000
		0.001						400000
								200000
		0.000						0
		0.0	0.2	0.4	0.6 TIME	0.8	1.0	1.2	0.0	0.2	0.4	0.6 TIME	0.8	1.0	1.2

  let us multiply the first equation in(3.22) by ((U I j ) n+1 ) T H and sum over Z, one obtains

	j∈Z

  1, still under the SKC, B u + g(εξ)B v is uniformly bounded away from 0 in εξ ∈ C + , and g(εξ) is uniformly bounded in εξ ∈ C + . Moreover one has |κ -(εξ)| < 1, we therefore obtain Similarly, by an application of the Plancherel's theorem for Z-transform, under the SKC, we have

	(3.33)	R -2n (U II 0 ) n 2	|f k | 2 .
		n≥0	k≥0
	(3.34)	R -2n (U II j ) n 2	|κ -(εξ)| 2j	|f k | 2 .
	n≥0 j≥0	j≥0	k≥0
	Following Lemma 2.6, since we assume that the condition (2.29) holds, one gets the following property
		j≥0

Substituting (3.17) into (3.12), the solution ( U j ) j∈N (z) to (3.2a)- (3.2c) is given by

Together with the definitions of Φ ± and Γ in (2.6) and (2.14), respectively, we have

Therefore, the solution ( U j ) j∈N (z) to (3.2a)-(3.2c) can be reformulated as

This ends the proof of Proposition 3.1.

Secondly, we can see that the solution ( U j ) j∈N (z) to (3.2a)-(3.2c) consists of three parts:

where

and 0.0515 0.0715 3 1.9 × 10 -5 0.0620 0.0894 Table 4. The energy terms E 1 (left) and E 2 (right).

• On the boundary x = 0, the value of E 2 for ε = 10 -2 increase slowly. This is due to the effect of incoming solution at the boundary when the initial relaxation of solution tends to the equilibrium system. In the case ε = 10 2 , its value increase fairly rapidly. It is not so much influenced by relaxation source term but more by the boundary dissipation. Clearly, in our numerical experiment, the numerical solution at the boundary x = 0 increase quickly as soon as the SKC (1.9) is not valid. The behavior of numerical solution is even worse if the UKC (1.4) is not satisfied. Indeed, it seems that the SKC (1.9) is also necessary condition to ensure the non-increase rapidly of the energy terms E 1 and E 2 under the effect of the relaxation source term and the boundary dissipation.

Proof. In the half plane {ζ ∈ C : Re ζ ≥ 0}, the complex function h(ζ) is analytic. As usual, we take √ ζ to be the principal branch with the branch cut along the negative real axis.

Let ζ = x + yi with x ≥ 0, y ∈ R and

Then,

Now, we observe that

Therefore,

This ends the proof of Lemma A.1.

Proof. Assume that ζ = x + yi with x < 0 and y ∈ R. Case 1: Consider first the easy case y = 0. Then We can see that the inequality (A.2) is equivalent to a 1 (x 2 + b 2 1 ) ≥ -xb 2 1 and since x < 0, the latter is now equivalent to its squared version, that reads

. By the definition of a 1 , b 1 above, the previous inequality is

But, for any p 1 , q 1 ∈ R, this is easy to see that p 1 -p 2 1 + q 2 1 ≤ 0 and thus any of the previous inequalities and so the expected one (A.2) follow. Now let us observe that the required inequality |ζ

that we prove now. According to the algebraic identities in (A.1), by eliminating the occurences of y, the previous formula is equivalent to

In addition, we observe that

x 2 , and thus the previous inequality is equivalent to

Since x < 0 and from the inequality (A.2), the formula (A.4) reads also

1 ) ≥ 0. This ends the proof of the inequality (A.3). Now since Re ζ < 0, the analysis of the first easy case again applies to get Re ζ + (Re ζ) 2 + 1 < 1. This ends the proof of Lemma A.2.

IRMAR (UMR CNRS 6625), Université de Rennes, Campus de Beaulieu, 35042 Rennes Cedex, France. Email address: benjamin.boutin@univ-rennes1.fr Email address: thi-hoai-thuong.nguyen@univ-rennes1.fr Email address: nicolas.seguin@univ-rennes1.fr