Juliette Chabassier

Sébastien Imperiale

Stability and dispersion analysis of improved time discretization for simply supported prestressed Timoshenko systems. Application to the stiff piano string

Keywords: Prestressed Timoshenko system, Theta schemes, Implicit time discretization, Dispersion analysis, Stability analysis, Energy techniques système de Timoshenko précontraint, theta schémas, discrétisation implicite en temps, analyse de dispersion, analyse de stabilité, techniques d'énergie

We study the implicit time discretization of piano strings governing equations within the Timoshenko prestressed beam model. Such model features two different waves, namely the flexural and shear waves, that propagate with very different velocities. We present a novel implicit time discretization that reduces the numerical dispersion while allowing the use of a large time step in the numerical computations. After analyzing the continuous system and the two branches of eigenfrequencies associated with the propagating modes, the classical θ-scheme is studied. We present complete new proofs of stability using energy-based approaches that provide uniform results with respect to the featured time step. A dispersion analysis confirms that θ = 1/12 reduces the numerical dispersion, but yields a severely constrained stability condition for the application considered. Therefore we propose a new θ-like scheme, which allows to reduce the numerical dispersion while relaxing this stability condition. Stability proofs are also provided for this new scheme. Theoretical results are illustrated with numerical experiments corresponding to the simulation of a realistic piano string.

Introduction

The piano strings can be modeled as linear vibrating rods governed by d'Alembert's equation. This model predicts that the eigenfrequencies of vibrating string are "harmonic", i.e. multiple of a fundamental frequency the so called musical pitch. However, measured spectra appear to be slightly shifted away from the ideal harmonic one. This phenomenon is referred to as inharmonicity and in fact contributes to the tone of the instrument (see [START_REF] Conklin | Design and tone in the mechanoacoustic piano. part III. Piano strings and scale design[END_REF][START_REF] Fletcher | Quality of piano tones[END_REF]). It appears that the diameters of realistic piano string are actually too large to neglect shear deformation and rotational inertia, hence it has been suggested in [START_REF] Chabassier | Modélisation et simulation numérique d'un piano par modèles physiques[END_REF] to model the piano string as a simply supported prestressed Timoshenko beam. This models describes the propagation of transversal displacements and shear angles along the string which are phenomena occurring with high velocities contrasts.

More generally, many physical models consist of coupled transient hyperbolic systems with different types of waves propagating at different velocities (for instance, elastodynamic propagation of shear and pressure waves in soft media or acoustic and elastic waves in poro-elastic media). When considering numerical schemes, a good accuracy is obtained when the discretization parameters are adapted to the time and space scales of these physical phenomena. When different scales are involved in the problem, a standard discretization may not be adapted to the entire range of dynamical phenomena. In [START_REF] Burel | Solving the Homogeneous Isotropic Linear Elastodynamics Equations Using Potentials and Finite Elements. The Case of the Rigid Boundary Condition[END_REF], which concerns the case of homogeneous isotropic elastodynamics, the authors take advantage of the fact that the elastic S and P waves can be decoupled up to the boundary of the considered domain, by adapting the spatial discretization, within the domain to each wave speed, while the time step is the same for both waves. When such decomposition into elementary modes is not possible, then alternative approaches must be considered. Although local time step procedures are an other interesting point of view, we choose to focus, in this article, on a specific class of implicit time discretizations. This approach is appealing since it allows to choose a large time step hence reducing computational costs. However, one of the main difficulties is then to limit numerical dispersion, which is inherent to most numerical schemes. This dispersion can significantly degrade the results, which is especially detrimental in the context of musical acoustics as the human ear can detect frequency variations of only a few percents.

In the global piano model, the string is only an elementary part of the problem which also includes the soundboard, modeled as a 2D Reissner Mindlin plate, and the sound radiation in the 3D domain. When considering a standard piano note, taking into account frequencies up to 10 kHz are necessary in order to represent the physics accurately. Hence, the spatial discretization of the 3D domain must represent wave lengths down to a few centimeters, and given the size of the problem, an explicit strategy must be adopted for the time discretization of the sound propagation problem. Unfortunately, the maximal time step allowed to ensure the stability of the 3D problem is too large to account for the physical phenomena occurring in the string (which would require a relatively to small time step if using a classical leap frog scheme). However it seems unreasonable to use such a small time step, since for the piano model, the resolution of the 3D problem is by far the most costly operation. Therefore, as mentioned before, the use of implicit methods to discretize the piano string equations is not only appealing but also necessary to allow the use of large time step for the global coupled problem.

After being discretized using continuous finite elements, the considered problem reduces to a set of hyperbolic coupled ODEs (see (23) for the exact definition). Newmark scheme may then be chosen to deal with the time discretization of this ODE system (see chapter XX of [START_REF] Dautray | Mathematical analysis and numerical methods for science and technology[END_REF]). This scheme, that is widely employed in mechanical applications, depends on two positive real num-Inria bers θ and δ. For δ < 1/2, the Newmark scheme is dissipative and first order accurate. When δ = 1/2, it is at least second-order accurate and preserves the energy of the discretized problem. In this case, it only depends on the parameter θ, it is therefore referred to as the θ-scheme. For θ = 0, this scheme, whose benefits are a less restrictive stability condition (if any) and an improved accuracy, leads to the inversion of a linear system at each time step, which is sometimes considered to be too expensive compared to explicit schemes (i.e. the case θ = 0). In this study and as discussed before, we are willing to use implicit schemes if it allows to obtain good accuracy with a large time step. Therefore the θ-scheme will be our starting point to derive an improved time discretization. This paper is divided into 3 main sections. First, in section 2, we study the continuous system of PDEs of the simply supported prestressed Timoshenko beam and give an energy identity that leads to a-priori estimates on the solution's norm. Then the eigenfrequencies are derived. These steps will be reproduced at a discrete level in sections 3 and 4 to construct and study numerical schemes : the numerical stability is analyzed using energy estimates and the consistency through the calculation of the discrete eigenfrequencies.

In section 3 we discuss the classical properties of the θ-schemes: they preserve an energy-like discrete quantity, which is positive for any ∆t when θ ≥ 1/4, or if ∆t is smaller than a maximal value when θ < 1/4. We provide a stability proof using standard arguments based on the proof in the continuous case for θ ≥ 1/4, and provide a new proof for θ < 1/4, which remains valid for ∆t at the stability limit, as opposed to the intuitive extension of the former classical proof. We show, through a discrete dispersion analysis of θ-schemes that the value θ = 1/12 reduces the numerical dispersion, but leads to a very restrictive CFL condition when applied to realistic piano strings.

Finally section 4 is concerned with the development and analysis of new θ-schemes, based on two different θ-approximations (defined by relation (28)) of the different wave propagating waves in the considered system. By adequately choosing two different values of θ we show that it is possible to construct stable schemes that reduce the numerical dispersion while allowing the use of a relatively large time-step.

All stability proofs of section 2, 3 and 4 will be done using energy techniques, which can easily be adjusted when dissipative terms, couplings or even non linear terms (see [START_REF] Chabassier | Energy preserving schemes for nonlinear Hamiltonian systems of wave equations: Application to the vibrating piano string[END_REF]) are added to the model (as opposed to Fourier techniques).

Continuous equation

The prestressed Timoshenko model considers two unknowns (u, ϕ): R + × [0, L] → R which stand respectively for the transversal displacement (m) and the shear angle (in radians) of the cross section of the the string. We denote L (m) its length, T 0 (N) the tension at rest, S (m 2 ) the cross section's area, ρ (kg•m -3 ) the density, I (m 4 ) the inertia momentum, E (Pa) the Young's modulus, G (Pa) the shear modulus, and κ ∈ [0, 1] the shear coefficient introduced in [START_REF] Timoshenko | On the correction for shear of the differential equation for transverse vibrations of bars of uniform cross-section[END_REF] and which value is discussed in [START_REF] Cowper | The shear coefficient in Timoshenko's beam theory[END_REF]. Finally we assume that these physical parameters are positive and that ES > T 0 (which is true in practice for piano strings). Moreover, we consider "simply supported" boundary conditions (zero displacement and zero torque).

The system reads: Find (u, ϕ) such that ∀

x ∈]0, L[, ∀ t > 0,        ρS ∂ 2 u ∂t 2 -T 0 ∂ 2 u ∂x 2 + SGκ ∂ ∂x ϕ - ∂u ∂x = σ, ρI ∂ 2 ϕ ∂t 2 -EI ∂ 2 ϕ ∂x 2 + SGκ ϕ - ∂u ∂x = 0, (1)
with boundary conditions

∂ x u(x = 0, t) = 0, ∂ x u(x = L, t) = 0, ∂ x ϕ(x = 0, t) = 0, ∂ x ϕ(x = L, t) = 0, (2)
and initial conditions

∂ t u(x, t = 0) = u 0 (x), ∂ t ϕ(x, t = 0) = ϕ 0 (x), ∂ t u(x, t = 0) = u 1 (x), ∂ t ϕ(x, t = 0) = ϕ 1 (x), (3)
where σ stands for a source term (which will come from the interaction with the hammer, in the context of the piano, see [START_REF] Chabassier | Modélisation et simulation numérique d'un piano par modèles physiques[END_REF], chapter I.2). In order to give existence and uniqueness properties we define the functional space naturally associated with simply supported boundary conditions

U 0 = {U ∈ H 1 0 ([0, L]) × H 1 ([0, L])}. (4)
System of equations ( 1), ( 2), (3) can be written as:

Find U ∈ U 0 such that ∀ x ∈ [0, L], ∀ t > 0,            M ∂ 2 t U -∂ x A∂ x U + BU + CU + t B∂ x U = M Σ, U (x, t = 0) = U 0 (x), ∂ t U (x, t = 0) = U 1 (x), ∂ x ϕ(x = 0, t) = 0, ∂ x ϕ(x = L, t) = 0,
where U = t (u, ϕ), Σ = M -1 t (σ, 0) and

M = ρS 0 0 ρI , A = T 0 + SGκ 0 0 EI , B = 0 -SGκ 0 0 , C = 0 0 0 SGκ . (5)
Lemma 1. The bilinear symmetric form defined for any U = (u, ϕ) ∈ U 0 and V = (v, ψ) ∈ U 0 by:

U, V U 0 = L 0 T 0 ∂ x u ∂ x v + L 0 EI ∂ x ϕ ∂ x ψ + L 0 SGκ (ϕ -∂ x u)(ψ -∂ x v), (6)
is a scalar product on U 0 . Hence U 0 is a Hilbert space with the induced norm denoted || • || U 0 .

Let us now introduce the operator associated with the prestressed Timoshenko system (1):

T (U ) = T (u, ϕ) =     T 0 ∂ 2 u ∂x 2 -SGκ ∂ ∂x ϕ - ∂u ∂x EI ∂ 2 ϕ ∂x 2 -SGκ ϕ - ∂u ∂x     (7)
and the functional space which is naturally induced:

D(T ) = U ∈ U 0 , such that T (U ) ∈ L 2 ([0, L]) 2 .
Theorem 2 (Unique strong solution).

If Σ ∈ C 1 0, T ; L 2 ([0, L]) , there exists a unique strong solution to (1)-(2)-(3) (u, ϕ) ∈ C 2 0, T ; (L 2 ([0, L]) 2 ) ∩ C 1 0, T ; U 0 ) ∩ C 0 0, T ; D(T ) . (8)
Proof. This result is a simple application case of Hille Yosida theorem, using lemma 1.

Inria

Energy and a priori estimates

Definition 2.1. We will use in the following the notation:

U ∈ U 0 , U 2
M := L 0 M U • U dx (9)
The next lemma shows the "energy preserving" nature of the system of wave equations ( 1), ( 2), (3). Our construction of a stable time discretization of the continuous equations will be based on this conservation property (as for instance in the work of [START_REF] Collino | A conservative space-time mesh refinement method for the 1-d wave equation. Part I: Construction[END_REF], [START_REF] Gilbert | Higher order time stepping for second order hyperbolic problems and optimal CFL conditions[END_REF], [START_REF] Chabassier | Introduction and study of fourth order theta schemes for linear wave equations[END_REF]).

Lemma 3 (Energy identity). Let U = t (u, ϕ) be the solution of (1)-( 2)-( 3), then

dE dt = L 0 M Σ • ∂ t U, with E(t) = E k (t) + E p (t). ( 10
)
where E k is the kinetic energy, defined by

E k (t) = 1 2 ∂ t U 2 M := 1 2 L 0 ρ S |∂ t u| 2 + 1 2 L 0 ρ I |∂ t ϕ| 2 ,
and E p is the potential energy, defined by

E p (t) = 1 2 U 2 U0 := 1 2 L 0 T 0 |∂ x u| 2 + 1 2 L 0 EI |∂ x ϕ| 2 + 1 2 L 0 SGκ |ϕ -∂ x u| 2 .
Proof. This identity is obtained by multiplying (1) by t (∂tu, ∂tϕ) and doing appropriate integrations by parts. Boundary terms vanish because of simply supported conditions (equations ( 2)).

It is important to understand what the implications of the relation [START_REF] Timoshenko | On the correction for shear of the differential equation for transverse vibrations of bars of uniform cross-section[END_REF] are at the continuous level, in order to later derive equivalent properties at a discrete level. More precisely it will allow us to establish stability conditions for the numerical schemes. For that reason and for the sake of completeness we give the proposition below that is a standard result in the literature (see [START_REF] Joly | Variational Methods for Time Dependent Wave Propagation[END_REF]).

Theorem 4 (A priori estimates). Let U = t (u, ϕ) be the solution of (1)-( 2)-( 3) then

E(t) ≤ E(0) + 1 √ 2 t 0 Σ(•, s) M 2 ,
where

E(0) = 1 2 U 1 2 M + 1 2 U 0 2 U0 . (11)
U (•, t) M ≤ U 0 M + t 2 E(0) + t 0 (t -s) Σ(•, s) M . (12)
Proof. Cauchy-Schwarz inequality on the energy identity gives

dE dt (t) = L 0 M Σ • ∂tU ≤ Σ M ∂tU M (13)
The potential energy being positive, we write

E(t) ≥ 1 2 ∂tU 2 M ⇒ ∂tU M ≤ 2 E(t). (14)
It is useful to notice that

dE dt = 2 √ E d √ E dt ⇒ 2 √ E d √ E dt ≤ Σ M ∂tU M ≤ Σ M √ 2 √ E ⇒ d √ E dt ≤ 1 √ 2 Σ M RR n°8088
thanks to [START_REF] Burel | Solving the Homogeneous Isotropic Linear Elastodynamics Equations Using Potentials and Finite Elements. The Case of the Rigid Boundary Condition[END_REF]. We now integrate over time between 0 and t:

E(t) -E(0) ≤ 1 √ 2 t 0 Σ(•, s) M ds, (15)
which gives [START_REF] Van Rensburg | Natural frequencies and modes of a Timoshenko beam[END_REF] with no difficulty. Moreover, we write

U (•, t) = U 0 (•) + t 0 ∂tU (•, s) ds ⇒ U M ≤ U 0 (•) M + t 0 ∂tU (•, s) M ≤ 2E(s)
ds which gives ( 12) using [START_REF] Van Rensburg | Natural frequencies and modes of a Timoshenko beam[END_REF] and an integration by parts.

Dispersion analysis

On finite domains, the dispersion analysis is the study of standing waves and their associated eigenfrequencies. In the case of the simply supported prestressed Timoshenko beam, these eigenfrequencies can be explicitly given. Our results supplement the work carried out in [START_REF] Van Rensburg | Natural frequencies and modes of a Timoshenko beam[END_REF][START_REF] Bhaskar | Elastic waves in Timoshenko beams: the "lost and found" of an eigenmode[END_REF] where implicit expressions are given for several boundary conditions (including the condition given by equation ( 2)).

Theorem 5 (Eigenfrequencies). If U (•, t) ∈ U 0 , under the form

U (x, t) = e -i 2πf t V (x), (16)
is solution to (1), (2), then there exists ℓ ∈ N * such that:

f = f ± ℓ = ω ± ℓ 2π , ω ± ℓ = Ω ± ℓ , (17)
where

Ω - ℓ = ℓ 2 π 2 L 2 (α + β) + αγ - √ ∆ ℓ 2αβ < Ω + ℓ = ℓ 2 π 2 L 2 (α + β) + αγ + √ ∆ ℓ 2αβ , (18)
with

∆ ℓ = ℓ 2 π 2 L 2 (α + β) + αγ 2 -4αβ γ(1 -δ) ℓ 2 π 2 L 2 + ℓ 4 π 4 L 4 ,
where we have defined

α = ρS T 0 + SGκ , β = ρ E , γ = SGκ EI , δ = SGκ T 0 + SGκ .
Proof. The proof of this theorem is based on a Fourier transform in time of the system (1), which leads to a condition under which eigenmodes satisfy the boundary conditions [START_REF] Bhaskar | Elastic waves in Timoshenko beams: the "lost and found" of an eigenmode[END_REF]. The associated eigenfrequencies are the researched real numbers f . The detailed proof is given in A.

Two branches of eigenfrequencies arise, which correspond to the two propagating waves in the system: flexural (corresponding to frequencies f - ℓ ) and shear (corresponding to frequencies f + ℓ ) waves. Explicit formula (17) can be developed for the first frequencies (ℓ small) of the flexural and shear branches: Inria Corollary 6. The first flexural eigenfrequencies of system (1), (2) satisfy, for small ℓ ≥ 1

f - ℓ = ℓf - 0 1 + ǫ ℓ 2 + O(ℓ 5 ), where f - 0 = 1 2L T 0 ρS , ǫ = π 2 2L 2 EI T 0 1 - T 0 ES . ( 19
)
The first shear eigenfrequencies satisfy, for small ℓ ≥ 1

f + ℓ = f + 0 1 + η ℓ 2 + O(ℓ 4 ), where f + 0 = 1 2π SGκ ρI , η = π 2 2L 2 EI + IGκ SGκ . ( 20
)
Remark 1 (High frequency behavior).

The high frequency behavior of the flexural branch is given by

f - ℓ ∼ ℓ 2L T 0 + SGκ SGκ , for ℓ → +∞. ( 21
)
On the other hand, The high frequency behavior of the shear branch reads:

f + ℓ ∼ ℓ 2L E ρ , for ℓ → +∞. ( 22
)
Figure 1 illustrates those two branches and their behaviors for the D♯1 string of the model D Steinway grand piano whose homogenized parameters are given in table 1. 1: Parameter values of the string D♯1. These parameters correspond to homogenized properties of experimentally measured piano strings (see [START_REF] Chabassier | Modélisation et simulation numérique d'un piano par modèles physiques[END_REF] for more details).

L S ρ T 0 m m 2 kg•m -3 N 1.
Remark 2. In most applications, only the flexural wave of Timoshenko beams is of interest: modeling the shear angle ϕ is a technical way of conferring inharmonicity to the flexural motion. In the case of the piano, the shear wave is not even transmitted to the structure (see [START_REF] Chabassier | Time domain simulation of a piano. Part I : model description[END_REF]). In consequence, it is important to develop a numerical scheme that reduces numerical dispersion mostly for flexural waves and this will be the objective of the scheme developed in section 4.

Classical theta-scheme : stability and dispersion analysis

We now recall some classical properties of the θ-scheme time discretization of (1), after a space discretization done by a variational method. Energy analysis is chosen to show stability because such analysis is easily adjustable when dissipative terms, couplings or even nonlinearities are added to the model. First we present a classical proof of stability for θ ≥ 1/4, then an innovative one for θ < 1/4 which is still valid when ∆t reaches its CFL upper bound. Finally, discrete eigenfrequencies will be calculated and we will see that the specific value θ = 1/12 reduces numerical dispersion, but leads to a severe CFL condition.

Let us consider a variational approximation of (1), with a subspace U h ⊂ U 0 of finite dimension N h , where h is a small parameter devoted to tend towards zero. The problem becomes a system of ODEs: Find

u h ∈ U h such that ∀ t > 0,    M h d 2 t U h (t) + K h U h (t) = M h Σ h (t), U h (0) = u 0 h , ∂ t U h (0) = u0 h , (23)
where U h is the vector of coordinates of u h in a chosen basis of U h , and for all

(U h , V h ) ∈ U 2 h , M h U h • V h := L 0 M u h • v h , (24)
K h U h • V h := L 0 A ∂ x u h • ∂ x v h + L 0 B u h • ∂ x v h + L 0 C u h • v h + L 0 t B ∂ x u h • v h . ( 25
)
where denotes the use of a quadrature formula. M h is a symmetric, positive definite matrix, while K h is a symmetric positive semi-definite matrix. In what follows we denote by {λ h,ℓ } (ℓ varies from 1 to N h ) the set of positive and increasing eigenvalues of of M -1 h K h .

Hypothesis 1. Let us consider the ordered sets of continuous and discrete eigenvalues:

• Λ = λ + i = 4π 2 (f + i ) 2 , i ≥ 1 ∪ λ - i = 4π 2 (f - i ) 2 , i ≥ 1 where f ± i are the eigenfrequencies of theorem 5, so that Λ = {λ i , i ≥ 1}, with λ i-1 ≤ λ i , • Λ h = {λ h,i , 1 ≤ i ≤ N h }, with λ h,i-1 ≤ λ h,i , the set of eigenvalues of M -1
h K h . We assume that for a fixed ℓ ∈ N * , the eigenvalue λ h,ℓ ∈ Λ h satisfies for h small enough

λ h,ℓ = λ ℓ + O(h 8 ). (26)
Remark 3 (Finite elements and eigenvalues).

If a spatial discretization is done with p th order finite elements on a regular mesh of size h, we can use a convergence theorem enunciated in [START_REF] Allaire | Numerical analysis and optimization: an introduction to mathematical modelling and numerical simulation[END_REF] to show that hypothesis 1 is true for p ≥ 4. In section 3.3 and 4.3 numerical results are obtained with fourth order finite elements.

In order to define the time discretization, we introduce a time step ∆t > 0 and define t n = n ∆t.

We denote the approximate unknown U n h :

U n h ≃ U h (t n ).
Applying a θ-scheme to problem (23) gives:

   M h U n+1 h -2U n h + U n-1 h ∆t 2 + K h {U h } n θ = M h Σ n h , U 0 h = u 0 h , U 1 h = u 1 h , (27)
Inria where the θ-approximation of U h (t n ) is a weighted average on three time steps:

{U h } n θ = θ U n+1 h + (1 -2θ) U n h + θ U n-1 h , (28)
which satisfies:

{U h } n θ = (θ - 1 4 )∆t 2 U n+1 h -2U n h + U n-1 h ∆t 2 + U n+1 h + 2U n h + U n-1 h 4 . (29)
We note that θ-scheme (27) can be written as an explicit scheme with a modified mass matrix:

M h,θ U n+1 h -2U n h + U n-1 h ∆t 2 + K h U n h = M h Σ n h where M h,θ = M h + θ ∆t 2 K h (30)
This modified mass matrix being positive definite, the numerical scheme admits a unique solution.

Definition 3.1. We will use in the following the notation, for all symmetric semi-definite matrix

A h ∀v ∈ R N h , v 2
A h := A h v • v (31)

Stability analysis

Stability of the numerical scheme ( 27) can be shown with energy techniques. First, we will show that any numerical solution satisfies an energy identity. If θ ≥ 1/4, this discrete energy is always positive, while if θ < 1/4, the time step ∆t must be lower than a maximal value ∆t θ . Then, we will show that the scheme is stable if the energy is positive. This last proof will be done in two steps : if θ ≥ 1/4, the proof is classical and the estimation is optimal compared to the estimation obtained at a continuous level, if θ < 1/4, we present an innovative proof based on a spectral decomposition, which is valid even if ∆t = ∆t θ (as opposed to the more classical proof given in [START_REF] Joly | Variational Methods for Time Dependent Wave Propagation[END_REF]).

Let us begin with the classical energy identity for the θ-scheme.

Lemma 7 (Discrete energy of θ-scheme). Any numerical solution to (27) satisfies

E n+1/2 θ -E n-1/2 θ ∆t = M h Σ n h • U n+1 h -U n-1 h 2∆t (32)
where

E n+1/2 θ
is the discrete energy defined by

E n+1/2 θ = E n+1/2 k,θ + E n+1/2 p,θ (33)
where

E n+1/2 k,θ and E n+1/2 p,θ
are respectively the discrete kinetic and potential energies at time

t n+1/2 : E n+1/2 p,θ = 1 2 U n+1 h + U n h 2 2 K h and E n+1/2 k,θ = 1 2 M h,θ U n+1 h -U n h ∆t • U n+1 h -U n h ∆t with M h,θ = M h + (θ - 1 4 )∆t 2 K h . (34)

RR n°8088

Proof. We take the scalar product of (27) with a centered approximation of the time derivative of

U h (t n ), which is (U n+1 h -U n-1 h
)/2∆t. The first term can directly be factorized, and we use (29) to factorize the second one. Adding the two results and taking into account the contribution of the right hand side gives (32).

Lemma 8 (θ-scheme's energy positivity). The discrete energy (33) is positive if and only if the matrix M h,θ is positive, which happens when:

• if θ ≥ 1/4,
the energy is positive for any ∆t > 0,

• if θ < 1/4, the energy is positive if and only if

∆t 2 ρ (M h ) -1 K h ≤ 4 1 -4θ where ρ (M h ) -1 K h = sup v =0 (M h ) -1 K h v • v v • v . ( 35
)
Proof. Since K h is positive semi-definite, the discrete potential energy is positive. The positivity of the discrete kinetic energy leads to the result.

Stability analysis for θ ≥ 1/4

Theorem 9 (A priori estimate for θ ≥ 1/4). We suppose that the discrete energy (33) is positive. Let U n h be the numerical solution of (27). For any n ≥ 1, we have:

E n+1/2 θ ≤ E 1/2 θ + ∆t √ 2 n ℓ=1 Σ ℓ h M h where E 1/2 θ = 1 2 U 1 h -U 0 h ∆t 2 M h,θ + 1 2 U 1 h + U 0 h 2 2 K h , (36)
U n+1 h M h ≤ U 0 h M h + t n+1 2E 1/2 θ + ∆t 2 n ℓ=1 ℓ k=1 Σ k h M h . ( 37
)
Proof. We use Cauchy-Schwarz inequality on (32), and triangular inequality on the right hand side:

1 ∆t E n+1/2 θ -E n-1/2 θ ≤ 1 2 Σ n h M h   U n+1 h -U n h ∆t M h + U n h -U n-1 h ∆t M h   . (38)
The discrete potential energy

E n+1/2 p,θ
being positive, we get:

U n+1 h -U n h ∆t M h,θ ≤ 2E n+1/2 θ . ( 39
) Since θ - 1 4 ≥ 0, we have for any X ∈ R N h , X M h ≤ X M h,θ , (40)
which implies, with (39), that we can estimate the right hand side of (38):

1 ∆t E n+1/2 θ -E n-1/2 θ ≤ 1 2 Σ n h M h 2E n+1/2 θ + 2E n-1/2 θ (41) which gives 1 ∆t E n+1/2 θ -E n-1/2 θ ≤ 1 √ 2 Σ n h M h (42)

Inria

and finally

E n+1/2 θ ≤ E n-1/2 θ + ∆t √ 2 Σ n h M h (43)
We sum this inequality from iterate n = 1, which gives a telescopic sum. We get:

E n+1/2 θ ≤ E 1/2 θ + ∆t √ 2 n ℓ=1 Σ ℓ h M h (44)
We bound above the norm of U n+1 h thanks to the triangular inequality and (39):

U n+1 h M h ≤ ∆t U n+1 h -U n h ∆t M h + U n h M h ≤ ∆t 2E n+1/2 θ + U n h M h ≤ √ 2∆t E 1/2 θ + ∆t √ 2 n ℓ=1 Σ ℓ h M h + U n h M h ,
we can again sum from n = 1 to obtain:

U n+1 h M h ≤ U 0 h M h + t n+1 2E 1/2 θ + ∆t 2 n ℓ=1 ℓ k=1 Σ k h M h .
Remark 4 (Optimality of theorem 9). This estimation is the exact discrete equivalent of the continuous estimations of lemma 3 and theorem 4.

Remark 5 (Transposition for θ < 1/4). It is still possible to adapt this proof when θ < 1/4, when the matrix M h,θ is not singular, which is the case when ∆t is chosen such that (35) is a strict inequality. In this case, the constant of a priori estimates depends on ∆t, more precisely one can show that instead of (40) we have

X M h ≤ ∆t θ (∆t θ ) 2 -∆t 2 1/2 X M h,θ ,
where ∆t θ is the maximal time step allowed by relation equation (35). This estimates blows up when ∆t reaches ∆t θ , which is not satisfying.

Remark 6. Similar estimates to those of lemma 9 can be obtained with a different technique (not based on the natural dicsrete energy

E n+1/2 θ
) as in [START_REF] Joly | Variational Methods for Time Dependent Wave Propagation[END_REF]. However, this technique cannot easily be extended if dissipative, coupling or nonlinear terms are added to the modeled problem, whereas the energy technique can. Moreover, it will later be possible to extend our proof to the new numerical scheme we will propose.

Stability analysis for θ < 1/4

As explained in remark 5, it is not possible to control the M h norm with the M h,θ norm when inequality (35) is an equality, i. e.:

ρ h := ρ (M h ) -1 K h = 4 (1 -4θ)∆t 2 . ( 45
)
If K h is singular (which can happen for other boundary conditions than (2)), there exist s eigenvalues of (M h ) -1 K h such that: λ h,1 = . . . = λ h,s-1 = λ h,s = 0. They prevent us from controlling the usual norm of associated eigenvectors with the K h induced semi-norm:

X M h ≤ C X K h , RR n°8088
for any strictly positive constant C. On the other hand, if M h,θ is singular (i.e. ∆t = ∆t θ ), there exists s such that λ

h,N h -s = . . . = λ h,N h -1 = λ h,N h = 4
(1-4θ)∆t 2 , which prevents us from controlling the usual norm of associated eigenvectors with the M h,θ induced semi-norm:

X M h ≤ C X M h,θ ,
for any strictly positive constant C. The original idea of our proof is then to divide the spectrum of (M h ) -1 K h in two parts (either sides of a given frequency), and to write the solution as a projection on the resulting high and low frequency subspaces. We then control the usual high frequency norm with the K h induced norm, and the usual low frequency norm with the M h,θ induced norm; which leads to an energy majoration in both cases since the discrete energy is the sum of the M h,θ semi-norm and the K h semi-norm of linear combinations of the solution at several time steps.

To this purpose we will use a high frequency projector P α h , defined below Definition 3.2. For any symmetric semi-definite matrix R h and any positive definite matrix M h let {W h,ℓ } N h ℓ=1 be the eigenvectors basis associated to the increasing ordered set of positive real eigenvalues

{λ h,ℓ } N h ℓ=1 such that    R h W h,ℓ = λ h,ℓ M h W h,ℓ , 1 ≤ ℓ ≤ N h , λ h,ℓ ≥ 0, (46a)
M h W h,ℓ • W h,m = δ ℓm , 1 ≤ ℓ, m ≤ N h . (46b)
Definition 3.3. For α ≥ 0, we define P α h associated to the family of eigenvectors and eigenvalues by:

∀ U h ∈ R N h , P α h U h = N h ℓ=Lα (M h W h,ℓ • U h ) W h,ℓ , (47)
where

L α ∈ [1, N h ] is the smallest integer such that λ Lα ≥ α.

Using definition 3.2 with

M h ≡ M h and R h ≡ K h ,
and applying the result of B (the reader will check that the hypotheses of definition 3.2 are satisfied) we can state the following upper bounds.

Lemma 10. Let 0 < α < ρ(M -1 h R h ) ≤ β, we have, for any U h ∈ R N h I h -P α h U h 2 M h ≤ β β -α (M h - 1 β R h )U h • U h and P α h U h 2 M h ≤ 1 α U h 2 R h ( 48
)
where I h is the identity matrix of size N h .

Proof. See the proof in B.

In addition to these non trivial identities, the high frequency projector P α h satisfies the classical properties

P α h P α h U h = P α h U h and ||U h || 2 M h = ||P α h U h || 2 M h + ||(I h -P α h ) U h || 2 M h . (49)
Let us now show the stability by energy techniques. The first result is an energy identity.

Inria

Lemma 11 (Energy estimate for θ < 1/4). We suppose that the discrete energy (33) is positive (i.e. condition (35) is satisfied). For all n ≥ 1, we have

E n+1/2 θ ≤ E 1/2 θ + √ 2 γ(θ) ∆t n ℓ=1 Σ ℓ h M h , ( 50
)
where E n+1/2 θ is given by (33), and

γ(θ) = 1 a(θ) + 1 4 -(1 -4θ)a(θ) 2 , a(θ) = 4 (1 -4θ) 2/3 + (1 -4θ) . ( 51
)
Proof. As in the other proof for θ ≥ 1/4, (32) gives:

1 ∆t E n+1/2 θ -E n-1/2 θ ≤ Σ n h M h U n+1 h -U n-1 h 2∆t M h I n (52) Let 0 < α ≤ ρ(M -1 h R h ),
we can write

I n ≤ P α h U n+1 h -U n-1 h 2∆t M h J n + (I h -P α h ) U n+1 h -U n-1 h 2∆t M h L n
The ideas is to use lemma (10) on J n and on L n . To do so we write

∆t J n ≤ P α h U n+1 h + U n h 2 M h + P α h U n h + U n-1 h 2 M h ≤ 1 √ α U n+1 h + U n h 2 R h + 1 √ α U n h + U n-1 h 2 R h and 2L n ≤ (I h -P α h ) U n+1 h -U n h ∆t M h + (I h -P α h ) U n h -U n-1 h ∆t M h ≤ β β -α U n+1 h -U n h ∆t M h -1 β R h + β β -α U n h -U n-1 h ∆t M h -1 β R h . We choose β = 4 (1 -4θ)∆t 2 ≥ ρ(M -1 h R h ) (53)
so that the latter term is the discrete kinetic energy. We also choose α = a(θ) 2 /∆t 2 with α < ρ(M -1 h R h ) and where a(θ) will be set later. We get

I n ≤ γ(θ) 2 E n+1/2 θ + 2 E n-1/2 θ ( 54
)
where γ(θ) was given in (51). Coming back to (52), we get

1 ∆t E n+1/2 θ -E n-1/2 ≤ √ 2 γ(θ) Σ n h M h ⇒ E n+1/2 θ ≤ E n-1/2 θ + ∆t √ 2 γ(θ) Σ n h M h .
Adding this telescopic sum from 1 to n, we get

E n+1/2 θ ≤ E 1/2 θ + ∆t √ 2 γ(θ) n ℓ=1 Σ ℓ h M h . ( 55
)
By inspection of the positivity condition (35) we see that this result is true for any

a(θ) ∈ ]0, 4 1 -4θ [ , RR n°8088
the best estimate is obtained when γ(θ) is minimal which happens for a(θ) given by (51). Unfortunately γ(θ) is still greater than 1/2, which will prevent us from reaching an "optimal" estimate (by comparison with the estimate obtained at the continuous level) as we did in proposition 9. Note that one can show that

γ(θ) ≤ √ 2,
which is obtained when θ = 0.

Thanks to lemma 11, it is now possible to establish a uniform a priori estimate on the numerical solution's M h norm, which results in the stability of the numerical scheme (27) provided that the discrete energy (33) is positive.

Theorem 12 (A priori estimate for θ < 1/4). We suppose that the discrete energy (33) is positive (i.e. condition (35) is satisfied). Let U n h be the numerical solution of (27). For any n ≥ 1, we have

U n+1 h M h ≤ √ 2 U 0 h M h + 2 γ(θ) t n+1 2E 1/2 θ + 4 γ(θ) 2 ∆t 2 n ℓ=1 ℓ k=1 Σ k h M h . ( 56
)
Proof. We will again use the frequency decomposition on the solution.

U n+1 h M h ≤ P α h U n+1 h M h M n + (I h -P α h ) U n+1 h M h N n . ( 57
)
To handle M n , we write:

P α h U n+1 h M h ≤ P α h U n h M h + 2 P α h U n+1 h + U n h 2 M h ≤ P α h U n h M h + 2 √ 2 √ α E n+1/2 θ ≤ P α h U 0 h M h + 2 √ 2 √ α n ℓ=1 E ℓ+1/2 θ .
To handle N n , we write:

(I h -P α h ) U n+1 h M h ≤ (I h -P α h ) U n h M h + ∆t (I h -P α h ) U n+1 h -U n h ∆t M h ≤ (I h -P α h ) U n h M h + ∆t 4 4 -∆t 2 (1 -4θ)α 2 E n+1/2 θ
, finally after successive substitutions we obtain

(I h -P α h ) U n+1 h M h ≤ (I h -P α h ) U 0 h M h + ∆t 4 4 -∆t 2 (1 -4θ)α n ℓ=1 2 E ℓ+1/2 θ . Combining M n and N n gives U n+1 h M h ≤ P α h U 0 h M h + (I h -P α h ) U 0 h M h + 2 √ 2 √ α + 2 ∆t √ 2 1 4 -∆t 2 (1 -4θ)α δ(θ) n ℓ=1 E ℓ+1/2 θ . We set α = a(θ) 2 /∆t 2 , so that δ(θ) = 2 √ 2 ∆t γ(θ)
. We now use [START_REF] Van Rensburg | Natural frequencies and modes of a Timoshenko beam[END_REF] to state that

U n+1 h M h ≤ P α h U 0 h M h + (I h -P α h ) U 0 h M h + 2 √ 2 ∆t γ(θ) n ℓ=1 E 1/2 θ + √ 2 γ(θ) ∆t ℓ k=1 Σ ℓ h M h .
The final result is obtained by noticing that, from property given by equation (49),

P α h U 0 h M h + (I h -P α h ) U 0 h M h ≤ √ 2 U 0 h M h .

Inria

The stability result of theorem 12 can be used to state a convergence theorem (at a semi discrete level). To do so, we denote the semi discrete error

e n h = U n h -U h (t n ) where U h (t n
) is the solution of (23) evaluated at time t n . Then, using a simple Taylor expansion, it is possible to show that e n h satisfies the scheme (27) where the source term Σ n h has been substituted by the consistency error (for simplicity we assume here that the initial conditions are zero for the semi-discrete problem). It is then easy to show (see [START_REF] Joly | Variational Methods for Time Dependent Wave Propagation[END_REF]) that, if the solution U h (t n ) is regular enough, the consistency error is of order 2 in ∆t for the θ-scheme (||Σ n h || M h ≤ C∆t 2 , for all n ≥ 1 and with C a generic constant depending only on U h ). Assuming that the restriction on the time step is satisfied, the conclusion of theorem 12 gives for any n ≥ 1:

e n h M h ≤ C (t n ) 2 ∆t 2 ,
where C is another generic constant depending only on U h .

In the specific case of the piano string, we are specifically interested in the low frequency behavior (for many reasons : the human ear cannot detect pitches higher than 20kHz, the hammer excitation barely exceeds 10kHz, and the spatial discretization with finite elements samples the solution down to a minimal wavelength which results in a maximal frequency). In the following we will focus on establishing accuracy estimates that depend on the frequency, by leading a discrete dispersion analysis.

Dispersion analysis

Numerical dispersion, which measures in finite domains the deviation of discrete eigenvalues from continuous ones, is an inherent fault in any numerical method (except for very special cases). In this section, we will quantify this deviation, which will give a good assessment on the quality of the approximation and is especially relevant in the context of musical acoustics.

Theorem 13. If the discrete problem (27) admits a solution of the form U n h = e i 2πf h n∆t V 0 h , and if the spatial discretization respects hypothesis 1, then there exists a positive integer ℓ such that f h = f h,ℓ with, for h and ∆t sufficiently small:

f h,ℓ = f ℓ + f 3 ℓ 2 1 12 -θ ∆t 2 + O(∆t 4 + h 4 ) (58)
where f ℓ = f ± ℓ is one of the eigenfrequencies of the continuous problem given in theorem 5. Proof. Although the proof is rather classical we present it for the sake of completeness. As we did in the continuous case (see section 2.2), we consider here a numerical plane wave:

U n h = e inω h ∆t V 0 h with V 0 h = 0.
Then, using this definition of U n h , we get

U n+1 h -2U n h + U n-1 h ∆t 2 = - 4 ∆t 2 sin 2 ω h ∆t 2 U n h .
We set

Ω 2 = 4 ∆t 2 sin 2 ω h ∆t 2 ∈ [0, 4 ∆t 2 ]
. Homogeneous θ-scheme (27) applied to this plane wave, gives:

-Ω 2 M h U n h + (1 -θ ∆t 2 Ω 2 )K h U n h = 0.

RR n°8088

If we denote λ h,ℓ the eigenvalues of M -1 h K h we find that there is only a finite number of eigenpulsations to the discrete system, which can be recovered by inverting the relation:

Ω 2 ℓ = (1 -θ ∆t 2 Ω 2 i )λ h,ℓ ⇔ ω h,ℓ = 2 ∆t Asin 1 2 ∆t 2 λ h,ℓ 1 + θ ∆t 2 λ h,ℓ .
Dividing by 2 π gives the classical expression of discrete eigenfrequencies:

f h,ℓ = 1 π∆t Asin 1 2 ∆t 2 λ h,ℓ 1 + θ ∆t 2 λ h,ℓ . (59)
It is possible to conduct a Taylor expansion for small enough ∆t to get:

f h,ℓ = λ h,ℓ 2π + λ 3/2 h,ℓ 4π 1 12 -θ ∆t 2 + O(∆t 4 ). ( 60
)
Since the continuous eigenfrequencies are f ℓ = √ λ ℓ /(2π), we get the expected result using the hypothesis 1.

Remark 7 (Value θ = 1/12). Equation (58) recalls the well known fact that θ = 1/12 plays a specific role: it sends the numerical dispersion back to fourth order. Unfortunately, this value is lower than 1/4, hence leads to a conditionally stable scheme (see proposition 8). We will see in section 3.3 that for realistic values of the piano string's coefficients, the CFL condition will be very severe.

Numerical illustration

Let us illustrate these results with numerical experiments. We use the D♯1 note of a model D Steinway grand piano, whose parameters were given above in table 1. Spatial discretization is done with fourth order finite elements on a 300 points regular mesh, so that hypothesis 1 is satisfied. We choose θ = 1/4 and ∆t = 10 -4 s. The right hand side is a pulse located 12.7 cm from one extremity, and we record the string motion along time at a point located 6 cm from the other extremity. Figure 2 shows the discrete Fourier transform obtained from the recorded signal during one second of the transversal displacement, for different frequency ranges : from 0 to 600 Hz in figure 2(a), from 1700 to 2700 Hz in figure 2(b), and from 3300 to 3700 Hz in figure 2(c). The first branch of continuous formula (17) is represented in diamonds (⋄), theoretical formula of θ-scheme approximation (59) is represented in circles (•). We can see that numerical dispersion causes a deviation of discrete eigenfrequencies from continuous eigenfrequencies, which is more and more pronounced as frequency increases. Moreover, theoretical formula (59) gives a very good assessment of the numerical behavior (blue spikes of the Fourier transform), both in low and high frequency ranges.

According to remark 7, the choice θ = 1/12 should reduce numerical dispersion. However, it leads to a conditionally stable scheme, which in our application case, is stable if ∆t ≤ 3.5 × 10 -7 s. If we had discretised the classical d'Alembert wave equation, the restriction would have been ∆t ≤ 5 × 10 -6 s. The shear wave, that we have modeled in order to account for the inharmonicity of the flexural wave, travels around 14 times faster, hence leads to a more severe CFL condition. This seems a great price to pay, especially as we are not particularly interested in a good approximation of the shear wave (see remark 2). In the remainder of this article we propose a new θ-scheme which allows us to reduce numerical dispersion for the flexural wave, while giving a less restrictive CFL condition.

Inria 4 New theta scheme : stability and dispersion analysis

The conclusion drawn in the previous paragraph is the same every time two waves propagate with coupling in the same system, but with very different velocities. In this context, using a conditionally stable scheme will constrain the time step to be adapted to the fastest wave (the shear wave in the Timoshenko system) whereas a larger time step would be sufficient for the slowest wave (the flexural wave in the Timoshenko system). Unconditionally stable schemes would be appealing if they did not induce so much numerical dispersion.

We propose a new time discretization where two different θ-approximations are done in the system, one being adapted to the slow wave and will in practice be done with θ = 1/12 to reduce numerical dispersion, and the other being adapted to the fast wave and will in practice be done with θ = 1/4 to avoid the stability condition. We will follow the same approach as in the previous section: after writing the scheme, we will write a discrete energy identity, a priori estimates on the solution, and lead a dispersion analysis before showing numerical results.

In order to write this scheme, we have to rewrite the continuous system by splitting the contribution to the flexural wave and shear wave. The choice of the splitting is done by focusing on the low frequency behavior of the flexural and shear waves. Corollary 6 shows that the flexural waves travel at speed T 0 /ρS at low frequency (i.e. small n in equation ( 19)), which suggests to decompose the matrix A into two sub-matrices:

A = T 0 + SGκ 0 0 EI = T 0 0 0 0 A + SGκ 0 0 EI A . (61)
We then define two matrices K h and K h as in (25), where A, B and C contribute for K h while only A contributes for K h . We consider the following scheme, with (θ, θ)

∈ [0, 1/2] 2 :      M h U n+1 h -2U n h + U n-1 h ∆t 2 + K h {U h } n θ + K h {U h } n θ = M h Σ n h U 0 h = u 0 h , U 1 h = u 1 h (62)
As we did in order to obtain (30), we can write this new scheme as an explicit scheme with a modified mass matrix:

(M h + θ ∆t 2 K h + θ ∆t 2 K h ) U n+1 h -2U n h + U n-1 h ∆t 2 + K h U n h = M h Σ n h
which shows the existence of a unique numerical solution since M h + θ ∆t 2 K h + θ ∆t 2 K h is always invertible.

Stability analysis

We present some energy techniques to show the stability of the new (θ, θ)-scheme (62) presented above.

Lemma 14 (Discrete energy of (θ, θ)-scheme). Any numerical solution to (62) satisfies

E n+1/2 θ,θ -E n-1/2 θ,θ ∆t = M h Σ n h • U n+1 h -U n-1 h 2∆t (63)
where the discrete energy is defined as

2 E n+1/2 θ,θ = U n+1 h -U n h ∆t 2 M h,θ,θ + U n+1 h + U n h 2 2 K h (64) with M h,θ,θ = M h - ∆t 2 4 (1 -4 θ)K h + (1 -4 θ)K h (65)
Proof. Both θ-approximations are written as in expression (29). We then take the scalar product of the scheme with the centered approximation of time derivative of U h , and the proof follows as for proposition 7.

By a simple rearrangement of the positivity of both terms of (64) we can derive the following sufficient conditions for the positivity of the energy E n+1/2 θ,θ independently of the solution:

Lemma 15 ((θ, θ)-scheme's energy positivity). The discrete energy (64) is positive if and only if the matrix M h,θ,θ is positive. We can give more precise sufficient conditions according to the values of (θ, θ):

• If θ ≥ 1/4 and θ ≥ 1/4, the energy is positive for any ∆t > 0.

• If θ < 1/4 and θ ≥ 1/4, the energy is positive if

∆t 2 ρ M -1 h K h ≤ 4 1 -4 θ (66)
• If θ < 1/4 and θ ≥ 1/4, the energy is positive if

∆t 2 ρ M -1 h K h ≤ 4 1 -4 θ (67)
• If θ ≤ 1/4 and θ ≤ 1/4, the energy is positive if and only if

∆t 2 ρ M -1 h (1 -4 θ)K h + (1 -4 θ)K h ≤ 4 (68)
Proof. Only the second and third case deserve some comments. Let us assume that θ < 1/4 and θ ≥ 1/4.

We want the matrix

M h -(∆t 2 /4) (1 -4 θ)K h + (1 -4 θ)K h to be positive. Since -(1 -4 θ) K h is positive, it is enough to verify that M h -(∆t 2 /4) (1 - 4
θ)K h is a positive matrix and we obtain the condition (66) with basic manipulations (in the same way we can derive (67)). This condition is however not optimal in the general case but easier to check in practice, indeed the sufficient and necessary condition reads

∆t 2 ρ M h + ∆t 2 1 -4 θ 4 K h -1 K h ≤ 4 1 -4 θ ,
which does not permit us to easily give an upperbound on ∆t .

We are now able to establish an energy identity, for which several cases arise according to the position of (θ, θ) compared to 1/4.

Inria

Lemma 16 (Energy estimate). We suppose that the discrete energy (64) is positive (i.e. the conditions given by lemma 15 are fulfilled). Let U n h be the numerical solution of (62). For any n ≥ 1, we have:

o If θ ≥ 1/4 and θ ≥ 1/4 E n+1/2 θ,θ ≤ E 1/2 θ,θ + ∆t √ 2 n ℓ=1 Σ ℓ h M h , o otherwise E n+1/2 θ,θ ≤ E 1/2 θ,θ + ∆t √ 2 γ min(θ, θ) n ℓ=1 Σ ℓ h M h ,
where γ(•) is defined by (51).

Proof. In this proof, in order to simplify the presentation and without any loss of generality, we assume that θ < θ. The proof for the three different cases will be similar, first to the proof of theorem 9, then to the proof of lemma 11.

First case: θ ≥ 1/4 and θ ≥ 1/4. The proof is similar to the proof of theorem 9 (formally replaced M h,θ by M h,θ,θ and

E n+1/2 θ by E n+1/2 θ,θ
) if one notices that, since θ ≥ 1/4:

U n+1 h -U n-1 h 2∆t M h ≤ U n+1 h -U n-1 h 2∆t M h,θ,θ
.

Then using the definition of the energy given by ( 64) we find

U n+1 h -U n-1 h 2∆t M h ≤ 1 √ 2 E n+1/2 θ,θ + E n-1/2 θ,θ
, the conclusion of the lemma is then a direct application of the steps (41)-(43) use in the proof of theorem 9.

Second case: θ < 1/4 and θ ≥ 1/4. By setting ). The key point being that

M h ≡ M h + ∆t 2 (4θ -
U n+1 h -U n-1 h 2∆t M h ≤ U n+1 h -U n-1 h 2∆t M h , since, by assumption 4θ -1 ≥ 0 and (U n+1 h + U n-1 h )/2 2 R h ≤ 2 E 1/2 θ,θ
.

Third case: θ < 1/4 and θ < 1/4. In a similar way, we introduce the matrices

M h ≡ M h and R h ≡ 1 -4 θ 1 -4 θ K h + K h (70)
and extend the proof of lemma 11. To do so we need to use the estimate

(U n+1 h + U n h )/2 2 R h ≤ 2 E n+1/2 θ,θ , valid since (1 -4 θ)/(1 -4 θ) < 1.
Theorem 17 (A priori estimate). We suppose that the discrete energy (64) is positive (i.e. the conditions given by lemma 15 are fulfilled). Let U n h be the numerical solution of (62). For any n ≥ 1, we have:

• If θ ≥ 1/4 and θ ≥ 1/4 U n+1 h M h ≤ U 0 h M h + t n+1 2E 1/2 θ,θ + ∆t 2 n ℓ=1 ℓ k=1 Σ k h M h , • otherwise U n+1 h M h ≤ √ 2 U 0 h M h +2 γ min(θ, θ) t n+1 2E 1/2 θ,θ +4 γ min(θ, θ) 2 ∆t 2 n ℓ=1 ℓ k=1 Σ k h M h .
Proof. Again, the proof for the three different cases are similar, first to the proof of theorem 9, then to the proof of lemma 12. When θ < 1/4 or θ < 1/4 the proof of lemma 12 must be adapted by choosing M h and R h as in ( 69) or (70) (depending of the values of θ and θ), the rest of the proof is identical after substituting formally

E n+1/2 θ by E n+1/2 θ,θ
Remark 8 (Stability when θ = 1/4). If the value θ = 1/4 is chosen to approximate the fast wave, we obtain the sufficient stability condition (66) which would be the same if we had applied a classical θ-scheme (with θ) on the equation described by the matrix K h corresponding to the slow wave (hence, less restrictive).

As before, the previous energy identity enabled us to deduce the stability estimates of theorem 17. For the (θ, θ)-scheme this estimate allows us to prove a convergence of order 2 in time.

However, as we will see in the next theorem, for specific values of (θ, θ) we can achieve higher order of accuracy for specific estimates (in terms of low frequency dispersion relation).

Dispersion analysis

Theorem 18. If the discrete problem (62) admits a solution of the form U n h = e i 2πf h n∆t V 0 h , and if the spatial discretization respects hypothesis 1, then there exists a fixed positive integer ℓ, such that f h = f h,ℓ where for h and ∆t sufficiently small:

f h,ℓ = f ℓ + O(∆t 2 + h 4 ) ( 71
)
where f ℓ are the eigenfrequencies of the continuous problem given in theorem 5. Moreover, for small ℓ, f h,ℓ = f + h,ℓ or f - h,ℓ , where

f - h,ℓ = ℓ f - 0 1 + ǫ ∆t ℓ 2 + O(ℓ 5 + ∆t 4 + h 4 ), (72a)
f + h,ℓ = f + 0,∆t 1 + η ∆t ℓ 2 + O(ℓ 3 + ∆t 4 + h 4 ), (72b)
with                    ǫ ∆t = ǫ + 2π 2 ∆t 2 1 12 -θ f - 0 2 , ( 73a
)
f + 0,∆t = f + 0 1 + (2πf + 0 ) 2 1 12 -θ ∆t 2 , ( 73b
)
η ∆t = η + π 2 (E + Gκ) 2ρL 2 θ - 1 12 ∆t 2 , ( 73c
)
where f - 0 , f + 0 , ǫ and η were defined in the corollary 6.

Inria

Proof. We prove this result by considering the semi discretization in time of the continuous system with our (θ, θ)-scheme. It is then very close to the continuous calculation of eigenfrequencies. More details are shown in C, including the explicit exact theoretical discrete eigenfrequencies in relation (103).

Remark 9 (Value θ = 1/12). We note that the specific value θ = 1/12 exactly provides, for discrete eigenfrequencies of the flexural wave, the same Taylor expansion as in the continuous case given up to O(ℓ 5 ) by (19) for small ℓ (see the corollary 6). We also notice that up to this term, this branch does not depend on the chosen value for θ.

Numerical illustration

Let us illustrate these results with numerical experiments. We use again D♯1 note of previous sections, whose parameters are summed up in table 1. Spatial discretization is done with fourth order finite elements on a 300 points regular mesh, as before, so that hypothesis 1 is satisfied. Time discretization is first done with the choice θ = 1/4 and θ = 1/2, in order to get an unconditionally stable scheme and to choose ∆t = 10 -4 s. As before, we record during one second the transversal displacement of a point located 6 cm from one extremity, and we represent its discrete Fourier transform in figure 3, for different frequency ranges. The first branch of continuous formula ( 17) is represented in red diamonds, theoretical formula of (θ, θ)-scheme approximation (103) is represented in magenta circles. We can see that the theoretical formula represent very well the numerical behavior, and that numerical dispersion again causes a deviation of discrete eigenfrequencies from continuous eigenfrequencies, which is more and more pronounced as frequency increases.

The main interest of this scheme was to choose, for the slow wave, a value of θ that diminishes numerical dispersion, and for the fast wave, a value of θ that ensures stability. Concretely, let us present a second numerical experiment where θ = 1/4, θ = 1/12 and ∆t ≤ 5×10 -6 s. Figure 4(a) shows the spectrum of the transversal displacement of the numerical solution as well as theoretical continuous (in red diamonds) and discrete (in magenta circles) eigenfrequencies, for a very high frequency range (5500 to 6500 Hz). We can see that numerical dispersion is very low, since the shift from continuous frequencies is only around 1 % at 5500 Hz, which is really good given the "large" chosen time step. Figure 4(b) shows the plot coming from the same experiment conducted with the usual θ-scheme with θ = 1/4, and we can clearly see that numerical dispersion is greater.

The last numerical illustration is presented in figure 5. We compare the explicit theoretical expressions of flexural eigenfrequencies for the continuous system (1)-( 2), for the classical θscheme (27) and for our new (θ, θ)-scheme (62). The explicit formulas (respectively (17), ( 59) and ( 103)) are used, since they have proven to very well reflect the numerical behavior of totally discrete schemes. Different time steps are considered ∆t = 10 -4 s , ∆t = 5 × 10 -6 s . As explained in remark 11 of C, relation (94) is no longer invertible after a certain rank, which depends on the time step (the maximal frequency for which the explicit formula is valid is shown in dashed black line). This is why the maximal considered rank changes between the subfigures. Figure 5 shows that the flexural dispersion deviation is always improved for the considered (θ, θ)-scheme. However this improvement is less pronounced when the time step ∆t decreases. This is explained by the fact that the (θ, θ)-scheme is designed to improve the approximation of the low frequency components of the solution. If ∆t is chosen relatively large then the (θ, θ)-scheme offers a real advantage because the classical θ-scheme is not even able to catch the low frequency behavior (see figure 5(a)). If ∆t is chosen very small then the low frequency components are well approximated by the two schemes and the benefit of the (θ, θ)-scheme is less spectacular even for the medium-frequency components of the solution (see figure 5(b)).

Conclusions and prospects

A simply supported prestressed Timoshenko beam can be used to model the motion of a stiff string as a piano string, for example. This system of PDEs describes the coupled propagation of flexural and shear waves, which have very different velocities. This article was concerned with finding a time discretization that reduced numerical dispersion while allowing quite a large time step. First, the classical θ-scheme were investigated, for which a new proof of stability was provided based on energy methods. This scheme is unconditionally stable when θ ≥ 1/4 and stable under a CFL condition when θ < 1/4. A dispersion analysis recalled that the value θ = 1/12 minimized numerical dispersion, but led to a very restrictive upper bound on the time step. We proposed a new time discretization based on two different θ-approximations according to the considered wave in the system. The stability analysis was done with energy methods, and the dispersion analysis was done on the semi-discrete scheme. In practice, the slow wave was evaluated with θ = 1/12 to reduce numerical dispersion, while the fast wave was evaluated with θ = 1/4 to avoid the severe stability condition. Numerical experiments showed that this new scheme gives very good results with no computational overcost.

The dissipative case is a very easy improvement of the work presented here. A centered term can be added to each damped equation, and all proofs based on energy identities still apply in this case. It would be interesting to see how this idea adapts to other wave systems where different velocities arise (as S and P waves for elastodynamics propagation in soft media, or acoustic and elastic waves in poro-elastic media, for instance). Another natural extension of this work would be to investigate the possibility to gain consistency orders. The authors have proposed in [START_REF] Chabassier | Introduction and study of fourth order theta schemes for linear wave equations[END_REF] new fourth order schemes based on the classical θ-scheme and modified equation technique, and it could be of interest to apply the same method on the new scheme presented above.

A Proof of Theorem 5

In order to obtain a dispersion relation, let us do a Fourier transform in time on the homogeneous system (1), with ω the Fourier variable:

       -ρS ω 2 û -(SGκ + T 0 ) ∂ 2 û ∂x 2 + SGκ ∂ φ ∂x = 0, -ρI ω 2 φ -EI ∂ 2 φ ∂x 2 -SGκ ∂ û ∂x + SGκ φ = 0. (74)

Inria

This system can be written as a first order formulation. Let Y = t (û, φ, ∂ x û, ∂ x φ). Then (74) is equivalent to the following system, with four unknowns:

∂ x Y + A(ω)Y = 0, where A(ω) =     0 0 -1 0 0 0 0 -1 a(ω) 0 0 d(ω) 0 b(ω) c(ω) 0     (75)
with

a(ω) = ρS T 0 + SGκ α ω 2 , b(ω) = ρI EI β ω 2 - SGκ EI γ , c = SGκ EI γ , d(ω) = - SGκ T 0 + SGκ δ . (76)
This system of coupled ODEs can be solved by studying the eigenvalues of the matrix A(ω), ie. the complex numbers λ(ω) such that

det(A -λ I 4 ) = 0 ⇔ λ 4 + λ 2 a + b -cd) + ab = 0. ( 77
)
This equation is a fourth degree equation that couples eigenvalues λ and Fourier variable ω:

λ 4 + (α + β) ω 2 λ 2 + αβ ω 4 -γ(1 -δ) λ 2 -αγ ω 2 = 0. (78)
Only even degrees are involved. We introduce Λ = λ 2 and Ω = ω 2 , giving the equation:

Λ 2 + (α + β) Ω -γ(1 -δ) Λ + αβ Ω 2 -αγ Ω = 0. (79)
This second order equation seen as a function of Λ can be solved by evaluating its discriminant:

∆(Ω) = (α+β) Ω-γ(1-δ) 2 -4 αΩ β Ω-γ = (α-β) 2 Ω 2 + 4αγ-2γ(α+β)(1-δ) Ω+γ 2 (1-δ) 2 .
This discriminant is itself a second order polynomial function in Ω whose roots of ∆ are negative. Indeed,

∆(Ω = 0) = γ 2 (1 -δ) 2 ≥ 0, ∆ ′ (Ω = 0) = 4αγ -2γ(α + β)(1 -δ) ≥ 0 since ES -T 0 ≥ 0, ∆ ′′ (Ω) = 2(α -β) 2 ≥ 0.
The previous inequalities imply that ∆(Ω) > 0 for Ω > 0. Since we are interested in ω ∈ R ⇒ Ω ≥ 0, equation (79) has two different real solutions Λ -≡ Λ -(Ω) and Λ + ≡ Λ + (Ω), as soon as Ω = 0, they read:

Λ -(Ω) = γ(1 -δ) -(α + β) Ω -∆(Ω) 2 , Λ + (Ω) = γ(1 -δ) -(α + β) Ω + ∆(Ω) 2 .
One can show that they satisfy:

Λ -(Ω) ≤ 0 ∀ Ω ≥ 0, Λ + (Ω) ≤ 0 ⇔ Ω ≥ γ β . (80)
Hence, matrix A has four eigenvalues which satisfy: Solutions of (75) can now be written in the basis of (e λ - -x , e λ - + x , e λ + -x , e λ + -x ) : there exists four functions (P, Q, R, S) of ω such that:

(λ - ± ) 2 = Λ -and (λ + ± ) 2 = Λ + . ( 81
û(x, ω) = P (ω)e λ - -x + Q(ω)e λ - + x + R(ω)e λ + -x + S(ω)e λ + + x . (82)
Any solution satisfies the boundary conditions ( 2) which can both be expressed on û:

û(x = 0, ω) = û(x = L, ω) = û′′ (x = 0, ω) = û′′ (x = L, ω) = 0
and give compatibility equations:

(L1) P + Q + R + S = 0, (L2) P (λ - -) 2 + Q(λ - + ) 2 + R(λ + -) 2 + S(λ + + ) 2 = 0, (L3) P e λ - -L + Qe λ - + L + Re λ + -L + Se λ + + L = 0, (L4) P (λ - -) 2 e λ - -L + Q(λ - + ) 2 e λ - + L + R(λ + -) 2 e λ + -L + S(λ + + ) 2 e λ + + L =
After algebraic manipulations, we get:

P e λ - -L -e λ - + L + R e λ + -L -e λ + + L = 0 P (Λ -) e λ - -L -e λ - + L + R(Λ + ) e λ + -L -e λ + + L = 0.
Finding a non-zero solution û is only possible if previous system has a non-zero solution. It must be singular, which gives:

det e λ - -L -e λ - + L e λ + -L -e λ + + L (Λ -) e λ - -L -e λ - + L (Λ + ) e λ + -L -e λ + + L = 0 ⇔ e λ - -L -e λ - + L (Λ + -Λ -) e λ + -L -e λ + + L = 0.
Since Λ + and Λ -are distinct, two cases arise: e λ - -L -e λ - + L = 0 or e λ + -L -e λ + + L = 0.this is equivalent to

sin( |Λ -|L) = 0 or sinh( √ Λ + L) = 0 if Ω < γ/β, sin( |Λ + |L) = 0 if Ω ≥ γ/β, which implies ∃ ℓ ∈ Z * , |Λ -|L = ℓ π, or √ Λ + = 0 if Ω < γ/β, ∃ ℓ ∈ Z * , |Λ -|L = ℓ π, or ∃ m ∈ Z * , |Λ + |L = m π if Ω ≥ γ/β. (83)
This illustrates the fact that since the domain is finite, only discrete wave numbers can arise in the string. In order to find the associated eigenfrequencies that satisfy (83), we must solve:

Λ -(Ω) = γ(1 -δ) -(α + β) Ω -∆(Ω) 2 = - ℓ 2 π 2 L 2 , (84)
Λ + (Ω) = γ(1 -δ) -(α + β) Ω + ∆(Ω) 2 = - m 2 π 2 L 2 . ( 85
)
One can easily show the uniqueness of Ω - ℓ and Ω + m such that

Λ -(Ω - ℓ ) = - ℓ 2 π 2 L 2 and Λ + (Ω + m ) = - m 2 π 2 L 2 .

Inria

The asymptotic behavior of these solutions are studied for Timoshenko beam (ie for T 0 = 0) in [START_REF] Van Rensburg | Natural frequencies and modes of a Timoshenko beam[END_REF][START_REF] Bhaskar | Elastic waves in Timoshenko beams: the "lost and found" of an eigenmode[END_REF]. For general boundary conditions, these equations must be solved numerically.

In this case (simply supported conditions), it is possible to explicitly express Ω - ℓ and Ω + m . Indeed, Λ -(Ω - ℓ ) is by definition solution of the second order equation (79) parametrized by Ω - ℓ , so it satisfies:

Λ -(Ω - ℓ ) 2 + (α + β) Ω - ℓ -γ(1 -δ) Λ -(Ω - ℓ ) + αβ (Ω - ℓ ) 2 -αγ Ω - ℓ = 0. (86)
We can substitute Λ -(Ω - ℓ ) with its value in this equation:

αβ (Ω - ℓ ) 2 - ℓ 2 π 2 L 2 (α + β) + αγ Ω - ℓ + γ(1 -δ) ℓ 2 π 2 L 2 + ℓ 4 π 4 L 4 = 0.
We (once again) have to solve a second order equation in Ω - ℓ , whose discriminant is:

∆ ℓ = ℓ 2 π 2 L 2 (α+β)+αγ 2 -4αβ γ(1-δ) ℓ 2 π 2 L 2 + ℓ 4 π 4 L 4 = α 2 γ 2 +2αγ ℓ 2 π 2 α + β -2β(1 -δ)] L 2 +(α-β) 2 ℓ 4 π 4 L 4 ≥ 0. (87)
This equation has two solutions, which are exactly the two solutions we would have obtained by inverting Λ -and Λ + with equations ( 84) and (85). They are positive (as roots of an upturned parabola being positive with negative slope at origin) and write, for a fixed ℓ ∈ N * , as (18).

B Definition of the high frequency projection operator and proof of lemma 10

Let us first re-introduce some notations.

Definition. For any symmetric semi-definite matrix R h and any positive definite matrix M h let {W h,ℓ } N h ℓ=1 be the eigenvectors basis associated to the increasing ordered set of positive real eigenvalues {λ h,ℓ } N h ℓ=1 such that

   R h W h,ℓ = λ h,ℓ M h W h,ℓ , 1 ≤ ℓ ≤ N h , λ h,ℓ ≥ 0, (88a)
M h W h,ℓ • W h,m = δ ℓm , 1 ≤ ℓ, m ≤ N h . (88b)
Definition. For α ≥ 0, we define the high frequency projector P α h associated to the family of eigenvectors and eigenvalues by:

∀ U h ∈ R N h , P α h U h = N h ℓ=Lα (M h W h,ℓ • U h ) W h,ℓ , (89)
where

L α ∈ [1, N h ] is the smallest integer such that λ Lα ≥ α.
As the vectors {W h,ℓ } N h ℓ=1 are orthonormal with respect to the scalar product induced by M h one can expand any vector U h ∈ R N h in the following way

U h = N h ℓ=1 (M h W h,ℓ • U h ) W h,ℓ . (90)
The first property we want to prove shows that we can bound the low frequency components of any U h ∈ R N h using a semi-norm induced by

M h -(1/β) R h : Lemma 19. Let 0 ≤ α < ρ(M -1 h R h ) ≤ β, we have I h -P α h U h 2 M h ≤ β β -α (M h - 1 β R h )U h • U h . (91)
Proof. To obtain such an estimate, we first expand the vector U h on the basis of the {W h,ℓ } N h ℓ=1 and use the definition of the eigenvectors. We have

(M h - 1 β R h )U h • U h = N h ℓ=1 (1 - λ h,ℓ β ) M h W h,ℓ • U h 2 .
By hypothesis, (1 -

β -1 λ h,ℓ ) is positive, hence N h ℓ=1 1 - λ h,ℓ β M h W h,ℓ • U h 2 ≥ Lα-1 ℓ=1 1 - λ h,ℓ β M h W h,ℓ • U h 2 ≥ 1 - α β Lα -1 ℓ=1 M h W h,ℓ • U h 2 ,
by definition,

Lα-1 i=1 M h W h,ℓ • U h 2 = (I h -P α h ) U h 2 M h ,
we get the wanted inequality by inverting (1β -1 α).

The other result we prove give a bound on the higher frequency part of any vector:

Lemma 20. For α > 0, we have, for any

U h ∈ R N h P α h U h 2 M h ≤ 1 α U h 2 R h . (92)
Proof. By definition of the projection operator P α h we find, using (88a) and (88b),

P α h U h 2 M h = N h ℓ=Lα M h W h,ℓ • U h 2 = N h ℓ=Lα λ h,ℓ λ h,ℓ M h W h,ℓ • U h M h W h,ℓ • U h .
Now, as, by assumptions λ h,ℓ ≥ α we obtain

P α h U h 2 M h ≤ 1 α N h ℓ=Lα λ h,ℓ M h W h,ℓ • U h M h W h,ℓ • U h ≤ 1 α N h ℓ=1 λ h,ℓ M h W h,ℓ • U h 2
where the last inequality is true as we add only positive terms. Using the decomposition (90) valid for all U h and omitting details, we find that,

N h ℓ=1 λ h,ℓ M h W h,ℓ • U h 2 = N h ℓ=1 λ h,i M h W h,ℓ • U h M h W h,ℓ • U h = N h ℓ=1 R h W h,ℓ • U h M h W h,ℓ • U h = N h ℓ=1 M h W h,ℓ • U h W h,ℓ • R h U h which gives N h ℓ=1 λ h,ℓ M h W h,ℓ • U h 2 = U h • R h U h ,
and so, we obtain the final result combining this inequality and the inequality (B) obtained above.

Inria

C Proof of Theorem 18

As before, we consider a numerical plane wave:

U n h = e inω h ∆t V 0 h ( 93
)
and want to determine the frequency f n = ω h /2 π, to do so we denote

Ω 2 h = 4 sin 2 ( ω h ∆t 2 ) ∆t 2 ∈ [0, 4 ∆t 2 ]. (94)
Homogeneous scheme (62) applied to this plane wave gives:

-

Ω 2 h M h V 0 h + (1 -θ∆t 2 Ω 2 h )K h V 0 h + (1 -θ∆t 2 Ω 2 h )K h V 0 h = 0,
which can be rewritten as

Ω 2 h (M h + θ ∆t 2 K h + θ ∆t 2 K h ) V 0 h = K h V 0 h . (95)
We propose to use the semi discrete system in order to perform a dispersion analysis of Timoshenko system where a (θ, θ)-scheme has been used, with a method similar to the continuous case. The equation (95) (which is fully discrete) can be used to state the following eigenvalues problem: Find every real number Φ, such that there exists V 0 ∈ U 0 that satisfies

Φ M + θ∆t 2 ∂ x A ∂ x + B + C + t B ∂ x + θ∆t 2 ∂ x A∂ x V 0 = ∂ x A ∂ x + B + C + t B ∂ x V 0 (96)
where M , B and C are Timoshenko system's matrices, while the matrix A is separated as in (61). From the positivity and symmetry properties of the operator in the left and right hand sides of the previous equations we know that there exists an increasing positive sequence of Φ ℓ (with associated V 0 ℓ ) such that (96) is satisfied. Then we assume that for any fixed ℓ the spatial discretization is sufficiently fine (as in hypothesis 1) so that

Ω 2 h,ℓ = Φ ℓ + O(h 8 ). (97)
We also assume that ∆t is sufficiently small, indeed from relation (94) we see that Ω 2 h,ℓ ≤ 4/∆t 2 , and so equation (97) can not be valid uniformly with respect to ℓ (this reflects the fact that the problem discretized in time but not in space is ill-posed).

To determine Φ we write the first order formulation of this problem, as in the proof of proposition 5:

∂ x Y + AY = 0 (98)
where the matrix A is however different and reflects the time discretization:

A =     0 0 -1 0 0 0 0 -1 a(Φ) 0 0 d(Φ) 0 b(Φ) c(Φ) 0     with a(Φ) = ρSΦ 2 c 1 (Φ) , b(Φ) = ρIΦ 2 c 2 (Φ) - SGκ c 2 (Φ) (1 -θ∆t 2 Φ 2 ), c(Φ) = SGκ EI , d(Φ) = - SGκ c 2 (Φ) EI c 1 (Φ) c 1 (Φ) = (1 -θ∆t 2 Φ 2 )SGκ + (1 -θ∆t 2 Φ 2 )T 0 , c 2 (Φ) = (1 -θ∆t 2 Φ 2 )EI.
This matrix admits complex eigenvalues such that

det(A -λI 4 ) = 0 ⇔ λ 4 + λ 2 a(Φ) + b(Φ) -c(Φ)d(Φ) + a(Φ)b(Φ) = 0.
We set Λ = λ 2 , and we want to nullify

f φ (Λ) = Λ 2 c 1 (Φ)c 2 (Φ)+Λ ρSΦc 2 (Φ) + ρIΦc 1 (Φ) -SGκ(1 -θ∆t 2 Φ)c 1 (Φ) + (SGκ) 2 (1 -θ∆t 2 Φ) 2 + ρSΦ ρIΦ -(1 -θ∆t 2 Φ)SGκ . ( 99
)
We are only interested in solutions Λ(Φ) ∈ R -, associated to imaginary λ which lead to propagating solutions of the system. Let us assume (see remark 10) that there exists one or two solutions Λ -(Φ) < 0 to (99). We can apply the same reasoning as in the continuous case : boundary conditions imply the existence of an integer ℓ ∈ N * such that

Λ -(Φ) = - ℓπ L 2 , f Φ Λ -(Φ) = 0 (100)
and we need to invert this relation in order to find associated Φ ℓ .

Remark 10 (Solutions to (99)). Only the situation where Φ > 0 deserves a remark: we see that when ∆t 2 Φ = 0, equation (99) is similar to equation (79) (with Φ = Ω 2 ) for which we know that there always exist two distinct real solutions with one of them always negative. This implies, by continuity arguments, that for ∆t 2 Φ small enough such negative solution exists.

As in the continuous case, we inject the resulting expression of Λ -(Φ) in the equation (99) and express it as an equation in Φ ℓ :

α ℓ Φ 2 ℓ + β ℓ Φ ℓ + γ ℓ = 0, (101)
where

                 α ℓ = ∆t 4 ℓ 4 π 4 L 4 θEI(θSGκ + θT 0 ) + ∆t 2 ℓ 2 π 2 L 2 ρI(θES + θSGκ + θT 0 ) + ∆t 2 SGκT 0 θθ + ρS ρI + ∆t 2 θSGκ , β ℓ = -∆t 2 ℓ 4 π 4 L 4 EI(2θSGκ + (θ + θ)T 0 ) - ℓ 2 π 2 L 2 ρI(ES + SGκ + T 0 ) + ∆t 2 (θ + θ)T 0 SGκ -ρS(SGκ), γ ℓ = ℓ 4 π 4 L 4 EI(SGκ + T 0 ) + SGκT 0 ℓ 2 π 2 L 2 .
We introduce ∆ ℓ = β 2 ℓ -4α ℓ γ ℓ . It is then possible to express ∆ ℓ as a polynom of ∆t parametrized by θ, θ and ℓ, so we introduce the notation ∆ ℓ ≡ ∆ ℓ (∆t ; θ, θ) This discriminant is not simple, but we can give some specific behaviors:

⊲ When θ = θ. Then, the discriminant simplifies to ∆ ℓ (∆t, θ, θ) = ℓ 4 ρ 4 π 4 I 4 L 4 (ES -T 0 -SGκ) 2 + ℓ 2 2π 2 ρ 2 A 2 GκI L 2 (SGκ + ES -T 0 ) + ρ 2 S 4 G 2 κ 2 .
It does not depend on ∆t nor θ and is always positive (we recall that by assumption ES -T 0 ≥ 0).

⊲ When θ = θ. It can be shown, by explicit computations, that the roots ∆t ± ℓ of ∆ ℓ (∆t ; θ, θ) are complex, and we deduce that ∆ ℓ (∆t ; θ, θ) always stays positive since ∆ ℓ (0 ; θ, θ) ≥ 0 (the value of ∆ ℓ (0 ; θ, θ) being given by equation (87) up to a positive multiplicative factor).

Inria

As seen above, ∆ ℓ is positive, hence there exist two real roots which read:

Φ + ℓ = -β ℓ + √ ∆ ℓ 2α ℓ > Φ - ℓ = -β ℓ - √ ∆ ℓ 2α ℓ . ( 102
)
When possible (see remark 11), we invert equation (94) for sufficiently small ∆t and we find

f ± h,ℓ = 1 π∆t Asin ∆t 2 Φ ± ℓ + O(h 8 ) 4 . (103)
It is then possible to conduct a Taylor expansion for small ∆t, to obtain

Φ ± ℓ = (2 π f ± ℓ ) 2 + O(∆t 2 ) ⇒ f ± h,ℓ = f ± ℓ + O(∆t 2 + h 4
). Finally, using a Taylor expansion with respect to ∆t and ℓ in expression (102, 103) we obtain the result of Theorem 18. Remark 11. A necessary condition for (94) to be invertible is that Φ ± ℓ ∈ [0, 4/∆t 2 ]. This leads to necessary conditions on ∆t and ℓ, which are not easy to comprehend (see [START_REF] Chabassier | Modélisation et simulation numérique d'un piano par modèles physiques[END_REF]), and become more and more restrictive on ∆t when ℓ grows.

)

  Their expression depends on the value ofΩ: if 0 ≥ Ω < γ/β, only λ - ± = ±i |Λ -(Ω)| corresponds to imaginary (hence propagative) eigenvalues, whereas if Ω ≥ γ/β, λ - ± = ±i |Λ -(Ω)| and λ + ± = ±i |Λ + (Ω)| correspond to imaginary eigenvalues.

  Two branches of eigenfrequencies : upper figure (in blue) for flexural frequencies f - ℓ , lower figure (in red) for shear frequencies f + ℓ . The high frequency behaviors are represented in dashed lines. Low frequency behavior of the flexural branch f - ℓ of equation (17) (in blue) compared to its approximate expression for small ℓ of equation (19) (in red).

Figure 1 :

 1 Figure 1: Eigenfrequencies of the simply supported prestressed Timoshenko system, for note D♯1.

Figure 2 :

 2 Figure 2: Transversal displacement's spectrum of string D♯1 using a θ-scheme with θ = 1/4 and ∆t = 10 -4 s.

Figure 3 :Figure 4 :

 34 Figure 3: Transversal displacement's spectrum of string D♯1 using a (θ, θ)-scheme with θ = 1/2, θ = 1/4 and ∆t = 10 -4 .

  ∆t = 5 × 10 -6 s

Figure 5 :

 5 Figure5: Comparison of theoretical eigenfrequencies of the continuous system (black circles •), new (θ, θ)-scheme with θ = 1/4 and θ = 1/12 (red dimonds ⋄) and usual θ-scheme with θ = 1/4 (blue plus sign +). The theoretical curves are plotted for the first eigenfrequencies after which numerical stability is no longer granted.

51 × 10 -5 1.99 × 10 5 3.95 × 10 -6 Table

	945 E	1.31 × 10 -6 I	44290 G	1328 κ
	Pa	m 4	Pa	-
	2.02 × 10 11 1.78 × 10 -14 1.00 × 10 10 0.95 f -0 ǫ f + η 0 Hz -Hz -
	38.89 3.

RR n°8088

Avenue de la Vieille Tour, 33405 Talence Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Chabassier & ImperialeInria