Gilles Chardon
email: gilles.chardon@m4x.org

Laurent Daudet

Low-complexity computation of plate eigenmodes with Vekua approximations and the Method of Particular Solutions

Keywords: Kirchhoff plate theory, biharmonic equation, numerical methods, algorithms, eigenvalues 1

This paper extends the Method of Particular Solutions (MPS) to the computation of eigenfrequencies and eigenmodes of thin plates, in the framework of the Kirchhoff-Love plate theory. Specific approximation schemes are developed, with plane waves (MPS-PW) or Fourier-Bessel functions (MPS-FB). This framework also requires a suitable formulation of the boundary conditions. Numerical tests, on two plates with various boundary conditions, demonstrate that the proposed approach provides competitive results with standard numerical schemes such as the Finite Element Method, at reduced complexity, and with large flexibility in the implementation choices.

Love plate theory. In this framework, an eigenmode is a non-zero solution to D∆ 2 u + T ∆w -ρhω 2 u = 0 [START_REF] Babuška | The h, p and h-p version of the finite element method: basis theory and applications[END_REF] with boundary conditions, where D is the rigidity of the plate, ρ the specific mass of its material, h its thickness, and T the normal tension applied at its edges, assumed to be uniform. The eigenfrequencies are the frequencies ω such that a non-zero solution exists. Apart from particular cases where these quantities can be analytically computed (e.g., circular plates with simple boundary conditions), they must be obtained by numerical methods.

The finite element method (FEM), which uses piecewise polynomials to approximate the solutions, can be used to compute these eigenmodes and eigenfrequencies. However, it can be computationally intensive at high frequencies, as the size of the numerical problem scales as the square of the spatial frequency, or in cases where singularities appear (e.g. polygonal plates). These singularities can be treated by adaptive methods, that require more complex theory and implementation [START_REF] Babuška | The h, p and h-p version of the finite element method: basis theory and applications[END_REF]. Because finite element methods look for solutions in the entire space (e.g. H 1 for the case of eigenmodes of the Laplace operator), the number of elements needed can be quite large, particularly at high frequencies.

To deal with this problem, alternative methods to the FEM have been developed, where only the solutions to the Helmholtz equation with a particular wavenumber are approximated, accelerating the convergence. Within a family (different for each wavenumber) of such solutions, eigenmodes are the linear combinations that also satisfy the boundary conditions. In the case of thin plates, one can cite the boundary element method (BEM) [START_REF] Chen | Mathematical analysis and numerical study of true and spurious eigenequations for free vibration of plates using realpart BEM[END_REF], the method of fundamental solution (MFS) [START_REF] Alves | The method of fundamental solutions applied to the calculation of eigensolutions for 2D plates[END_REF], or its variant proposed by Kang et al. [START_REF] Kang | Free vibration analysis of arbitrarily shaped plates with clamped edges using wave-type functions[END_REF], the Non Dimensional Influence Function (NDIF). Eigenmodes are identified by computing the determinant of the operator mapping the vectors of size n containing the coefficients of the expansion to the value of the function at n points on the border. A determinant equal (or near) to zero indicates a solution to Eqn. [START_REF] Babuška | The h, p and h-p version of the finite element method: basis theory and applications[END_REF] that approximately satisfies the boundary conditions, with non-zero coefficients. In practice, local minima of the function mapping the wavenumber to the determinant are assumed to be eigenfrequencies. Implementation of these methods can however be delicate, because of the appearance of spurious modes, the difficult choice of the charge points for the MFS and NDIF, or severe ill-conditioning of the numerical problems at high frequencies or with large approximation orders.

Here, we investigate a new computational method derived from the Method of Particular Solutions, proposed by Fox, Henrici and Moller (FHM) [START_REF] Fox | Approximations and bounds for eigenvalues of elliptic operators[END_REF], analyzed by Eisenstat [START_REF] Eisenstat | On the rate of convergence of the Bergman-Vekua method for the numerical solution of elliptic boundary value problems[END_REF], and further improved by Betcke and Trefethen [START_REF] Betcke | Reviving the method of particular solutions[END_REF], for the computation of eigenmodes of the Laplace operator. While the fundamental ideas used in this method are similar to the previously cited methods, it has several advantages:

the stability of the numerical problems is improved, multiple eigenvalues are easier to determine, and basis of the associated eigenspaces are readily available.

These improvements are obtained by a change of criteria for an eigenmode being non-zero : we here look at its L 2 -norm, instead of its expansion coefficients. This means, implementation-wise, sampling not only the border of the domain, but also its interior, and replacing the computation of a determinant by the search of the largest eigenvalue of a generalized eigenvalue problem.

The rest of this paper is constructed as follows. Section 2 presents the MPS framework as introduced by Fox, Henrici and Moller, and some of its more recent developments. We then extend this theory to the computation of eigenmodes and eigenvalues of plates, assumed to be star-shaped plates with smooth boundaries -due to fundamental limitations of the Vekua theory. Our first contribution is the analysis of an approximation scheme based on the Vekua theory, given in Section 3. This provides some bounds on the approximation error of a solution of Eqn. (1) by sums of Fourier-Bessel functions in Sobolev norms, based on similar results for the simpler case of the Helmholtz equation. Our second contribution, described in Section 4, is the formulation of the problem in a way compatible with the MPS, and its numerical evaluation presented in Section 5. We dis-cuss the extension to more general shapes, and various implementation matters, in Section 6.

The method of particular solutions

The method of particular solutions was introduced by Fox, Henrici and Moller (FHM) [START_REF] Fox | Approximations and bounds for eigenvalues of elliptic operators[END_REF] for the case of eigenmodes of the Laplace operator in a L-shaped domain with a singular corner.

The basic idea of this method is, instead of considering the entire space in which the eigenmodes are searched (e.g. H 1 for the case of the Laplace operator, approximated by finite element spaces), to consider separately the spaces of solutions of the Helmholtz equation for different wave numbers. Then, in each of these spaces, we can look for a nonzero function which also satisfies the boundary conditions, i.e., an eigenmode. While building approximations for a lot of different spaces seem to be counter-productive compared to the unique approximation needed for a Galerkin approximation, this alternative scheme is interesting as efficient approximations can be obtained for these spaces of solutions to the Helmholtz equation.

The method developed by FHM, for the computation of eigenmodes of the Laplace operator with Dirichlet boundary conditions, is as follows. For each wavenumber k, we consider N points x j on the border of the domain, and a family of N functions φ i spanning a subspace which approximates the set of solutions to the Helmholtz equation. The considered family, Fourier-Bessel functions of fractional orders, was specifically constructed to take into account the singularity arising in the reentrant corner of the L-shaped domain. In order to find the eigenfrequencies, one has to construct a square matrix M (k) that contains the values of the functions φ j at the N points of the border, and to compute its determinant d(k):

d(k) = det M (k) = φ 1 (x 1 ) • • • φ N (x 1 ) . . . . . . φ 1 (x N ) • • • φ N (x N ) (2)
If k is an eigenfrequency, there is a non-zero solution to the Helmholtz equation with values zero on the boundary, and its approximation α i φ i is thus close to zero at the sampling points, with nonzero coefficients α i . The image of the vector of coefficients (α i ) by the matrix M (k) is precisely the vector of the values of α i φ i on the points of the border. This means that the determinant of the matrix M (k) is close to zero. Therefore, the eigenfrequencies are obtained as local minima of d(k).

Numerous variants of this method, using different approximation schemes, have been developed since the original article. They mostly differ by the functions used to approximate the solutions: the MFS uses fundamental solutions, the NDIF [START_REF] Kang | Free vibration analysis of arbitrarily shaped plates with clamped edges using wave-type functions[END_REF] uses Bessel functions of the first kind of order 0, etc.

As pointed out in [START_REF] Betcke | Reviving the method of particular solutions[END_REF], this simple method has known limitations. The discretization of the space of solutions and the sampling of the boundaries must be of the same size, and more importantly, the matrix M (k) gets illconditioned as the number N of functions grows. Indeed, using a larger family of functions φ i , in order to have better approximations of the modes, makes the problem numerically unstable. An interpretation of this fact is that the algorithm does not search for a nonzero function inside the domain with value zero on its boundary, but actually for a function with non-zero coefficients of its expansion, with zero value on the boundary. The properties of the approximating families are such that having non-zero expansion coefficients does not ensure significant values of the function inside the domain.

To avoid these problems, Betcke and Trefethen [START_REF] Betcke | Reviving the method of particular solutions[END_REF] suggest to solve, for each k, the following optimization problem :

τ (k) = min u T u 2 L2(∂Ω) such that u 2 L2(Ω) = 1 and ∆u + k 2 u = 0 ( 3
)
where T is the trace operator on the boundary ∂Ω.

Then, k is an eigenfrequency if and only if τ (k), called the tension, is zero. This problem can be discretized as follows. A family of functions (φ i ) is chosen for the approximation of the solutions of the Helmholtz equation. For a function with expansion coefficients u = (u 1 , . . . , u N ),

T u 2 L2(∂Ω) = u Fu u 2 L2(Ω) = u
Gu where the coefficients of the matrices F and G are

F ij = T φ i , T φ j L2(∂Ω) , G ij = φ i , φ j L2(Ω) (4)
These scalar products can be estimated by sampling the domain and its border, and using a Monte-Carlo approximation of the integrals. The optimization problem (3) can be replaced by the generalized eigenvalue problem

λFu = Gu (5)
for which the largest eigenvalue is the inverse of τ (k). Note that here, the size of the matrices is the size of the approximating family, and does not depend on the number of samples used in the domain and on its boundary. The eigenfrequencies are found as the local minima of τ (k). This method offers additional advantages. The coefficients of the expansion of the eigenmodes are readily available as the first eigenvector of Eqn. ( 5), and nmultiple eigenfrequencies are characterized by the fact that the n first eigenvalues of Eqn. ( 5) exhibit a local minimum. Here, a basis of the eigenspace is obtained using the n first eigenvectors of Eqn. [START_REF] Fox | Approximations and bounds for eigenvalues of elliptic operators[END_REF]. Note that this basis has no reason to be orthogonal, as it is the vectors of the coefficients of the expansions of the basis functions that are orthogonal, not the functions themselves.

The same method, with a slightly different implementation, was used by Barnett and Berry [8] to compute high frequency modes of quantum cavities. As they considered only domains with smooth boundaries, plane wave families were sufficient for the application of the method.

To generalize this method to plates, two adaptations are needed:

an approximation scheme for solutions of Eqn. [START_REF] Babuška | The h, p and h-p version of the finite element method: basis theory and applications[END_REF] has to be developed the boundary conditions encountered in plate problems have to be modeled in a way compatible with formulation Eqn. (3).

These two points are the topics of the next two sections.

Approximation of plate eigenmodes

In this section, we prove that solutions of Eqn. (1), with arbitrary boundary conditions, can be approximated by sums of Fourier-Bessel functions and modified Fourier-Bessel functions. We first give a short account of the Vekua theory for the Laplace operator, and then extend these results to the bi-Laplace operator.

The domain where the functions are defined is assumed to be star-shaped, and to contain the ball centered on the origin of radius ρh, where h is the diameter of the domain. Furthermore, we assume that the domain satisfies the exterior cone condition with angle βπ. This means that each point of the border is the vertex of a cone of angle βπ which does not intersect the interior of the domain. Such a domain is pictured Figure 1.

In the following, we give the approximation bounds in weighted Sobolev norms defined by

u 2 m,k = m p=0 1 k 2p p1+p2=p Ω ∂ p ∂x p1 ∂y p2 u 2 dxdy (6)
Fig. 1 A star-shaped domain satisfying the exterior cone condition with angle βπ

Vekua theory for the Laplace operator

A simple example of an approximation of a solution to a differential equation is the case of holomorphic functions. These functions, solutions to the Cauchy-Riemann equations, can be, on a simply connected domain, approximated by polynomials of the complex variable. By taking the real part of a holomorphic function and of its approximations, it is shown that we can approximate a harmonic function, solution to ∆u = 0, by harmonic polynomials in R 2 . The Vekua theory, exposed in [START_REF] Vekua | New methods for solving elliptic equations[END_REF] and summarized in [START_REF] Henrici | A survey of I. N. Vekua's theory of elliptic partial differential equations with analytic coefficients[END_REF], gives similar results for solutions of elliptic PDEs by generalizing the operation "taking the real part", which allows us to map holomorphic functions to harmonic functions, to solutions of these PDEs. Using these operators, which are continuous and continuously invertible, approximation of holomorphic functions by polynomials of the complex variable is translated to approximation of solutions to the PDEs by the images of polynomials. Further details on these operators in the case of the Helmholtz equation ∆u + λu = 0 [START_REF] Betcke | Reviving the method of particular solutions[END_REF] and their properties in Sobolev spaces can be found in [START_REF] Moiola | Vekua theory for the Helmholtz operator[END_REF].

The main result of Moiola et al. is that the solutions of Eqn. [START_REF] Betcke | Reviving the method of particular solutions[END_REF] can be approximated by generalized harmonic polynomials, i.e. functions of the form

u L = L n=-L α n J n (kr)e inθ (8)
where (r, θ) are the polar coordinates and k = √ λ. When λ is strictly positive, Moiola et al. [START_REF] Moiola | Vekua theory for the Helmholtz operator[END_REF] have shown that there exists a generalized harmonic polyno-mial of order at most L such that

u -Q L j,k ≤ C(1 + (kh) j+6 )e 3(1-ρ)kh/4 log(L + 2) L + 2 β(K+1-j) u K+1,k (9)
In the case where λ < 0, the bound is multiplied by e 3kh/2 . When the function to be approximated is infinitely differentiable in an open domain containing Ω, the convergence is exponential in L [START_REF] Melenk | Operator adapted spectral element methods I: harmonic and generalized harmonic polynomials[END_REF].

Moiola et al. used these results to analyze the approximation of solutions of Eqn. ( 7) by sums of plane waves [START_REF] Moiola | Plane wave approximation of homogeneous Helmholtz solutions[END_REF]. The orders of convergence are identical to the generalized harmonic polynomials case.

Extension to plates eigenmodes

Eigenmodes of a homogenous plate of rigidity D, specific mass ρ and thickness h, with in-plane tension T , are solutions of

D∆ 2 w + T ∆w -ρhω 2 w = 0. (10)
In order to approximate such functions, we reduce this problem to the approximation of two solutions of the Helmholtz equation. Indeed, solutions of equation [START_REF] Henrici | A survey of I. N. Vekua's theory of elliptic partial differential equations with analytic coefficients[END_REF] can be decomposed as a sum of two solutions of equation [START_REF] Betcke | Reviving the method of particular solutions[END_REF], with parameters deduced from the properties of the plate.

Lemma 1 Let w be a solution to [START_REF] Henrici | A survey of I. N. Vekua's theory of elliptic partial differential equations with analytic coefficients[END_REF] in the sense of distributions. Then w can be decomposed as the sum of w 1 solution of ∆w 1 -λ 1 w 1 = 0, and w 2 solution of ∆w 2 -λ 2 w 2 = 0, where λ 1 and λ 2 are the zeros of

Dλ 2 + T λ -ρhω 2 . ( 11
)
Furthermore, if w ∈ H K+2 , then

||w 1 || K,k1 ≤ 2k 2 + δ λ ||w|| K+2,k1 (12)
||w 2 || K,k2 ≤ 2k 2 + δ λ ||w|| K+2,k2 . ( 13
)
where k 1 = |λ 1 | and k 2 = |λ 2 |, k + being the largest, and δ λ = T 2 + 4Dρhω 2 /D the difference between λ 1 and λ 2 .

Proof Analysis : assuming the decomposition, we find

∆w -λ 2 w = (∆w 1 -λ 2 w 1 ) + (∆w 2 -λ 2 w 2 ) (14) = (λ 1 -λ 2 )w 1 + 0 (15)
A similar computation for w 2 gives

w 1 = 1 λ 1 -λ 2 (∆w -λ 2 w), w 2 = 1 λ 2 -λ 1 (∆w -λ 1 w).
Note that λ 1 and λ 2 are always distinct as δ λ is always strictly positive. Synthesis : we check that w = w 1 + w 2 :

w 1 + w 2 = 1 λ 1 -λ 2 (∆w -λ 2 w) -(∆w -λ 1 w) (16) = w (17)
Then that ∆w 1 -λ 1 w 1 = 0 :

∆w 1 -λ 1 w 1 = (∆ 2 w -λ 2 ∆w) -(λ 1 ∆w -λ 1 λ 2 w) λ 1 -λ 2 (18) = ∆ 2 w -(λ 2 + λ 1 )∆w + λ 1 λ 2 w λ 1 -λ 2 (19)
= 0 (20)

The last equality comes from the fact that λ 1 and λ 2 are the zeros of the polynomial Dλ 2 + T λ -ρhω 2 . We also find that ∆w 2 -λ 2 w 2 = 0. Finally, if w ∈ H K+2 , then

w 1 K,k1 ≤ 1 |λ 1 -λ 2 | ( ∆w K,k1 + λ 2 w K,k1 ) (21)
≤ 1 δ λ k 2 1 w K+2,k1 + k 2 2 w K+2,k1 (22)
≤ 2k 2 + δ λ w K+2,k1 (23)
The result is identical for w 2 .

Remark 2 In the case where no tension is applied to the plate, i.e. T = 0, we have 2 and the results can be simplified as :

k 1 = k 2 = (ρh/D) 1/4 √ ω = k, δ λ = 2k
||w 1 || K,k ≤ ||w|| K+2,k , ||w 2 || K,k ≤ ||w|| K+2,k .
Remark 3 Two orders are lost in the majorizations (12) and (13) : the norm of order K of both components of w are bounded by their norm of order K + 2. It is impossible, in the general case, to have a better bound. Let us consider, on a disk sector centered at the origin, the function f defined by f = e iθ/2 (J 1/2 (r) -I 1/2 (r)) in polar coordinates. This function is solution of ∆ 2 w -w = 0, and can be decomposed as the sum of f 1 = e iθ/2 J 1/2 (r), solution of ∆f 1 + f 1 = 0, and f 2 = -e iθ/2 I 1/2 (r), solution of ∆f 2 -f 2 = 0. These two functions have a radial behavior at the origin similar to r 1/2 , and thus are not in H 3 . Their sum however behaves like r 5/2 , and is in H 4 .

If λ 1 , or λ 2 , is negative, the corresponding component of w can be readily approximated by generalized harmonic polynomials or plane waves using the results of Moiola et al. [START_REF] Moiola | Plane wave approximation of homogeneous Helmholtz solutions[END_REF]. If λ 1 or λ 2 is positive, the associated component can be approximated in a similar way. In that case, we use either a family of modified Fourier-Bessel functions; where the Bessel functions J n are replaced by modified Bessel functions I n , or a family of exponential functions e k•x instead of plane waves. The bounds on the approximation error are similar.

Theorem 4

Let Ω be a domain satisfying the assumptions of this section, K ≥ 1 integer, and w ∈ H K+2 verifying conditions of lemma 1. Then for all L > K, there exist two generalized harmonic polynomials P L and Q L with parameters λ 1 and λ 2 , of degree at most L such that for all j ≤ K, w -

(P L + Q L ) j,k+ ≤ C k 2 + δ λ (1 + (k + h) j+6 )e 3 4 (3-ρ)k+h ln(L + 2) L + 2 β(K-j) (k + h) K-j w K+2,k-, (24)
where k -is the smallest of k 1 and k 2 , and k + the largest, and C is a constant which depends only on the shape of Ω, j and K. If λ 1 and λ 2 are both negative, the bound can be divided by e We have

w -(P L + Q L ) j,k+ ≤ w 1 j,k+ + w 2 j,k+ ≤ w 1 j,k1 + w 2 j,k2 ≤ C(1 + (k 1 h) j+6 )e 3 4 (1-ρ)k1h ln(L + 2) L + 2 β(K-j) (k 1 h) K-j w 1 K,k1 + C(1 + (k 2 h) j+6 )e 3 4 (3-ρ)k2h ln(L + 2) L + 2 β(K-j) (k 2 h) K-j w 1 K,k2 ≤ C k 2 + δ λ (1 + (k + h) j+6 )e 3 4 (3-ρ)k+h ln(L + 2) L + 2 β(K-j) (k + h) K-j w K+2,k- Remark 5
For plates without in-plane tension, the result is slightly simpler. In that case, the roots of Eqn. [START_REF] Moiola | Vekua theory for the Helmholtz operator[END_REF] have same absolute values and opposite signs. We thus have

w -(P L + Q L ) j,k ≤ C(1 + (kh) j+6 )e 3 4 (3-ρ)kh ln(L + 2) L + 2 β(K-j) (kh) K-j w K+2,k . (25)
Remark 6 We prove the approximation result for sums of Fourier-Bessel functions. However, as Fourier-Bessel functions can be approximated by sums of plane waves, modified Fourier-Bessel functions can be approximated by sums of exponential function of the type e k•x with constant norm of k. This allows an approximation of solutions of Eqn. (10) by sums of plane waves and exponential functions.

Boundary conditions

Boundary conditions usually encountered in plate problems (clamped edges, simply supported edges and free edges) are not as such readily usable for the numerical scheme proposed here, and thus need to be modeled differently.

For clamped edges, the displacement and its normal derivative are zero. The tension for clamped boundary conditions is then:

t c = Γ |w| 2 + 1 k 2 Γ ∂ w ∂n 2 .
A simply supported edge has zero displacement and torsion moment

M n t s = Γ |w| 2 + 1 D 2 k 4 Γ |M n | 2 .
Finally, free edges have zero torsion moment and Kelvin-Kirchhoff edge reaction

K n t l = 1 D 2 k 4 Γ |M n | 2 + 1 D 2 k 6 Γ |K n | 2 .
The torsion moment and the Kelvin-Kirchhoff edge reaction writes [START_REF] Géradin | Mechanical vibrations: theory and application to structural dynamics[END_REF] M

n = -D ∂ 2 w ∂n 2 + ν ∂ 2 w ∂t 2 K n = -D ∂ 3 w ∂n 3 + (2 -ν) ∂ 3 w ∂n∂t 2 + 1 -ν R ∂ 2 w ∂n 2 - ∂ 2 w ∂t 2
where ν is the Poisson ratio of the material, and R the curvature radius of the boundary.

The constants 1/k 2 , 1/D 2 k 4 and 1/D 2 k 6 in front of some integrals not only make the quantities homogeneous, but also improve the numerical stability, as the contribution of the two integrals are rescaled to have the same order of magnitude. For instance, in the case of t c estimated using plane waves, the first term contains products of planes waves, while the second term contains products of derivatives of plane waves, which are the plane waves themselves multiplied by the scalar product of the wave vector and a unit vector normal to the boundary. With no rescaling, the second term would be more and more influent as the frequency increases. This rescaling is similar to the weights used to define the norms • m,k in Eqn. [START_REF] Eisenstat | On the rate of convergence of the Bergman-Vekua method for the numerical solution of elliptic boundary value problems[END_REF].

Note that for a plate with various boundary conditions along the border, the corresponding tension is the sum of the tension for the different boundary conditions, integrated on their respective domain.

Numerical results

We now compare the methods with analytical results for simple cases, and with numerical results obtained with Cast3M [16], a widely-used FEM simulation program. To avoid technicalities, the simulated plates are starshaped with smooth boundaries. The treatment of more general shapes is discussed in the next section. Since these numerical tests are done without in-plane tension, we can use the wavenumber k = (ρ/D) 1/4 √ ω to express the eigenfrequencies. In the case of in-plane tension, the wavenumbers used to generate the Fourier-Bessel functions (resp. plane waves) and modified Fourier-Bessel functions (resp. exponential functions) are computed according to Lemma 1.

We first test the method on a circular plate of radius 1, with various boundary conditions. Here, we use plane waves, as eigenmodes of circular plates are sums of a Fourier-Bessel function and a modified Fourier-Bessel function, making the numerical problem trivial. Table 1 gives the computed eigenfrequencies and the values given in Leissa [START_REF] Leissa | Vibration of Plates[END_REF] for low-frequency modes. Figure 2 shows the tension as a function of the frequency in the clamped case, and the next three eigenvalues of the numerical problem [START_REF] Fox | Approximations and bounds for eigenvalues of elliptic operators[END_REF]. These can be used to identify multiple eigenvalues and to compute a basis of the associated eigenspaces. Examples of eigenmodes are given for the three boundary conditions in Figure 3. In these numerical experiments, the boundary of the disk was discretized with 2048 points, its interior with 1024 points drawn using the uniform distribution on Ω, and the number of plane waves (which is also the size of the numerical problems to solve) was taken as 8k, which in this case ranges from 28 to 80. We now compute eigenfrequencies and eigenmodes of a second plate, of more complex shape, with bound-aries defined by the parametric equations

x = cos t y = sin t + 1 3 sin 2t t ∈ [0, 2π) (26)
Table 2 gives the first ten eigenfrequencies of this plate with the three types of boundary conditions : clamped, simply supported and free. The proposed method, MPS with plane waves (MPS-FB), uses 2048 points on the border and 1024 points inside. The number of plane waves is 10k, which here ranges from 20 to 80. For the clamped conditions, the results using Fourier-Bessel functions are also given (method MPS-FB). We compare our results with those obtained by the FEM, as implemented in the Cast3M package using a triangular mesh and finite elements of order 3, with 6926 elements and the border discretized by 180 segments.

For two eigenmodes of the clamped plate, we give in table 3 the estimated eigenfrequencies with varying size of the discrete problems for FEM and the MPS-PW, and compare them to the values obtained in [START_REF] Alves | The method of fundamental solutions applied to the calculation of eigensolutions for 2D plates[END_REF] with the MFS (with an unspecified size). These two modes, computed with the MPS-PW, are shown on Figure 4.

Discussion

We now discuss some implementation details, relative to the treatment of more complex shapes, or the acceleration of the computations, and the relative merits of the MPS and other methods.

Fourier-Bessel vs. plane waves

As shown by the numerical experiments, both Fourier-Bessel functions and plane waves can be used to approximate eigenmodes. They have similar approximation properties, but differ implementation-wise. Fourier-Bessel functions are orthogonal on a disc, ensuring better stability, while plane waves are more and more ill-conditioned as their number increases. However, this can be treated by pre-conditioning the plane waves family with a discrete Fourier Transform, mapping the plane waves to approximations of the Fourier-Bessel functions.

The main advantage of the plane waves is the straightforward computation of their derivatives : differentiating a plane wave along a certain direction amounts to multiplying it with the scalar product of its wave vector with a unit vector. This makes the construction of the matrices both easy to implement and fast.

Shapes

The proposed method relies on an approximation scheme for solutions to the equation ( 1), and therefore is limited to cases where such approximations are available. In particular, the approximation scheme assumes a starshaped domain. In the case of a simply connected, but non-star convex, domain, the approximation is not guaranteed to succeed. A possible way to overcome the problem is to cut the domain into star-convex subdomains, to approximate the solutions of (1) in these subdomain, and to add terms in the tension ensuring that the displacement, normal derivative of the displacement, bending moment and strain have the same value at both sides of the internal boundaries. A similar method as already been applied to the particular case of polygonal membranes [START_REF] Descloux | An accurate algorithm for computing the eigenvalues of a polygonal membrane[END_REF].

Other cases of non-convex domains are domains with holes, but star-convex if the holes are filled. In that case, following Vekua, one can approximate solutions of (1) by adding to the family of planes waves or Fourier-Bessel functions, the sets of Fourier-Bessel functions of second kind e inθ Y n (kr) and e inθ K n (kr), one for each hole, centered on a point chosen in each of them. This type of approximation can be compared to the approximation of holomorphic functions given by the Runge theorem. An application to membranes can be found in [START_REF] Betcke | Numerical computation of eigenfunctions of planar regions[END_REF].

Singularities

The domains considered here have smooth boundaries. This guarantees that the eigenmodes of the plates are 3 Eigenfrequencies for two modes of the clamped second plate, with various discretization sizes, computed by FEM (as implemented in Cast3M), MFS (results of [START_REF] Alves | The method of fundamental solutions applied to the calculation of eigensolutions for 2D plates[END_REF] on the same plate), and MPS-PW, with various number of plane waves.

smooth, and that the convergence of the approximations is fast. However, as shown in Theorem 4, singularities, which can appear at corners of a domain with non-smooth boundaries, slow down the convergence of the approximations. The size of the discretization for such cases is then larger than for smooth boundaries, possibly too large to guarantee the numerical stability of the computations. To accelerate the convergence, FHM, Eisenstat and BT use fractional Fourier-Bessel functions centered on the singular corner of the domain (in their numerical experiments, an L-shaped polygon). This idea has been used for plates in a different setting by De Smet et al. [START_REF] Vanmaele | An efficient wave based prediction technique for dynamic plate bending problems with corner stress singularities[END_REF].

Note that, for a plate with polygonal holes, at least a corner of the hole is singular. In that case, it is impossible to use fractional Fourier-Bessel functions, as it is impossible to define them on a domain containing a path around the origin. In that case, one can combine fractional Fourier-Bessel function with the method described in the previous section [START_REF] Betcke | Numerical computation of eigenfunctions of planar regions[END_REF].

Numerical considerations

Although the numerical stability is improved compared to the determinant based MPS, the improved version is still prone to instabilities when a large approximation order is used. These instabilities are further amplified in the case of plates, where modified Fourier-Bessel functions, or exponential functions, are included in the approximating families. The behavior of these functions are such that they are non-negligible only on a small region near the boundary of the domain. Using the Monte-Carlo approximation with uniform density to estimate the coefficients of the matrices is thus unstable. Using a non-uniform density of samples, with more samples near the border would be a way to improve the stability of the estimations, in a way similar to what is used in [START_REF] Chardon | Approximation of solutions to the Helmholtz equation from scattered data[END_REF], where the reconstruction of a solution of the Helmholtz equation on a disc is improved by placing a fraction of the samples on the border of the disc.

Speeding up the eigenvalue search

In order to locate the minima of the tension, the proposed algorithm simply computed it on linearly spaced values in the interval we were interested in. However, the particular behavior of the tension (and more generally, of the eigenvalues of problem ( 5)) could be used to accelerate the search of these minima. Indeed, the tension is a sequence of branches which behave more or less as parabolas. Newton iterations, along with the computation of the derivatives of the eigenvalues of problem (5), could therefore be used to quickly locate the minima.

Comparison with other methods

The simulations (Table 3) show that the MPS method needs a significantly lower order of approximation than the FEM to obtain accurate results. In particular, the MPS using 50 Fourier-Bessel functions achieve a much better estimation of the eigenfrequencies that the FEM with 56 elements. Indeed, while the FEM needs elements of size comparable with the tenth of the wavelength (and thus an order scaling like the square of the wavenumber), the MPS needs an order proportional to the wavenumber [START_REF] Barnett | Dissipation in Deforming Chaotic Billiards[END_REF], because only the solutions to Eqn. (1) are approximated.

The computation times for the FEM and the MPS to obtain Table 2 were similar. However the MPS was here implemented using Octave and Matlab, without the improvements given in the previous subsection, and significants gain in term of complexity can be expected. The implementation of the MPS is also simpler, as no meshing is required.

The MFS, like the MPS, needs low orders to achieve correct results. However, to use the MFS (or the NDIF), one has to choose a set of charge points. The accuracy of the method can vary with this choice, and while methods have been suggested for this task [START_REF] Alves | The method of fundamental solutions applied to the calculation of eigensolutions for 2D plates[END_REF], an optimal and efficient way is not yet known. Here, in the case of the plane wave approximation, no choice is needed. Furthermore, the use of Bessel functions makes the implementation delicate, as pointed out in Subsection 6.1.

Singularities can be treated with FEM methods: using adaptive methods, such as hp-FEM [START_REF] Babuška | The h, p and h-p version of the finite element method: basis theory and applications[END_REF], a fast convergence in presence of singularities can be achieved. The treatment of singularities is however simpler using in the MPS, as the singular functions to be used can be determined a priori directly from the geometry of the problem.

The solutions given by the MPS are smooth, as they are finite sums of smooth functions. Galerkin methods, such as the Element Free Galerkin method, have also been developed to improve the regularity of the solutions, which is limited in the FEM by the regularity of the elements. However, as pointed out in the case of the FEM, the fact that these methods approximate a functional space not limited to the solutions to Eqn. (1) leads to a large order of approximation to achieve good results.

While the MPS has significant advantages compared to other methods, the need of a family approximating the solutions to Eqn. (1) restricts its use to cases where such a family is known. Problems with plates of heterogeneous materials or thickness can be difficult or impossible to treat using this method.

Conclusion

This paper has described the extension of the Method of Particular Solutions for the computation of eigenmodes of plates. This method has numerous advantages. It can be used with any approximation scheme for the solutions of the studied equation ; in this paper, we used Fourier-Bessel functions and plane waves, but the method could be extended with fractional Fourier-Bessel functions in order to treat singularities. Its formulation offers a large flexibility in the discretization of the domain, and independently, in the size of the numerical problem. The determination of multiple eigenvalues and eigenmodes is also straightforward. Finally, the so-called tension, that has to be minimized to find the eigenfrequencies, has a specific shape that can be used to speed up the search. Future improvements include sampling schemes yielding better stability of the numerical problems, and more efficient computational methods for high frequency eigenmodes.

Reproducible research

The Matlab/Octave code used to compute eigenfrequencies and eigenmodes of Figures 3 and4 is available online [22].

Fig. 2 Fig. 4

 24 Fig. 2 Inverses of the four largest eigenvalues (the first one being the tension) of discrete problem (5), for the clamped circular plate

Table 2

 2 First ten eigenfrequencies, expressed in wavenumber k, of the second plate with various boundary conditions, for FEM, MPS with plane waves (MPS-PW) and MPS with Fourier-Bessel functions (MPS-FB, clamped conditions only)

			Clamped		Simply supp.	Free
		FEM	MPS-PW MPS-FB	FEM	MPS-PW	FEM	MPS-PW
	1	3.3786	3.3782	3.3562	2.4104	2.4061	2.1933	2.2000
	2	4.5937	4.5944	4.5864	3.6727	3.6733	2.3516	2.4001
	3	4.6267	4.9275	4.8815	4.0065	3.9993	2.7468	2.7481
	4	5.8359	5.8376	5.8366	4.9500	4.9505	3.4225	3.4222
	5	6.2048	6.2067	6.1987	5.3203	5.3216	3.7624	3.7623
	6	6.4551	6.4567	6.4527	5.5451	5.5456	4.0996	4.0993
	7	7.0716	7.0738	7.0738	6.2134	6.2147	4.2007	4.2004
	8	7.4914	7.4949	7.4949	6.6169	6.6188	4.7986	4.7985
	9	7.7709	7.7740	7.7590	6.9041	6.9058	4.9847	4.9845
	10 7.9676	7.9300	7.9660	7.0629	7.0308	5.4319	5.4316
			FEM			MFS	MPS-PW
		number of elements		number of plane waves
	56	242	1048	6926	28171	50	60	70
	8.7505 8.9527 9.0714 9.1066 9.1111 9.11259 9.1068 9.1121 9.1126
	9.0159 9.1555 9.2438 9.2735 9.2777 9.27903 9.2651 9.2722 9.2787
	Table

Acknowledgement

The authors acknowledge partial support from Agence Nationale de la Recherche (ANR), project ECHANGE (ANR-08-EMER-006), project LABEX WIFI (ANR-10-IDEX-0001-02 PSL*), and Austrian Science Fund (FWF) START-project FLAME (Y 551-N1).