Maria Inês

TOUFIC Maria Inês M Copetti
email: mimcopetti@ufsm.br

AND E L Arwadi
email: t.elarwadi@bau.edu.lb

Wael Youssef
email: youssef.wael@ul.edu.lb

Wael Youssef On

Lana

Toufic El

ON THE THEORETICAL AND NUMERICAL STABILITY OF THE THERMOVISCOELASTIC BRESSE SYSTEM

Keywords:

whether they are published or not. The documents may come

Introduction

In this work, we consider a new model of thermoviscoelastic Bresse beam which is described by the following equations (1.1)

                                                   ρ 1 φ tt -k(φ x + ψ + lw) x -γ 1 (φ x + ψ + lw) xt -k 0 l(ω x -lφ)
-γ 0 l(w x -lφ) t + lαθ 1 = 0 in (0, L) × (0, ∞), ρ 2 ψ tt -bψ xx -γ 2 ψ xxt + k(φ x + ψ + lω) +γ 1 (φ x + lω + ψ) t + αθ 2x = 0 in (0, L) × (0, ∞), ρ 1 ω tt -k 0 (ω x -lφ) x -γ 0 (ω x -lφ) xt + kl(φ x + ψ + lω) +γ 1 l(φ x + ψ + lω) t + αθ 1x = 0 in (0, L) × (0, ∞),

θ 1t -θ 1xx + α(w x -lφ) t = 0 in (0, L) × (0, ∞), θ 2t -θ 2xx + αψ xt = 0 in (0, L) × (0, ∞)
φ(t, x) = ψ(t, x) = ω(t, x) = θ 1 (t, x) = θ 2 (t, x) = 0, for x = 0, L, (φ, φ t , ψ, ψ t , w, w t , θ 1 , θ 2 )(x, 0) = (φ 0 , φ 1 , ψ, ψ 1 , w 0 , w 1 , θ 10 , θ 20 )(x).

with t > 0 and 0 < x < L, where L represents the distance between the ends of the center line of the beam, see figure. The functions φ, ψ and ω denote, respectively, the transverse displacement of a curved beam, the rotation angle of the filament and the longitudinal displacement, while θ 1 and θ 2 are the temperature deviations along the longitudinal and vertical directions (see [START_REF] Liu | Energy Decay of the thermoelastic Bresse System[END_REF]). Here, k 0 = EH, k = GH, b = EI and ρ 1 , ρ 2 , l, G, E, H and α denote positive constants characterizing physical properties of the beam and the filament. In addition, l = 1 R , where R is the radius of curvature, and γ 0 , γ 1 and γ 2 are viscosity coefficients. We denote by u t and u x the derivatives of u with respect to t and x.

The Bresse system is also known as the circular arch problem and is given by the following equations:

(1.2)

           ρ 1 φ tt = -Q x + lN , ρ 2 ψ tt = M x + Q, ρ 1 ω tt = -N x -lQ.
where

N = k 0 (ω x lφ), Q = k(φ x + ψ + lω), M = bψ x .
N , Q and M to denote the axial force, the shear force and the bending moment. The energy of solutions of the system (1.1) is defined by

(1.3)    E(t) = 1 2 ∫ L 0 ( ρ 1 |φ t | 2 + ρ 2 |ψ t | 2 + ρ 1 |ω t | 2 + b|ψ x | 2 +k|φ x + ψ + lω| 2 + k 0 |ω x -lφ| 2 + θ 2 1 + θ 2 2 ) dx.
This energy is decreasing. Indeed, a straightforward calculation gives (1.4)

     d dt E(t) = - ∫ L 0 ( γ 1 |(φ x + ψ + lω) t | 2 + γ 2 |ψ xt | 2 + γ 0 |ω xt -lφ t | 2 +|θ 1x | 2 + |θ 2x | 2 ) dx ≤ 0.
In the study of the thermoviscoelasticity of a beam, we take into account the effects of the viscoelastic and thermal dissipation energies upon the corresponding vibrations.

Many authors studied theoretically the stabilization of the thermoelastic Bresse beam and thermo(visco)elastic Timoshenko beam. In [START_REF] Rivera | Mildly dissipative nonlinear Timoshenko systems -global existence and exponential stability[END_REF], Racke and Rivera considered the thermoelastic Timoshenko system and they established the exponential stability when E = G. Moreover, they also proved nonnexponential stability when E ̸ = G. In [START_REF] Liu | Energy Decay of the thermoelastic Bresse System[END_REF], Liu and Rao studied the energy decay rate for the thermoelastic Bresse system formed by three wave equations coupled with two heat equations. The two wave equations about the longitudinal displacement and shear angle displacement are damped by dissipation from the two heat equations. However, the wave equation about the vertical displacement is only weakly damped. The authors showed that, under conditions on the wave speeds, the energy decays exponentially. Recently, in [START_REF] Gallego | Decay rates for solutions to thermoelastic Bresse systems of types I and III[END_REF], Gallego and Rivera studied the energy decay for the thermoelastic Bresse system in the whole line with two different dissipative mechanism, given by heat conduction (types I and III). They proved that the decay rate of the solutions are very slow.

In [START_REF] Arwadi | Theoretical and Numerical Observability of the Bresse Beam[END_REF], Arwadi, Wehbeh and Youssef studied the observability of the elastic Bresse system in one-dimensional bounded domain. They introduced a new notion of observability which is the discrete observability by using an implicit Euler scheme in time and Lagrange finite elements in space.

On an other side, in [START_REF] Copetti | A dynamic contact problem involving a Timoshenko beam model[END_REF], Copetti and Fernández considered the Timoshenko model for thick beam with viscosity and obtained some a priori error estimates for the numerical solution. In [START_REF] Bernardi | Discretization of a nonlinear dynamic thermoviscoelastic Timoshenko beam model[END_REF] Bernardi and Copetti considered a nonlinear model for a thermoviscoelastic Timoshenko beam that can enter in contact with obstacles. They proposed a discretization based on a combination of backward Euler and Crank-Nicolson schemes in time and finite elements in space and performed the a priori analysis of the discrete problem. They presented also some numerical experiments.

To our knowledge, this is the first paper to consider the thermoviscoelastic Bresse system (1.1) from the theoretical and numerical point of view. This system is a generalization of the thermoviscoelastic Timoshenko beam when l and w are neglected. The difference between this work and the others that study the Bresse system is the consideration of the viscoelasticity and thermoelasticity, which are physically important, and the performing of a priori analysis of a proposed discrete problem. In particular, exponential stability is derived without any restrictions on the wave velocities. Moreover, error estimates are obtained without using the standard technique which compares the solution of the discrete problem to some projection of the exact solution.

Well-posedness

The energy space associated to problem (1.1) is defined by:

H = ( H 1 0 (0, L) ) 3 × ( L 2 (0, L) ) 3 × ( H 1 0 (0, L) ) 2 .
The space H is equipped with the inner product that induces the energy norm. For

U = (φ, ψ, ω, , u, v, z, θ 1 , θ 2 ) ∈ H, we define (2.1) { ∥U ∥ 2 H = ρ 1 ∥u∥ 2 + ρ 2 ∥v∥ 2 + ρ 1 ∥z∥ 2 + b∥ψ x ∥ 2 + k∥φ x + ψ + lω∥ 2 +k 0 ∥ω x -lφ∥ 2 + ∥θ 1 ∥ 2 + ∥θ 2 ∥ 2 ,
where ∥ • ∥ is the L 2 (0, L) norm. Using (3.1), we remark that the norm ∥U ∥ H is equivalent to the usual norm of H which is

(2.2) ∥U ∥ 2 = ∥φ x ∥ 2 + ∥ψ x ∥ 2 + ∥ω x ∥ 2 + ∥u∥ 2 + ∥v∥ 2 + ∥z∥ 2 + ∥θ 1 ∥ 2 + ∥θ 2 ∥ 2 .
Now, we introduce the linear unbounded operator A in H as follows:

AU =                                     u v z k ρ 1 (φ xx + ψ x + lω x ) + lk 0 ρ 1 (ω x -lφ) + γ 1 ρ 1 (φ x + ψ + lw) xt + k 0 ρ 1 l(ω x -lφ) + lγ 0 ρ 1 (w x -lφ) t l ρ 1 αθ 1 b ρ 2 ψ xx - k ρ 2 (φ x + ψ + lω) + γ 2 ρ 2 ψ xxt - k ρ 2 (φ x + ψ + lω) - γ 1 ρ 2 (φ x + lω + ψ) t - α ρ 2 θ 2x k 0 ρ 1 (ω xx -lφ x ) - lk ρ 1 (φ x + ψ + lω) + γ 0 ρ 1 (ω x -lφ) xt - lγ 1 ρ 1 (φ x + ψ + lω) t -     , D(A) = { (φ, ψ, ω, u, v, z, θ 1 , θ 2 ) ∈ H; φ, ψ, ω, θ 1 , θ 2 ∈ H 1 0 (0, L) ∩ H 2 (0, L), u, v, z ∈ H 1 0 (0, L) } .
Then, by substituting u, v and z by u t .v t and z t the thermoviscoelastic Bresse beam system can be written in the form of a first order evolution on the Hilbert space H (2.3)

{ U ′ = AU, U (0) = U 0 .
Theorem 2.1. A generates a C 0 semigroup S(t) of contractions on H.

Proof. Firstly, it is clear that D(A) is dense in H and

⟨AU, U ⟩ H = - ( γ 1 ∥(φ x +ψ+lω) t ∥ 2 +γ 2 ∥ψ xt ∥ 2 +γ 0 ∥ω xt -lφ t ∥ 2 +∥θ 1x ∥ 2 +∥θ 2x ∥ 2 ) dx ≤ 0,
where ⟨ . , . ⟩ H is the scalar product associated to the energy norm defined in (2.1). It remains to prove that

(2.4) AU = f, ∀f ∈ D(A),
where f = (f 1 , f 2 , ..., f 8 ), has a unique solution U ∈ D(A). Indeed, the first three equations in (2.4) give

u = f 1 , v = f 2 , and z = f 3 .
Substituting u, v, and z into the last two equation in (2.4) and using the standard elliptic PDE theory yield to the existence of the unique solution

θ 1 , θ 2 ∈ H 1 0 ∩ H 2 .
On the other hand, if we consider the bilinear form defined by

B ( (φ, ψ, ω), ( φ), ψ, ω ) := ∫ L 0 [ k(φ x + ψ + lω)( φ x + ψ + l ω) + bψ x ψ x + k 0 (ω x -lφ)( ω x -l φ) ] dx,
and if we use the Lax-Milgram theorem we conclude that there exists a unique (φ, ψ, ω) solution of the fourth, fifth, and sixth equations of (2.4). Hence, 0 ∈ ρ(A) and so by the resolvent identity, for small λ > 0 we have R(λ -A) = H. Consequently, the Lumer-Phillip theorem implies that A generates a C 0 semigroup S(t) of contractions on H. (For details see theorem 1.2.4 in [START_REF] Liu | Semigroups Associated with Dissipative Systems[END_REF]) □

The following existence and uniqueness result follows from the latter theorem.

Theorem 2.2.

If U 0 = (φ 0 , ψ 0 , ω 0 , φ 1 , ψ 1 , ω 1 , θ 10 , θ 20 ) ∈ D(A), problem (1.1) ad- mits a unique classical solution U = (φ, ψ, ω, φ t , ψ t , ω t , θ 1 , θ 2 ) such that U ∈ C ( [0, +∞); D(A) ) ∩ C 1 ([0, +∞); H ) .

Dissipation of the energy

In this section, we shall prove the exponential stabilization of the Bresse system which is based on the following inequality used in Haraux [START_REF] Haraux | Semi-groupes lineéaires et équations d'évolution linéaires périodiques[END_REF] and Lagnese [START_REF] Lagnese | Boundary Stabilization of Thin Plates[END_REF] (see [START_REF] Komornik | Exact controllability and stabilization. The Multiplier Method[END_REF] for the proof). Lemma 3.1. Let E : R + → R + be a non-increasing function and assume that there exists a constant T > 0 such that

∫ ∞ t E(s)ds ≤ T E(t), ∀t ∈ R + .
Then,

E(t) ≤ E(0)e 1-t T , ∀t ≥ T.
The next result will be needed in the proof of the exponential stabilization (see [START_REF] Youssef | Contrôle et stabilisation de systèmes élastiques couplés[END_REF]). Lemma 3.2. There exists a positive constant C such that for every (φ, ψ, ω) ∈

V := ( H 1 0 (0, L)
) 3 , we have

(3.1) ∫ L 0 ( |φ x | 2 +|ψ x | 2 +|ω x | 2 ) dx ≤ C ∫ L 0 ( |ψ x | 2 +|φ x +ψ+lω| 2 + | ω x -lφ | 2 ) dx.
We obtain the same inequality if we replace the functions φ, ψ and ω by φ t , ψ t and

ω t i.e. (3.2) ∫ L 0 ( |φ xt | 2 +|ψ xt | 2 +|ω xt | 2 ) dx ≤ C ∫ L 0 ( |ψ xt | 2 +|(φ x +ψ+lω) t | 2 + | (ω x -lφ) t | 2 ) dx.
Now, we are ready to show that the energy will decrease to zero exponentially when t tends to +∞. Theorem 3.1. There exist two positive constants C 1 and η such that the energy of the solution of (1.1) satisfies

(3.3) E(t) ≤ C 1 E(0)e -ηt ∀t ≥ 0.
Proof. Several steps are required for the proof of this theorem. For simplicity, we shall denote by

∫ T S ∫ L 0 f dx dt = ∫ T S ∫ L 0 f .
Step 1. Multiplying the first equation of (1.1) by φ, the second by ψ, the third by ω, the fourth by θ 1 , the fifth by θ 2 and integrating over [S, T ] × (0, L), we obtain:

(3.4)                - ∫ T S ∫ L 0 ρ 1 |φ t | 2 + ρ 1 ∫ L 0 [ φ t φ ] T S dx + k ∫ T S ∫ L 0 (φ x + ψ + lω)φ x -γ 1 ∫ T S ∫ L 0 (φ x + ψ + lw) xt φ -k 0 l ∫ T S ∫ L 0 (ω x -lφ)φ -lγ 0 ∫ T S ∫ L 0 (w x -lφ) t φ + lα ∫ T S ∫ L 0 θ 1 φ = 0, (3.5)
               - ∫ T S ∫ L 0 ρ 2 |ψ t | 2 + ρ 2 ∫ L 0 [ ψ t ψ ] T S dx + b ∫ T S ∫ L 0 |ψ x | 2 +k ∫ T S ∫ L 0 (φ x + ψ + lω)ψ -γ 2 ∫ T S ∫ L 0 ψ xxt ψ +γ 1 ∫ T S ∫ L 0 (φ x + lω + ψ) t ψ + α ∫ T S ∫ L 0 θ 2x ψ = 0, (3.6)
               - ∫ T S ∫ L 0 ρ 1 |ω t | 2 + ρ 1 ∫ L 0 [ ω t ω ] T S dx + k 0 ∫ T S ∫ L 0 (ω x -lφ)ω x +lk ∫ T S ∫ L 0 (φ x + ψ + lω)ω -γ 0 ∫ T S ∫ L 0 (ω x -lφ) xt ω +γ 1 l ∫ T S ∫ L 0 (φ x + ψ + lω) t ω + α ∫ T S ∫ L 0 θ 1x ω = 0. (3.7)        1 2 ∫ L 0 [ |θ 1 (x, t)| 2 ] T S dx + ∫ T S ∫ L 0 |θ 1x | 2 +α ∫ T S ∫ L 0 (w x -lφ) t θ 1 = 0, and (3.8)
       1 2 ∫ L 0 [ |θ 2 (x, t)| 2 ] T S dx + ∫ T S ∫ L 0 |θ 2x | 2 +α ∫ T S ∫ L 0 ψ xt θ 2 = 0.
For simplicity we shall denote by I 1 -

Step 2. Adding (3.4), (3.5), (3.6), (3.7), and (3.8) leads to (3.9)

                                                         - ∫ T S ∫ L 0 ρ 1 |φ t | 2 - ∫ T S ∫ L 0 ρ 2 |ψ t | 2 - ∫ T S ∫ L 0 ρ 1 |ω t | 2 +k ∫ T S ∫ L 0 | φ x + ψ + lω | 2 +k 0 ∫ T S ∫ L 0 | ω x -lφ | 2 +b ∫ T S ∫ L 0 |ψ x | 2 +ρ 1 ∫ L 0 [ φ t φ ] T S dx + ρ 2 ∫ L 0 [ ψ t ψ ] T S dx + ρ 1 ∫ L 0 [ ω t ω ] T S dx + γ 1 2 ∫ T S ∫ L 0 ( (φ x + ψ + lω) 2 ) t + γ 2 2 ∫ T S ∫ L 0 (ψ 2 x ) t + γ 0 2 ∫ T S ∫ L 0 ( (ω x -lφ) 2 ) t + 1 2 ∫ L 0 [ |θ 1 (x, t)| 2 ] T S dx + ∫ T S ∫ L 0 |θ 1x | 2 +α ∫ T S ∫ L 0 (w x -lφ) t θ 1 + 1 2 ∫ L 0 [ |θ 2 (x, t)| 2 ] T S dx + ∫ T S ∫ L 0 |θ 2x | 2 +α ∫ T S ∫ L 0 ψ xt θ 2 + lα ∫ T S ∫ L 0 θ 1 φ + α ∫ T S ∫ L 0 θ 2x ψ +α ∫ T S ∫ L 0 θ 1x ω = 0.
From the definition of E we deduce

(3.10)                k ∫ L 0 | φ x + ψ + lω | 2 dx + k 0 ∫ L 0 | ω x -lφ | 2 dx + b ∫ L 0 |ψ x | 2 dx = 2E(t) - ∫ L 0 ρ 1 |φ t | 2 dx - ∫ L 0 ρ 2 |ψ t | 2 dx - ∫ L 0 ρ 1 |ω t | 2 dx - ∫ L 0 |θ 1 | 2 dx - ∫ L 0 |θ 2 | 2 dx.
By combining (3.9) and (3.10) we get

(3.11)                                                    2 ∫ T S E(t)dt = -ρ 1 ∫ L 0 [ φ t φ ] T S dx + ρ 2 ∫ L 0 [ ψ t ψ ] T S dx + ρ 1 ∫ L 0 [ ω t ω ] T S dx +2 ( ∫ T S ∫ L 0 ρ 1 |φ t | 2 + ∫ T S ∫ L 0 ρ 2 |ψ t | 2 + ∫ T S ∫ L 0 ρ 1 |ω t | 2 ) + ∫ T S ∫ L 0 |θ 1 | 2 + ∫ T S ∫ L 0 |θ 2 | 2 - ∫ T S ∫ L 0 |θ 1x | 2 - ∫ T S ∫ L 0 |θ 2x | 2 - γ 1 2 ∫ T S ∫ L 0 ( (φ x + ψ + lω) 2 ) t - γ 2 2 ∫ T S ∫ L 0 (ψ 2 x ) t - γ 0 2 ∫ T S ∫ L 0 ( (ω x -lφ) 2 ) t - 1 2 ∫ L 0 [ |θ 1 (x, t)| 2 ] T S dx - 1 2 ∫ L 0 [ |θ 2 (x, t)| 2 ] T S dx -α ∫ T S ∫ L 0 (w x -lφ) t θ 1 -α ∫ T S ∫ L 0 ψ xt θ 2 -lα ∫ T S ∫ L 0 θ 1 φ -α ∫ T S ∫ L 0 θ 2x ψ -α ∫ T S ∫ L 0 θ 1x ω. The first term - γ 1 2 ∫ T S ∫ L 0 ( (φ x + ψ + lω) 2 ) t is given by - γ 1 2 ∫ T S ∫ L 0 ( (φ x + ψ + lω) 2 ) t = - γ 1 2 ∫ L 0 [ (φ x + ψ + lω) 2 ] T S dx ≤ γ 1 2 ∫ L 0 [ (φ x + ψ + lω) 2 ] T S ≤ γ 1 2 ∫ L 0 [ (φ x + ψ + lω) 2 (T ) ≤ γ 1 k E(T ) (by the definition of E) ≤ γ 1 k E(S) (since E is decreasing).
Hence,

(3.12) - γ 1 2 ∫ T S ∫ L 0 ( (φ x + ψ + lω) 2 ) t ≤ γ 1 k E(S).
Similarly, we have

(3.13) - γ 2 2 ∫ T S ∫ L 0 (ψ 2 x ) t ≤ γ 2 b E(S)
and

(3.14) - γ 0 2 ∫ T S ∫ L 0 ( (ω x -lφ) 2 ) t ≤ γ 0 k 0 E(S).
Therefore, using (3.12), (3.13), and (3.14) in (3.11) yields

(3.15)                                          2 ∫ T S E(t)dt ≤ -ρ 1 ∫ L 0 [ φ t φ ] T S dx + ρ 2 ∫ L 0 [ ψ t ψ ] T S dx + ρ 1 ∫ L 0 [ ω t ω ] T S dx +2 ( ∫ T S ∫ L 0 ρ 1 |φ t | 2 + ∫ T S ∫ L 0 ρ 2 |ψ t | 2 + ∫ T S ∫ L 0 ρ 1 |ω t | 2 ) + ∫ T S ∫ L 0 |θ 1 | 2 + ∫ T S ∫ L 0 |θ 2 | 2 - 1 2 ∫ L 0 [ |θ 1 (x, t)| 2 ] T S dx - 1 2 ∫ L 0 [ |θ 2 (x, t)| 2 ] T S dx -α ∫ T S ∫ L 0 (w x -lφ) t θ 1 -α ∫ T S ∫ L 0 ψ xt θ 2 -lα ∫ T S ∫ L 0 θ 1 φ -α ∫ T S ∫ L 0 θ 2x ψ -α ∫ T S ∫ L 0 θ 1x ω +( γ 0 k 0 + γ 1 k + γ 2 b )E(S) ∀S ∈ (0, T ).
Step 3. In this step, we shall estimate the terms of the right member of (3.15). By using the Young inequality, the Poincaré inequality, and the estimation (3.1) we deduce

(3.16) -ρ 1 ∫ L 0 [ φ t φ ] T S dx ≤ cE(S) ∀S ∈ (0, T ), (3.17) -ρ 2 ∫ L 0 [ ψ t ψ ] T S dx ≤ cE(S) ∀S ∈ (0, T ), and
(3.18) -ρ 1 ∫ L 0 [ ω t ω ] T S dx ≤ cE(S) ∀S ∈ (0, T ),
where c is a generic constant that depends on ρ 1 , ρ 2 , k 0 , k, b, and the Poincaré constant.

On the other hand, the Poincaré inequality and the relation (3.2) lead to

2 ( ∫ T S ∫ L 0 (ρ 1 |φ t | 2 + ρ 2 |ψ t | 2 + ρ 1 |ω t | 2 ) ) ≤ c ( ∫ T S ∫ L 0 (|φ xt | 2 + |ψ xt | 2 + |ω xt | 2 ) ) ≤ c ( ∫ T S ∫ L 0 (|(φ x + ψ + lω) t | 2 + |ψ xt | 2 + |(ω x -lφ) t | 2 ) ) ≤ cE(S) ∀S ∈ (0, T ). But, ( ∫ T S ∫ L 0 (|(φ x + ψ + lω) t | 2 + |ψ xt | 2 + |(ω x -lφ) t | 2 ) ) ≤ c ∫ T S (-E ′ )dt ≤ cE(S) ∀S ∈ (0, T ).
Therefore,

(3.19) 2 ( ∫ T S ∫ L 0 (ρ 1 |φ t | 2 + ρ 2 |ψ t | 2 + ρ 1 |ω t | 2 ) ) ≤ cE(S) ∀S ∈ (0, T ).
On the other hand, from the Poincaré inequality, easily we conclude the following estimation

(3.20)          ∫ T S ∫ L 0 |θ 1 | 2 + ∫ T S ∫ L 0 |θ 2 | 2 ≤ c ( ∫ T S ∫ L 0 |θ 1x | 2 + ∫ T S ∫ L 0 |θ 2x | 2 ) c ∫ T S (-E ′ )dt ≤ cE(S) ∀S ∈ (0, T ).
Now, the Cauchy-Schwarz inequality gives us

(3.21) ∫ T S ∫ L 0 (w x -lφ) t θ 1 ≤ ∫ T S ( ∫ L 0 |(w x -lφ) t | 2 dx ) 1 2 ( ∫ L 0 |θ 1 | 2 dx ) 1 2 dt.
But, from (1.4) and the definition of E, we have

∫ L 0 |(w x -lφ) t | 2 dx ≤ c ∫ T S (-E ′ )dt, ∫ L 0 |θ 1 | 2 dx ≤ 2E(t) ∀t > 0.
Hence, (3.21) implies

∫ T S ∫ L 0 (w x -lφ) t θ 1 ≤ c ∫ T S (-E ′ ) 1 2 E 1 2 dt.
Now, for all ε > 0, Young's inequality gives

∫ T S ∫ L 0 (w x -lφ) t θ 1 ≤ c ε ∫ T S (-E ′ )dt + ε ∫ T S Edt
and so

(3.22) ∫ T S ∫ L 0 (w x -lφ) t θ 1 ≤ cE(S) + ε ∫ T S Edt ∀S ∈ (0, T ).
Repeating the same argument to

∫ T S ∫ L 0 ψ xt θ 2 , we conclude that (3.23) ∫ T S ∫ L 0 ψ xt θ 2 ≤ cE(S) + ε ∫ T S E(t)dt ∀S ∈ (0, T ).
Concerning the term -1 2

∫ L 0 [ |θ 1 | 2 ]
T S dx, we have:

- 1 2
∫ L 0 [ |θ 1 | 2 ] T S dx = 1 2 ∫ L 0 |θ 1 (S)| 2 dx- 1 2
∫ L 0 |θ 1 (T )| 2 dx ≤ 1 2 ∫ L 0 |θ 1 (T )| 2 dx ≤ E(T ) ≤ E(S),
since E is decreasing. Therefore,

(3.24) - 1 2
∫ L 0 [ |θ 1 | 2 ] T S dx ≤ E(S)
and by the same calculation we have

(3.25) - 1 2
∫ L 0 [ |θ 2 | 2 ] T S dx ≤ E(S).
The estimation of the terms in the fifth row of (3.15) can be obtained by using Poincaré's inequality, (3.1), and (1.4). So,

(3.26)    -lα ∫ T S ∫ L 0 θ 1 φ -α ∫ T S ∫ L 0 θ 2x ψ -α ∫ T S ∫ L 0 θ 1x ω ≤ cE(S) ∀S ∈ (0, T ).
Step 4. Using (3.16), (3.17

∫ T S E(t)dt ≤ cE(S), ∀S ∈ (0, T ).
Thus, if we fix S and make T -→ +∞, the proof is complete by applying lemma 3.1. □

Numerical approximation

In this section, we propose a numerical approximation to the solution of the continuous problem (1.1). Let us introduce the functions φ = φ t , ψ = ψ t , ω = ω t and the weak form associated to system (1.1) obtained by multiplying the equations by test functions ζ, χ, ξ, η 1 , η 2 ∈ H 1 0 (0, L) and integrating by parts:

(4.1)

                                           ρ 1 ( φt , ζ) + k(φ x + ψ + lω, ζ x ) + γ 1 ( φx + ψ + lω, ζ x ) -k 0 l(ω x -lφ, ζ) -γ 0 l(ω x -l φ, ζ) + lα(θ 1 , ζ) = 0, ρ 2 ( ψt , χ) + b(ψ x , χ x ) + γ 2 ( ψx , χ x ) + k(φ x + ψ + lω, χ) +γ 1 ( φx + ψ + lω, χ) -α(θ 2 , χ x ) = 0, ρ 1 (ω t , ξ) + k 0 (ω x -lφ, ξ x ) + γ 0 (ω x -l φ, ξ x ) + kl(φ x + ψ + lω, ξ) +γ 1 l( φx + ψ + lω, ξ) -α(θ 1 , ξ x ) = 0, (θ 1t , η 1 ) + (θ 1x , η 1x ) + α(ω x -l φ, η 1 ) = 0, (θ 2t , η 2 ) + (θ 2x , η 2x ) + α( ψx , η 2 ) = 0.
Partition the interval Ω into subintervals

I j = (x j-1 , x j ) of length h = L/s, with 0 = x 0 < x 1 < • • • < x s = L,
and denote by S h ⊂ H 1 0 (0, L) the space of continuous piecewise linear functions defined on this partition. For a given final time T > 0 and a positive integer N, let ∆t = T /N be the time step.

The finite element method using the backward Euler scheme is to find φn

h , ψn h , ωn h , θ n 1h , θ n 2h ∈ S h such that, for n = 1, . . . , N and for all ζ h , χ h , ξ h , η 1h , η 2h ∈ S h (4.2)                                                      ρ1 ∆t ( φn h -φn-1 h , ζ h ) + k(φ n hx + ψ n h + lω n h , ζ hx ) +γ 1 ( φn hx + ψn h + lω n h , ζ hx ) -k 0 l(ω n hx -lφ n h , ζ h ) -γ 0 l(ω n hx -l φn h , ζ h ) + lα(θ n 1h , ζ h ) = 0, ρ2 ∆t ( ψn h -ψn-1 h , χ h ) + b(ψ n hx , χ hx ) + γ 2 ( ψn hx , χ hx ) +k(φ n hx + ψ n h + lω n h , χ h ) + γ 1 ( φn hx + ψn h + lω n h , χ h ) -α(θ n 2h , χ hx ) = 0, ρ1 ∆t (ω n h -ωn-1 h , ξ h ) + k 0 (ω n hx -lφ n h , ξ hx ) + γ 0 (ω n hx -l φn h , ξ hx ) +kl(φ n hx + ψ n h + lω n h , ξ h ) + γ 1 l( φn hx + ψn h + lω n h , ξ h ) -α(θ n 1h , ξ hx ) = 0, 1 ∆t (θ n 1h -θ n-1 1h , η 1h ) + (θ n 1hx , η 1hx ) + α(ω n hx -l φn h , η 1h ) = 0, 1 ∆t (θ n 2h -θ n-1 2h , η 2h ) + (θ n 2hx , η 2hx ) + α( ψn hx , η 2h ) = 0, where (4.3) φn h = φ n h -φ n-1 h ∆t , ψn h = ψ n h -ψ n-1 h ∆t , ωn h = ω n h -ω n-1 h ∆t are approximations to φ t (t n ), ψ t (t n ) and ω t (t n ), respectively. Here, φ 0 h , φ0 h , ψ 0 h , ψ0 h , ω 0 h , ω0
h , θ 0 1h , θ 0 2h are given approximations to the initial conditions φ 0 , φ 1 , ψ 0 , ψ 1 , ω 0 , ω 1 , θ 10 , θ 20 , respectively.

The following equality will often be used:

(4.4) (a -b, a) = 1 2 ( ∥a -b∥ 2 + ∥a∥ 2 -∥b∥ 2 )
The next result is a discrete version of the energy decay property satisfied by the solution of system (1.1) Theorem 4.1. Let the discrete energy be

E n h = 1 2 ( ρ 1 (∥ φn h ∥ 2 + ∥ω n h ∥ 2 ) + ρ 2 ∥ ψn h ∥ 2 + b∥ψ n hx ∥ 2 + k∥φ n hx + ψ n h + lω n h ∥ 2 +k 0 ∥ω n hx -lφ n h ∥ 2 + ∥θ n 1h ∥ 2 + ∥θ n 2h ∥ 2
) .

Then, the decay property

E n h -E n-1 h ∆t ≤ 0 holds for n = 1, 2, . . . , N. Proof. Taking ζ h = φn h , χ h = ψn h , ξ h = ωn h , η 1h = θ n 1h and η 2h = θ n 2h it results that (4.5)      ρ1 2∆t ( ∥ φn h -φn-1 h ∥ 2 + ∥ φn h ∥ 2 -∥ φn-1 h ∥ 2 ) + k(φ n hx + ψ n h + lω n h , φn hx ) +γ 1 ( φn hx + ψn h + lω n h , φn hx ) -k 0 l(ω n hx -lφ n h , φn h ) -γ 0 l(ω n hx -l φn h , φn h ) + lα(θ n 1h , φn h ) = 0, (4.6)
       ρ2 2∆t ( ∥ ψn h -ψn-1 h ∥ 2 + ∥ ψn h ∥ 2 -∥ ψn-1 h ∥ 2 ) + k(φ n hx + ψ n h + lω n h , ψn h ) + b 2∆t ( ∥ψ n hx -ψ n-1 hx ∥ 2 + ∥ψ n hx ∥ 2 -∥ψ n-1 hx ∥ 2 ) + γ 2 ∥ ψn hx ∥ 2 +γ 1 ( φn hx + ψn h + lω n h , ψn h ) -α(θ n 2h , ψn hx ) = 0, (4.7)
     ρ1 2∆t ( ∥ω n h -ωn-1 h ∥ 2 + ∥ω n h ∥ 2 -∥ω n-1 h ∥ 2 ) + kl(φ n hx + ψ n h + lω n h , ωn h ) +γ 1 l( φn hx + ψn h + lω n h , ωn h ) + k 0 (ω n hx -lφ n h , ωn hx ) +γ 0 (ω n hx -l φn h , ωn hx ) -α(θ n 1h , ωn hx ) = 0, (4.8)
{ 1 2∆t ( ∥θ n 1h -θ n-1 1h ∥ 2 + ∥θ n 1h ∥ 2 -∥θ n-1 1h ∥ 2 ) + ∥θ n 1hx ∥ 2 +α(ω n hx -l φn h , θ n 1h ) = 0,
and (4.9)

{ 1 2∆t ( ∥θ n 2h -θ n-1 2h ∥ 2 + ∥θ n 2h ∥ 2 -∥θ n-1 2h ∥ 2 ) + ∥θ n 2hx ∥ 2 +α( ψn hx , θ n 2h ) = 0.
Recalling (4.3) and (4.4) we deduce that

k(φ n hx + ψ n h + lω n h , φn hx + ψn h + lω n h ) ≥ k 2∆t ( ∥φ n hx + ψ n h + lω n h ∥ 2 -∥φ n-1 hx + ψ n-1 h + lω n-1 h ∥ 2
) .

and that

k 0 (ω n hx -lφ n h , ωn hx -l φn h ) ≥ k 0 2∆t ( ∥ω n hx -lφ n h ∥ 2 -∥ω n-1 hx -lφ n-1 h ∥ 2
) .

Thus, after summing equations (4.5)-(4.9), we arrive at ρ 1 2∆t

( ∥ φn h ∥ 2 -∥ φn-1 h ∥ 2 + ∥ω n h ∥ 2 -∥ω n-1 h ∥ 2 ) + ρ 2 2∆t ( ∥ ψn h ∥ 2 -∥ ψn-1 h ∥ 2 ) + b 2∆t ( ∥ψ n hx ∥ 2 -∥ψ n-1 hx ∥ 2 ) + γ 2 ∥ ψn hx ∥ 2 + 1 2∆t ( ∥θ n 1h ∥ 2 -∥θ n-1 1h ∥ 2 + ∥θ n 2h ∥ 2 -∥θ n-1 2h ∥ 2 ) + ∥θ n 1hx ∥ 2 + ∥θ n 2hx ∥ 2 + k 2∆t ( ∥φ n hx + ψ n h + lω n h ∥ 2 -∥φ n-1 hx + ψ n-1 h + lω n-1 h ∥ 2 ) + k 0 2∆t ( ∥ω n hx -lφ n h ∥ 2 -∥ω n-1 hx -lφ n-1 h ∥ 2
)

+γ 1 ∥ φn hx + ψn h + lω n h ∥ 2 + γ 0 ∥ω n hx -l φn h ∥ 2 ≤

Error estimate

Now, we will obtain some a priori estimates for the difference between the exact solution and the numerical solution.

Theorem 5.1. There exists a positive constant C, independent of the discretization parameters h and ∆t, such that for all

{ζ i h , χ i h , ξ i h , η i 1h , η i 2h } n i=0 ⊂ S h , ∥ φn -φn h ∥ 2 + ∥ ψn -ψn h ∥ 2 + ∥ω n -ωn h ∥ 2 + ∥ψ n x -ψ n hx ∥ 2 +∥φ n x + ψ n + lω n -(φ n hx + ψ n h + lω n h )∥ 2 +∥ω n x -lφ n -(ω n hx -lφ n h )∥ 2 + ∥θ n 1 -θ n 1h ∥ 2 + ∥θ n 2 -θ n 2h ∥ 2 ≤ C∆t n ∑ i=1 ( ∥ φi t - 1 ∆t ( φi -φi-1 )∥ 2 + ∥ ψi t - 1 ∆t ( ψi -ψi-1 )∥ 2 +∥ω i t - 1 ∆t (ω i -ωi-1 )∥ 2 + ∥θ i 1t - 1 ∆t (θ i 1 -θ i-1 1 )∥ 2 +∥θ i 2t - 1 ∆t (θ i 2 -θ i-1 2 )∥ 2 + ∥ φi -ζ i h ∥ 2 + ∥ φi x -ζ i hx ∥ 2 +∥ ψi -χ i h ∥ 2 + ∥ ψi x -χ i hx ∥ 2 + ∥ω i -ξ i h ∥ 2 + ∥ω i hx -ξ i hx ∥ 2 +∥θ i 1 -η i 1h ∥ 2 + ∥θ i 1x -η i 1hx ∥ 2 + ∥θ i 2 -η i 2h ∥ 2 + ∥θ i 2x -η i 2hx ∥ 2 ) + C ∆t n-1 ∑ i=1 ( ∥ φi -ζ i h -( φi+1 -ζ i+1 h )∥ 2 + ∥ ψi -χ i h -( ψi+1 -χ i+1 h )∥ 2 +∥ω i -ξ i h -(ω i+1 -ξ i+1 h )∥ 2 + ∥θ i 1 -θ i 1h -(θ i+1 1 -η i+1 1h )∥ 2 +∥θ i 2 -θ i 2h -(θ i+1 2 -η i+1 2h )∥ 2 ) +C ( ∥φ 1 -φ0 h ∥ 2 + ∥ψ 1 -ψ0 h ∥ 2 + ∥ω 1 -ω0 h ∥ 2 + ∥ψ 0 -ψ 0 h ∥ 2 +∥φ 0 x + ψ 0 + lω 0 -(φ 0 hx + ψ 0 h + lω 0 h )∥ 2 +∥ω 0 x -lφ 0 -(ω 0 hx -lφ 0 h )∥ 2 + ∥θ 0 1 -θ 0 1h ∥ 2 + ∥θ 0 2 -θ 0 2h ∥ 2
) .

Proof. The proof is divided into three steps to help the reader follow the arguments and the calculations.

Step 1. For a continuous function

v(t), let v n = v(t n ). Subtracting variational equation (4.1) 1 , at time t n for ζ = ζ h ∈ S h , and the discrete variational equation (4.2) 1 results in ρ 1 ( φn t - 1 ∆t ( φn h -φn-1 h ), ζ h ) + k(φ n x + ψ n + lω n -(φ n hx + ψ n h + lω n h ), ζ hx ) +γ 1 ( φn x + ψn + lω n -( φn hx + ψn h + lω n h ), ζ hx ) -k 0 l(ω n x -lφ n -(ω n hx -lφ n h ), ζ h ) -γ 0 l(ω n x -l φn -(ω n hx -l φn h ), ζ h ) +lα(θ n 1 -θ n 1h , ζ h ) = 0.
Thus,

(5.1)

                                   ρ 1 ( φn t -1 ∆t ( φn h -φn-1 h ), φn -φn h ) +k(φ n x + ψ n + lω n -(φ n hx + ψ n h + lω n h ), φn x -φn hx ) +γ 1 ( φn x + ψn + lω n -( φn hx + ψn h + lω n h ), φn x -φn hx ) -k 0 l(ω n x -lφ n -(ω n hx -lφ n h ), φn -φn h ) -γ 0 l(ω n x -l φn -(ω n hx -l φn h ), φn -φn h ) + lα(θ n 1 -θ n 1h , φn -φn h ) = ρ 1 ( φn t -1 ∆t ( φn h -φn-1 h ), φn -ζ h ) +k(φ n x + ψ n + lω n -(φ n hx + ψ n h + lω n h ), φn x -ζ hx ) +γ 1 ( φn x + ψn + lω n -( φn hx + ψn h + lω n h ), φn x -ζ hx ) -k 0 l(ω n x -lφ n -(ω n hx -lφ n h ), φn -ζ h ) -γ 0 l(ω n x -l φn -(ω n hx -l φn h ), φn -ζ h ) + lα(θ n 1 -θ n 1h , φn -ζ h ).
Now, we subtract variational equation (4.1) 2 , at time t n for χ = χ h ∈ S h , and the discrete variational equation (4.2) 2 to find

ρ 2 ( ψn t - 1 ∆t ( ψn h -ψn-1 h ), χ h ) + b(ψ n x -ψ n hx , χ hx ) + γ 2 ( ψn x -ψn hx , χ hx ) +k(φ n x + ψ n + lω n -(φ n hx + ψ n h + lω n h ), χ h ) +γ 1 ( φn x + ψn + lω n -( φn hx + ψn h + lω n h ), χ h ) -α(θ n 2 -θ n 2h , χ hx ) = 0. Therefore, (5.2)
                                     ρ 2 ( ψn t -1 ∆t ( ψn h -ψn-1 h ), ψn -ψn h ) + b(ψ n x -ψ n hx , ψn x -ψn hx ) +γ 2 ( ψn x -ψn hx , ψn x -ψn hx ) +k(φ n x + ψ n + lω n -(φ n hx + ψ n h + lω n h ), ψn -ψn h ) +γ 1 ( φn x + ψn + lω n -( φn hx + ψn h + lω n h ), ψn -ψn h ) -α(θ n 2 -θ n 2h , ψn x -ψn hx ) = ρ 2 ( ψn t -1 ∆t ( ψn h -ψn-1 h ), ψn -χ h ) + b(ψ n x -ψ n hx , ψn x -χ hx ) +γ 2 ( ψn x -ψn hx , ψn x -χ hx ) +k(φ n x + ψ n + lω n -(φ n hx + ψ n h + lω n h ), ψn -χ h ) +γ 1 ( φn x + ψn + lω n -( φn hx + ψn h + lω n h ), ψn -χ h ) -α(θ n 2 -θ n 2h , ψn x -χ hx ).
Next, we turn to the longitudinal displacement. We subtract (4.1) 3 , at time t n for ξ = ξ h ∈ S h , and the discrete variational equation (4.2) 3 to get

ρ 1 ( ωn t - 1 ∆t (ω n h -ωn-1 h ), ξ h ) + γ 0 (ω n x -l φn -(ω n hx -l φn h ), ξ hx ) +kl(φ n x + ψ n + lω n -(φ n hx + ψ n h + lω n h ), ξ h ) +γ 1 l( φn x + ψn + lω n -( φn hx + ψn h + lω n h ), ξ h ) +k 0 (ω n x -lφ n -(ω n hx -lφ n h ), ξ hx ) -α(θ n 1 -θ n 1h , ξ hx ) = 0.
As a consequence, we have

(5.3)                                    ρ 1 ( ωn t -1 ∆t (ω n h -ωn-1 h ), ωn -ωn h ) +γ 0 (ω n x -l φn -(ω n hx -l φn h ), ωn x -ωn hx ) +kl(φ n x + ψ n + lω n -(φ n hx + ψ n h + lω n h ), ωn -ωn h ) +γ 1 l( φn x + ψn + lω n -( φn hx + ψn h + lω n h ), ωn -ωn h ) +k 0 (ω n x -lφ n -(ω n hx -lφ n h ), ωn x -ωn hx ) -α(θ n 1 -θ n 1h , ωn x -ωn hx ) = ρ 1 ( ωn t -1 ∆t (ω n h -ωn-1 h ), ωn -ξ h ) +γ 0 (ω n x -l φn -(ω n hx -l φn h ), ωn x -ξ hx ) +kl(φ n x + ψ n + lω n -(φ n hx + ψ n h + lω n h ), ωn -ξ h ) +γ 1 l( φn x + ψn + lω n -( φn hx + ψn h + lω n h ), ωn -ξ h ) +k 0 (ω n x -lφ n -(ω n hx -lφ n h ), ωn x -ξ hx ) -α(θ n 1 -θ n 1h , ωn x -ξ hx )
. Similar ideas applied to the variational equations to the temperatures yield (5.4)

               ( θ n 1t -1 ∆t (θ n 1h -θ n-1 1h ), θ n 1 -θ n 1h ) + (θ n 1x -θ n 1hx , θ n 1x -θ n 1hx ) +α(ω n x -l φn -(ω n hx -l φn h ), θ n 1 -θ n 1h ) = ( θ n 1t -1 ∆t (θ n 1h -θ n-1 1h ), θ n 1 -η 1h ) + (θ n 1x -θ n 1hx , θ n 1x -η 1hx ) +α(ω n x -l φn -(ω n hx -l φn h ), θ n 1 -η 1h
) ,

               ( θ n 2t -1 ∆t (θ n 2h -θ n-1 2h ), θ n 2 -θ n 2h ) + (θ n 2x -θ n 2hx , θ n 2x -θ n 2hx ) +α( ψn x -ψn hx , θ n 2 -θ n 2h ) = ( θ n 2t -1 ∆t (θ n 2h -θ n-1 2h ), θ n 2 -η 2h ) + (θ n 2x -θ n 2hx , θ n 2x -η 2hx ) +α( ψn x -ψn hx , θ n 2 -η 2h and (5.5)
) .

Step 2. To handle the first term in equation (5.1) we observe that

( φn t - 1 ∆t ( φn h -φn-1 h ), φn -φn h ) = ( φn t - 1 ∆t ( φn -φn-1 ), φn -φn h ) + 1 ∆t ( φn -φn-1 -( φn h -φn-1 h ), φn -φn h ) = ( φn t - 1 ∆t ( φn -φn-1 ), φn -φn h ) + 1 2∆t ∥ φn -φn h -( φn-1 -φn-1 h )∥ 2 + 1 2∆t ( ∥ φn -φn h ∥ 2 -∥ φn-1 -φn-1 h )∥ 2
) .

where (4.4) was used. We deal with similar terms in (5.2)-(5.5) in the same way.

Then, using again (4.3) and (4.4) and adding (5.1)-(5.5) we end up with ρ 1 2∆t

( ∥ φn -φn h ∥ 2 -∥ φn-1 -φn-1 h ∥ 2 + ∥ω n -ωn h ∥ 2 -∥ω n-1 -ωn-1 h ∥ 2 ) + ρ 2 2∆t ( ∥ ψn -ψn h ∥ 2 -∥ ψn-1 -ψn-1 h ∥ 2 ) + k 2∆t ( ∥φ n x + ψ n + lω n -(φ n hx + ψ n h + lω n h )∥ 2 -∥φ n-1 x + ψ n-1 + lω n-1 -(φ n-1 hx + ψ n-1 h + lω n-1 h )∥ 2 ) + k 0 ∆t ( ∥ω n x -lφ n -(ω n hx -lφ n h )∥ 2 -∥ω n-1 x -lφ n-1 -(ω n-1 hx -lφ n-1 h )∥ 2 ) + b 2∆t ( ∥ψ n x -ψ n hx ∥ 2 -∥ψ n-1 x -ψ n-1 hx ∥ 2 ) + 1 2∆t ( ∥θ n 1 -θ n 1h ∥ 2 -∥θ n-1 1 -θ n-1 1h ∥ 2 + ∥θ n 2 -θ n 2h ∥ 2 -∥θ n-1 2 -θ n-1 2h ∥ 2 ) +C(∥ ψn x -ψn hx ∥ 2 + ∥ φn x + ψn + lω n -( φn hx + ψn h + lω n h )∥ 2 ) +C(∥ω n x -l φn -(ω n hx -l φn h )∥ 2 + ∥θ n 1x -θ n 1hx ∥ 2 + ∥θ n 2x -θ n 2hx ∥ 2 ) ≤ C ( ∥ φn t - 1 ∆t ( φn -φn-1 )∥ 2 + 1 ∆t ( φn -φn-1 -( φn h -φn-1 h ), φn -ζ h ) ) +C(∥ φn -φn h ∥ 2 + ∥ φn -ζ h ∥ 2 + ∥ φn x -ζ 2 hx ∥ 2 ) +C(∥φ n x + ψ n + lω n -(φ n hx + ψ n h + lω n h )∥ 2 ) + ∥ω n x -lφ n -(ω n hx -lφ n h )∥ 2 +C ( ∥θ n 1 -θ n 1h ∥ 2 + ∥ ψn t - 1 ∆t ( ψn -ψn-1 )∥ 2 + ∥ ψn -ψn h ∥ 2 ) +C ( 1 ∆t ( ψn -ψn-1 -( ψn h -ψn-1 h ), ψn -χ h ) + ∥ ψn -χ h ∥ 2 ) +C(∥ψ n x -ψ n hx ∥ 2 + ∥ ψn x -χ hx ∥ 2 + ∥θ n 2 -θ n 2h ∥ 2 ) +C ( ∥ω n t - 1 ∆t (ω n -ωn-1 )∥ 2 + 1 ∆t (ω n -ωn-1 -(ω n h -ωn-1 h ), ωn -ξ h ) ) +C(∥ω n -ωn h ∥ 2 + ∥ω n -ξ h ∥ 2 + ∥ω n hx -ξ hx ∥ 2 ) +C ( ∥θ n 1t - 1 ∆t (θ n 1 -θ n-1 1 )∥ 2 + 1 ∆t (θ n 1 -θ n-1 1 -(θ n 1h -θ n-1 1h ), θ n 1 -η 1h ) ) +C ( ∥θ n 2t - 1 ∆t (θ n 2 -θ n-1 2 )∥ 2 + 1 ∆t (θ n 2 -θ n-1 2 -(θ n 2h -θ n-1 2h ), θ n 1 -η 2h ) ) +C(∥θ n 1 -η 1h ∥ 2 + ∥θ n 1x -η 1hx ∥ 2 + ∥θ n 2 -η 2h ∥ 2 + ∥θ n 2x -η 2hx ∥ 2 ).
Step 3. Multiplying the latter inequality by ∆t and summing over n we obtain, for all

{ζ i h , χ i h , ξ i h , η i 1h , η i 2h } n i=0 ⊂ S h , ∥ φn -φn h ∥ 2 + ∥ ψn -ψn h ∥ 2 + ∥ω n -ωn h ∥ 2 + ∥ψ n x -ψ n hx ∥ 2 +∥φ n x + ψ n + lω n -(φ n hx + ψ n h + lω n h )∥ 2 +∥ω n x -lφ n -(ω n hx -lφ n h )∥ 2 + ∥θ n 1 -θ n 1h ∥ 2 + ∥θ n 2 -θ n 2h ∥ 2 ≤ C∆t n ∑ i=1 ( ∥ φi -φi h ∥ 2 + ∥ ψi -ψi h ∥ 2 + ∥ω i -ωi h ∥ 2 + ∥ψ i x -ψ i hx ∥ 2 +∥φ i x + ψ i + lω i -(φ i hx + ψ i h + lω i h )∥ 2 +∥ω i x -lφ i -(ω i hx -lφ i h )∥ 2 + ∥θ i 1 -θ i 1h ∥ 2 + ∥θ i 2 -θ i 2h ∥ 2 +∥ φi t - 1 ∆t ( φi -φi-1 )∥ 2 + 1 ∆t ( φi -φi-1 -( φi h -φi-1 h ), φi -ζ i h ) +∥ ψi t - 1 ∆t ( ψi -ψi-1 )∥ 2 + 1 ∆t ( ψi -ψi-1 -( ψi h -ψi-1 h ), ψi -χ i h ) +∥ω i t - 1 ∆t (ω i -ωi-1 )∥ 2 + 1 ∆t (ω i -ωi-1 -(ω i h -ωi-1 h ), ωi -ξ i h ) +∥θ i 1t - 1 ∆t (θ i 1 -θ i-1 1 )∥ 2 + 1 ∆t (θ i 1 -θ i-1 1 -(θ i 1h -θ i-1 1h ), θ i 1 -η i 1h ) +∥θ i 2t - 1 ∆t (θ i 2 -θ i-1 2 )∥ 2 + 1 ∆t (θ i 2 -θ i-1 2 -(θ i 2h -θ i-1 2h ), θ i 1 -η i 2h ) +∥ φi -ζ i h ∥ 2 + ∥ φi x -ζ i hx ∥ 2 + ∥ ψi -χ i h ∥ 2 + ∥ ψi x -χ i hx ∥ 2 +∥ω i -ξ i h ∥ 2 + ∥ω i hx -ξ i hx ∥ 2 +∥θ i 1 -η i 1h ∥ 2 + ∥θ i 1x -η i 1hx ∥ 2 + ∥θ i 2 -η i 2h ∥ 2 + ∥θ i 2x -η i 2hx ∥ 2 ) +C ( ∥φ 1 -φ0 h ∥ 2 + ∥ψ 1 -ψ0 h ∥ 2 + ∥ω 1 -ω0 h ∥ 2 + ∥ψ 0 x -ψ 0 hx ∥ 2 +∥φ 0 x + ψ 0 + lω 0 -(φ 0 hx + ψ 0 h + lω 0 h )∥ 2 +∥ω 0 x -lφ 0 -(ω 0 hx -lφ 0 h )∥ 2 + ∥θ 0 1 -θ 0 1h ∥ 2 + ∥θ 0 2 -θ 0 2h ∥ 2
) .

Taking into account that (with an equivalent result for similar terms)

n ∑ i=1 ( φi -φi-1 -( φi h -φi-1 h ), φi -ζ i h ) = ( φn -φn h , φn -ζ n h ) +( φ0 h -φ0 , φ1 -ζ 1 h ) + n-1 ∑ i=1 ( φi -φi h , φi -ζ i h -( φi+1 -ζ i+1 h )) ≤ C(∥ φn -φn h ∥ 2 + ∥ φn -ζ n h )∥ 2 + ∥ φ0 h -φ0 ∥ 2 + ∥ φ1 -ζ 1 h ∥ 2 ) +C∆t n-1 ∑ i=1 ∥ φi -φi h ∥ 2 + C ∆t n-1 ∑ i=1 ∥ φi -ζ i h -( φi+1 -ζ i+1 h )∥ 2 .
and applying a discrete version of Gronwall's inequality, the result follows. □

The linear convergence of the numerical method is summarized in the following corollary.

Corollary 5.1. Suppose that the solution to the continuous problem is sufficiently regular, that is

φ, ψ, ω ∈ H 3 (0, T ; L 2 (I)) ∩ W 1,∞ (0, T ; H 2 (I)) ∩ H 2 (0, T ; H 1 (I)), θ 1 , θ 2 ∈ H 2 (0, T ; L 2 (I)) ∩ L ∞ (0, T ; H 2 (I)) ∩ H 1 (0, T ; H 1 (I)).
Then, there exists a positive constant C, independent of the discretization parameters h and ∆t, such that

∥ φ(t n ) -φn h ∥ 2 + ∥ ψ(t n ) -ψn h ∥ 2 + ∥ω(t n ) -ωn h ∥ 2 + ∥ψ(t n ) -ψ n h ∥ 2 +∥φ x (t n ) + ψ(t n ) + lω(t n ) -(φ n hx + ψ n h + lω n h )∥ 2 +∥ω x (t n ) -lφ(t n ) -(ω n hx -lφ n h )∥ 2 + ∥θ 1 (t n ) -θ n 1h ∥ 2 + ∥θ 2 (t n ) -θ n 2h ∥ 2 ≤ C(h 2 + (∆t) 2 ).
Proof. The result is a consequence of estimates like

1 ∆t n-1 ∑ i=1 ∥ φi -ζ i h -( φi+1 -ζ i+1 h )∥ 2 ≤ Ch 2 ∥ φt ∥ 2 L 2 (0,T ;H 1 (I)) ,
as proved in the work of Han, Shillor and Sofonea [START_REF] Han | Variational and numerical analysis of a quasi-static viscoelastic problem with normal compliance, friction and damage[END_REF], and classical results on the approximation properties of the finite element spaces [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]. □

Simulations

We describe in this section the results of some numerical experiments performed to demonstrate the accuracy of the approximation and the behavior of the solution.

To find the solution of system (4.2) an iterative procedure was used. Assuming that φn-1

h , ψn-1 h , ωn-1 h , θ n-1 1h , θ n-1 2h
are known and setting

φ n,0 h = φ n-1 h , φn,0 h = φn-1 h , ψ n,0 h = ψ n-1 h , ψn,0 h = ψn-1 h , ω n,0 h = ω n-1 h , ωn,0 h = ωn-1 h , θ n,0 1h = θ n-1 1h , θ n,0 2h = θ n-1 2h
. we solved the problem (6.1)

                                                       1 ∆t (θ n,j 1h -θ n-1 1h , η 1h ) + (θ n,j 1hx , η 1hx ) + α(ω n,j-1 hx -l φn,j-1 h , η 1h ) = 0, 1 ∆t (θ n,j 2h -θ n-1 2h , η 2h ) + (θ n,j 2hx , η 2hx ) + α( ψn,j-1 hx , η 2h ) = 0, ρ2 ∆t ( ψn,j h -ψn-1 h , χ h ) + b(ψ n,j hx , χ hx ) + γ 2 ( ψn,j hx , χ hx ) +k(φ n,j-1 hx + ψ n,j h + lω n,j-1 h , χ h ) +γ 1 ( φn,j-1 hx + ψn,j h + lω n,j-1 h , χ h ) -α(θ n,j 2h , χ hx ) = 0, ρ1 ∆t ( φn,j h -φn-1 h , ζ h ) + k(φ n,j hx + ψ n,j h + lω n,j-1 h , ζ hx ) +γ 1 ( φn,j hx + ψn,j h + lω n,j-1 h , ζ hx ) -k 0 l(ω n,j-1 hx -lφ n,j h , ζ h ) -γ 0 l(ω n,j-1 hx -l φn,j h , ζ h ) + lα(θ n,j 1h , ζ h ) = 0, ρ1 ∆t (ω n,j h -ωn-1 h , ξ h ) + k 0 (ω n,j hx -lφ n,j h , ξ hx ) + γ 0 (ω n,j hx -l φn,j h , ξ hx ) +kl(φ n,j hx + ψ n,j h + lω n,j h , ξ h ) + γ 1 l( φn,j hx + ψn,j h + lω n,j h , ξ h ) -α(θ n,j 1h , ξ hx ) = 0, where, for j = 1, 2 • • • , φ n,j h = φ n-1 h + ∆t φn,j h , ψ n,j h = ψ n-1 h + ∆t ψn,j h , ω n,j h = ω n-1 h + ∆tω n,j h .
Problem (6.1) consists of five, uncoupled, linear systems of algebraic equations, with tri-diagonal matrices, that have unique solution.

A curved beam with radius of curvature R = 1 and length L = 1 was considered with ρ 1 = 1, ρ 2 = 2, k = 1, k 0 = 0.5, and b = 1. The viscosity parameters are γ 0 = γ 1 = γ 2 = 0.1 and the coupling constant is α = 0.017. The iterative procedure was stopped when the difference between two successive iterations became smaller than a given tolerance of T OL = 10 -7 .

In the first experiment, the time evolution of the system and the energy decay were investigated. The discretization parameters are h = 0.01 and ∆t = 10 -4 and the initial data φ 0 (x) = 0.01x(x -1), φ 1 (x) = 20x(x -1) 2 , ψ 0 (x) = ψ 1 (x) = ω 0 (x) = 0, ω 1 (x) = x 3 -x 2 , θ 10 (x) = 20x cos(0.5πx), θ 20 (x) = 10 sin(πx). As expected, the system evolved towards the zero steady-state with the energy decaying to zero very quickly. In our simulations, the temperatures θ 1 and θ 2 vanished faster than φ, ψ and ω. See figures 1-3.

Next, we performed a simulation to test, numerically, the error estimate. We solved the modified problem The computed errors at time t = 1.2 are shown in table 1 where Error is defined by

                                       ρ 1 φ tt -k(φ x + ψ + lw) x -γ 1 (φ x + ψ + lw) xt -k 0 l(ω x -lφ) -γ 0 l(w x -lφ) t + lαθ 1 = f 1 , ρ 2 ψ tt -
Error = ( ∥ φ(t n ) -φn h ∥ 2 + ∥ ψ(t n ) -ψn h ∥ 2 + ∥ω(t n ) -ωn h ∥ 2 + ∥ψ(t n ) -ψ n h ∥ 2 +∥φ x (t n ) + ψ(t n ) + lω(t n ) -(φ n hx + ψ n h + lω n h )∥ 2 +∥ω x (t n ) -lφ(t n ) -(ω n hx -lφ n h )∥ 2 + ∥θ 1 (t n ) -θ n 1h ∥ 2 + ∥θ 2 (t n ) -θ n 2h ∥ 2 ) 1 2 .
We note that the error decreased by a factor of 2. The linear convergence rate can also be observed in the curves in figure 4.

ϕ(x,t)

Acknowledgements

The work of M.I.M. Copetti was partially supported by the Brazilian institution CNPq (grant 304709/2017-4).

MFigure :

 : Figure: After deformation the particle M 0 of the beam take the position M

  ), (3.18), (3.19), (3.20), (3.22), (3.23), (3.24), (3.25), and (3.26) in (3.15) yields 2 ∫ T S E(t)dt ≤ cE(S) + ε ∫ T S E(t)dt ∀S ∈ (0, T )

Figure 1 .

 1 Figure 1. The evolution in time of φ, ψ, and ω.

  0. which proves the result.

□

Remark 4.1. Note that to find the solution of (4.2) a square linear system of algebraic equations needs to be solved. It follows from the above proof that when all data are zero, the solution { φn h , ψn h , ωn h , θ n 1h , θ n 2h } is zero. Then, (4.2) admits a unique solution.

  bψ xx -γ 2 ψ xxt + k(φ x + ψ + lω) +γ 1 (φ x + lω + ψ) t + αθ 2x = f 2 , ρ 1 ω tt -k 0 (ω x -lφ) x -γ 0 (ω x -lφ) xt + kl(φ x + ψ + lω) +γ 1 l(φ x + ψ + lω) t + αθ 1x = f 3 , ρ 1 θ 1t -θ 1xx + α(w x -lφ) t = f 4 , ρ 2 θ 2t -θ 2xx + αψ xt = f 5 ,where f 1 , f 2 , f 3 , f 4 , f 5 , and the initial data are calculated from the exact solution

φ(x, t) = -0.1(x 2 -x)t 2 / √ 2, ψ(x, t) = 0.5t(x 2 -x), ω(x, t) = t(x 3 -x 2 ), θ 1 (x, t) = t 2 sin(πx), θ 2 (x, t) = e t (x 2 -x).

Table 1 .

 1 Figure 2. The evolution in time of θ 1 , θ 2 , and the energy.Figure 3. The evolution in time of φ, ψ, and ω at point x = 0.6. Computed errors when T = 1.2.

		8
		6				t=0
	1 (x,t)	4
	θ
		2
					t=0.2
		0		t=0.6	t=1.0
		-2
		0	.2	.4	.6	.8		1
					t
		10
						t=0
	2 (x,t)	5
	θ
					t=0.2
		0		t=0.6	t=1.0
		0	.2	.4	.6	.8		1
					t
		40
		30
	ε(t)	20
		10
		0
		0	.1	.2	.3	.4	.5	.6
					t


Ribas Dayana

Vivolab

Aragon

Emmanuel Vincent

AN IMPROVED UNCERTAINTY PROPAGATION METHOD FOR ROBUST I-VECTOR BASED SPEAKER RECOGNITION

Keywords: Uncertainty propagation, speaker verification, data distortion, robustness, i-vector

The performance of automatic speaker recognition systems degrades when facing distorted speech data containing additive noise and/or reverberation. Statistical uncertainty propagation has been introduced as a promising paradigm to address this challenge. So far, different uncertainty propagation methods have been proposed to compensate noise and reverberation in i-vectors in the context of speaker recognition. They have achieved promising results on small datasets such as YOHO and Wall Street Journal, but little or no improvement on the larger, highly variable NIST Speaker Recognition Evaluation (SRE) corpus. In this paper, we propose a complete uncertainty propagation method, whereby we model the effect of uncertainty both in the computation of unbiased Baum-Welch statistics and in the derivation of the posterior expectation of the i-vector. We conduct experiments on the NIST-SRE corpus mixed with real domestic noise and reverberation from the CHiME-2 corpus and preprocessed by multichannel speech enhancement. The proposed method improves the equal error rate (EER) by 4% relative compared to a conventional i-vector based speaker verification baseline. This is to be compared with previous methods which degrade performance.

INTRODUCTION

Uncertainty propagation has emerged as a paradigm for robust signal processing whereby the data are not treated as point estimates anymore, but as a parametric posterior distribution, typically approximated as a Gaussian. This approach provides a principled framework to deal with the loss of information due to signal distortion -epistemic uncertainty -or to the finite number of data points -aleatoric uncertainty -. The uncertainty is represented by a set of scalar variances or covariance matrices, which are first estimated on the data, and then propagated through the subsequent processing steps in order to compensate for the effect of uncertainty on the computed quantities and ultimately improve the system performance [START_REF] Arrowood | Using observation uncertainty in HMM decoding[END_REF][START_REF] Becerra | Speaker verification in noise using a stochastic version of the weighted Viterbi algorithm[END_REF][START_REF] Deng | Dynamic compensation of HMM variances using the feature enhancement uncertainty computed from a parametric model of speech distortion[END_REF][START_REF] Liao | Joint uncertainty decoding for noise robust speech recognition[END_REF].

In speaker recognition, the uncertainty propagation approach has gained traction motivated by the uncertain nature of the system pipeline when the application scenario tends to more real-world situations. The necessity of improving the system robustness in noisy environments has inspired the development of epistemic uncertainty approaches for speaker modeling based on Gaussian mixture models [START_REF] Zhao | Variational Bayesian joint factor analysis for speaker verification[END_REF][START_REF] Ozerov | Uncertainty-based learning of acoustic models from noisy data[END_REF] and, more recently, i-vectors [START_REF] Yu | Uncertainty propagation in front end factor analysis for noise robust speaker recognition[END_REF][START_REF] Ribas | Uncertainty propagation for noise robust speaker recognition: the case of NIST-SRE[END_REF]. Other work based on the aleatoric uncertainty concept has also given rise to uncertainty propagation approaches for speaker recognition. However, these approaches focused on the issue of computing representations with in-sufficient data, caused by utterances with different, possibly short durations [START_REF] Hautamäki | Minimax i-vector extractor for short duration speaker verification[END_REF][START_REF] Stafylakis | Text dependent speaker recognition using PLDA with uncertainty propagation[END_REF][START_REF] Kenny | PLDA for speaker verification with utterances of arbitrary duration[END_REF][START_REF] Cumani | On the use of ivector posterior distributions in Probabilistic Linear Discriminant Analysis[END_REF][START_REF] Saeidi | Accounting for uncertainty of i-vectors in speaker recognition using uncertainty propagation and modified imputation[END_REF][START_REF] Kenny | Uncertainty modeling without subspace methods for text-dependent speaker recognition[END_REF][START_REF] Lin | Fast scoring for plda with uncertainty propagation via i-vector grouping[END_REF]. This paper pursues the same line as our preliminary study [START_REF] Ribas | Uncertainty propagation for noise robust speaker recognition: the case of NIST-SRE[END_REF], that considered an epistemic uncertainty propagation approach for noise-robust text-independent speaker verification using a system based on i-vectors [START_REF] Dehak | Front-end factor analysis for speaker verification[END_REF] and probabilistic linear discriminant analysis (PLDA) [START_REF] Garcia | Analysis of i-vector length normalization in speaker recognition systems[END_REF]. Despite the recent introduction of deep learning based modules in the speaker recognition pipeline, reports of the last NIST Speaker Recognition Evaluation (SRE) in 2016 [START_REF] Lee | The I4U submission to the 2016 NIST speaker recognition evaluation[END_REF] and the experience in the recent campaign NIST-SRE 20181 evidenced that the i-vector-PLDA approach still performs among the best systems of the state-of-the-art. It is also the focus of current challenges in the field such as the 2018 Multi-target Speaker Detection and Identification Challenge Evaluation2 .

In [START_REF] Ribas | Uncertainty propagation for noise robust speaker recognition: the case of NIST-SRE[END_REF], we proposed a method to estimate and propagate the residual uncertainty after multichannel speech enhancement from the acoustic features to the i-vectors. Specifically, we modified the posterior probability of each Gaussian mixture component to obtain unbiased Baum-Welch (BW) statistics. Preliminary experiments yielded good results on the Wall Street Journal (WSJ) corpus, but little or no improvement on a subset of the NIST-SRE 2008 corpus, similarly to the findings of Yu et al. [START_REF] Yu | Uncertainty propagation in front end factor analysis for noise robust speaker recognition[END_REF] on the YOHO and NIST-SRE 2010 datasets. We studied the causes of this under-performance on NIST-SRE and found out that the high variability of NIST-SRE makes it intrinsically harder to obtain accurate BW statistics.

In this paper, we propose a new uncertainty propagation method that models the effect of epistemic uncertainty both in the computation of the BW statistics and in the derivation of the i-vector. This method provides a more complete strategy towards compensating for background noise in all steps of the i-vector computation process. Furthermore, this study contributes to clarifying how the i-vector speaker representations are affected by the environmental distortion. The results are evaluated on a subset of the NIST-SRE 2008 corpus mixed with real domestic noise and reverberation from the CHiME-2 corpus [START_REF] Vincent | The second 'CHIME' speech separation and recognition challenge: Datasets, tasks and baselines[END_REF].

Section 2 recalls the i-vector computation process. Section 3 provides a novel analysis of the limitations of related previous works and Section 4 introduces the proposed method. Speaker verification experiments are conducted in Section 5. The results are reported and discussed in Section 6. Finally, the conclusions of the study and future work are presented in Section 7.

I-VECTOR COMPUTATION

Front-end factor analysis [START_REF] Dehak | Front-end factor analysis for speaker verification[END_REF] relies on a universal background model (UBM) that is a mixture of C Gaussian components indexed by c. Denoting by F the feature dimension, the CF × 1 supervector M (u) for one utterance u is expressed as

M (u) = m + T w(u) + (u) (1)
where m consists of the means mc of all UBM components, T is the CF × D low-rank total variability matrix, w(u) is the D × 1 vector of total factors or i-vector, and (u) represents the residual data variability not captured by T . The i-vector is modeled as a zero-mean standard Gaussian random vector. It is obtained by computing the posterior expectation of w(u) over the feature sequence {y1, . . . , yL}, with L the number of time frames:

E[w(u)] = (I + T V -1 N (u)T ) -1 T V -1 F (u). (2)
In this equation, N (u) is a CF × CF diagonal matrix with diagonal blocks Nc(u)I where Nc(u) are the zeroth-order BW statistics for all components c, F (u) is a CF × 1 supervector obtained by concatenating the centralized first-order BW statistics Fc(u), V is the diagonal CF × CF covariance matrix of (u), and denotes matrix transposition. The BW statistics are given by

Nc(u) = L t=1 γt(c) (3)
Fc(u) = L t=1 γt(c)(yt -mc) (4)
where

γt(c) = πc N (yt|µc, Σc) C i=1 πi N (yt|µi, Σi) (5)
is the posterior probability of the c-th UBM component, as obtained from its mean mc, covariance Σc and weight πc.

Note that (2) can be equivalently rewritten in terms of the "normalized" statistics Ñ (u) and F (u) as

E[w(u)] = (I + T Ñ (u)T ) -1 T F (u). (6)
where

Ñ (u) = V -1 N (u) is a diagonal matrix with diagonal blocks Ñc(u) = Nc(u)V -1 c , F (u) = V -1 F (u) is obtained by concate- nating Fc(u) = V -1 c
Fc(u), and Vc is the c-th diagonal block of V . The multiplication by V -1 c can be distributed at each time t of the summations in ( 3) and ( 4). The purpose of this rewriting will become clear in the next section.

UNCERTAINTY PROPAGATION TO THE I-VECTOR

Let us assume that we now observe a corrupted speech signal xt involving noise and/or reverberation. Using a speech enhancement algorithm together with an uncertainty estimation technique, the posterior probability of the clean speech features yt can be modeled as [START_REF] Kolossa | Robust speech recognition of uncertain or missing data[END_REF] p(yt|xt) = N (yt|ȳt, Σt)

with ȳt the enhanced features and Σt the uncertainty covariance matrix at time t. In other words, Σt is the covariance of the estimation error between the enhanced features and the (unknown) clean features at a given time.

To the best of our knowledge, the studies in [START_REF] Yu | Uncertainty propagation in front end factor analysis for noise robust speaker recognition[END_REF][START_REF] Ribas | Uncertainty propagation for noise robust speaker recognition: the case of NIST-SRE[END_REF] are the only ones exploiting this model for noise and reverberation robustness in i-vector based speaker recognition systems, while other works focused on earlier, now deprecated systems. They proposed two different ways to propagate the uncertainty from the enhanced features to the i-vectors.

Uncertainty propagation through the front-end factor analysis model

The authors in [START_REF] Yu | Uncertainty propagation in front end factor analysis for noise robust speaker recognition[END_REF] considered the generative data model corresponding to [START_REF] Arrowood | Using observation uncertainty in HMM decoding[END_REF]. By integrating over the unknown clean features, they accounted for the impact of uncertainty on the expression of the joint posterior probability. They derived the posterior expectation of w(u) over the feature sequence in a similar way to (6) as3

E[wunc(u)] = (I + T Ñunc(u)T ) -1 T Func(u). (8)
Ñunc(u) becomes the CF × CF block-diagonal matrix with diagonal blocks

Ñunc,c(u) = L t=1 γt(c)V -1 unc,c,t , (9)
Func(u) is the CF × 1 supervector obtained by concatenating

Func,c(u) = L t=1 γt(c)V -1 unc,c,t (ȳt -mc), (10)
and Vunc,c,t is the total covariance of the residual variability and the uncertainty:

Vunc,c,t = Vc + Σt. (11)
In order to compute the posterior probability of the c-th UBM component, the authors substituted the clean features yt in (5) by the enhanced features ȳt:

γt(c) = πc N (ȳt|µc, Σc) C i=1 πi N (ȳt|µi, Σi) (12)
This expression does not account for the difference between the enhanced features and the clean features, hence the BW statistics are biased. Indeed, it is well known from the field of speech recognition that using enhanced data as inputs to an acoustic model trained on clean data often results in poor recognition performance due to the residual distortions in the enhanced data [START_REF] Deng | Dynamic compensation of HMM variances using the feature enhancement uncertainty computed from a parametric model of speech distortion[END_REF][START_REF] Kolossa | Robust speech recognition of uncertain or missing data[END_REF]. We conclude that the uncertainty is not fully propagated to the i-vector domain, and therefore the obtained i-vectors remain affected by the data distortion. This may be the reason for the poor results achieved by this method on the NIST-SRE 2010 corpus, even in the situation when oracle (ideal) uncertainty estimates are used [START_REF] Yu | Uncertainty propagation in front end factor analysis for noise robust speaker recognition[END_REF].

Uncertainty propagation through the UBM

The study in [START_REF] Ribas | Uncertainty propagation for noise robust speaker recognition: the case of NIST-SRE[END_REF] presented an algorithm to compute unbiased BW statistics instead. By considering the generative data model associated with the UBM and integrating over the unknown clean features, the likelihood p(yt|c) = N (yt|µc, Σc) of the c-th UBM component is classically substituted by [START_REF] Deng | Dynamic compensation of HMM variances using the feature enhancement uncertainty computed from a parametric model of speech distortion[END_REF] p(xt|c) = N (ȳt|µc, Σunc,c,t)

where Σunc,c,t = Σc + Σt.

The authors considered the i-vector computation in (2) as a deterministic operation and consequently derived the i-vector as [START_REF] Lin | Fast scoring for plda with uncertainty propagation via i-vector grouping[END_REF] where Nunc(u) is a diagonal matrix with blocks Nunc,c(u)I, Func(u) is the concatenation of Func,c(u), and the unbiased BW statistics [START_REF] Ozerov | Uncertainty-based learning of acoustic models from noisy data[END_REF] are given by γunc,t(c) = πc N (ȳt|µc, Σunc,c,t)

E[wunc(u)] = (I + T V -1 Nunc(u)T ) -1 T V -1 Func(u)
C i=1 πi N (ȳt|µi, Σunc,i,t) (16
)

Nunc,c(u) = L t=1 γunc,t(c) (17)
Func,c(u) = L t=1 γunc,t(c)Wc,t(ȳt -mc) (18)
with Wc,t the Wiener filter [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF] defined as

Wc,t = ΣcΣ -1 unc,c,t . (19)
Despite the use of unbiased BW statistics, the resulting i-vectors are not fully compensated for the uncertainty at the input. Since the ivector computation is not according to the modifications introduced to the speaker model due to the propagation of the uncertainty, it raises a mismatch that negatively impacts the system performance. Indeed, i-vector computation is not a deterministic operation and it must be modified to account for this uncertainty. This may partly explain the poor results reported on the NIST-SRE 2008 corpus in [START_REF] Ribas | Uncertainty propagation for noise robust speaker recognition: the case of NIST-SRE[END_REF], in addition to the high variability of this corpus which makes it harder to obtain meaningful BW statistics (see study in [START_REF] Ribas | Uncertainty propagation for noise robust speaker recognition: the case of NIST-SRE[END_REF]).

PROPOSED METHOD

The analysis of previous works on uncertainty propagation for ivector based speaker recognition has uncovered the fact that the uncertainty is propagated either through the front-end factor analysis model or through the UBM. In order to fully compensate i-vectors for the distortion of the input data, we propose to propagate the uncertainty through both models. I-vectors are then computed in the same way as in [START_REF] Ribas | Uncertainty propagation for noise robust speaker recognition: the case of NIST-SRE[END_REF], albeit with unbiased normalized statistics. Specifically, Ñunc(u) is now the block-diagonal matrix with diagonal blocks

Ñunc,c(u) = L t=1 γunc,t(c)V -1 unc,c,t (20)
and Func(u) is obtained by concatenating

Func,c(u) = L t=1 γunc,t(c)V -1 unc,c,t Wc,t(ȳt -mc), (21)
where γunc,t(c), Vunc,c,t, and Wc,t are defined in ( 16), [START_REF] Kenny | PLDA for speaker verification with utterances of arbitrary duration[END_REF], and [START_REF] Vincent | The second 'CHIME' speech separation and recognition challenge: Datasets, tasks and baselines[END_REF], respectively.

This time, the enhanced data ȳt and the uncertainty Σt are both used to compute the unbiased BW statistics. Furthermore, the i-vector obtained from the computation of the posterior expectation of w(u) also considers the influence of Σt.

In order to assess the impact of unbiased BW statistics independently of the i-vector computation process, we compared them to biased BW statistics as follows. We chose the unnormalized F statistics as a rough representation of the speaker characteristics before i-vector computation. We computed the cosine distance between the F vectors of each trial of the dataset used in the experimental setup (see section 5.1). The results displayed in Fig. 1 show higher cosine distances for the unbiased statistics in non-target trials. The average cosine distance over all non-target trials is 0.7939 for the unbiased F statistics vs. 0.7422 for the biased F statistics. This suggests that non-target utterances can be better separated using unbiased statistics. The poor results reported on the NIST-SRE 2008 corpus in [START_REF] Ribas | Uncertainty propagation for noise robust speaker recognition: the case of NIST-SRE[END_REF] using i-vector based cosine distances suggest that the benefit of unbiased BW statistics observed when considering the statistics themselves was lost when considering i-vectors instead. Overall, this suggests that unbiased BW statistics and uncertainty-aware i-vector computation are both necessary to achieve good results.

EXPERIMENTAL SETUP

Dataset description

We compare the proposed method with [START_REF] Yu | Uncertainty propagation in front end factor analysis for noise robust speaker recognition[END_REF] and [START_REF] Ribas | Uncertainty propagation for noise robust speaker recognition: the case of NIST-SRE[END_REF] using the same experimental setup as in [START_REF] Ribas | Uncertainty propagation for noise robust speaker recognition: the case of NIST-SRE[END_REF]. The speech signals are male conversations in English from the NIST-SRE 2004 and 2005 corpora for the training stage (3285 speech signals from 262 speakers), and from the NIST-SRE 2008 short2 and short3 corpora for the evaluation stage (enrollment: 470 utterances, test: 671 utterances). A total of 6615 verifications were performed on the det7 condition of NIST-SRE. We used two-channel noise and reverberation from Track 1 of the 2nd CHiME Challenge [START_REF] Vincent | The second 'CHIME' speech separation and recognition challenge: Datasets, tasks and baselines[END_REF], recorded in a real domestic environment. Each training or evaluation signal was convolved with one of 121 two-microphone room impulse responses with a reverberation time of 0.3 seconds. Moreover, the enrollment and test signals were mixed with a random segment of real background noise including, e.g., voices, TV, game console, cutlery sounds, and footsteps. This mixing process resulted in noisy speech signals with different SNRs ranging from about -10 to +20 dB with an average of 6.1 dB.

Speech enhancement and uncertainty estimation

For speech enhancement, we applied the CHiME-2 recipe [START_REF] Salaün | The flexible audio source separation toolbox version 2.0[END_REF] of the Flexible Audio Source Separation Toolbox (FASST) [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF]. This recipe was designed to reduce noise, but not reverberation. Hence we consider yt to be the reverberated noiseless data. For uncertainty estimation, in order to assess the potential of each method, we used the oracle diagonal uncertainty covariance matrices Σt = diag(yt -ȳt)

2 [START_REF] Kolossa | Robust speech recognition of uncertain or missing data[END_REF]. However, for a real application there are several uncertainty estimators [START_REF] Dung T Tran | Nonparametric uncertainty estimation and propagation for noise robust ASR[END_REF].

Speaker recognition system

The features consist of 19 Mel frequency cepstral coefficients (MFCCs), the log-energy, and their first-and second-order derivatives, followed by voice activity detection (VAD) and cepstral mean and variance normalization (CMVN) [START_REF] Benesty | Springer Handbook of Speech Processing[END_REF]. The UBM (C = and the T matrix (D = 400) were trained on either clean or reverberated noiseless data [START_REF] Ribas | Full multicondition training for robust i-vector based speaker recognition[END_REF]. Since the chosen speech enhancement method reduces noise but not reverberation, training on reverberated speech provides better results than training on clean speech. The i-vectors are centered, whitened and length-normalized, and subsequently projected with 330-dimensional LDA. Classification relies on Gaussian PLDA [START_REF] Garcia | Analysis of i-vector length normalization in speaker recognition systems[END_REF].

RESULTS AND DISCUSSION

The speaker verification results expressed in terms of equal error rate (EER) are presented in Table 1. For comparison purposes, the first two rows show the results obtained when training and testing on clean (original NIST) signals or when training and testing on reverberated noiseless signals. The following rows show the results obtained when training on reverberated noiseless signals and testing on noisy or enhanced signals. The last three rows correspond to the two previous uncertainty propagation methods reviewed in Section 3 and to the proposed one in Section 4. Yu et al.'s method [START_REF] Yu | Uncertainty propagation in front end factor analysis for noise robust speaker recognition[END_REF] was tested based on an implementation from the authors. The best achievable EER on reverberated data (without performing dereverberation) is 4.33%. Background noise strongly degrades it up to 31.85%. Multichannel speech enhancement improves it down to 10.48%. Previous uncertainty propagation methods do not outperform this baseline; actually, the method in [START_REF] Yu | Uncertainty propagation in front end factor analysis for noise robust speaker recognition[END_REF] significantly degrades it. Only the proposed method manages to improve performance down to an EER of 10.02%. This corresponds to a 4% relative improvement with respect to the speech enhancement baseline without uncertainty propagation. Although this might seem limited, it is to the best of our knowledge the first time an improvement in terms of noise robustness is reported for i-vector based systems by means of uncertainty propagation on a large NIST-SRE corpus.

Figure 2 shows the histograms of scores for the three uncertainty propagation algorithms in Table 1. The scores are normalized using the maximum score value across the three sets. The scores for nontarget trials obtained with the proposed method appear to be smallest among the three tested systems. Indeed the corresponding curve is the leftmost one in the zoomed view. This behavior is consistent with our findings regarding unbiased vs. biased statistics in Section 4. This indicates there is a wider separation between non-target and target trials in the system using the proposed uncertainty propagation algorithm, hence the performance improvement obtained in Table 1.

CONCLUSIONS

In this paper, we proposed a new method for the propagation of uncertainty due to speech distortion from the input data to the i-vectors in a speaker recognition framework. We analyzed the limitations of previously proposed methods, and introduced a more complete propagation method that operates both on the UBM and on the front-end factor analysis model. We evaluated the speaker verification performance on the NIST-SRE 2008 corpus mixed with real domestic noise and reverberation and obtained a 4% relative EER improvement compared to a speech enhancement baseline without uncertainty propagation, while previous methods degraded performance. Future work will target the extension of the uncertainty propagation approach to the full speaker verification pipeline including the PLDA and the joint handling of epistemic uncertainty due to signal distortion and aleatoric uncertainty due to finite signal size. Also we plan to extend the experimental setup for including female speakers.

ACKNOWLEDGMENTS

Authors would like to thank Dr. Chengzhu Yu for sharing the implementation of the UP referenced approach in [START_REF] Yu | Uncertainty propagation in front end factor analysis for noise robust speaker recognition[END_REF]. We are also grateful for the useful discussions with Dr. Antonio Miguel and Dr. Alfonso Ortega from ViVoLab.

Fig. 1 .

 1 Fig. 1. Cosine distance for the unbiased and biased F statistics in non-target trials.

Fig. 2 .

 2 Fig. 2. Histogram of scores for non-target (left) and target trials (right). The top right corner is a zoomed view of the dashed area.

Table 1 .

 1 Speaker recognition results on NIST-SRE.

	Training set	Evaluation set	EER (%)
	Clean	Clean	3.19
	Reverberated	Reverberated	4.33
	Reverberated	Noisy	31.85
	Reverberated	Enhanced	10.48
	Reverberated	Uncertainty propagation [7]	11.39
	Reverberated	Uncertainty propagation [8]	10.69
	Reverberated Proposed uncertainty propagation	10.02

https://www.nist.gov/itl/iad/mig/nist-2018-speaker-recognitionevaluation

2 http://mce.csail.mit.edu/

For the sake of clarity, we modified the notations in[START_REF] Yu | Uncertainty propagation in front end factor analysis for noise robust speaker recognition[END_REF] for consistency with the above.