Yacine Bouzidi
email: bouzidi.yacine@gmail.com

Alban Quadrat
email: alban.quadrat@inria.fr

Fabrice Rouillier
email: fabrice.rouillier@inria.fr

Grace Younes
email: grace.younes@inria.fr

Computation of the L ∞ -norm of finite-dimensional linear systems

Keywords: L∞-norm computation, Real roots, Symbolic computation, Complexity computation, Implementation, Control theory

In this paper, we study the problem of computing the L∞norm of finite-dimensional linear time-invariant systems. This problem is first reduced to the computation of the maximal x-projection of the real solutions (x, y) of a bivariate polynomial system Σ ={P, ∂P ∂y }, with P ∈ Z[x, y]. Then, we use standard computer algebra methods to solve the problem. In this paper, we alternatively study a method based on rational univariate representations, a method based on root separation, and finally a method first based on the sign variation of the leading coefficients of the signed subresultant sequence and then based on the identification of an isolating interval for the maximal x-projection of the real solutions of Σ.

Introduction

An important issue in robust control theory is the computation of the L ∞norm of linear systems [START_REF] Doyle | Feedback Control Theory[END_REF][START_REF] Zhou | Robust and Optimal Control[END_REF]. Contrary to the L 2 -norm, no tractable formula is known for the characterization of the L ∞ -norm of finite-dimensional systems (i.e., systems defined either by linear ordinary differential equations or by linear recurrence relations) [START_REF] Doyle | Feedback Control Theory[END_REF][START_REF] Zhou | Robust and Optimal Control[END_REF]. Hence, the standard methods for the L ∞ -norm computation are numerical (e.g., bisection algorithms, eigenvalues computation of Hamiltonian matrices) [START_REF] Boyd | A bisection method for computing the h∞ norm of a transfer matrix and related problems[END_REF][START_REF] Chen | A fast algorithm to compute compute the h∞-norm of a transfer function matrix[END_REF]. In their paper [START_REF] Kanno | Validated numerical computation of the l∞-norm for linear dynamical systems[END_REF], Kano and Smith develop a validated numerical algorithm for the L ∞ -norm computation. They reduce the problem to the localization of the real solutions of a bivariate polynomial and then use Sturm chain tests to guarantee the accuracy of their algorithm. In [START_REF] Chen | Computing the supremum of the real roots of a parametric univariate polynomial[END_REF], Chen, Mazza and Xie provide an equivalent study using the theory of border polynomials, which makes the presentation of their solution simpler.

When numerical methods are used, it is worth mentioning that the result is usually obtained within a short time but with a slight error up to a precise accuracy. In contrast, when using symbolic methods, the result usually takes more time to be computed but is exact. In this paper, to compute the L ∞ -norm, we try to develop the right balance between these two approaches.

In this paper, following the approach developed in [START_REF] Kanno | Validated numerical computation of the l∞-norm for linear dynamical systems[END_REF][START_REF] Chen | Computing the supremum of the real roots of a parametric univariate polynomial[END_REF], we shall study three different certified symbolic-numeric algorithms for the computation of the L ∞ -norm with the goal of minimizing the drawback of the symbolic part of the computation. This symbolic part consists in computing an isolating interval of the maximal projection of the real solutions of a system of bivariate polynomials. We develop the complexity analysis of each algorithm. Finally, we compare the theoretical complexities of the algorithms and then their performances using an implementation in the computer algebra system Maple.

Given two coprime polynomials P and Q in Z[x, y] of degree bounded by d and of coefficient bitsize bounded by τ , the solving of the system Σ = {P, Q} can be studied using numerous methods. Typically, isolating boxes of the solutions can be computed either directly from the input system using numerical methods (such as subdivision or homotopy methods) or indirectly by first computing intermediate symbolic representations such as triangular sets, Gröbner bases, or rational parameterizations [START_REF] Aubry | Triangular sets for solving polynomial systems: a comparative implementation of four methods[END_REF][START_REF] Boztas | Applied Algebra, Algebraic Algorithms and Error-Correcting Codes: 14th International Symposium, AAECC-14[END_REF].

Two methods used in the paper require putting the system in a generic position, i.e., require to finding a separating linear form x + a y that defines a shear of the coordinate system (x, y), i.e., (x, y) -→ (t -a y, y), so that no two distinct solutions of the sheared system Σ a = {P (t -a y, y), Q(t -a y, y)} are vertically aligned. This approach has long been used in the literature. For instance, a separating linear form x + a y with a ∈ {0, . . . , 2 d 4 } can be computed as shown in [START_REF] Bouzidi | Separating linear forms and rational univariate representations of bivariate systems[END_REF][START_REF] Bouzidi | Solving bivariate systems using rational univariate representations[END_REF]. We can then use a Rational Univariate Representation (RUR) for the polynomial system Σ a followed by the computation of isolating boxes for its real solutions. For more details, see [START_REF] Bouzidi | Solving bivariate systems using rational univariate representations[END_REF]. We simply apply this approach (i.e., the so-called RUR method ) to the polynomial system associated with the L ∞norm computation problem and then choose the maximal x-projection of the real solutions of the system. The complexity analysis shows that this algorithm performs ÕB (d x d 3 y (d 2 x + d x d y + d y τ )) bit operations in the worst case, where

d x = max(deg x (P ), deg x (Q)), d y = max(deg y (P ), deg y (Q)), (1)
and τ is the maximal coefficient bitsize of the polynomials P and Q. Alternatively, we can also localize the maximal x-projection of the real solutions of the polynomial system Σ by simply applying a linear separating form on the system Σ. The linear separating form t = x + s y proposed in [START_REF] Cheng | Root isolation for bivariate polynomial systems with local generic position method[END_REF] preserve the order of the solutions of the sheared system Σ s = {P (t -s y, y), Q(t -s y, y)} with respect to the x-projection of the real solutions of the original system Σ. Thus, the projection of the solutions of Σ s onto the new separating axis t can be done so that we can simply choose the x-projection corresponding to the maximal t-projection of the real solutions of Σ s . The drawback of this method lies on the growth of the size of the coefficients of the sheared system for the linear separating form t = x + s y due to the large size of s. The complexity analysis shows that this algorithm performs ÕB (d 4

x d 5 y τ ) bit operations in the worst case. The third method developed in this paper localizes the maximal x-projection of the system real solutions -denoted by x -by first isolating the real roots of the univariate resultant polynomial Res(P, ∂P ∂y , y) and then verifying the ex-istence of a real root of the univariate polynomial P (x, y) ∈ R[y] as done in [START_REF] Chen | Computing the supremum of the real roots of a parametric univariate polynomial[END_REF].

A key point is that we can compute a Sturm-Habicht sequence [START_REF] González-Vega | Sturm-habicht sequences, determinants and real roots of univariate polynomials[END_REF] of P (x, y) without any consequent overhead. As P (x, y) = 0 is bounded in the x-direction, it is then possible to compute the number of real solutions of P (x, y) with a good complexity in the worst case, which gives an efficient algorithm as soon as the curve P (x, y) = 0 has no isolated real singular points. The complexity analysis shows that this algorithm performs ÕB (d 2 x d 4 y (d x + τ )) bit operations in the worst case and ÕB (d 2

x d 4 y τ ) when the plane curve P (x, y) = 0 has no isolated real singular points.

Finally, we conclude the paper by comparing the bit complexity of those three algorithms and then the experimental time obtained by the implementation of each of these algorithms in Maple.

Problem description

Before stating the problem studied in this paper, we first introduce a few standard notations and definitions. If k is a field and P ∈ k[x, y], then Lc var (P ) is the leading coefficient of P with respect to the variable var ∈ {x, y} and deg var (P ) the degree of P in the variable var ∈ {x, y}. We also denote by deg(P ) the total degree of P . Moreover, let π x : R 2 -→ R be the projection map from the real plane R 2 onto the x-axis, i.e., π x (x, y) = x for all (x, y) ∈ R 2 . For P, Q ∈ k[x, y], let gcd(P, Q) be the greatest common divisor of P and Q, I := P, Q the ideal of k[x, y] generated by P and Q, and Definition 1 ( [START_REF] Doyle | Feedback Control Theory[END_REF][START_REF] Zhou | Robust and Optimal Control[END_REF]). Let RH ∞ be the R-algebra of all the proper and stable rational functions with real coefficients, namely:

V K (I) := {(x, y) ∈ K 2 | ∀ R ∈ I : R(x, y) = 0},
RH ∞ := n d | n, d ∈ R[s], gcd(n, d) = 1, deg s (n) ≤ deg s (d), V C ( d ) ∩ C + = ∅ .
An element g of RH ∞ is holomorphic and bounded on C + , i.e.,

g ∞ := sup s∈C+ |g(s)| < +∞,
RH ∞ is a sub-algebra of the Hardy algebra H ∞ (C + ) of bounded holomorphic functions on C + . The maximum modulus principle of complex analysis yields:

g ∞ = sup ω∈R |g(i ω)|.
Note that this equality shows that the function g |i R : i ω ∈ i R -→ g(i ω) belongs to the Lebesgue space L ∞ (i R) or, more precisely, to the following R-algebra

RL ∞ := n(i ω) d(i ω) | n, d ∈ R[i ω], gcd(n, d) = 1, deg ω (n) ≤ deg ω (d), V R ( d ) = ∅},
i.e., the algebra of real rational functions on the imaginary axis i R which are proper and have no poles on i R, or simply, the algebra of real rational functions with no poles on i P 1 (R), where P 1 (R) := R ∪ ∞.

We can extend the above L ∞ -norms defined on functions of RH ∞ (resp., RL ∞ ) to matrices as follows. Let G ∈ RH u×v ∞ (resp., G ∈ RL u×v ∞ , R(s) u×v ), i.e., G is a u × v matrix with entries in RH ∞ (resp., RL ∞ , R(s)) and let σ (•) denote the largest singular value of a complex matrix. Then, we can define:

G ∞ := sup s∈C+ σ (G(s)) resp., G ∞ := sup ω∈R σ (G(i ω)) .
If G ∈ RH u×v ∞ , then, as above, we have G ∞ = sup ω∈R σ (G(i ω)). The paper aims at developing certified symbolic-numeric algorithms for the computation of

G ∞ for G ∈ R(s) u×v satisfying G |i R ∈ RL u×v ∞ .
This problem plays a fundamental role in H ∞ -control theory [START_REF] Doyle | Feedback Control Theory[END_REF][START_REF] Zhou | Robust and Optimal Control[END_REF].

The conjugate G of G ∈ R(s) u×v is defined by G(s) := G T (-s).

The next proposition gives a first characterization of G ∞ .

Proposition 1 ([13]). Let γ > 0, G ∈ R(s) u×v be such that G |i R ∈ RL u×v ∞ and let us consider Φ γ (s) = γ 2 I v -G(s) G(s). Then, γ > G ∞ if and only if γ > σ (G(i ∞)) and det(Φ γ (i ω)) = 0 for all ω ∈ R.
Let n(γ, ω) and d(ω) be two coprime polynomials over R[γ, ω] satisfying:

det(Φ γ (i ω)) = n(γ, ω) d(ω) . (2)
Note that det(Φ γ (s)) is a real function in s 2 and γ 2 , and thus, det(Φ γ (i ω)) is a real function in ω 2 and γ 2 . A consequence of Proposition 1 is the next result.

Proposition 2. Let G ∈ RL u×v ∞
and n ∈ R[γ, ω] be defined by [START_REF] Basu | of algorithms and computation in mathematics[END_REF]. We denote by n ∈ R[γ, ω] the square free part of n. Then, we have:

G ∞ = max π γ V R n, ∂ n ∂ω ∪ V R ( Lc ω (n) ) . Example 1. If G ∈ RL ∞ then, by definition, G ∞ is the supremum of the continuous function ω ∈ P 1 (R) := R ∪ {∞} -→ |G(i ω)|,
and thus, we have

G ∞ = max ω∈P 1 (R) |G(i ω)|, i.e., G ∞ = max {|G(i ∞)|, γ max }
, where:

γ max := max ω∈R |G(i ω)| = max γ ∈ R | ∃ ω ∈ R : γ 2 = |G(i ω)| 2 .
We find again Proposition 2, i.e., γ > G

∞ iff Φ γ (i ω) = γ 2 -|G(i ω)| 2 = 0 for all ω ∈ R and γ > |G(i ∞)|. Using |G(i ω)| 2 = G(-i ω) G(i ω) ∈ R(ω 2
), a computation of G ∞ amounts to first computing the zeros of the numerator of d|G(i ω)| 2 dω , then evaluating |G(i ω)| on these zeros and finally choosing the maximal occurring value, that to say γ, and (iii)

G ∞ = max{|G(i ∞)|, γ}.
More explicitly, if we write G as G(s) = a(s)/b(s), where a, b ∈ R[s], gcd(a, b) = 1, q = deg s (a) ≤ r = deg s (b), and b does not vanish on i R, then G(i ∞) = 0 if q < r (i.e., G is strictly proper) or G(i ∞) = a r /b r if q = r (i.e., G is proper), where a r = Lc s (a) and b r = Lc s (b). Moreover, we have

|G(i ω)| 2 = N (ω)/D(ω), where N (ω) = |a(i ω)| 2 and D(ω) = |b(i ω)| 2 ∈ R[ω 2 ]. Since b(i ω) has not real roots, D(ω) = 0 for all ω ∈ R. Hence, if we note Z := {ω ∈ R | N (ω) D(ω) -N (ω) D (ω) = 0}, then we obtain: G ∞ = max{|G(i ∞)|, γ}, γ := max ω∈Z (N (ω)/D(ω)) 1/2 . Note that if Z ∩ V R ( D (ω) ) = V R ( N (ω) D(ω), D (ω) ) = ∅, then we also have γ = max ω∈Z (N (ω)/D (ω)) 1/2 . For instance, if G(s) = (2 s + 1)/(s + 1), then N (ω) = 4 ω 2 + 1, D(ω) = ω 2 + 1, Z = {0},Z ∩ V R ( D (ω) ) = {0}, γ = (N (0)/D(0)) 1/2 = 1, |G(i ∞)| = 2, and G ∞ = max{2, γ} = 2.
Finally, according to Proposition 2, we have

n(γ, ω) = D(ω) γ 2 -N (ω) and d(ω) = D(ω). Now, using gcd(N, D) = 1, n = n, Lc ω (n) = b 2 r γ 2 if q < r or Lc ω (n) = (b 2 r γ 2 -a 2 r ) if q = r, which yields: V R ( Lc ω (n) ) =    0, if q < r, ± a r b r , if q = r,
and using the fact that D(ω) = 0 for all ω ∈ R, we have

π γ V R n, ∂ n ∂ω = π γ V R ( D(ω) γ 2 -N (ω), D (ω) γ 2 -N (ω) ) = γ ∈ R | γ 2 = N (ω) D(ω) , ω ∈ Z , and thus, G ∞ = max{γ, 0} = γ if q < r and G ∞ = max{γ, a r /b r } if q = r. Corollary 1. Let G ∈ RL u×v ∞ and n ∈ R[γ, ω] be the numerator of det γ 2 I v -G(i ω) G(i ω)
defined by [START_REF] Basu | of algorithms and computation in mathematics[END_REF]. Then, the real γ-projection

π γ (V R ( n )) of V R ( n ) is bounded by G ∞ .
According to Proposition 2, given G ∈ RL u×v ∞ , the problem of computing G ∞ can be reduced to the computation of the maximal γ-projection of the real solutions of the following bivariate polynomial system:

Σ := n(γ, ω), ∂ n(γ, ω) ∂ω . (3)
For studying this problem, we propose three different symbolic-numeric methods -Rational Univariate Representation method, Roots Separation Method and Sturm-Habicht method -and compare them. Without loss of generality, we shall suppose that n is squarefree in Z[γ, ω].

Rational Univariate Representation method

In this section, we briefly state a straightforward algorithm which computes the maximal γ-projection of the real solutions of (3) based on a Rational Univariate Representation Method [START_REF] Rouillier | Solving zero-dimensional systems through the rational univariate representation[END_REF][START_REF] Rouillier | Efficient isolation of polynomial's real roots[END_REF][START_REF] Bouzidi | Separating linear forms and rational univariate representations of bivariate systems[END_REF][START_REF] Bouzidi | Solving bivariate systems using rational univariate representations[END_REF]. This algorithm consists in first computing a rational parametrization (RUR) of the solutions of (3), then isolating the roots of a univariate polynomial p defining the associated field extension, using the intervals obtained to compute isolating boxes for the real solutions of (3) and finally selecting the real solution of (3) with the maximal γ-projection.

If P, Q ∈ Q[x, y] are two coprime polynomials, then the computation of the RUR of V K ( P, Q ), K = R, C, consists in finding s ∈ N such that x+s y separates the K-zeros of {P, Q} and four polynomials p, q, p 0 , q 0 ∈ Q[T ] which define a 1-1 correspondence between V K ( Σ ) and V K ( p ), i.e., the following bijection:

V K ( P, Q ) -→ V K ( p ) (x, y) -→ ξ = x + s y, p 0 (ξ) q(ξ) , p 1 (ξ) q(ξ) ←-ξ
Roughly speaking, using the RUR of V K ( P, Q ), we can transform the study of problems on V K ( P, Q ) into corresponding problems on V K ( p ) [START_REF] Rouillier | Solving zero-dimensional systems through the rational univariate representation[END_REF][START_REF] Bouzidi | Solving bivariate systems using rational univariate representations[END_REF].

For the L ∞ -norm computation, the polynomials of Σ ⊂ Z[γ, ω], defined by (3), are coprime. Hence, to compute G ∞ , we first use the RUR method to obtain isolating boxes for the real solutions (γ, ω) of Σ, choose the maximal γ-projection γ 1 , then compute an isolation box γ 2 for the maximal real root of the univariate polynomial Lc ω (n), and finally compute G ∞ = max{γ 1 , γ 2 }.

Algorithm 1 RUR method

Input: A zero dimensional polynomial system {n, ∂n ∂ω } ⊂ Z[γ, ω]. Output: An isolating interval of max πγ VR n, ∂n ∂ω ∪ VR (Lcω(n)) .
1. Apply the RUR function (Isolate) for solving the polynomial system {n, ∂n ∂ω }. 2. Let γ1 be the maximal γ-projection of the system's real solutions. 3. Let γ2 be the maximal real root of Lcω(n). 4. Return the isolating interval of max{γ1, γ2}.

Remark 1 In step 1 of Algorithm 1, we obtain isolating boxes

[a i , b i ] × [c i , d i ]
of the system's real solutions (γ i , ω i ). To compare the real values in step 2 and then in step 4, we can apply a straightforward strategy consisting in computing the resultant polynomial R γ = Res(n, ∂n ∂ω , ω) and then refining the boxes until each interval [a i , b i ] is included in an isolating interval of R γ . Even with this naive approach, the asymptotic complexity of these operations does not exceed the algorithm's overall worst case bit complexity.

In what follows, O B denotes the bit complexity and ÕB means that logarithmic factors have been omitted. Given two coprime polynomials P, Q ∈ Z[x, y] of degree bounded by d and coefficient bitsize bounded by τ , an algorithm for computing linear separating forms, RUR representations and isolating boxes of the solutions can be obtained in the worst case bit complexity ÕB (d 6 + d 5 τ ) [START_REF] Bouzidi | Solving bivariate systems using rational univariate representations[END_REF].

Let us compute the complexity of Algorithm 1. We first need the next result.

Lemma 1. Let G ∈ RL u×v ∞ , Φ γ (i ω) = γ 2 I v -G(i ω) G(i ω), n ∈ R[γ, ω] be defined by (2), d γ = deg γ (n), d ω = deg ω (n), and τ n the coefficient bitsize of n. Moreover, let α = max{u, v}, N = max 1≤i≤u,1≤j≤v {deg ω (Q i,j )}, where G jk = P jk Q jk denotes the (j, k) th entry of G and P jk , Q jk ∈ R[i ω] are coprime, and τ G the maximal coefficient bitsize of {P jk , Q jk } 1≤j≤u,1≤k≤v
. Then, we have:

d γ = O(α), d ω = O(N α 2 ), τ n = Õ(τ G α 2 ). Proof. Let G jk = P jk /Q jk be the (j, k) th entry of G, where P jk , Q jk ∈ R[i ω] are coprime. Since G ∈ RL u×v
∞ , G jk is a proper rational function, and thus, deg ω (P jk ) ≤ deg ω (Q jk ) ≤ N , which shows that the degrees in ω of the numerators and the denominator of the entries of Φ γ (i ω) are bounded by 2 N α, and thus, d ω is bounded by 2 N α 2 . Similarly, the maximal coefficient bitsize of the entries of Φ γ (i ω) is 2 τ n α, which yields τ n is bounded by 2 τ G α 2 . Finally, d γ is clearly bounded by 2 α. Theorem 1. With the above notations, the complexity of Algorithm 1 for the computation of G ∞ , where G ∈ RL u×v ∞ , is given by:

ÕB (d γ d 3 ω (d 2 γ + d γ d ω + d ω τ n )) = ÕB (α 9 N 4 (α + τ n )).
Proof. According to [START_REF] Bouzidi | Solving bivariate systems using rational univariate representations[END_REF], using the RUR method, the complexity of the resolution of a zero dimensional bivariate polynomial system comes first from the computation of the triangular decomposition of the system after shearing -using a separating linear form γ + s ω -then from the root isolation of the univariate polynomial defining the associated field extension, and finally from the computation of the isolating boxes for the solutions.

In the present case, the degrees in ω and γ are not of the same order. Hence, the results of [START_REF] Bouzidi | Solving bivariate systems using rational univariate representations[END_REF] must be adapted.

First, we determine the size and the degree of the sheared system up to the method used in [START_REF] Bouzidi | Solving bivariate systems using rational univariate representations[END_REF]: the degree with respect to the variable ω is Õ(d γ + d ω ) and d γ with respect to the variable t = γ + s ω. The size of the sheared system is Õ(τ n + d γ ). From [START_REF] Lazard | Bivariate triangular decompositions in the presence of asymptotes[END_REF], the complexity of the computation of a triangular decomposition of a system over Z[x, y] costs ÕB (d

3 x d 3 y + (d 2 x d 3 y + d x d 4 y ) τ )
, where d x and d y are defined by (1) and τ is its maximal coefficient bitsize. Thus, using Lemma 1, we obtain that the complexity of the computation of the triangular decomposition of the sheared system in Z[t, ω] is given by:

ÕB (d 3 γ (d γ + d ω ) 3 + (d 2 γ (d γ + d ω ) 3 + d γ (d γ + d ω ) 4 ) (τ n + d γ )) = ÕB (d γ d 3 ω (d 2 γ + d γ d ω + d ω τ n )).
Moreover, by Lemma 1, we obtain ÕB (α 9 N 4 (α + τ n )). This triangular decomposition yields (RUR) polynomials of degree that sum up to Õ((

d γ + d ω ) d γ ) = Õ(α 3 N ) with coefficients of bitsize Õ((d γ + d ω )(d γ + τ n )) = Õ(α 2 N (α + τ n )).
Finally, according to [START_REF] Bouzidi | Solving bivariate systems using rational univariate representations[END_REF], it is known that the computation of isolating boxes of all the roots of the system can be done in ÕB ((

α 3 N ) 3 +(α 3 N ) 2 (α 2 N (α+τ n ))) = ÕB (α 8 N 3 (α + τ n )) bit operations.

Roots Separation Method

In this section, we localize the maximal γ-projection of the real solutions of the polynomial system Σ by only shearing the system Σ using a special linear separating form [START_REF] Cheng | Root isolation for bivariate polynomial systems with local generic position method[END_REF]. With this linear separating form t = γ + s ω, we obtain:

t 1 = γ 1 + s ω 1 < t 2 = γ 2 + s ω 2 =⇒ γ 1 ≤ γ 2 . Let P, Q ∈ Z[x, y] be coprime and R x = Res(P, Q, y) ∈ Z[x] be their resul- tant. Moreover, let x 1 ≤ . . . ≤ x m be the real roots of R x with isolating intervals [c 1 , d 1 ], . . . , [c m , d m ]
. Moreover, let the real numbers δ, M and s be defined by:

δ < 1 2 min i=1,...,m-1 (x i+1 -x i ), M > max{y | (x, y) ∈ V R ( P, Q )}, 0 < s < δ M .
(4) We can use the general root bounds for zero dimensional systems to estimate δ and M . In fact, M is the measure of the univariate polynomial Res(P, Q, x). Note that the resultant computation can be avoided by using the concept of sleeve functions studied in [START_REF] Cheng | Complete numerical isolation of real roots in zero-dimensional triangular systems[END_REF] and [START_REF] Cheng | Root isolation for bivariate polynomial systems with local generic position method[END_REF]Lemma 3.3].

Let us consider an invertible linear map (a shear) of R 2 to R 2 defined by Ψ s : (x, y) -→ (t, y) = (x + s y, y). Let us also note Ψ s (P ) = P (t -s y, y), Ψ s (Q) = Q(t -s y, y), R t = Res(Ψ s (P ), Ψ s (Q), y) and let t 1 ≤ . . . ≤ t m = t max be the real roots of R t .

To get a 1-1 correspondence between the zeros of {P, Q} and the roots of R t , Lc y (Ψ s (P )) and Lc y (Ψ s (Q)) must not both vanish. It is always possible to choose s such that this condition is satisfied. In what follows, we shall consider this case. For more details for the computation of s up to this condition, see [START_REF] Cheng | Root isolation for bivariate polynomial systems with local generic position method[END_REF]. Remark 2 In Fig. 1, we only draw a part of the plot. In fact, since n ∈ Z[γ 2 , ω 2 ], Σ = 0 is symmetric with respect to the γ and ω axes. Proposition 3. With the above notations, let x m be a real root of R x with an isolating interval [c m , d m ] and t max the maximal real root of R t . If

t max ∈ [c m -δ, d m + δ], then the maximal x-projection of V R ( P, Q ) is equal to x m .
Proof. For each real root x i of R x with an isolating interval [c i , d i ], let us denote by P i,j = (x i , y i,j ) the real solutions of {P, Q} which project onto x i . Then, Ψ s (P i,j ) = (x i + s y i,j , y i,j ), where x i + s y i,j is the first coordinate of a real solution of {Ψ s (P ), Ψ s (Q)}. Using (4), we obtain that different Ψ s (P i,j )

ω γ γ1 γ2 γ3 M δ
Fig. 1: The blue dots represent the real solutions of Σ = 0 and the blue crosses are the complex ones; the red dots are the solutions of Ψ s (Σ) = {Ψ s (n), Ψ s ( ∂n ∂ω )}; the orange dots on the γ-axis are the roots of univariate polynomial R t .

have different first coordinates, and thus, {Ψ s (P, Ψ s (Q)} is in a generic position. Furthermore, we have elif D ≤ 2 , let 1 = D/2, δ = D -1 and compute s up to the required conditions ; 3. Expand {Ψs(P ), Ψs(Q)} and compute Rt = Res(Ψs(P ), Ψs(Q), y). 4. Isolate Rt up to an accuracy less than 1 and set [pt, qt] to be the isolating interval of its maximal real root tmax. 5. for j from 1 to m do:

|x i + s y i,j -x i | = |s y i,j | < ( δ M ) M = δ. Consequently Ψ s (P i,j ) ∈ I i = [x i -δ, x i + δ] × [-M, M ]. In addition, since δ < 1 2 (x i+1 -x i ), I i are disjoint for different i. Hence, a real solution (x, y) of {P, Q} is mapped to (η, y), where η ∈ [x -δ, x + δ]. Now, since x i ∈ [c i , d i ], the real roots of R t associated with x i are in the interval [c i -δ, d i + δ].

Algorithm 2 Roots Separation method

if [pt, qt] ⊂ [cj -D, dj + D], then X1 = xj. 6. Let X2 be the maximal real root of Lcy(P ). 7. Return the isolating interval of max{X1, X2}. Proof. The proof is a direct consequence of the proof of [4, Lemma 7] by taking into account the bitsize τ s of s.

Theorem 2. We consider a zero dimensional system {P, Q} ⊂ Z[x, y], where d x = max(deg x (P ), deg x (Q)), d y = max(deg y (P ), deg y (Q)) and τ the maximal coefficients bitsize of the polynomials. We can compute an isolating interval for the maximal x-projection of the real solutions of {P, Q} in ÕB (d 4

x d 5 y τ ) bit operations, using Algorithm 2.

Proof. For each root t j of R t defined in Algorithm 2 with an isolating interval [p j , q j ], there exists a unique i ∈ {1, . . . , m} such that [p j , q j ] ⊂ [c i -D, d i + D]: we know that there exists a unique [c i , d i ] such that t j ∈ [c i -δ, d i + δ]. From step 4, q j -p j < 1 and D = δ + 1 . Hence, q j < t j + 1 < d i +δ + 1 < d i +D. And similarly p j > c i -D. Consequently, based on Proposition 3, Algorithm 2 outputs an isolating interval for the maximal x-projection of the real solutions of {P, Q}. and M = 2 O(dx τ ) [START_REF] Bouzidi | Solving bivariate systems using rational univariate representations[END_REF]. The bitsize of s is then equal to Õ(d x d 2 y τ ). Consequently, the coefficient bitsize of the sheared system is Õ(d 2

x d 2 y τ ), and the worst case bit complexity of step 5 is Õ(d 3

x d 2 y (d x + d y ) 3 τ ), as computed in Lemma 2. In step 6, we isolate the resultant of the sheared system. Considering the size and degree of the sheared polynomials computed using Lemma 2, the size and degree of the resultant of the sheared system are Õ(d 2

x d 3 y τ ) and Õ(d x (d x + d y )) respectively. Then, knowing the complexity of the isolation mentioned in [START_REF] Bouzidi | Solving bivariate systems using rational univariate representations[END_REF], we can say that the worst case bit complexity in this line is equal to ÕB ((d

x d y ) 3 + (d x d 2 y ) 2 (d 2 x d 3 y τ )) = ÕB (d 4 x d 5 y τ ).
Finally, in the step 7, we simply compare two rational numbers. The maximal coefficients bitsize of these rationals is in Õ(d 3

x d 3 y (d x + d y ) τ ) and the computation in this step is done in Õ(d 3 x d 3 y (d x + d y ) τ ) bit operations. Hence, the overall bit complexity is given by ÕB (d 4

x d 5 y τ ).

Corollary 2. With the notations of Lemma 1, the worst case bit complexity for the computation of G ∞ with the separation method (Algorithm 2) is given by ÕB (α 14 N 5 τ n ).

From this section, we can conclude that trying to concentrate only on the solution with the maximal γ-projection, after putting the system in a generic position, costs much more than computing isolating boxes for all the real solutions, due to the large size of the separating bound that we must use. Hence, in the next section, using another strategy than shearing the system, we shall try to find the maximal γ-projection of the polynomial solutions without computing isolating boxes for all the real solutions.

Sturm-Habicht method

In this section, as in Section 4, we shall concentrate only on the maximal γprojection γ of the real solutions of the polynomial system. But instead of shearing the system, we shall verify the existence of a real root of the polynomial system over γ by studying the sign variation of the leading coefficients of subresulant polynomials over γ. Hence, before explaining the third proposed method, we first state again a few standard preliminaries on subresultants and Sturm-Habicht sequences.

We denote by K the unique factorization domain Q[x] and we consider P, Q ∈ K [y], where p = deg y (P ) and q = deg y (Q). We assume that p ≥ q. For 0 ≤ i ≤ min(q, p -1), the i th subresultant polynomial of P and Q is denoted by Sres y,i (P, p, Q, q). When there is no ambiguity on the degrees of the polynomials P and Q, we simply denote it by Sres y,i (P, Q). It has degree at most i in y and the coefficient of y i is denoted by sres y,i (P, Q). It is called the i th principal subresultant coefficient. We recall that sres y,i (P, Q) = 0 implies that Sres y,i (P, Q) vanishes identically. Note that Sres y,0 (P, Q) = sres y,0 (P, Q) is the resultant of P and Q with respect to y, also denoted by Res(P, Q, y). The greatest common divisor gcd(P, Q) of the polynomials P and Q (uniquely defined up to units of K ) is the first non-zero subresultant polynomial Sres y,i (P, Q) for increasing i.

Letting v = p + q -1 and δ k = (-1)

k(k+1) 2
for k ∈ Z ≥0 , the j th polynomial in the Sturm-Habicht sequence associated to (P, Q), denoted by StHa j (P, Q), is then defined by δ v-j Sres y,j (P, v + 1, P Q, v), where P denotes the derivative of P with respect to y. The principal j th Sturm-Habicht coefficient is denoted by stha j (P, Q) for j = 0, . . . , v + 1. We also denote by SignVar the function which maps {sign(stha j (P, 1))} j=0,...,v+1 to the number of real roots of P . For more details on the function SignVar, see [12, Definition 4.1, Theorem 4.1].

As stated above, we aim at computing:

γ = max π γ V R n, ∂ n ∂ω ∪ V R ( Lc ω (n) ) .
Hence, γ is either the maximal real root of Lc ω (n) or an algebraic value over which gcd(n(γ, ω), ∂n ∂ω (γ, ω)) ∈ R[ω] has at least one real root. We recall that gcd(n(γ, ω), ∂n ∂ω (γ, ω)) is proportional to the first subresultant polynomial Sres ω,i (n, ∂n ∂ω ) (for i increasing) that does not identically vanish for γ = γ. If γ is not a real root of Lc ω (n), then we can compute the Sturm-Habicht sequence of the univariate polynomial n(γ, ω) ∈ R[ω] to check the existence of a real root for gcd(n(γ, ω), ∂n ∂ω (γ, ω)). In what follows, we shall need the next result.

Lemma 3. Let P ∈ Z[x, y] and x be a root of Res(P, ∂P ∂y , y). Moreover, let G = gcd P (x, y), ∂P ∂y (x, y) ∈ R[y]. If the x-projection of the points of P is bounded by x, then we have V R ( P (x, y)

) = V R ( G(x, y) ). Proof. If V R ( G(x, y) )
V R ( P (x, y) ), then there exists y 0 ∈ R such that P (x, y 0 ) = 0 and G(x, y 0 ) = 0. This is equivalent to saying that P (x, y 0 ) = 0 and ∂P ∂y (x, y 0 ) = 0. Hence, based on the theorem of implicit functions, there exists a real function ϕ of class C p (p > 0), defined on an open interval V ⊂ R, containing x, and an open neighborhood Ω of (x, y 0 ) in R 2 such that for all (x, y) in

R 2 , {(x, y) ∈ Ω | P (x, y) = 0} is equivalent to {x ∈ V | y = ϕ(x)}.
This cannot be true since the x-projection of the points of the curve P = 0 is bounded by x, and thus, an open interval containing x, such as V , does not exist. Consequently, we obtain V R ( P (x, y) ) = V R ( G(x, y) ).

Algorithm 3 Sturm-Habicht method

Input: A bivariate polynomial P ∈ Z[x, y] such that P = 0 is bounded in the xdirection. Output: Isolating interval of max πx VR P, ∂P ∂y ∪ VR ( Lcy(P ) ) . Lemma 4. Let P ∈ Z[x, y], d x = deg x (P ), d y = deg y (P ) and τ be the maximal coefficients bitsize of P . Let {StHa j (P (x, y), 1)} j=0,...,dy be the Sturm-Habicht sequence and x j a real root of sres y,0 P, ∂P ∂y . Then, {sign(stha k (x j ))} k=dy,...,1 can be computed in ÕB (d 2

x d 4 y (d x + τ )) bit operations.

Proof. We denote by sres 0 (resp., sres i ) sres y,0 P, ∂P ∂y (resp., sres y,i P, ∂P ∂y ), where sres i ∈ Z[x]. We first recall that stha i (x j ) = δ dy-1-i sres i (x j ). Based on [2, Proposition 8.46], sres i is of degree d x d y and of coefficients bitsize d y τ . Thus, the square free part of sres 0 is of coefficients bitsize O(d y (d x + τ )) and, based on [4, Lemma 5], can be computed in ÕB (d 2

x d 3 y τ ). By following the proof of [START_REF] Strzebonski | Univariate real root isolation in an extension field and applications[END_REF]Proposition 6], the overall cost for obtaining the list {sign(stha dy (x j )), . . . , sign(stha 1 (x j ))} is ÕB (d 2

x d 4 y (d x + τ )). Proof. The maximal x-projection of the real solutions of {P, ∂P ∂y } is the maximal real root of sres 0 P, ∂P ∂y , say x m , such that gcd P (x m , y), ∂P ∂y (x m , y) has at least one real root. If the x-projection of the points of P is bounded by x m , then, by Lemma 3, the real roots of gcd P (x m , y), ∂P ∂y (x m , y) are the real roots of P (x m , y). Consequently, we can compute an isolating interval of x m using Algorithm 3. According to [ In Algorithm 4, we suppose that there are no real isolated points, and thus, we replace the computation of signs of polynomials at real algebraic numbers by signs of polynomials at rational numbers. Syntactically, these are small modifications but the effect on the computations is consequent in practice, as well as in theory, since the evaluation of signs of polynomials at real algebraic numbers carries the theoretical worst case complexity of Algorithm 3. We first compute V R (Lc ω (n)) = {0}. Then, applying the RUR method, we obtain the following rational univariate representation:

       p = (t 2 + t -1)(t 2 -t -1), γ = 3 t 2 -2 t (2 t 2 -3) , ω = 0.
Thus, the system's real solutions (γ, ω) are:

- √ 5 2 - 1 2 , 0 , - √ 5 2 + 1 2 , 0 , √ 5 2 - 1 2 , 0 , √ 5 2 + 1 2 , 0 .
Thus, we simply pick their maximal γ-projection to obtain √ 5 2 + 1 2 , which yields:

G ∞ = max 0, √ 5 2 + 1 2 = √ 5 2 + 1 2 .
Following the second approach, which consists in directly focusing on the maximal γ-projection of the system's real solutions, we first compute Res(n, ∂ n ∂ω , ω) and denote by R = γ (γ 2 + γ -1) (γ 2 -γ -1) ∈ Z[γ] its square free part. Then, the maximal real root of R has the following isolating interval: Following the Root Separation method, we obtain:

       s = 12060328540887 281474976710656 , δ = 43490275647441 140737488355328 . We have R t = Res(Ψ s (n), Ψ s ( ∂ n ∂ω ), ω) = α (t 2 + t -1)(t 2 -t -1)
, where α is a rational of size about 3000 bits. We denote by t max the maximal real root of R t . An isolating interval of t max is then given by:

Experiments

The three proposed methods can be implemented in a few lines of Maple but we then have to use implementations at different levels that do not give valuable information about the intrinsic efficiency. For instance, the RUR is implemented in C but for general zero dimensional polynomial systems: a variant for bivariate polynomials, the one used for the complexity analysis, is not part of Maple and is much more efficient for bivariate systems.

In order to have fair comparisons, we extract the dominating operations and compare them using exactly the same implementations. Namely, resultant computations of sheared/non sheared systems and Root Isolation carry the largest percentage of the computation time. For instance, Algorithm 4 saves time on the resultant computation since it does not perform any shear while it loses time on the root isolation.

For the three methods, the principle subresultant sequence is computed using the routine SubResultantChain of the Maple package RegularChain.

Isolating the real roots of univariate polynomials is another common basic block shared between the three algorithms for which we use Isolate provided by the Maple routine package RootFinding.

In left table of Table 1, we list the main steps of the three algorithms. The check marks mean that the step makes part of the method and the double check marks indicate that this step is the bottleneck of the method. Note that Res1 stands for the resultant of the original system and Res2 for the resultant of the sheared system. Keep in mind that the shear done in Hinf RUR is different than the one done in Hinf Sep. Finally, Iso means Isolate.

In the right by random proper rational functions of degree N (degree of the denominators) 3 . It corresponds to a fixed input coefficient bitsize τ = 2, i.e., the rational functions involved in the entries of the matrices have coefficients of magnitude O(2 τ ). We finally mention that with these experiments, our goal is not to illustrate the theoretical complexity, but, on the contrary, to show that on practical examples, the results in practice are different than in theory. In theory, the RUR algorithm might asymptotically be the fastest while in practice the Sturm method performs better.

Conclusion

In this paper, we have presented three different algorithms for the computation of the L ∞ -norm of the transfer matrix of a finite-dimensional linear control system. By reformulating this problem as the search for the maximal projection of the real solutions of a zero dimensional polynomial system, we have used existing methods such as the rational univariate representation (RUR method). As for the second algorithm, we have only used a special separating linear transformation to shear the polynomial system and put it in a generic position. The last method (Sturm-Habicht method) was based on verifying the existence of a real root for a univariate polynomial having real coefficients.

The complexity analysis has showed that the RUR method has the best theoretical efficiency in comparison with other algorithms. Practically, as we can notice in the tables given in Section 6, the practical efficiency is nearly matching with the theoretical efficiency but with a slight advantage for the Sturm-Habicht method probably due to the fact that it is the most adaptive one.

  where K is a field containing k. Finally, let C + := {s ∈ C | Re(s) > 0} be the open right-half plane of C.

  Input: A zero dimensional system {P, Q} ⊂ Z[x, y], where Q = ∂P ∂y . Output: An isolating interval of max {πx (VR ( P, Q )) ∪ VR ( Lcy(P ) )} 1. Isolate Rx = Res(P, Q, y) up to an accuracy and let RI := {[c1, d1], . . . , [cm, dm]} be the isolating intervals of the real roots {x1, . . . , xm} of Rx. 2. Compute M and D = 1 2 mini=1,...,m-1 |ci+1 -di|: if D > 2 , let 1 = , δ = D -1 and compute s up to the required conditions.

Lemma 2 .

 2 Let P ∈ Z[x, y], d x = deg x (P ), d y = deg y (P ) and τ be the maximal coefficient bitsize of P . The sheared polynomial P (t -s y, y) satisfies deg y (P (t -s y, y)) = d x + d y , deg t (P (t -s y, y)) = d x , and it can be expanded in ÕB (d x d 2 y (τ + d x (1 + τ s )). The maximal bitsize of the coefficients of P (t -s y, y) is equal to Õ(τ + d x (1 + τ s )), where τ s denotes the bitsize of s.

  As for the complexity, step 1 has worst-case bit complexity ÕB (d x d 3 y τ ) based on [2, Proposition 8.46]. Step 2 is of worst case bit complexity ÕB (d 3 + d 2 τ ), where d = deg(R x ) and the coefficient size of R x is equal to τ [3, Lemma 54]. From [2, Proposition 8.46], d = O(d x d y ) and τ = Õ(d y τ ). Consequently, step 2 is of worst case bit complexity ÕB (d 2 x d 3 y (d x + τ ). In steps 3 and 4, we get δ = 2 -Õ(dx d 2 y τ )

1 .

 1 Compute {Sresj(P, ∂P ∂y )} j=0,...,deg y (P ) . 2. Compute x1 < . . . < xm the real roots of sres0. 3. for i from 1 to m do: if x1-i+m ∈ VR( Lcy(P ) ) then return the isolating interval of x1-i+m; -elif SignVar({sign(stha dy (x1-i+m)), . . . , sign(stha1(x1-i+m))}) > 0, then return the isolating interval of x1-i+m. 4. end if end do.

Theorem 3 .

 3 Let P ∈ Z[x, y] be such that d x = deg x (P ), d y = deg y (P ) and of maximal coefficients bitsize τ . Then, we can compute an isolating interval of the maximal x-projection of the real solutions of {P, ∂P ∂y } (Algorithm 3) in ÕB (d 2 x d 4 y (d x + τ )) bit operations in the worst case.

Corollary 3 .

 3 Based on Theorem 3, G ∞ can be computed by Sturm-Habicht method in the worst case bit complexity ÕB (α 10 N 4 (α + τ n )).

Theorem 4 .Fig. 2 :

 42 Fig. 2: Plot of n(γ, ω) = 0, where ω/γ is in the horizontal/vertical axis.

  In this case, [c, d] ⊂ [a -δ, b + δ], which shows that G ∞ is equal to the maximal real root of R of isolating interval [a, b]. Following the Sturm-Habicht method, we have to check the existence of a real root for the univariate polynomial n([a, b], ω). To do that, we first compute L = [sres ω,i (n, ∂ n ∂ω )] i=1,...,deg ω (n)=4 = [5 γ 14 (2 γ 2 -3), 2 γ 10 (2 γ 2 -3), 4 γ 4 , γ 4 ]. We then compute the list of signs of the elements of L over [a, b]. We obtain the list L s = [-, -, +, +]. Then, SignVar(L s ) = 1 and we conclude that n([a, b], ω) admits one real root. Hence, G ∞ is equal to the maximal real root of R of isolating interval [a, b].

  2, Proposition 8.46], we can compute the set of principal subresultants in Õ(d x d 3 y τ ) bit operations and each subresultant polynomial is of degree O(d x d y ) and of coefficient bit size Õ(d y τ ). Thus, step 2, which performs real root isolation of sres 0 , is of complexity Õ((d x , d y ) 3 + (d x d y ) 2 d y τ ) [3, Lemma 54]. Using Lemma 4, step 3 can be done in ÕB (d 2 x d 4 y (d x +τ )) operations since its first step is of complexity ÕB (d 3 x + d 2 x τ ). Hence, the overall complexity of this algorithm is ÕB (d 2 x d 4 y (d x + τ )). Considering the notations of Lemma 1, the following result is an immediate consequence of Corollary 1 and Theorem 3.

Table 1 :

 1 table of Table 1, we report the average running time in CPU seconds of the marked steps listed in the table on the left of Table1for the three proposed algorithms run on square matrices of size α, with entries given Left: main steps considered in the implementation of the proposed method. Right: timings for L ∞ -norm for random matrices with τ G = 2.

			α	N Hinf RUR Hinf Sep Hinf Sres
				2	0.2	3	0.2
	Hinf RUR Hinf Srep Hinf Sres	Res1 + Iso Res2 + Iso List of signs	2	3 4 5 6 7	0.5 2.5 10 37 50	7 25 83 96 186	0.5 2 6 10 47.5
				8	133.5	353	59
				9	236	394	130

The experiments were conducted on Intel(R) Core(TM) i7-7500U CPU @ 2.70GHz 2.90 GHz, Installed RAM 8.00 GB under a Windows platform

Algorithm 4 Sturm-Habicht method -equidimensional

Input: A bivariate polynomial P ∈ Z[x, y] such that the curve P = 0 is bounded in the x-direction and has not real isolated singular points. Output: An isolating interval of max πx VR P, ∂P ∂y ∪ VR (Lcy(P )) 1. Compute {StHaj(P, 1)} j=0,...,deg y (P ) . 2. Let x1 < . . . < xm be the roots of sres0. 3. for i from 1 to m do: From the above complexity analysis, we can conclude that RUR method and the Sturm-Habicht method have comparable theoretical complexities since, in our case, we have α N .

6 Experiments

Practical example

We consider the following transfer matrix:

. We study the real solutions of the polynomial system Σ = {n, ∂ n ∂ω }, where: n(γ, ω) = γ 4 ω 4 + γ 2 (2 γ 2 -3) ω 2 + (γ 2 + γ -1)(γ 2 -γ -1).


Alexis Thibault
email: alexis.thibault@inria.fr

Juliette Chabassier

Dissipative time-domain 1D model for viscothermal acoustic propagation in wind instruments

Dissipative time-domain one-dimensional model for viscothermal acoustic propagation in wind instruments.

I. INTRODUCTION

Time-domain simulations of wind musical instruments require the choice of a balance between efficiency and accuracy. Real-time synthesis techniques, such as digital waveguide synthesis with lumped wall loss filter (Abel et al., 2003), [START_REF] Mignot | Digital Waveguide Modeling for Wind Instruments: Building a State-Space Representation Based on the Webster-Lokshin Model[END_REF], modal decomposition [START_REF] Silva | MoReeSC: a framework for the simulation and analysis of sound production in reed and brass instru-ments[END_REF], and Finite Difference -Time Domain [START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF], favor their performance objective at the expense of model and discretization errors. By contrast, computer aided instrument prototyping requires a fine assessment of the physical phenomena occurring in the instrument, including dissipation and dispersion effects on propagative waves inside air filled pipes, caused by viscous and thermal boundary layers. The Navier-Stokes equations of fluid dynamics indeed do account for all these effects, but a practical understanding of the macroscopic phenomena from direct numerical simulation requires a tremendous computational effort. In the different characteristic regimes found in musical instruments, several approximate 3D or 1D models exhibit good accuracy and performance ; they are briefly reviewed thereafter. More details about the qualitative and quantitative fulfilment of their underlying assumptions can be found in [START_REF] Thibault | Viscothermal models for wind musical instruments[END_REF]. Most of them, including the prevalent model initially derived by (Zwikker and Kosten, 1949) (ZK), are naturally expressed in the frequency domain, which hampers their use for synthesising the sound of musical instruments by coupling the pipe to a nonlinear embouchure and using time-varying commands. This work focuses on the time domain formulation of a stable and accurate 1D model, and its numerical implementation based on a variational approach, aiming at computer aided instrument prototyping.

A. Viscous and thermal boundary layers in 3D

The Linearized Navier Stokes (LNS) equations, also called thermoviscoacoustic equations, describe the evolution of the acoustic density, pressure, temperature and velocity, which are small perturbations of a background medium assumed steady with no mean flow. They feature acoustic wave propagation, as well as viscous and thermal diffusion [START_REF] Tijdeman | On the propagation of sound waves in cylindrical tubes[END_REF].

For a domain filled with air bounded with rigid, non-porous isothermal walls, and when the acoustic wavelength is greater than the viscous and thermal characteristic lengths, the acoustic velocity and temperature are functions of the acoustic pressure and two potentials [START_REF] Kampinga | An efficient finite element model for viscothermal acoustics[END_REF], which are solution to non homogeneous heat equations with Dirich-let boundary conditions. The potentials are close to one in the bulk of the domain, and drop to zero near the walls. Then the acoustic pressure satisfies a Helmholtz equation with variable coefficients that explicitly depend on these potentials. This leads to the emergence of viscous and thermal boundary layers affecting the acoustic pressure.

When these boundary layers are thin with respect to the radius of the pipe, an asymptotic analysis justifies to model them as an effective wall impedance. In the interior domain, the acoustic pressure is then solution to a uniform 3D Helmholtz equation [START_REF] Cremer | On the acoustic boundary layer outside a rigid wall[END_REF], [START_REF] Bruneau | General formulation of the dispersion equation in bounded visco-thermal fluid, and application to some simple geometries[END_REF], [START_REF] Jith | Boundary layer impedance model to analyse the viscothermal acousto-elastic interactions in centrifugal compressors[END_REF], [START_REF] Berggren | Acoustic boundary layers as boundary conditions[END_REF].

B. Model reduction to 1D

In the case of a cylindrical pipe, [START_REF] Kirchhoff | Ueber den einfluss der wärmeleitung in einem gase auf die schallbewegung[END_REF] derives an implicit analytical dispersion relation to the LNS equations that must be solved iteratively, which does not offer a practical modeling framework. The reduced 1D model derived from an approximation of LNS by (Zwikker and Kosten, 1949) has been shown to be valid for audible frequencies and assuming the reduced frequency k = ωR/c is small, where ω is the angular frequency, R the radius of the pipe and c the celerity of the wave [START_REF] Tijdeman | On the propagation of sound waves in cylindrical tubes[END_REF]. More precisely, the average pressure along a pipe section p(x, ω), and the air flow v(x, ω), along the longitudinal spatial variable x ∈ [0, L], are solution to the so-called Zwikker and Kosten (ZK) model:

(ZK)        dp dx + jω ρ S 1 1 -F v (ω) v = 0, ( 1a
)
dv dx + jω S ρc 2 [1 + (γ -1)F θ (ω)] p = 0, ( 1b
)
where S is the pipe section, L its length, ρ the static density, γ the ideal gas constant, and coefficients F v and F θ describe the contribution of respectively viscous and thermal dissipation, as:

F v (ω) = φ R - jω ρ µ , F θ (ω) = φ R - jω ρ C P κ , with φ(α) = 2J 1 (α) αJ 0 (α) , ( 2
)
where µ is the gas viscosity, κ the thermal conductivity, C p the specific heat with constant pressure (see Table I), and where J 0 and J 1 are zeroth-and first-order Bessel functions.

Although it has been derived in the context of cylindrical geometries, model (ZK) has been used intensively for varying geometries, namely for S depending on the longitudinal variable

x. The classical horn equations describing plane wave propagation in an axisymmetric lossless pipe can be retrieved from an asymptotic analysis relying on Euler's equations in a pipe with varying section [START_REF] Rienstra | Webster's horn equation revisited[END_REF]. Model (ZK) can be seen as a perturbation of these horn equations, and has been employed for dissipative pipes with varying section for instance in [START_REF] Chaigne | Acoustics of musical instruments[END_REF], [START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF], (Tournemenne Viscosity: µ = 1.708 × 10 -5 (1 + 0.0029 t) kg m -1 s -1

Thermal conductivity: κ = 5.77 × 10 -3 (1 + 0.0033 t) Cal/(ms and [START_REF] Tournemenne | A comparison of a one-dimensional finite element method and the transfer matrix method for the computation of wind music instrument impedance[END_REF] in the harmonic regime. Curvature of the wave fronts can occur in varying geometries and especially in the instrument bell, which can be modeled by an equation similar to (1) [START_REF] Hélie | One-dimensional acoustic models of horns and comparison with measurements[END_REF]. The present work focuses on viscothermal effects, and will neglect the curvature effects which may be included in future work.

The highly nonlinear dependency of F v and F θ with ω induces a nonlocal formulation in the time domain leading to mathematical and numerical intricacy, which motivates the present work.

C. Model approximations and time domain representations

Model (ZK) has been approximated for different frequency regimes and/or pipe sizes, especially when the Stokes number s = R ρω/µ is large [START_REF] Tijdeman | On the propagation of sound waves in cylindrical tubes[END_REF], [START_REF] Keefe | Acoustical wave propagation in cylindrical ducts: Transmission line parameter approximations for isothermal and nonisothermal boundary conditions[END_REF], [START_REF] Stinson | The propagation of plane sound waves in narrow and wide circular tubes, and generalization to uniform tubes of arbitrary cross-sectional shape[END_REF], [START_REF] Scheichl | On the calculation of the transmission line parameters for long tubes using the method of multiple scales[END_REF]. Terms in √ jω arise in the derived equations, leading in the time domain to fractional derivatives. Note that similar terms arise in the Webster-Lokshin 1D model which models the acoustic pressure close to the boundary layers using Cremer 3D effective wall impedance [START_REF] Hélie | One-dimensional acoustic models of horns and comparison with measurements[END_REF]. These terms can be treated numerically with approximations of diffusive representations [START_REF] Hélie | Diffusive representations for the analysis and simulation of flared acoustic pipes with visco-thermal losses[END_REF], [START_REF] Berjamin | Time-domain numerical modeling of brass instruments including nonlinear wave propagation, viscothermal losses, and lips vibration[END_REF]. Other approaches are based on direct diffusive representations of model (ZK) [START_REF] Thompson | Analog model for thermoviscous propagation in a cylindrical tube[END_REF], [START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF], [START_REF] Schmutzhard | Parameter optimisation of a viscothermal time-domain model for wind instruments[END_REF].

D. Contribution and outline

One important aspect of space and time discretization of a model is to control the trade-off between accuracy and efficiency, in all targeted applications and configurations.

In the context of computer aided instrument prototyping, our purpose is to ensure numerical stability and quantify accuracy with respect to model (ZK). Energy-based methods have proven especially efficient to discretize wave equations in time [START_REF] Van Der Schaft | Port-hamiltonian systems: an introductory survey[END_REF], [START_REF] Bilbao | Direct simulation of reed wind instruments[END_REF], [START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF], [START_REF] Hélie | Self-oscillations of a vocal apparatus: a port-hamiltonian formulation[END_REF], [START_REF] Chabassier | Construction and analysis of fourth order, energy consistent, family of explicit time discretizations for dissipative linear wave equations[END_REF] [START_REF] Chatziioannou | Investigating clarinet articulation using a physical model and an artificial blowing machine[END_REF]. The model of [START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF], which satisfies an energy balance identity, is particularly suitable for physics-based sound synthesis. However, synthetizing sounds of a specific instrument with this model requires to run an optimization algorithm for every different value taken by the pipe radius and temperature before running time iterations. The present work makes this optimization step geometry-independent, by re-writing the model of [START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF] with a new expression of the coefficients. Optimized constants are given in Table II and are usable directly in synthesis algorithms. The new coefficients of the model explicitly depend on the non-constant geometrical and physical parameters, therefore no optimization must be implemented nor launched before time iterations.

The model, presented in Sec. II along with its energy balance, is suitable for variational approximation, therefore naturally compatible for coupling with other energy-based models.

Its mathematical derivation is developed in Sec. III and a practical implementation method is proposed in Sec. IV along with a space and time discretization which guarantees a discrete energy identity and an algorithmic strategy for explicit update of the unknowns. A numerical scheme using 1D mixed spectral finite elements is proposed, allowing high-order space discretization leading to arbitrary accuracy. A comparison with other models is proposed in Sec. V.

Boundary conditions model the way waves are introduced at one boundary, or how they radiate in exterior air at the bell [START_REF] Rabiner | Digital processing of speech signals[END_REF], [START_REF] Dalmont | Radiation impedance of tubes with different flanges: numerical and experimental investigations[END_REF], [START_REF] Silva | Approximation formulae for the acoustic radiation impedance of a cylindrical pipe[END_REF], [START_REF] Monteghetti | Energy analysis and discretization of nonlinear impedance boundary conditions for the time-domain linearized euler equations[END_REF]. Their energy-based time domain formulation can be intricate, and since it is not the scope of the present work, they are considered elementary (closed pipe, induced flow, or open pipe), and abusively omitted when non necessary.

II. MAIN RESULT

The acoustic pressure averaged on a pipe section and the volume flow can be modeled by the solutions p(x, t) and v(x, t) to the following system for x ∈]0, L[, t > 0, where N is an integer and v i (x, t), p 0 (x, t), p i (x, t) are 2N + 1 auxiliary variables:

                                                         ρ S ∂v ∂t + R 0 v + N i=1 R i (v -v i ) + ∂p ∂x = 0, (3a)
S ρc 2 ∂p ∂t + G 0 (p -p 0 ) + N i=1 G i (p -p 0 -p i ) + ∂v ∂x = 0, (3b)
L i dv i dt = R i (v -v i ), for 1 ≤ i ≤ N, (3c)
C 0 dp 0 dt = G 0 (p -p 0 ) + N i=1 G i (p -p 0 -p i ), (3d)
C i dp i dt = G i (p -p 0 -p i ) for 1 ≤ i ≤ N. (3e)
The coefficients of this system are defined as

                         R 0 (x) = πµ S(x) 2 a 0 , (4a)
L i (x) = ρ S(x) a i , R i (x) = πµ S(x) 2 a i b i , (4b)
C 0 (x) = S(x)(γ -1) ρc 2 , G 0 = πκ(γ -1) ρ 2 c 2 C P a 0 , (4c)
C i (x) = S(x)(γ -1) ρc 2 a i , G i = πκ(γ -1) ρ 2 c 2 C P a i b i , (4d)
where the coefficients a i and b i are dimensionless constants obtained from an optimization procedure described in Section III, see Table II for N ∈ {2, 4, 8} and supplementary material SuppPub1.txt up to N=16 1 . We expect that increasing the value of N will make the solution closer to the original model (1), although it will also increase the number of auxiliary unknowns. The form of (3) is the same as the system used in [START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF], but the values of the coefficients change, because the variable v represents here the acoustic flow. To conform to their notations, our formulas for R i and L i should be multiplied by S, and C i and G i should be divided by S.

This model is shown in Section III to satisfy the following energy balance identity

d dt E = -Q -p(L)v(L) + p(0)v(0), (5)
E = 1 2 ˆΩ ρ S v 2 + S ρc 2 p 2 + C 0 p 2 0 + N i=1 C i p 2 i + N i=1 L i v 2 i , (6)
Q = ˆΩ R 0 v 2 + N i=1 R i (v -v i ) 2 + G 0 (p -p 0 ) 2 + N i=1 G i (p -p 0 -p i ) 2 ≥ 0 (7)
The model hence dissipates energy through the negative term -Q, and can exchange work with other systems via its boundaries x = 0 and x = L.

III. MATHEMATICAL DERIVATION

System (1) can be written as 023 15 × 10 -1 1.031 48 × 10 -3 2 6.452 52 × 10 -3 4.096 97 × 10 -6 Coefficients for N = 2 (17 it.) where the characteristic viscous and thermal times are given by

               dp dx + ρ S jω + G(τ v ω) τ v v = 0, ( 8a
)
dv dx + S ρc 2     jω + γ -1 1 jω + τ θ G(τ θ ω)     p = 0, (8b) i a i b i 0 8 1 1.
i a i b i 0 8 1 2.101 57 × 10 -1 1.046 29 × 10 -2
τ v (x) := R(x) 2 ρ µ and τ θ (x) := R(x) 2 ρ C P κ , (9)
and the loss coefficients depend on

G(ζ) := jζ φ( √ -jζ) 1 -φ( √ -jζ) , ∀ ζ ∈ R, ( 10
)
where φ is given by ( 2). [START_REF] Hélie | Diffusive representations for the analysis and simulation of flared acoustic pipes with visco-thermal losses[END_REF].

Instead, in the same spirit as [START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF], function G is approximated in the whole range of interest with a family of functions G N of the following form where (a i , b i ) i=1...N are real coefficients, that must be positive in order to ensure the passivity of the resulting time domain model. In contrast to what is used in [START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF][START_REF] Schmutzhard | Parameter optimisation of a viscothermal time-domain model for wind instruments[END_REF], given a number of oscillators N , a single optimization must be performed to approximate G once and for all, rather than one optimization for every possible radius.

G N (ζ) = a 0 + N i=1 a i jζ b i jζ + 1 , (11)

A. Optimization of the model coefficients

The objective function is chosen to be The positivity of (a i , b i ) i=1...N is enforced by expressing them as exp(a i ), exp(b i ) . This reparametrization warrants the use of unconstrained optimization algorithms, and is suitable to control the model's behavior on frequencies spanning several orders of magnitude. The BFGS algorithm [START_REF] Nocedal | Numerical optimization[END_REF] is then used to find a minimizer of function II gives values of coefficients (a i , b i ) that minimize this 2 error, for different values of N . Note that these coefficients do not depend on the geometry of the instrument, and can be readily used as given. The choice of N is a trade-off between the better precision enabled by using more additional variables, and the higher computing power required for simulation (see Sec. IV D).

E = M k=1 G N (ζ k ) G(ζ k ) -1 2 , ( 12
)
E (a i , b i ) i=1...N . Table

B. Resulting family of models in the frequency domain

Replacing G with G N yields the following family of approximate models:

                     dp dx + ρ S jω + a 0 τ v + N i=1 a i jω b i τ v jω + 1 v = 0, dv dx + S ρc 2      jω + γ -1 1 jω + 1 a 0 τ θ + N i=1 a i jω b i τ θ jω+1      p = 0,
These equations can be written and represented using an equivalent electronic circuit with Foster structure at each abscissa x, as is done in [START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF], by defining the coefficients R i , L i , C i and G i as in equations (4). Such a representation is useful, not only to obtain an explicit choice of state-space representation for the rational functions in (13), but also for the derivation of the energy balance.

C. Resulting family of models in the time domain

This family of models can be written in the time domain as follows

(M N )                                                    ρ S ∂v ∂t + ∂p ∂x + ∆ = 0, (14a)
S ρc 2 ∂p ∂t + ∂v ∂x + m = 0, (14b)
∆ = ∆ 0 + N i=1 ∆ i , m = m 0 + N i=1 m i , (14c)
∆ 0 = R 0 v, m 0 = G 0 q 0 , (14d)
∆ i = R i w i = L i d t v i , v = w i + v i , (14e)
m = C 0 d t p 0 , p = p 0 + q 0 , (14f)
m i = C i d t p i = G i q i , q 0 = p i + q i . ( 14g
)
Eliminating the unknowns ∆, m, w i and q i directly leads to (3), where only 2N + 1 auxiliary variables are necessary.

D. Energy balance

The energy identity comes from the circuit representation [START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF], and the associated equations ( 14). Let us multiply the first equation with v and the second equation with p, and integrate both in space over [0, L]. The terms v∆ and p m that arise can be interpreted as the power brought into the electronic circuit. They can be written, using the rest of system ( 14), as

             v∆ = R 0 v 2 + N i=1 R i w 2 i + L i 2 dv 2 i dt p m = C 0 2 dp 2 0 dt + G 0 q 2 0 + N i=1 G i q 2 i + C i 2 dp 2 i dt
Eliminating the same unknowns as before establishes the energy balance identity (5) where the energy is defined as ( 6) and the losses are given by (7).

IV. PRACTICAL IMPLEMENTATION

Space and time discretization of system (3) can be done with various numerical methods depending on the situation. In the context of sound synthesis, it is essential to design stable, accurate and efficient numerical schemes that couple the pipe with other elements such as the radiation at the bell, junctions with tone holes, or the embouchure, that can behave nonlinearly.

The current article proposes to use one-dimensional finite elements in space followed by an energy-consistent time discretization, in order to ensure numerical stability via an energy technique and to provide a numerical method that will easily extend to the presence of couplings. For the sake of simplicity, but without loss of generality, it is assumed that the outwards pipe flow is equal to given values λ -(t) and λ + (t) respectively at the pipe entrance and bell. In the presence of realistic coupling terms at the entrance and bell, these values will be unknowns and will require additional equations to be evaluated.

A. One-dimensional finite elements for space discretization

The finite element method (FEM) relies on a variational formulation (weak form) of the entire system in usual infinite dimensional functional spaces [START_REF] Brezis | Functional analysis, Sobolev spaces and partial differential equations[END_REF] to which belong p, v, p 0 , p i and v i , followed by the definition of finite dimensional spaces in which we seek numerically the approximate solutions p h , v h , p h,0 , p h,i and v h,i . One possible choice, called Mixed Spectral FEM, is described in (Tournemenne and Chabassier, 2019) and is followed here. It consists in using as finite dimensional spaces the set of piecewise polynomial functions of the spatial variable x, element by element, where jumps across element edges are authorized for v h and v h,i but not for p h , p h,0 and p h,i . These polynomial functions are chosen as the Lagrange interpolation polynomials on the Gauss-Lobatto points of each element. The order of the polynomial functions is called the order of the FEM and will be noted r thereafter. Finally, the integral terms coming from the variational formulation are evaluated using a quadrature formula based on the same Gauss-Lobatto points, so that mass matrices are diagonal [START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF].

This procedure results in the definition of matrices that replace all the spatial operators in (3). The following "semi-discrete system" is obtained. Let V h , V h,i , P h , P h,0 , P h,i be the vectors of coordinates of resp. v h , v h,i , p h , p h,0 and p h,i in the spanning basis of the finite elements, and see [START_REF] Tournemenne | A comparison of a one-dimensional finite element method and the transfer matrix method for the computation of wind music instrument impedance[END_REF] for the definition of the matrices

M V h , M P h , B h , E h . The diagonal matrices R h,0 , R h,i , G h,0 , G h,i , C h,0 , C h,i and L h,i , are
obtained via a similar procedure to the diagonal mass matrices M V h and M P h .

                                                       M V h ∂V h ∂t + R h,0 V h + N i=1 R h,i (V h -V h,i ) -B h P h = 0, (16a)
M P h ∂P h ∂t + G h,0 (P h -P h,0 ) + N i=1 G h,i (P h -P h,0 -P h,i ) (16b) +B * h V h + λ + E + h + λ -E - h = 0, (16c) L h,i ∂V h,i ∂t = R h,i (V h -V h,i ), for 1 ≤ i ≤ N, (16d)
C h,0 ∂P h,0 ∂t = G h,0 (P h -P h,0 ) (16e) + N i=1 G h,i (P h -P h,0 -P h,i ), (16f) C h,i ∂P h,i ∂t = G h,i (P h -P h,0 -P h,i ), for 1 ≤ i ≤ N. ( 16g
)
One advantage of this formulation is the natural treatment of the boundary conditions λ ± , which can become Lagrange multipliers for coupling with other systems. This semi-discrete system satisfies an analogue of the previous energy identity (5) where the continuous spatial norms are replaced with their semi discrete counterparts.

Note that the semi-discrete system ( 16) can be straightforwardly adapted for computation in the harmonic regime. This leads to one totally discrete system per value of ω, where the operator ∂ t is replaced with a multiplication with jω, which requires a sparse matrix inversion per ω.

B. Energy consistent time discretization

For time discretization, an interleaved scheme similar to [START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF] is used. The time step is denoted ∆t, and the step number n, so that t = n∆t. Unknowns P h , P h,0 and P h,i are evaluated at integer times {0, 1, 2, . . .}, whereas V h and V h,i are evaluated on a staggered time grid { 1 /2, 3 /2, 5 /2, . . .}. Let us define the discrete operators δ and µ acting on any vector sequence {X n } n∈N as

δX n+ 1 2 = X n+1 -X n ∆t , µX n+ 1 2 = X n+1 + X n 2 . ( 17
)
Conversely, if {Y n+ 1 2 } n is a sequence with indices on the staggered grid, {δY n } n is defined at integer times. The proposed totally discrete scheme reads

                                                                         M V h δV n h + R h,0 µV n h + N i=1 R h,i µ(V h -V h,i ) n -B h P n h = 0, ( 18a
)
M P h δP n+ 1 2 h + G h,0 µ(P h -P h,0 ) n+ 1 2 + N i=1 G h,i µ(P h -P h,0 -P h,i ) n+ 1 2 + B * h V n+ 1 2 h +λ n+ 1 2 + E + h + λ n+ 1 2 -E - h = 0, (18b) L h,i δV n h,i = R h,i µ(V h -V h,i ) n , ∀i, (18c)
C h,0 δP

n+ 1 2 h,0 = G h,0 µ(P h -P h,0 ) n+ 1 2 + N i=1 G h,i µ(P h -P h,0 -P h,i ) n+ 1 2 , (18d) C h,i δP n+ 1 2 h,i = G h,i µ(P h -P h,0 -P h,i ) n+ 1 2 , ∀i. (18e)
This system can be interpreted as an interleaved leap-frog scheme for pressure and flow, combined with an implicit midpoint resolution of the electronic circuit identities. It satisfies the following discrete equivalent of the energy balance (5) (see Appendix):

δE n+ 1 2 h = -µQ n+ 1 2 h,visc + Q n+ 1 2 h,therm + S n+ 1 2 h ( 19
)
where the discrete energy is defined as

E n h = E n h,V + E n h,P + µE n h,visc + E n h,therm + e n h . (20)
The different terms of energy are defined as

E n h,V := 1 2 µV n h 2 M V h , E n h,P := 1 2 P n h 2 M P h E n+ 1 2 h,visc := N i=1 1 2 V n+ 1 2 h,i 2 L h,i E n h,therm := 1 2 P n h,0 2 C h,0 + N i=1 1 2 P n h,i 2 C h,i e n h := ∆t 2 8 µV n h + N i=1 R -1 h,0 R h,i µ(V h -V h,i ) n 2 R h,0 Q n h,visc := µV n h,0 2 R h,0 + N i=1 µ(V h -V h,i ) n 2 R h,i Q n+ 1 2 h,therm := µ(P h -P h,0 ) n+ 1 2 2 G h,0 + N i=1 µ(P h -P h,0 -P h,i ) n+ 1 2 2 G h,i S n+ 1 2 h := -λ n+ 1 2 + (E + h ) * µP n+ 1 2 h -λ n+ 1 2 -(E - h ) * µP n+ 1 2 h where R h,0 = R h,0 (M V h ) -1 R h,0 ,
and for any time series of vectors {X n } n , {Y n } n , their scalar product is defined as

(X n , Y m ) = k X n k Y m
k and the weighted norm, for any non negative matrix A, is

X n A = (AX n , X n ) 1 /2
. This energy is positive as soon as the modified mass matrix M P h , defined as

M P h = M P h - ∆t 2 4 B * h (M V h ) -1
B h is a positive quadratic form. This leads to the following stability condition:

∆t ≤ 2 ρ (M P h ) -1 B * h (M V h ) -1 B h
which is the same condition as for solving the lossless wave equation with the classical leapfrog scheme [START_REF] Chabassier | Introduction and study of fourth order theta schemes for linear wave equations[END_REF]. Showing the convergence of the discrete scheme, as ∆t and h tend to zero, is out of the scope of the present work but it is expected that the dissipation of a positive energy will enable such a result [START_REF] Chabassier | Space/time convergence analysis of a class of conservative schemes for linear wave equations[END_REF].

C. Explicit update of the unknowns

One of the advantages of the of spectral finite elements is that

M P h , M V h , L h,i , C h,0 and C h,i
are diagonal matrices, making them trivial to invert. Thanks to an algorithmic elimination strategy, which can also be interpreted as a Schur complement, an explicit update of the unknowns is possible and is given in Appendix B. This approach generalizes the one given in Table I of [START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF] for finite differences, to finite elements of arbitrary order in matrix notations.

D. Accuracy and efficiency

Model ( 18) is applied to a simplified natural trumpet of total length 1.335 m. The pipe comprises two parts: a cylinder, followed by a "Bessel horn":

R(x) =            6 × 10 -3 m if 0 < x < 0.716, γ(x -x P ) -α if 0.716 < x < 1.335,
where α = 0.7 is the parameter of the Bessel horn. Constants γ = 4.404 37 × 10 -3 m and

x P = 1.358 97 m are chosen such that the radius is continuous at the junction, and that the bell radius is 60 mm, see Fig. 3. For simplicity, we assume a zero-pressure condition at the bell. An impulse response is computed by setting the input boundary term to a prescribed flow -λ

n+ 1 2 - = v 0 (n + 1 2 )∆t , with v 0 (t) =            8V 0 3t 1 sin 4 π t t 1 if 0 < t < t 1 , 0 otherwise, (22)
where the duration of the impulse is t 1 = 4 × 10 -4 s and the total injected volume is V 0 = 1 × 10 -7 m 3 . At the bell an elementary open condition p = 0 is applied: the generalization to more realistic radiation impedances is a natural possible extension of the present work.

The numerical scheme ( 18) is used to compute the impulse response at the temperature t = 20 • C with 34 elements of order 10 and a time step ∆t = 3.185 × 10 -6 s, which is the largest value satisfying the stability condition (21) . The final time is set to T = 0.2 s. Finally Fig. 7 shows the consecutive relative L ∞ -error in time of the impulse response with respect to the number of finite elements. As the spatial discretization of the pipe is refined, the time step is decreased by choosing the largest possible value given by the stability condition (21). The obtained impulse response is interpolated on a fixed time grid and compared. The resulting space-time convergence curve displays a second order rate of convergence.

V. COMPARISON WITH OTHER MODELS

Model (3) is compared with several models of the literature, in time domain and in frequency domain. The comparison includes two reduced 1D models: model (ZK) (1) and Webster-Loskhin model (WL) developed in [START_REF] Hélie | One-dimensional acoustic models of horns and comparison with measurements[END_REF], without the hypothesis of curved wavefronts, which is out of the scope of the present article. They are compared to 3D linear acoustic equations in the air column, associated with effective wall impedance boundary conditions of two types: from [START_REF] Cremer | On the acoustic boundary layer outside a rigid wall[END_REF] with an incident angle of π/2 noted (Cr), and from (Berggren et al., 2018) noted (BBN). All 3D solutions are obtained using finite elements for the spatial discretization using the software Montjoie 2 , in an axisymmetric configuration. The numerical solution relies on a curved mesh of the simplified natural closer to (WL) which indeed relies on Cremer wall impedances, and (BBN) to (ZK). It can be noted that the difference in resonance frequency (up to 1% i.e. 17 cents) between models is of the same order as the differential pitch sensitivity of the ear (around 4 cents) [START_REF] Micheyl | Influence of musical and psychoacoustical training on pitch discrimination[END_REF]. Fig. 10 quantitatively compares the models. With model (ZK) taken as a reference, the largest error of the other models on the impulse response is calculated. This error amounts to 1.5% for (WL). For the numerical schemes, the discrepancy results from two kinds of

G h = G h,0 + i G h,i , R h = R h,0 + i R h,i , C h,0 = C h,0 + ∆t 2 G h .
The next iterates P n+1 h , P n+1 h,0 , P n+1 h,i , V

n+ 3 /2 h , V n+ 3 /2 h,i
can be computed as follows, where all the matrices to invert are diagonal, beginning with P n+1 h :

M P h + ∆t 2 C -1 h,0 C h,0 G h P n+1 h = M P h -∆t 2 C -1 h,0 C h,0 G h P n h + ∆t C -1 h,0 C h,0 G h P n h,0 + ∆t i C -1 h,0 (G -1 h,i + ∆t 2 C -1 h,i ) -1 C h,0 P n h,i -∆t B * h V n+ 1 2 + λ n+ 1 2 + E + h + λ n+ 1 2 -E - h .

Now that

M V h + ∆t 2 R h V n+ 3 2 h = M V h -∆t 2 R h V n+ 1 2 h + ∆t i (R -1 h,i + ∆t 2 L -1 h,i ) -1 V n+ 1 2 h,i + ∆tBP n+1 h Finally, V n+ 3 2 h,i
is obtained as

L h,i + ∆t 2 R h,i V n+ 3 2 h,i = L h,i - ∆t 2 R h,i V n+ 1 2 h,i + ∆t 2 R h,i (V n+ 1 2 h + V n+ 3 2 h
).

1 See supplementary material at https://asa.scitation.org/doi/suppl/10.1121/10.0005537 for the tables of coefficients (a i , b i ) for N = 0 to 16.

  Sound velocity: c = 331.45 T /T 0 m s -1 Density: ρ = 1.2929 T 0 /T kg m -3

  2 4.075 43 × 10 -2 4.020 92 × 10 -4 3 8.148 25 × 10 -3 1.622 09 × 10 -5 4 1.961 59 × 10 -3 5.688 60 × 10 -7 Coefficients for N = 4 (23 it.) 864 11 × 10 -1 3.168 42 × 10 -2 2 8.063 38 × 10 -2 5.883 91 × 10 -3 3 3.520 99 × 10 -2 1.112 01 × 10 -3 4 1.533 51 × 10 -2 2.116 66 × 10 -4 5 6.695 83 × 10 -3 4.045 03 × 10 -5 6 2.932 51 × 10 -3 7.735 96 × 10 -6 7 1.328 25 × 10 -3 1.444 92 × 10 -6 8 9.403 66 × 10 -4 1.483 83 × 10 -7 Coefficients for N = 8 (114 it.) TABLE II. Coefficients (a i , b i ), optimized for M = 100 values of ζ ranging from 8 to 2 × 10 6 , covering radii of 1 × 10 -3 m to 0.1 m and frequencies of 20 Hz to 2 × 10 4 Hz. In parentheses: number of iterations to reach the stagnation threshold 10 -8 .

  FIG. 1. Dimensionless loss coefficient G(ζ), w.r.t. dimensionless frequency ζ. The range of interest

  FIG. 2. Relative error between G and G N , as a function of dimensionless frequency ζ, for several

FIG. 3 .

 3 FIG. 3. The simplified natural trumpet : radius w.r.t. position along bore

Fig. 4 FIG. 4 .

 44 Fig. 4 displays the computed impulse response, i.e. the evolution of the pressure at the

Fig 5 FIG

 5 Fig 5 displays the energy distribution with respect to time, according to the definitions

  FIG. 8. Impulse response: evolution of the pressure at the entrance of the simplified trumpet

  error: model error, and discretization error. The former is due to approximation (11) and can be reduced by increasing N ; the latter is due to the numerical scheme (18) and can be reduced by decreasing the time step ∆t, as done in Fig 10.For N = 8 the model error compared to (ZK) is observed to be about 0.14%, and for N = 16 it is less than 0.02%.VI. CONCLUSIONS AND PROSPECTSThis work presents a 1D model for viscothermal wave propagation suitable for timedomain simulation, which makes use of an adjustable number of auxiliary unknowns. It is derived from an approximation of the loss coefficients of the Zwikker-Kosten model with rational functions over the whole range of frequencies and radii of musical acoustics. The proposed improvement upon previous work is that the coefficient optimization procedure has been done only once, as the objective function does not depend on temperature or pipe radius. It leads to a closed-form for the resulting model, where the user only needs to specify the instrument geometry and physical constants. Numerical comparison with 3D models show that the approximation error is smaller than the discrepancies between different models. The use of auxiliary variables induces a numerical burden which was to be expected for accounting for viscothermal effects. The model satisfies an energy identity and is therefore suitable for time-domain coupling with other models as sound radiation, or reed evolution, which are a possible extension of this work. This work could be extended to include additional forms of acoustic losses, such as wall admittance due to porosity, or to better justify its use in tubes with variable cross-sectional area.Taking the scalar product of µ(18a) with µµV yields M V µδV + µ∆, µµV -(BµP, µµV ) = 0. The scalar product of (18b) with µP gives M P δP + m, µP + (B * V, µP ) + λ -E -+ λ + E + , µP B and C are treated separately.Using (A3) and the values of A , B and C , and replacing ∆ with its value, leads to the expected relation (19

2

  http://montjoie.gforge.inria.fr 3 http://openwind.gitlabpages.inria.fr/web Abel, J., Smyth, T., and O. Smith, J. (2003). "A simple, accurate wall loss filter for acoustic tubes," International Conference on Digital Audio Effects 2003 Proceedings, London, UK 2003, 53-57.

TABLE I

 I

. Numerical values

[START_REF] Chaigne | Acoustics of musical instruments[END_REF]

of air constants used in the model. t is the temperature in Celsius, and T the absolute temperature with T 0 = 273.15K.

  ).

	APPENDIX B: EXPLICIT UPDATE OF THE UNKNOWNS
	Assume P n h , P n h,0 , P n h,i , V h n+ 1 /2	, V h,i n+ 1 /2	are known. Define

 (21)

APPENDIX A: DISCRETE ENERGY BALANCE 1. Useful identities Using discrete-time operators µ and δ defined in (17), the following identities, implicitly centered in t n+ 1 2 , hold

Notations

In this appendix, for the sake of brevity, subscripts h and time integers n and n + 1 /2 are omitted. We denote

≡ C 0 δP 0 . Moreover we use the notations defined in section IV B.

Now, the discrete work (∆, µV ) writes

Moreover, the discrete work (m, µP ) writes

, µ(P -P 0 ))

Finally, using (A1b) and commutativity of µ and δ,