Xavier Antoine
email: xavier.antoine@mines.inpl-nancy.fr

Chokri Chniti

Karim Ramdani
email: ramdani@iecn.u-nancy.fr

On the numerical approximation of high-frequency acoustic multiple scattering problems by circular cylinders

The aim of this paper is to propose a numerical strategy for computing the solution of two-dimensional time-harmonic acoustic multiple scattering problems at high-frequency. The scatterers are assumed to be circular, leading therefore to semi-analytical representation formulae of the scattered field through the solution of a large linear system of equations. Taking advantage of the special block Toeplitz structure of the matrix of the linear system, a fast iterative and preconditioned numerical method yielding large memory savings is proposed. Several numerical experiments for general configurations are presented to show the efficiency of the numerical method.

Introduction

Multiple scattering problems find their origins in many applications related to different areas of applied sciences: acoustics, electromagnetism, elasticity and water waves. For such problems, the scattered field appears as the superposition of elementary scattered fields resulting from the interaction between the incident wave and the scatterers on one hand, and between the scatterers on the other hand. A better understanding of such problems requires a precise knowledge of the influence of the different physical and geometrical parameters of the problem on these interactions.

Due to the complexity of this problem, computing a numerical solution requires a special care, especially in the case where the number of objects is large and/or for high frequencies. For a recent and complete overview of these issues with extended bibliographical references, we refer the reader to the textbook [START_REF] Martin | Multiple Scattering, Interaction of Time-Harmonic Waves with N obstacles[END_REF] by P.A. Martin. In this paper, we investigate the two-dimensional Helmholtz equation for circular cylinders.

Following Záviska and next Row [START_REF] Záviska | Über die Beugeng elektromagnetischer Wellen an parallelen, unendlich langen Kreiszylindern[END_REF][START_REF] Row | Theoretical and experimental study of electromagnetic scattering by two-dimensional conducting cylinders[END_REF][START_REF] Martin | Multiple Scattering, Interaction of Time-Harmonic Waves with N obstacles[END_REF], a multipole expansion formulation of the solution can be obtained. Unlike the case of one object where the Mie series solution expressed by the Fourier coefficients is explicitly obtained through the inversion of a diagonal matrix, the case of multiple scattering is more complex and requires the use of numerical methods. Indeed, it can be shown (see Section 2 or [START_REF] Martin | Multiple Scattering, Interaction of Time-Harmonic Waves with N obstacles[END_REF] for more details) that the scattered field admits an expansion as a superposition of Fourier series. The set of Fourier coefficients is then solution of a nondiagonal and dense complex linear system which can be extremely large for many scatterers at high-frequency. The numerical solution of this system has been already obtained by some authors [START_REF] Linton | The interaction of waves with arrays of vertical circular cylinders[END_REF][START_REF] Linton | Handbook of Mathematical Techniques for Wave/Structure Interactions[END_REF][START_REF] Tsang | Scattering of Electromagnetic Waves, Numerical Simulation[END_REF][START_REF] Yilmaz | An iterative procedure for the diffraction problem of water waves by multiple cylinders[END_REF][START_REF] Yilmaz | Analytical solutions of the diffraction problem of a group of truncated vertical cylinders[END_REF][START_REF] Yilmaz | Wave enhancement due to blockage in semi-submersible and TLP structures[END_REF][START_REF] Zhang | Fast multipole scattering matrix method for multiple scattering of a large number of cylinders[END_REF] (see also [START_REF] Gumerov | Multiple scattering from N spheres using multipole reexpansion[END_REF] for scattering by spheres). However, this remains restricted to low-moderate wavenumbers. Alternative solutions using numerical methods based on Dirichletto-Neumann non-reflecting boundary conditions have been recently developped by Grote and Kirsch [START_REF] Grote | Dirichlet-to-Neumann boundary conditions for multiple scattering problems[END_REF]. This technique has the advantage of allowing general shaped scatterers but, on the other hand, it requires to mesh the finite computational domain. In particular, this leads to strong limitations (as in the single scattering problem) as the wavenumber increases since e.g. small element sizes are required.

To the best of the authors' knowledge, the only papers dealing with high frequencies for multiple scattering problems are the ones by Bruno, Reitich and their co-authors [START_REF] Anand | Analysis of multiple scattering iterations for high-frequency scattering problems. II: The three-dimensional scalar case[END_REF][START_REF] Bruno | Prescribed error tolerances within fixed computational times for scattering problems of arbitrarily high-frequency: the convex case[END_REF][START_REF] Bruno | A new high-order high-frequency integral equation method for the solution of scattering problems. II: Multiple-scattering configurations[END_REF][START_REF] Ecevit | Analysis of multiple scattering iterations for high-frequency scattering problems. I: The two-dimensional case[END_REF] where scattering by three convex and non-convex objects is analyzed for two-and three-dimensional problems. The authors consider an integral equation formulation of the Dirichlet problem and derive a fast highorder method for computing the multiple scattered field. This robust solver therefore allows to consider general shaped obstacles. While the technique is efficient and robust, it is very technical and requires specific developments in the background of integral equations. The present paper proposes to come back to the scattering problem by circular cylinders using the multipole expansion and to develop a reliable and robust accurate numerical solution for large frequencies and many scatterers.

The outline of the paper is the following. In Section 2, we recall some classical results on the problem of multiple scattering by circular cylinders and on its resolution using multipole expansion formulation. In particular, we will see that the Fourier coefficients of the local scattered field solve an infinite linear system (see [START_REF] Martin | Multiple Scattering, Interaction of Time-Harmonic Waves with N obstacles[END_REF]). Section 3 explains how to suitably truncate the above infinite dimensional system to get an accurate numerical solution. Section 4 deals with the numerical resolution of the complex, dense and large size linear system obtained. We take advantage of the special block Toeplitz structure of the linear system to reduce the memory storage and accelerate the computation of the matrix. Next, we propose a fast iterative solution of the linear system using a Krylov solver in conjunction with an original geometrical preconditioner. Several examples, corresponding to many geometrical configurations, are provided to understand the characteristics of the algorithm for high wavenumbers and for a large number of scatterers. Finally, in Section 5, we conclude by summing up our results and outlining some future directions of work.

Multiple scattering by M circular cylinders

Consider a homogeneous acoustic medium filling the whole space R 2 and containing M disjoint scatterers B 1 , . . . , B M . We assume that each scatterer B p , p = 1, . . . , M , is a bounded subdomain of R 2 of boundary ∂B p . We denote by Ω -= M p=1 B p the domain occupied by the obstacles. We consider the scattering problem of an incident plane wave u inc (r) = e ikβ•r of direction β by Ω - (where the time dependence is assumed to be of the form e -iωt and where the wavenumber k is real). In other words, we want to determine the scattered wave u solving the exterior boundary value problem

(E)        ∆u + k 2 u = 0 (R 2 \ Ω -) Λu = -Λu inc (∂Ω -) lim |x|→+∞ |x| 1/2 ∇u • x |x| -iku = 0.
To model sound-soft and sound-hard obstacles, the boundary operator Λ appearing in the above relations denotes either the trace or the normal trace operator on ∂Ω -. A natural idea to solve numerically the multiple scattering problem (E) is to reduce it to a family of single scattering problems. By linearity, this can be achieved by introducing M fictitious scattered waves u 1 , . . . , u M , where each field u p corresponds to the wave reflected only by the scatterer p when it is illuminated by the incident wave and the scattered waves u q , for q = 1, . . . , M with q = p. More precisely, we have the following result (see [START_REF] Balabane | Boundary decomposition for Helmholtz and Maxwell equations. I. Disjoint subscatterers[END_REF] for the proof).

Theorem 1. Let u be the solution of the multiple scattering problem (E). Then, the family of M coupled single scattering problems for p = 1, . . . , M :

(E p )                ∆u p + k 2 u p = 0 (R 2 \ B p ) Λu p = -Λ   u inc + M q=1,q =p u q   (∂B p ) lim |x|→+∞ |x| 1/2 ∇u • x |x| -iku = 0.
admits a unique solution (u 1 , . . . , u p ). Furthermore, the following decomposition holds

u = M p=1 u p . (1)
In this paper, we focus on the case where the scattering obstacles B p are circular cylinders of radius a p centered at the points O p = (x p , y p ) of a given orthonormal system of coordinates

( -→ Ox, -→ Oy). O p b p α O p α pq b pq O q x y a q Figure 1: A view of two typical cylinders.
We set for all p = 1, . . . , M (see Figure 1):

b p = --→ OO p b p = |b p | α p = Angle(Ox, b p )
and for all q = 1, . . . , M , with q = p:

b pq = ---→ O q O p b pq = |b pq | α pq = Angle(Ox, b pq ).
Any point M of the plane will be described by its cartesian coordinates (x, y) or by its polar coordinates:

r = --→ OM r = |r| θ = Angle(Ox, r).
We will also use in the sequel the local polar coordinates of the point M in the orthonormal system of coordinates associated to the scatterer p:

r p = ---→ O p M r p = |r p | θ p = Angle(Ox, r p ).
Let us introduce for all m ∈ Z the following cylindrical wavefunctions, which are particular solutions of Helmholtz equations for r > 0:

ψ m (r) = H (1) m (kr)e imθ , ψ m (r) = J m (kr)e imθ , (2)
where J n is the nth order Bessel function and H

n is the nth order Hankel function of the first kind.

Note that for all m ∈ Z, the function ψ m satisfies the outgoing Sommerfeld radiation condition.

We also define for all m ∈ Z the local cylindrical wavefunctions associated with the scatterer p, for p = 1, . . . , M , by setting:

ψ p m (r) = ψ m (r p ) = H (1) m (kr p )e imθp , ψ p m (r) = ψ m (r p ) = J m (kr p )e imθp , ∀ m ∈ Z. (3)
Since each field u p is an outgoing solution of a single scattering problem outside a disk, it admits the following modal decomposition in the local cylindrical outgoing wavefunctions:

u p (r) = m∈Z c p m ψ p m (r), ∀ p = 1, . . . , M, ∀r p > a p . (4)
where the complex coefficients (c p m ) m∈Z are determined by imposing the boundary condition on the boundary of the scatterer p:

Λu p = -Λu inc - M q=1,q =p Λu q on ∂B p . (5)
In order to explicit this equation, we have to express the incident field u inc and the fictitious scattered fields u q , for q = p, in the local system of coordinates of the scatterer p. For the incident plane wave of direction β = (cos β, sin β) we recall that (cf. [19, p. 125])

u inc (r) = m∈Z d p m ψ p m (r) (6)
where d p m = e ikβ•bp e im( π 2 -β) . Concerning the fields u q , we make use of the separation theorem (see for instance [START_REF] Martin | Multiple Scattering, Interaction of Time-Harmonic Waves with N obstacles[END_REF]Theorem 2.12]).

Theorem 2. Let 1 ≤ p, q ≤ M , with p = q. Then, we have the following relations:

ψ q m (r) =        n∈Z S mn (b pq ) ψ p n (r) for r p < b pq , n∈Z S mn (b pq )ψ p n (r) for r p > b pq , ∀ m ∈ Z, (7)
where we have set

S mn (b pq ) = ψ m-n (b pq ), S mn (b pq ) = ψ m-n (b pq ). ( 8
)
The infinite matrices S p,q = (S mn (b pq )) m,n∈Z and S p,q = S mn (b pq ) m,n∈Z are called separation (or transfer) matrices.

Using relations (1), ( 4), [START_REF] Balabane | Boundary decomposition for Helmholtz and Maxwell equations. I. Disjoint subscatterers[END_REF] and the first equation in [START_REF] Bruno | Prescribed error tolerances within fixed computational times for scattering problems of arbitrarily high-frequency: the convex case[END_REF], straightforward computations show that the unknown Fourier coefficients solve the following equations

c p m + J m (ka p ) H (1) m (ka p ) M q=1,q =p n∈Z S nm (b pq )c q n = - J m (ka p ) H (1)
m (ka p )

d p m , ∀ m ∈ Z, ∀p = 1, . . . , M. (9)
in the sound-soft case and

c p m + J m (ka p ) H (1) m (ka p ) M q=1,q =p n∈Z S nm (b pq )c q n = - J m (ka p ) H (1)
m (ka p )

d p m , ∀ m ∈ Z, ∀p = 1, . . . , M. ( 10
)
in the sound-hard case (note that H

n (ka p ) and H

n (ka p ) never vanish for real wavenumbers).

The infinite linear systems ( 9) and ( 10) can be written in the more compact vector form:

C p + D p M q=1,q =p (S p,q ) T C q = B p ∀ p = 1, . . . , M, (11) where
• C p = (c p n ) n∈Z is the infinite vector containing the coefficients of the cylindrical decomposition (4) of u p ,

• (S p,q ) T denotes the transpose of the separation matrix S p,q between the obstacles B p and B q defined by S p,q = (S p,q mn ) m∈Z,n∈Z S p,q mn = ψ m-n (b pq ),

• D p = (D p mn ) mn∈Z is the diagonal infinite matrix, with diagonal terms

D p m,m =          J n (ka p ) H (1)
n (ka p )

for sound-soft obstacles,

J n (ka p ) H (1)
n (ka p )

for sound-hard obstacles,

• B p = -D p d p , where d p = (d p m ) m∈Z is the infinite vector containing the coefficients of the cylindrical decomposition (6) of the incident wave.

The M infinite linear systems [START_REF] Chew | On the spectrum of the electric field integral equation and the convergence of the moment method[END_REF] can equivalently be written in the abstract form

AC = B (13)
where

A =       I D 1 S 1,2 T . . . D 1 S 1,M T D 2 S 2,1 T I . . . D 2 S 2,M T . . . . . . D M S M,1 T D M S M,2 T . . . I       , (14)
C =      C 1 C 2 . . . C M      B =      B 1 B 2 . . . B M     
and where I denotes the identity operator on 2 (C).

Let us emphasize that many physical quantities of interest can be deduced from the solution of [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]. The scattered field u can be obtained in the neighborhood of B p (namely, for r p < min

1≤q≤M q =p b pq )
by the relation

u(r) = m∈Z c p m ψ p m (r) + m∈Z   M q=1,q =p n∈Z S nm (b pq )c q n   ψ p m (r). ( 15
)
Using

r p = r -b p cos(θ -α p ) + O(1/r), the scattering amplitude a(θ) defined by u(r) = e ikr √ r a(θ) + O 1 r , as r → +∞.
is obtained by

a(θ) = e -iπ/4 2 πk p=M p=1 e -ibpk cos(θ-αp) n∈Z e in(θ-π 2 ) c p n . (16)
For the Dirichlet multiple scattering problem, one can be interested in computing the normal derivative ∂u(r) ∂r p on the boundary of the scatterer p. This can be achieved by taking the normal derivative of [START_REF] Grote | Dirichlet-to-Neumann boundary conditions for multiple scattering problems[END_REF]. Following [START_REF] Linton | The interaction of waves with arrays of vertical circular cylinders[END_REF], let us notice that the resulting formula can be simplified (to eliminate the double sum) using ( 13) and the Wronskian relation for Bessel functions. The final expression reads

∂ n U p = -∂ n U inc + 2i πa p J -p C p , p = 1, ..., M, (17)
where ∂ n U p (respectively ∂ n U inc ) is the infinite dimensional vector containing the Fourier coefficients of the normal derivative of the scattered (respectively incident) field on ∂B p and J -p is the infinite diagonal matrix with diagonal terms 1/J m (ka p ). A similar formula for the trace can be derived in the case of Neumann boundary conditions. Note that equation ( 17) is particularly well adapted for the numerical computations since it only involves diagonal matrices.

Finite dimensional approximation

Obviously, the infinite linear system (11) (or ( 13)) must be truncated to be solved numerically. In particular, only a finite number of modes can be used to describe the solution u p . Let 2N p +1 be the corresponding number of modes used for the scatterer p, i.e. assume that we keep only the modes

ψ p m such that -N p ≤ m ≤ N p .
Note that we have chosen a possibly different number of modes for each scatterer, in order to take into account geometrical configurations where the obstacles have different radii. With this notation, [START_REF] Chew | On the spectrum of the electric field integral equation and the convergence of the moment method[END_REF] should be truncated to the finite system

C p + D p M q=1,q =p (S p,q ) T C q = B p ∀ p = 1, . . . , M, (18)
where

• C p = (c p n ) n=-Np,.
..,Np is the finite vector containing approximations of the first 2N p + 1 modal coefficients of the cylindrical decomposition (6) of u p , that we still denote c p n for the sake of clarity.

• S p,q is the (2N p + 1) × (2N q + 1) finite dimensional separation matrix taking into account only the interactions between the first modes of the obstacles B p and B q :

S p,q = (S p,q mn ) -Np≤m≤Np,-Nq≤n≤Nq S p,q mn = ψ m-n (b pq ),

• D p = (D p mn ) -Np≤m≤Np,-Nq≤n≤Nq is the diagonal finite matrix, with diagonal terms

D p m,m =          J m (ka p ) H (1)
m (ka p ) for sound-soft obstacles,

J m (ka p ) H (1)
m (ka p ) for sound-hard obstacles,

• B p = -D p d p , where d p = (d p m ) -Np≤m≤Np is the finite vector containing the 2N p + 1 first coefficients of the cylindrical decomposition (4) of the incident wave.

The M coupled finite dimensional systems [START_REF] Linton | Handbook of Mathematical Techniques for Wave/Structure Interactions[END_REF] can equivalently be written as

AC = B ( 19
)
where A ∈ C N,N is the full complex square matrix of size N = M p=1 (2N p + 1) defined by

A =       I 1 D 1 S 1,2 T . . . D 1 S 1,M T D 2 S 2,1 T I 2 . . . D 2 S 2,M T . . . . . . D M S M,1 T D M S M,2 T . . . I M       (20)
where I p denotes the identity matrix of C 2Np+1 and

C =      C 1 C 2 . . . C M      B =      B 1 B 2 . . . B M     
are to complex-valued vector fields from C N .

The number of modes N p to retain in the approximation must be fixed to get an accurate solution. On one hand, N p must be large enough to capture both the propagating and grazing parts of the solution (typically, N p ≥ ka p ). On the other hand, taking too many modes for approximating the solution makes the matrix A ill-conditioned. Therefore, a numerical stagnation occurs when using an iterative algorithm for solving [START_REF] Martin | Multiple Scattering, Interaction of Time-Harmonic Waves with N obstacles[END_REF] (see Section 4.3.1). This phenomenon is due to the fact that high order spatial modes |m| correspond to the evanescent part of the field and therefore, computing their Fourier coefficients is definitively out of reach using an iterative solver (with a fixed tolerance). For our simulations, we will use the following empirical formula

N p = ka p + 1 2 √ 2 ln(2 √ 2πka p ε -1 ) 2 3 (ka p ) 1/3 + 1 , (21)
where [x] denotes the integer part of a real number x, and ε is the desired error bound on the Fourier coefficients. The above formula has been proposed in the literature in the contexts of single scattering [START_REF] Chew | Fast and Efficient Algorithms in Computational Electromagnetics[END_REF] and multipole methods [START_REF] Carayol | Error estimates in the Fast Multipole Method for scattering problems. Part 1: Truncation of the Jacobi-Anger series[END_REF]. Nevertheless, according to our numerical results

(see Section 4.3.1), it turns out that it can also be successfully used in the framework of multiple scattering investigated in this paper.

Let us give here a formal derivation of the above formula in the case of simple scattering (see the proof given in [START_REF] Chew | Fast and Efficient Algorithms in Computational Electromagnetics[END_REF] for more details). First of all, the convergence of the Fourier series

m∈Z c p m H (1)
m (ka p )e imθp in L 2 (0, 2π) implies that for |m| large enough, we have

|c p m | ≤ 1 |H (
|c p m | ≤ 2 √ 2πka p 2 exp - √ 2 ka p ζ 3/2
Therefore, we can ensure an error bound of ε on the unknown Fourier coefficients by imposing the right-hand side of the above relation to be less than ε. This gives

ζ = 1 2 √ 2 ln(2 √ 2πka p ε -1 ) 2 3 (ka p ) -2/3
and thus an order of truncation N p = ka p (1 + ζ) which is exactly formula [START_REF] Row | Theoretical and experimental study of electromagnetic scattering by two-dimensional conducting cylinders[END_REF].

4 Implementation and validation of a fast preconditioned iterative numerical method

We want to solve the linear system [START_REF] Martin | Multiple Scattering, Interaction of Time-Harmonic Waves with N obstacles[END_REF] for realistic complex configurations, including numerous obstacles M and high frequencies ka p . Therefore, according to [START_REF] Row | Theoretical and experimental study of electromagnetic scattering by two-dimensional conducting cylinders[END_REF], the number of modes 2N p + 1 needed to approximate the solution u p with a reasonable precision must then be large. This implies that a huge memory storage is a priori required and large computational times are needed to build the dense matrix A. Moreover, we are led to solve a large scale complex-valued linear system with size N × N , with N = M p=1 (2N p + 1). Using a direct linear solver would yield a prohibitive computational time, especially for high frequencies. Our solution consists in developing a strategy based on an iterative solver leading to problems related to fast evaluations of dense matrix-vector products as well as convergence questions. We carefully analyze these delicate problems in the present section after presenting the different geometrical test configurations.

Geometrical configurations

In the numerical computations, we consider three kinds of geometrical configurations:

• The single-row configuration: this structure is composed of M x (= M ) equally spaced obstacles aligned along the x-axis, the distance between two successive scatterers being denoted b x = b 12 (see Figure 2). Moreover, the row is centered at x = 0.

Storage

Although being full, the matrix A has a particular structure. Indeed, each one of its off-diagonal blocks is obtained (see [START_REF] Nédélec | Acoustic and electromagnetic equations[END_REF]) by multiplying the diagonal matrix D p ∈ C 2Np+1,2Np+1 by the matrix (S p,q ) T ∈ C 2Np+1,2Nq+1 which has a Toeplitz structure [START_REF] Chen | Matrix Preconditioning Techniques and Applications[END_REF] since S p,q mn = ψ m-n (b pq ).

Consequently, using the notations from [START_REF] Chen | Matrix Preconditioning Techniques and Applications[END_REF], the storage of (S p,q ) T can be optimized using a compressed version based on the root vector σ p,q = (S p,q Nq,-N p , ..., S p,q -Nq+1,-Np , S p,q -Nq,-Np , ..., S p,q -Nq,Np ) T .

To take advantage in our future algorithms of the special structure of each off-diagonal block of A, we store both the Toeplitz matrix (S p,q ) T (through the root-vector ( 22)) and the diagonal matrix D p . According to our notations, the compressed storage needs 2(2N p + N q + 1) entries instead of the (2N p + 1)(2N q + 2) complex coefficients required for the full version. For A, this must be repeated for the M (M -1) off-diagonal blocks by summing over p and q. This results in a global storage equal to 3N (M -1) entries which must be compared to the N 2 -N M cost of the full storage. Furthermore, the computational time involved in the construction of the global matrix is also reduced according to the memory storage. In the case where we have a p = a for any p, the vector root version of A leads to a memory storage and a CPU time of the order of O(6kaM 2 ) while it is O(4k 2 a 2 M 2 ) for the full version. This is a crucial point for solving a multiple scattering problem for a large wavenumber. To show the improvement induced by the compressed storage version using [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] over the full version, we represent on Figure 5 the logarithm of the CPU time1

scaled by the computational memory requirement with respect to the wave number ka p needed for building the global matrix A. We consider a single-row configuration for M = 2, a p = 1, for p = 1, 2, and b x = 3. As expected, the CPU time for the compressed version is linear according to ka p while it is quadratic for the full version.

Iterative solution

Since we only have access now to the compressed Toeplitz format [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] of the matrix A, the linear system [START_REF] Martin | Multiple Scattering, Interaction of Time-Harmonic Waves with N obstacles[END_REF] can only be solved through an iterative linear solver. We consider both the GMRES, possibly with a restart parameter η (then denoted by GMRES(η)), and the BICGStab algorithms.

The tolerance error of the iterative solver is set to tol and the number of iterations to get this tolerance is denoted by n iter . Let us recall that the main CPU cost of the GMRES is due to one Matrix-Vector Product (MVP) per iteration. Moreover, a large amount of memory can be necessary to store the Krylov basis if too many iterations are needed to reach the tolerance. We will see that it leads us to consider a restart parameter if the number of obstacles is large and the frequency high. Two MVPs are required for each iteration of the BICGStab. This finally gives the total cost which is related to each algorithm. We first study numerically the influence of the order of truncation N p on the convergence of the iterative solver. In particular, this validates the choice of N p given by formula [START_REF] Row | Theoretical and experimental study of electromagnetic scattering by two-dimensional conducting cylinders[END_REF]. Next, we explain how to reduce the computational cost of a global MVP using a fast algorithm for each Toeplitz subblock. We investigate numerically the dependence of the rate of convergence of the iterative solver with respect to the geometrical and physical parameters. Finally, we propose and study a geometrically-based preconditioner.

From now on, we consider an angle of incidence β = 0 and we restrict the presentation to the Dirichlet case (which appears to be harder to solve than the Neumann problem in terms of iterations).

Influence of the order of truncation N p on the convergence

As noticed in Section 3, the order of truncation N p must be fixed carefully e.g. through formula [START_REF] Row | Theoretical and experimental study of electromagnetic scattering by two-dimensional conducting cylinders[END_REF]: N p must be large enough to compute accurately the solution, but not too large to avoid the stagnation of the iterative solver. To make this statement precise, let us consider a uniform square lattice with M x = M y = 2 (M = 4) and b x = b y = 3 for a radius a p = 1, 1 ≤ p ≤ M , and k = 100.

The linear system ( 19) is solved by the GMRES with tol = 10 -8 . For a given value of tol and all along the paper, we fix ε = tol in formula [START_REF] Row | Theoretical and experimental study of electromagnetic scattering by two-dimensional conducting cylinders[END_REF]. This gives here N p = 120 (represented by a red dot on Figures 6 and8). We report n iter versus N p on Figure 6. One can observe three distinct zones. First, from N p = 1 to N p < ka p + 2, n iter increases. This means that the computation of a correct solution requires more harmonics. This is obtained in the second stable zone (for ka p + 2 ≤ N p ≤ 138). However, if we include too many harmonics (third zone), typically N p ≥ 139 in our example, then we obtain a break down of the GMRES as it can be remarked on Figures 678. Indeed, stagnation occurs e.g. for N p = 150 while it does not for N p = 120. In particular, the relative error corresponding to the stagnation at N p = 150 is equal to 10 -2 while it is tol for N p = 120.

Fast MVPs for Toeplitz matrices

The direct computation of a MVP y = Ax can be computed by blocks. Let us set: x = (x 1 , ..., x M ), with x p ∈ C 2Np+1 , and y = (y 1 , ..., y M ), with y p ∈ C 2Np+1 . Then, we directly have from the structure of A given by ( 20)

y = x + D p 1≤p = ≤M z , z = (S ,p ) T x p , (23)
with 1 ≤ ≤ M . The main cost of the above evaluation is linked to the computation of z which is quadratic according to (2N p + 1). Moreover, this must be repeated for each sub-block and each component of y . This is very expensive when the frequency is large since the size of (S ,p ) T is (2N + 1) × (2N p + 1). Another way of computing a Toeplitz MVP for a matrix of size n × n is to use the fast algorithm explained in ([10, pp. 95-96]) for MVPs involving Toeplitz matrices. The idea consists in building an associated circulant matrix using the Toeplitz matrix and next applying an FFT-based MVP algorithm for circulant matrices. This algorithm is coded using Matlab FFT function. The resulting total cost in terms of real operations for computing y is composed from: one complex-valued Toeplitz MVP in 15n log 2 (n) + n operations (see [START_REF] Boyd | Chebychev and Fourier Spectral Methods[END_REF] page 193), with n = 2N p + 2N + 2, summing up next on 1 ≤ p = ≤ M , the computation of the diagonal matrix D p which is 6N + 3, and finally adding x . Again, summing up on = 1, ..., M gives the total cost. If a p = a, then, this requires asymptotically O(60(M -1) 2 ka log 2 (4ka)) operations compared to O(4(M -1) 2 ka 2 ) for a direct MVP. An example is given on Figure 9 showing the CPU time reduction with respect to ka (with a = a p , 1 ≤ p ≤ M ) using the fast MVP algorithm

Numerical study of the convergence rate

Iterative solver. As already mentioned, the linear system (19) must be solved iteratively using either the BICGStab or the GMRES [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. From a large set of numerical simulations, the GMRES provides the fastest convergence rate (see Figure 10). However, a more important memory storage is required which can significantly limit the possibility of prospecting high-frequencies for example. We present on Figure 10 the behaviour in terms of MVPs of the GMRES, BICGStab and GMRES(50) according to ka in the single-row configuration. We can see that the GMRES breaks down at ka = 140 due to memory limitations. This is not the case of BICGStab and GMRES(50) which keep on working. Moreover, it appears that the restarted GMRES generally leads to similar or better convergence results than the BICGStab. For these reasons, we choose the GMRES(50) in the sequel. Another way of improving the convergence rate consists in preconditioning the linear system. This solution is analyzed in Section 4.3.4 where we build a geometrical preconditioner.

Single-row configuration. Figure 10 shows the dependence of the number of MVPs according to the wavenumber ka. As seen, for this complex configuration, the number of MVPs increases with ka. For a fixed wavenumber, we also observe a dependence with respect to the number of obstacles M and the distance δ between two obstacles: δ = b x -a. We remark on Figure 11 that the number of MVPs increases linearly with M . Moreover, the slope of the line is more important as δ tends toward zero, meaning that we have two closer successive scatterers. This is observed again on Figure 12 where the number of MVPs is given according to δ in logarithmic scale for ka = 100 and M = 10. We see that the number of MVPs strongly decreases as the separation distance δ tends to infinity, i.e. δ λ. This corresponds to a weaker coupling between the obstacles in the multiple This generally creates a deterioration of the condition number at these frequencies, due to small eigenvalues in the matrix of the system. A similar problem also arises in the context of integral equations [START_REF] Chew | On the spectrum of the electric field integral equation and the convergence of the moment method[END_REF][START_REF] Chew | Fast and Efficient Algorithms in Computational Electromagnetics[END_REF][START_REF] Antoine | Analytic preconditioners for the electric field integral equation[END_REF][START_REF] Antoine | Analytic preconditioners for the boundary integral solution of the scattering of acoustic waves by open surfaces[END_REF].

Rectangular and triangular lattice configurations. We now analyze the rectangular lattice configuration. A first test-case is given on Figure 13. We consider a rectangular lattice with M x = 8 and increase the size of layers according to M y = 2, 4, 8. We fix b x = 8 and b y = 13. The number of MVPs required by the GMRES(50) is represented as a function of ka. We observe a stabilization with the frequency but the number of MVPs is larger with M y , and thus, with M . This is consistent with the previous observations in the single-row case. This situation is, in some sense, not extreme because the distance between the obstacles is sufficiently large.

A more difficult problem is considered on Figure 14 where we represent the number of GMRES(50) iterations of an 5 × M y lattice, for different values of M y , according to ka (for b x = b y = 3). We observe that the number of MVPs is again slightly dependent of ka but strongly varies with the number of layers, characterized by M y . Moreover, some peaks appear at some frequencies, and in particular at ka ≈ 20. This is more clearly visible on Figure 15 where we increase the lattice size according to M = M 2 x at ka = 50. We notice that the number of MVPs increases strongly with M x . In particular, on this example, it can be shown that, on the first values, the number of MVPs behaves like 3M 3.3 and not quadratically with M . Another way of considering particular frequencies where peaks occur consists in modifying δ. Some numerical computations, not reported here, confirm this property.

Concerning the triangular lattice, similar conclusions can be drawn. However, it appears that this situation is less dramatic in terms of MVPs compared to the rectangular lattice case.

A geometrically-based preconditioner

Since many iterations may be necessary in some situations, one way to improve the rate of convergence is to precondition the linear system. In our context, wa cannot directly apply an algebraic strategy like the incomplete LU or SPAI preconditioners [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF][START_REF] Chen | Matrix Preconditioning Techniques and Applications[END_REF] which would require to reconstruct the full version of the matrix A. An alternative direction is to build a geometrically-based preconditioner. We propose here a simple procedure in two steps. In fact, there is no reason to assume that F satisfies ρ(F) < 1, where ρ(F) stands for the spectral radius of F. However, relation [START_REF] Yilmaz | An iterative procedure for the diffraction problem of water waves by multiple cylinders[END_REF] must be viewed as a formal way of building an approximation of the inverse of A and so has also a subjacent limitation range. It could be possible to choose more terms in the approximation. However, extensive numerical computations show that this is not a good strategy. It can even lead to the divergence of the method.

Since P is still a matrix taking all the interactions between the obstacles, it is interesting to reduce its application cost by considering only the closest interactions. This can be done through a second approximation by introducing a parameter d > 0 representing a maximal coupling interaction distance. Then, the preconditioner, denoted by P d , only considers the interactions between obstacles with indices 1 ≤ p, q ≤ M satisfying: b pq < d. We must notice that the construction of To show the improvement induced by the proposed preconditioner, we present on Figure 16 the number of MVPs for a single-row configuration with M = 75 obstacles. We observe the improvement in terms of convergence rate if we compare these results to the ones obtained on ). We see an interesting gain in terms of reduction of iterations. However, it appears that for situations where the scatterers are close (b x -2a ≈ b y -2a ≤ a = min{a p } 1≤p≤M ), the preconditioner is not efficient and can even lead to a deterioration of the convergence. This means that, for this kind of configuration, more efforts must be done for building a suitable preconditioner.

We conclude this Section by analyzing the performance of our numerical method for an un- Figure 18 shows the convergence rate of our algortihm without preconditioner and with the preconditioner P d for d = 1.5 b min . Once again, we note a faster convergence for the preconditioned algorithm.

Conclusion and future work

In this paper, we investigated the numerical simulation of high frequency multiple scattering by circular cylinders. The main difficulty that arises in this context is due to the fact that the complex dense linear system to be solved is very large and ill-conditioned. This is in particular true when the number of scatterers is large and/or for high frequencies. Taking advantage of the particular block Toeplitz structure of the matrix of the linear system, we proposed an adapted storage of the system and an iterative algorithm of resolution, based on a fast MVP computation. We realized a thorough numerical study of the convergence rate with respect to different geometrical parameters of the problem (the reduced wavenumber, the distance between the scatterers, the number of scatterers). Finally, we proposed and tested the efficiency of a geometrically-based preconditioner, obtained by taking into account close interactions.

Considering the results obtained in this paper, we plan in the future to tackle the following questions.

• Series truncation : The heuristic formulae [START_REF] Row | Theoretical and experimental study of electromagnetic scattering by two-dimensional conducting cylinders[END_REF] defining the number of harmonics to be taken into account to obtain an accurate approximation needs to be proved rigorously. • Preconditioner : The efficiency of the preconditioner we proposed seems from our numerical experiments to be limited in some situations (very close scatterers for instance). This is due to the fact that the interactions become stronger in this case.

• Extensions : A natural continuation of our work is to investigate other scatterers with simple shapes, the three dimensional Helmholtz equation and Maxwell's equations.

1 )m

 1 (ka p )| . Taking in the above relation: |m| = ka p (1 + ζ), with ζ << 1, and using the asymptotics of Hankel functions for large arguments and large orders, one obtains that

Figure 2 :Figure 3 :

 23 Figure 2: Regular line with M = M x disks.

Figure 4 :

 4 Figure 4: Triangular lattice with M x disks on the first row, (M x -1) disks on the second one and a total number of M y odd rows.

Figure 5 :

 5 Figure 5: Behaviour of the CPU time according to the wavenumber ka p for building the matrix A in the case of the scattering by two circular cylinders fixing N p by formula (21).

Figure 6 :

 6 Figure 6: Number of iterations with respect to N p .

Figure 7 :

 7 Figure 7: Evolution of the residuals with respect to n iter for N p = 120 and N p = 150.

Figure 8 :

 8 Figure 8: Relative error with respect to N p .

Figure 9 :

 9 Figure 9: Comparison of the CPU time for computing one MVP using the direct and fast algorithms. We take: ka = 100, b x = 3 for the single-row configuration (Dirichlet problem).

Figure 10 :

 10 Figure 10: Number of MVPs with respect to the wavenumber ka: single-row configuration with 75 obstacles and b x = 3. We use different iterative solvers: GMRES, GMRES(50) and BiCGStab. The tolerance is tol = 10 -10 (Dirichlet problem).

1 Figure 11 :

 111 Figure 11: Number of MVPs with respect to the number of obstacles M : single-row configuration with ka = 100 and two different values of b x . The GMRES(50) solver is used fixing tol = 10 -8 (Dirichlet problem).

Figure 12 :

 12 Figure 12: Number of MVPs with respect to the distance δ between two obstacles. We fix: ka = 100, M = 10 obstacles. The solution is obtained with GMRES(50) for tol = 10 -8 (Dirichlet problem).

First, let us

  introduce the notation: A = I + F, where I = diag((I j ) 1≤j≤M ) is the identity diagonal block of A and F = A -I is its complement off-diagonal part. Using the first-order approximation of the Neumann series of A -1 gives A -1 ≈ I -F = P.

Figure 13 :

 13 Figure 13: Number of MVPs with respect to the wavenumber ka for the rectangular lattice configuration with b x = 8 and b y = 13. We fix M x = 8 and M y = 2, 4, 8. The solution is obtained without and with the preconditioned GMRES(50) for tol = 10 -8 (Dirichlet problem).

Figure 14 :

 14 Figure 14: Number of MVPs with respect to the wavenumber ka: rectangular lattice configurations with M x = 5, M y = 2, ..., 6. We consider b x = b y = 3 and solve the linear system using the GMRES(50) solver for tol = 10 -8 (Dirichlet problem).

Figure 15 :

 15 Figure 15: Number of MVPs with respect to M x = √ M (setting M y = M x ): rectangular lattice configuration with ka = 50 and b x = b y = 3. The GMRES(50) solver is used fixing tol = 10 -8 (Dirichlet problem).

Figure 10 .

 10 Figure 10. The preconditioner always improves the convergence even for very close scatterers. The situation is more delicate when considering the rectangular or triangular lattices. We report on Figure 13 the results obtained for different layers M y = 2, 4, 8, setting M x = 8. The distance between the obstacles is b x = 8 and b y = 13 (and so d = 13). We see an interesting gain in terms

Figure 16 :

 16 Figure 16: Number of MVPs with respect to the wavenumber ka for the single-row configuration with M = 75 and b x = 3. The solution is obtained with the preconditioned GMRES and GMRES(50) for tol = 10 -10 (Dirichlet problem).

Figure 17 :

 17 Figure 17: The considered unstructured configuration.

Figure 18 :

 18 Figure 18: Number of MVPs with respect to the wavenumber ka for a unstructured configuration (with and without preconditioner).

All the computations were performed on a Power Mac G4 1.67GHz with 1 Go DDR SDRM. The algorithms are developed under Matlab.