F Castella

P Chartier

S Descombes
email: sdescomb@unice.fr

G Vilmart
email: gilles.vilmart@unige.ch

E Hansen

A Ostermann

] F Castella

Splitting methods with complex times for parabolic equations

Keywords: splitting, complex time steps, composition method, higher order, parabolic equations Mathematics Subject Classification (2000) 65M12, 47D06

Using composition procedures, we build up high order splitting methods to solve evolution equations posed in finite or infinite dimensional spaces. Since high-order splitting methods with real time are known to involve large and/or negative time steps, which destabilizes the overall procedure, the key point of our analysis is, we develop splitting methods that use complex time steps having positive real part: going to the complex plane allows to considerably increase the accuracy, while keeping small time steps; on the other hand, restricting our attention to time steps with positive real part makes our methods more stable, and in particular well adapted in the case when the considered evolution equation involves unbounded operators in infinite dimensional spaces, like parabolic (diffusion) equations.

We provide a thorough analysis in the case of linear equations posed in general Banach spaces. We also numerically investigate the nonlinear situation. We illustrate our results in the case of (linear and nonlinear) parabolic equations.

Introduction

The goal of the present text is to derive high-order splitting methods, obtained by using complex time steps. These methods are obtained through composition procedures. For stability purposes, the retained methods only involve time steps that have positive real part: our motivation is to recover methods which can be used in the case when unbounded operators are involved, associated with propagators that are C 0 semi-groups only (instead of C 0 groups). Our paradigm is the case of diffusion equations.

Let us make our statement precise. Consider a linear evolution equation of the form d dt u(t) = A u(t) + B u(t), (1.1) where the right-hand-side involves the sum of two well identified operators A and B. Here, the unknown u(t) is assumed to belong to some finite or infinite dimensional Banach space X, while A and B are linear, possibly unbounded, operators. In the case when A and B are unbounded, we assume A, B, and A + B generate C 0 semi-groups of propagators over X, denoted by e tA , e tB and e t (A+B) , respectively, whenever t ≥ 0. The prototype we have in mind is the linear heat equation with potential

∂ t u(t, x) = ∆ u(t, x) + V (x)u(t, x),
where t ≥ 0 is time and x ∈ R d (or x ∈ T d ) is the space variable, ∆ denotes the standard Laplacian in x, while the potential V (x) is assumed bounded. In that case one may choose X = L 2 (R d ) and define the unbounded operator (A u)(x) = ∆ u(x) with domain D(A) = H 2 (R d ), together with the bounded operator (B u)(x) = V (x)u(x) (other choices are obviously possible for the Banach space X at this level, depending on the assumed smoothness of the potential: X may as well measure Sobolev regularity, or Hölder regularity of the solution u).

It is well-known that a possible approach to numerically solve (1.1) is to use a splitting method, i.e. to approximate the propagator of the full operator u → A u + B u by using an appropriate combination of the propagators u → A u and u → B u, both assumed to be numerically cheaper to evaluate. In the finite dimensional setting, splitting methods basically rely on the identity e h(A+B) = e hA e hB + O h 2 , where h is some small time-step. Higher order approximations may be obtained by writing e h(A+B) = e b 1 hA e a 1 hB e b 2 hA e a 2 hB . . . e b s hA e a s hB + O h r+1 , where a 1 , . . ., a s and b 1 , . . ., b s are (to be chosen) real or complex numbers, and s is usually referred to as the number of steps of the method. The exponent r depends on the chosen values of the a i 's and b i 's.

The above procedure immediately extends to the case when the operators A and B become nonlinear. In this case indeed, the above formulae remain unchanged, provided the factors e h (A+B) , e a i hA , and e b i hB are replaced by the true flows Φ A+B (h), Φ A (a i h), and Φ B (b i h) respectively, or by appropriate approximations of them. Here we have defined, for any y 0 , the flow Φ A (t)(y 0 ) as the solution to the differential equation

d dt (Φ A (t)(y 0 )) = A (Φ A (t)(y 0 )) ,
supplemented with the initial condition Φ A (0)(y 0 ) = y 0 (and the similar definition is used to prescribe the flows Φ A+B (t) and Φ B (t)). In this context however, it needs to be assumed that the vector fields A and B possess enough smoothness to have well-defined flows Φ A+B (h), Φ A (h), and Φ B (h) for small values of h. Note also that if a i is complex, the definition of Φ A (a i h) requires, say, that for any y 0 the differential equation

d dt (Φ A (t)(y 0 )) = A (Φ A (t)(y 0 )) be solvable along the complex line t = a i t ′ (t ′ ∈ R) whenever 0 ≤ t ′ ≤ h
and h is small enough (and similarly for B).

The formal extension of all above formulae in the infinite dimensional setting is easy as well, keeping in mind that the existence and well definiteness of all involved propagators over the retained Banach space X should then be carefully checked. In the paradigmatic case when A = ∆ , for instance, we recall that the propagator e z∆ (z ∈ C) is well-defined, in any reasonable distribution sense, if and only if Re(z) ≥ 0. Naturally, another key difficulty in the infinite dimensional situation is to check that the remainder terms indeed have size O(h 2 ) resp. O(h r+1 ) in the correct norm. Now, the derivation of high-order splitting methods is not straightforward in general, even in the finite dimensional case. The simplest, high-order splitting methods involve large negative time steps alternating with large positive time steps (i.e. large positive values of the a i 's or b i 's alternating with large negative values of the same coefficients), which eventually leads to poor accuracy in practice. Even more, a disappointing result shows that all splitting methods (or composition methods -see below for the definition) with real coefficients must have some negative coefficients a i and b i in order to achieve order 3 or more.

The existence of at least one negative coefficient was shown in [START_REF] Sheng | Solving linear partial differential equations by exponential splitting[END_REF][START_REF] Sussman | Chaotic evolution of the solar system[END_REF], and the existence of a negative coefficient for both operators was proved in [START_REF] Goldman | Nth-order operator splitting schemes and nonreversible systems[END_REF]. An elegant geometric proof can also be found in [START_REF] Blanes | On the necessity of negative coefficients for operator splitting schemes of order higher than two[END_REF]. As a consequence, such high-order splitting methods cannot be used in general when one operator A or B has large negative spectrum, or when it only generates a C 0 semi-group of propagators -and not a group (like the Laplacian).

In order to circumvent this order-barrier, there are two possibilities. One can use a linear, convex combination (see [START_REF] Gegechkori | High degree precision decomposition method for the evolution problem with an operator under a split form[END_REF][START_REF] Gegechkori | The fourth order accuracy decomposition scheme for an evolution problem[END_REF][START_REF] Bandrauk | Complex integration steps in decomposition of quantum exponential evolution operators[END_REF] for methods of orders 3 and 4) or non-convex combination (see [START_REF] Schatzman | Numerical integration of reaction-diffusion systems, Numerical methods for ordinary differential equations (Auckland[END_REF][START_REF] Descombes | Convergence of a splitting method of high order for reaction-diffusion systems[END_REF] where an extrapolation procedure is exploited), of elementary, low order splitting methods (some of the above mentioned works use elementary methods which involve one or two complex time steps). Another possibility is to systematically consider splitting methods with complex coefficients a i and b i , having yet positive real parts (see [START_REF] Chambers | Symplectic integrators with complex time steps[END_REF] in celestial mechanics). In 1962/1963, Rosenbrock [START_REF] Rosenbrock | Some general implicit processes for the numerical solution of differential equations[END_REF] considered complex coefficients in a similar context. We may also quote the text [START_REF] Suzuki | General theory of fractal path integrals with applications to many-body theories and statistical physics[END_REF] -see also [START_REF] Bandrauk | Complex integration steps in decomposition of quantum exponential evolution operators[END_REF] -, where some low order methods with complex coefficients are derived (one can find here an alternative proof of the existence of negative coefficients when only real time steps are allowed). This is the route we chose here.

In this article, we consider splitting methods of the form

e h(A+B) = s ∏ i=1 e b i hA e a i hB + O h r+1 ,
and we derive new high-order splitting methods (up to order fourteen), which involve complex time steps having positive real part. We state and prove error estimates that are valid both in the finite and in the infinite dimensional setting. We last investigate numerically the behaviour of the retained methods both in the case of the linear heat equation with bounded potential (the setting is one-dimensional with periodic boundary conditions), and in the case of nonlinear versions of the heat equation (in the similar setting).

Our derivation uses composition techniques that were originally developed for the geometric numerical integration of ordinary differential equations [START_REF] Hairer | Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF].

The main advantages of our approach are the following:

-the splitting method inherits the stability property of exponential operators; -in the retained methods, we can always replace the costly exponentials of the operators e hA etc. by cheap low order approximations of the latter, without altering the overall order of accuracy; -using complex coefficients allows to reduce the number of compositions needed to achieve any given order;

This paper is organized as follows. In Sect. 2, we derive new high-order splitting methods. In Sect. 3 we give a rigorous order estimate in the linear case, obtained as a direct consequence of the recent results by Hansen & Ostermann [START_REF] Hansen | Exponential splitting for unbounded operators[END_REF]. Sect. 4 presents several numerical simulations, confirming the formally expected order of accuracy in the non-linear case.

Deriving high order splitting/composition methods

Composition methods in the finite dimensional case

Composition methods were mainly developed in the 90's in the papers of Suzuki [START_REF] Suzuki | Fractal decomposition of exponential operators with applications to many-body theories and Monte Carlo simulations[END_REF], Yoshida [START_REF] Yoshida | Construction of higher order symplectic integrators[END_REF] and McLachlan [START_REF] Mclachlan | Composition methods in the presence of small parameters[END_REF] in the context of ordinary differential equations. They rely on the following observation.

Consider a (linear or nonlinear) ordinary differential equation of the form

d dt u(t) = f (u(t))
where u(t) ∈ R d belongs to some finite dimensional space. Denote by φ (t) the flow of the above equation, namely φ (t) :

R d → R d satisfies u(t) = φ (t) (u(0)) .
On the other hand take a time step h > 0, and consider any approximation φ h of φ (h) at order p, for some value of p ∈ N. In other words, φ h : R d → R d is assumed to satisfy

φ h = φ (h) + O(h p+1 ).
The above identity is assumed to hold between mappings on R d . It means that whenever K ⊂ R d is a compact set, there is a constant C > 0 and a small h 0 > 0, such that for any u ∈ K, and any 0 < h < h 0 , we have φ h (u) -φ (h)(u) ≤ C h p+1 , Lastly, take an integer s and choose (real or complex) coefficients γ 1 , . . . , γ s (in the classical theory, only real coefficients were considered).

Under these circumstances, a composition method is defined as the operator

ψ h = φ γ s h • . . . • φ γ 1 h , (2.1)
i.e. as the composition of the method φ h , successively used with time steps γ 1 h, γ 2 h, . . . γ s h. Naturally, if the γ i 's are complex, we implicitly assume here that the operators φ γ i h are well-defined for small values of h > 0, in the following sense: for any compact set K ⊂ R d , there is an h 0 > 0 such that the operator φ γ i h is welldefined over K whenever 0 < h < h 0 . Similarly, we also assume that the operators φ (γ i h) are well-defined for small values of h > 0, and that the estimates φ γ i h = φ (γ i h) + O(h p+1 ) hold for small h > 0 as well (in the above sense).

The following classical result in numerical integration illustrates that the composition procedure allows to transform a method φ h of order p, into a higher-order method ψ h of order p + 1.

Theorem 2.1 (see [START_REF] Hairer | Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF]Theorem II.4.1]) With the above notation and assumptions, let φ h be an approxi- mation of φ (h) of order p, namely

φ h = φ (h) + O(h p+1 ).
If the γ i 's satisfy γ 1 + . . . + γ s = 1 and γ p+1 1

+ . . . + γ p+1 s = 0, (2.2)
then the composition method ψ h = φ γ s h • . . . • φ γ 1 h approximates φ (h) at order p + 1, i.e.

ψ h = φ (h) + O(h p+2 ).
Remark 2.1 Whenever p is even and the composition is symmetric (i.e. γ s-i+1 = γ i for any i), then ψ h is of order p + 2.

Proof.

The idea of proof is to show that if the basic method has order p, i.e.

φ h (y) = φ (h)(y) + C(y)h p+1 + O(h p+2 ),
where φ (h) denotes the exact flow, then, using the fact that the sum of the γ i 's is one, we have

φ γ s h • . . . • φ γ 1 h (y) = φ (h)(y) + C(y)(γ p+1 1 + . . . + γ p+1 s )h p+1 + O(h p+2 ).
Here, the constant C(y) denotes a quantity that remains bounded whenever y belongs to a given compact set. The result follows.

⊓ ⊔

Given the above theorem, a classical idea is the following. Starting from a low order method φ h , we may increase the order by one, by appropriately choosing the γ i 's ; iterating the process, and choosing therefore possibly different γ i 's at each stage of the whole procedure, we may eventually come up with a high-order method.

This is the program we intend to follow, in the very case of splitting methods. The point is, such a program fails past order 2 when the γ i 's are restricted to only take real values: past order 2 indeed, negative γ i 's, as well as large γ i 's come up in the analysis, which makes the so-obtained methods have poor accuracy in practice.

We therefore rely on the use of complex γ i 's. In that perspective, our main constraint is to obtain high- order method for which the γ i 's all have positive real part: our goal is to eventually apply the methods in the case of diffusion equations. Secondarily, we try to keep the number of stages (the integer s in formula (2.1)) reasonably small (to reduce the computational cost), and the moduli |γ i | as small as possible as well (to reduce the size of the time steps). Lastly, we also try to keep the quantities |arg(γ i )| as small as possible. In the context of splitting methods, all these constraints are fairly natural, since the reader should keep in mind that one step of the full evolution equation ẏ = A(y) + B(y), along a time step h, is here approximated by s steps of either ẏ = A(y) or ẏ = B(y), along the successive time steps γ 1 h, . . ., γ s h (see e.g. figure 2.1 below).

A last, important remark is in order. Throughout this text, we will restrict our attention to the case of symmetric methods. The reason for this choice is, as stated in Remark 2.1, that symmetric composition methods applied to symmetric procedures allow to gain two orders of accuracy each time one applies Theorem 2.1. We nevertheless stress that this choice is arbitrary, and considering non-symmetric methods is relevant as well. For instance, we may quote the work of Hansen and Ostermann [START_REF] Hansen | High order splitting methods for analytic semigroups exist[END_REF] where non-symmetric methods are considered.

2.2 Buliding up high-order splitting/composition methods with complex coefficients -the linear, finite-dimensional case Throughout this paragraph, we take fixed matrices A and B acting on R d , and introduce new splitting/composition methods to solve the linear ODE ẏ = Ay + By.

Though A and B are finite dimensional matrices here, the reader may keep in mind that we will eventually consider the infinite dimensional situation where A is the Laplace operator, while B denotes the multiplication by a bounded potential V (see introduction). Hence A may be typically thought of as a matrix with 'large' negative eigenvalues.

Following the general methodology described in the previous paragraph, we first need to choose some low-order approximation of the true propagator exp (h(A + B)). We retain the simplest symmetric splitting algorithm, namely the Strang splitting operator, and we set

Φ h = exp hB 2 exp(hA) exp hB 2 , (2.3)
a symmetric second order approximation of exp (h(A + B)). While the methods we propose below are all based on this particular choice of a basic low order method, we readily mention that the analysis we provide remains unchanged when starting from any other symmetric second order method. For instance, the following

Φ P h = Id - h 2 B -1 implicit Euler Id - h 2 A -1 Id + h 2 A implicit midpoint Id + h 2 B .
explicit Euler (2.4) would provide such a basic choice. In the infinite dimensional setting and when formally choosing A as the Laplace operator while B is the multiplication by the bounded potential V , the method Φ P h coincides with the Peaceman-Rachford formula [START_REF] Peaceman | The numerical solution of parabolic and elliptic differential equations[END_REF] originally developed for the heat equation, and extended to reactiondiffusion problems in [START_REF] Descombes | Convergence of the Peaceman-Rachford approximation for reaction-diffusion systems[END_REF]. Note that the use of an implicit midpoint approximation for the operator A corresponds to a standard Crank-Nicolson scheme when A is the Laplace operator, a standard choice.

Triple Jump procedures

Starting from the basic, second order, Strang splitting algorithm Φ h , we wish to derive various higher order symmetric methods by applying Theorem 2.1. Since the symmetry requirement anyhow prevents the choice s = 2 in Theorem 2.1 (system (2.2) for s = 2 and symmetric γ i 's imposes γ 1 = γ 2 = 0), the simplest choice is to set s = 3 in Theorem 2.1 and to look for a three steps, symmetric composition procedure. Such a method is usually called 'Triple jump composition procedure'. Note in passing that double jump composition procedures (with complex time steps) have been considered in [START_REF] Hansen | Exponential splitting for unbounded operators[END_REF], where methods of order 6 are derived.

In the case when s = 3, and starting from an arbitrary, symmetric method φ h of order p (p is an arbitrary even integer here), a triple jump composition procedure is provided by a triple (γ 1,p , γ 2,p , γ 3,p ) ∈ C 3 such that γ 1,p = γ 3,p (symmetry),

γ 1,p + γ 2,p + γ 3,p = 1, γ p+1 1,p + γ p+1 2,p + γ p+1 3,p = 0 (Theorem 2.1), (2.5)
and the associated improved method is ψ h = φ γ 3,p h • φ γ 2,p h • φ γ 1,p h . The set of all complex solutions to (2.5) is given by the p + 1 values

γ (k) 1,p = γ (k) 3,p = 1 2 -2 1/(p+1) e 2ikπ/(p+1) , γ (k) 2,p = - 2 1/(p+1) e 2ikπ/(p+1) 2 -2 1/(p+1) e 2ikπ/(p+1) (k = 0, . . . , p) . (2.6)
Setting k = 0 above provides the unique real solution to (2.5), namely γ (0)

1,p = γ (0) 3,p = 1 2 -2 1/(p+1) , γ (0)
2,p = -

2 1/(p+1) 2 -2 1/(p+1) .
(2.7)

Unsurprisingly, the coefficient γ

(0)
2,p is negative. In any circumstance, if φ h has order 2, then the method ψ (4)

h := φ γ (0) 3,2 h • φ γ (0) 2,2 h • φ γ (0) 1,2 h
reaches order 4 using three steps of φ h , while, repeating the method at any order, the method ψ (6)

h := ψ (4) γ (0) 3,4 h • ψ (4) γ (0) 2,4 h • ψ (4) γ (0) 1,4 h
reaches order 6 with 9 steps of φ h , and so on. These methods are originally due to Creutz & Gocksch [START_REF] Creutz | Higher-order hybrid Monte Carlo algorithms[END_REF], Forest [START_REF] Forest | Canonical integrators as tracking codes[END_REF], Suzuki [START_REF] Suzuki | Fractal decomposition of exponential operators with applications to many-body theories and Monte Carlo simulations[END_REF], Yoshida [START_REF] Yoshida | Construction of higher order symplectic integrators[END_REF], the name 'Triple Jump composition methods' was given in [START_REF] Hairer | Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF]Example II.4.2]. However, since γ (0) 2,p < 0, these methods cannot be extended, in the infinite dimensional situation, to problems where A only generates a C 0 semi-group of propagators, and where the basic method of choice φ h coincides with the Strang splitting method Φ h , or with any low order splitting method. On top of that, and even in the finite dimensional setting, the estimate γ (0) j,p > 1, valid for any j = 1, 2, 3, implies a terrible zig-zag in the coefficients of the methods. Thus, the above technique is not even very efficient in the finite dimensional case.

Another choice of k is therefore in order in (2.6). Setting k = ±p/2 (recall that p is even), provides the two conjugate solutions of (2.5) which minimize the quantity γ (k)

1,p + γ (k) 2,p + γ (k)
3,p . These two conjugate solutions also minimize the quantity max i=1,2,3 arg(γ

(k) i,p ) .
In order to keep notation simple, we drop the upper index "±p/2" for the associated coefficients, and simply define (γ 1,p , γ 2,p , γ 3,p ) as γ 1,p = γ 3,p = e iπ/(p+1) 2e iπ/(p+1) + 2 1/(p+1) , γ 2,p = 2 1/(p+1) 2e iπ/(p+1) + 2 1/(p+1) .

(2.8)

Needless to say, these γ i 's have positive real part.

A triple jump composition strategy: reaching order 8 Symmetric composition methods Φ (p) h of order p (p even) can be constructed by induction, setting

Φ (2) h = Φ h , Φ (p+2) h = Φ (p) γ 3,p h • Φ (p) γ 2,p h • Φ (p) γ 1,p h for p ≥ 2, (2.9)
where γ 1,p , γ 2,p , γ 3,p are given in (2.8). The method Φ (p) h requires s = 3 p/2-1 compositions of the basic method Φ h . Taking care of the non-commuting products (product signs should be read from the right to the left), we may write

Φ (2) h = Φ h = exp hB 2 exp(hA) exp hB 2 ,
together with

Φ (4) h = 3 ∏ j=1 Φ γ j,2 h , = Φ γ 3,2 h • Φ γ 2,2 h • Φ γ 1,2 h , Φ (6)
h = 3 ∏ k=1 3 ∏ j=1 Φ γ k,4 γ j,2 h , Φ (8)
h = 3 ∏ ℓ=1 3 ∏ k=1 3 ∏ j=1 Φ γ ℓ,6 γ k,4 γ j,2 h
, and so on. In summary, each method Φ (p)

h reads Φ (p) h = 3 (p/2)-1 ∏ j=1 Φ α j,p h ,
for some coefficients α j,p that are obtained as products of the γ k,p-2 's, γ k,p-4 's, ..., γ k,2 's. This defines the coefficients α j,p .

Remark 2.2 Coming back to the value of Φ h in turn allows to write

Φ (p) h = 3 (p/2)-1 ∏ j=1 exp (b j,p hA) exp (a j,p hB) ,
where a 1,p = α 1,p /2, b 1,p = α 1,p , while a j,p = (α j,p + α j-1,p )/2, b j,p = α j,p whenever 2 ≤ j ≤ 3 (p/2)-1 -1, and lastly a j,p = α j,p /2, b j,p = 0 in the particular case when j = 3 (p/2)-1 . Now, the important point is Proposition 2.1 The above defined method Φ (p) h has order p. Besides, for p = 2, 4, 6, 8, the coefficients α j,p ( j = 1, . . . , 3 (p/2)-1 ) satisfy Re (α j,p ) > 0.

This property ceases to hold whenever p ≥ 10.

Remark 2.3

The fact that the method Φ (p) h reaches order p is here presented in the case of a linear equation ẏ = Ay + By, and when the basic method Φ h is the Strang splitting algorithm. Naturally the same result holds when the chosen basic method is any symmetric method of order 2. Mutatis mutandis (see introduction), and provided the appropriate assumptions described in the introduction are met, the same result holds in the nonlinear setting as well, where y → A(y) and y → B(y) become nonlinear operators.

As a consequence, starting from a second order symmetric method φ h (be it the Strang splitting algorithm Φ h as above, or any symmetric second order algorithm), the present composition technique can only improve numerical accuracy up to order 8, while preserving the use of time steps α j,p h that all have non-negative real part. This property may be somewhat precised. We observe in Figure 2.2 that the quantity max j arg(α j,p ) increases with the order p of the composition methods in (2.9). For the method (2.9) of order p = 10 this quantity is greater than π/2: it involves 81 factors Φ α j,10 and the middle coefficient α 41,10 has a negative real part, namely Re (α 41,10 ) ≈ -5 • 10 -5 < 0. Thus, this method cannot be used, in general, when the operator A or B has large negative eigenvalues, nor can it be extended, in the infinite dimensional case, when the operator A coincides with the Laplacian.

Another triple jump composition strategy: reaching order 14

Before concluding this paragraph about triple jump procedures, we mention an improvement of the above method. To reduce the quantity max i=1...s arg(α i,p ) , an idea is to alternate the coefficients (γ 1,p , γ 2,p , γ 3,p ) by (γ 1,p , γ 2,p , γ 3,p ) in (2.9). In other words, we propose here to set

Φ h (2) = Φ h , Φ h (p+2) = Φ γ 3,p h (p) • Φ γ 2,p h (p) • Φ γ 1,p h (p) if p/2 odd, (2.10) Φ h (p+2) = Φ γ 3,p h (p) • Φ γ 2,p h (p) • Φ γ 1,p h (p)
else.

This yields a family of composition methods which again reads

Φ h (p) = 3 (p/2)-1 ∏ j=1 Φ α j,p h ,
for some coefficients α j,p that can be explicitely computed. In that situation we have Proposition 2. [START_REF] Blanes | On the necessity of negative coefficients for operator splitting schemes of order higher than two[END_REF] The above defined method Φ h (p) has order p. Besides, for p = 2, 4, 6, 8, 10, 12, 14, the coefficients α j,p ( j = 1, . . . , 3 (p/2)-1 ) satisfy Re α j,p > 0.

This property ceases to hold whenever p ≥ 16.

Remark 2.4 Surprisingly, the sum of the moduli of coefficients

3 (p/2)-1 ∑ j=1 α j,p and 3 (p/2)-1 ∑ j=1 α j,p in the con- sidered composition methods Φ (p) h or Φ h (p)
is bounded as the order p goes to infinity. It is estimated by

∞ ∏ k=1 γ 1,2k + γ 2,2k + γ 3,2k ≤ ∞ ∏ k=1 2 + 2 1/(2k+1) |2e iπ/(2k+1) + 2 1/(2k+1) | = ∞ ∏ k=1 1 + π 2 36k 2 + O 1 k 3 < +∞
This means that the length of the family of polygons in figure 2.1 is bounded (this limit is ≈ 1.315).

Quadruple Jump composition methods

In the similar way we have derived symmetric, triple jump composition methods, we investigate here the symmetric quadruple jump case. To do so, for any even integer p, we need to find complex quadruples (γ 1,p , γ 2,p , γ 3,p , γ 4,p ), still denoted by the letters γ j,p not to overweight notation, such that γ 1,p = γ 4,p , γ 2,p = γ 3,p (symmetry),

γ 1,p + γ 2,p + γ 3,p + γ 4,p = 1, γ p+1 1,p + γ p+1 2,p + γ p+1 3,p + γ p+1 4,p = 0 (Theorem 2.

1). (2.11)

Starting from any basic symmetric method φ h of order p, the symmetric method φ γ 4,p h • φ γ 3,p h • φ γ 2,p h • φ γ 1,p h automatically has order p + 2. Now, the set of all complex solutions to (2.11) is given by the p values γ (k)

1,p = γ (k) 2,p = γ (k) 3,p = γ (k) 4,p = 1 2 -2e 2kiπ/(p+1) (k = 1, . . . , p).
(2.12)

The two complex conjugate solutions with minimal sum of moduli are obtained with k = ±p/2. These solutions also have minimal value of max i=1,...,4 arg γ

(k) i,p
. Therefore, we set γ 1,p = γ 2,p = γ 3,p = γ 4,p = 1 4 + i sin(π/(p + 1)) 4 + 4 cos(π/(p + 1))

.

With this notation we define the quadruple jump procedure

Ψ (2) h = Φ h , Ψ (p+2) h = Φ γ 4,p h (p) • Φ γ 3,p h (p) • Φ γ 2,p h (p) • Φ γ 1,p h (p) , (2.13)
where the method Φ (p) h has been defined before in (2.10). Naturally, each method

Ψ (p) h reads Ψ (p) h = 4×3 (p/2)-2 ∏ j=1 Φ β j,p h ,
for some coefficients β j,p that can be explicitly computed. We draw the reader's attention to the fact that Ψ (p+2) h is here defined recursively using Φ h

(p)
as a building block (instead of Ψ (p) h

-in other words, we

do not define Ψ (p+2) h as the composition Ψ (p) γ 4,p h • Ψ (p) γ 3,p h • Ψ (p) γ 2,p h • Ψ (p) γ 1,p h
). This choice is made to reduce the total number of compositions of the basic method Φ h needed to build up the method Ψ (p) h : in our case, Ψ (p) h requires 4 × 3 (p/2)-2 compositions of Φ h (instead of 4 (p/2)-1 ).

We have the Proposition 2. [START_REF] Chambers | Symplectic integrators with complex time steps[END_REF] The above defined method Ψ (p) h has order p. Besides, for p = 2, 4, 6, 8, 10, 12, the coefficients β j,p ( j = 1, . . . , 4 × 3 (p/2)-2 ) satisfy Re (β j,p ) > 0.

This property ceases to hold whenever p ≥ 14.

Another advantage of this composition procedure is that we obtain an accurate approximation of the solution at intermediate time steps as well. More precisely, we have the Proposition 2.4 Take an initial condition y 0 and an even integer p. Define for any integer n the sequence of vectors y n (n ≥ 0) by the recursion y n+1 = Ψ (p) h (y n ). Lastly, define y(t) as the solution to the ODE ẏ(t) = A y(t) + B y(t) with initial condition y(0) = y 0 .

Then, y n approximates y(nh) to within O(h p ). Moreover, writing

Ψ (p) h = Φ γ 4,p-2 h (p-2) • Φ γ 3,p-2 h (p-2) • Φ γ 2,p-2 h (p-2) • Φ γ 1,p-2 h (p-2)
, and setting

y n+1/2 = Φ γ 2,p-2 h (p-2) • Φ γ 1,p-2 h (p-2) (y n ),
we also have that y n+1/2 approximates y ((n + 1/2)h) to within O(h p ).

More generally, the same result holds, with the obvious change of notation, taking any symmetric method φ (p-2) h of order p -2 and defining the improved symmetric method of order p as ψ

(p) h = φ (p-2) γ 4,p-2 h • φ (p-2) γ 3,p-2 h • φ (p-2) γ 2,p-2 h • φ (p-2) γ 1,p-2 h .
The proof of this fact simply comes from observing that γ 1,p + γ 2,p = 1/2 for any p, so that the two couples (2γ 1,p , 2γ 2,p ) and (2γ 2,p , 2γ 1,p ) satisfy the order equations (2.2) with s = 2. Hence the recursion

y n+1/2 = Φ γ 2,p-2 h (p-2) • Φ γ 1,p-2 h (p-2)
(y n ) yields an approximation of the true solution at time t = nh + h/2 with local error O(h p ). Since this error is not propagated (it is only an inner stage), we obtain an approximation of order p both for y n+1 at time (n + 1)h and for y n+1/2 at time (n + 1/2)h.

Convergence analysis for unbounded operators

In this section, we extend the analysis of the previous paragraph, valid in the finite dimensional case, to the infinite dimensional situation. We first give a general statement, next specify the assertions in the case of splitting algorithms applied to linear diffusion equations.

A general statement

Hansen and Ostermann in [START_REF] Hansen | Exponential splitting for unbounded operators[END_REF] have provided an elegant and general framework, valid for linear equations, and which allows to assert that if a splitting method is p-th order accurate in the finite dimensional case, then the same method is p-th order accurate in the infinite dimensional case as well, provided some natural functional analytic assumptions are met.

The Hansen and Ostermann result states the following.

Theorem 3.1 (see Hansen and Ostermann [14]

) Let X be an arbitrary complex Banach space with norm • . Denote by the same symbol . the norm on the space of bounded linear operators over X. Consider s linear unbounded operators A j ( j = 1, . . . , s). Lastly, take a time T ≥ 0, an integer p, and an initial datum u 0 ∈ X. Assume that the following assumptions are met:

(i) (semi-group property).
The linear operators A j ( j = 1, . . . , s) and A 1 + • • • + A s generate C 0 semigroups on X. Moreover there exist a real ω and s real numbers ω j ( j = 1, . . . , s) such that ∀t ≥ 0, e t(A 1 +•••+A s ) ≤ e ωt and e tA j ≤ e ω j t .

(3.1)

(ii) (smoothness assumption).

For any operator E p+1 that is obtained as the product of exactly p + 1 factors chosen amongst the A j 's, there is a constant C > 0 such that

∀0 ≤ t ≤ T, E p+1 e t(A 1 +•••+A s ) u 0 ≤ C. ( 3.2)
(iii) (splitting method). Take a splitting method S of the form S = m ∏ j=1 e γ j,1 hA 1 e γ j,2 hA 2 . . . e γ j,s hA s , where γ j , 1 ≤ j ≤ m, are nonnegative reals, while m is an integer. Assume this splitting method is a p-th order approximation of e h(A 1 +•••+A s ) . This means that whenever the A j 's are replaced by finite dimensional matrices M j , say, we have ∏ s j=1 e γ j,1 hM 1 . . . e γ j,s M s = e h(M 1 +•••+M s ) + O(h p+1 ), in the sense of matrix norms. Under all these assumptions, the following holds. There exists a constant C > 0 such that for any integer n ≥ 0 and any time step h > 0 satisfying nh ≤ T , we have,

S n -e nh(A 1 +•••+A s ) u 0 ≤ Ch p .
Using this Theorem, the following is easily deduced

Corollary 3.1 (Banach space formulation)

Let X be an arbitrary complex Banach space with norm • . Denote by the same symbol . the norm on the space of bounded linear operators over X. Take two linear unbounded operators A and B. Lastly, take a time T ≥ 0, an integer p, and an initial datum u 0 ∈ X. Assume that the following assumptions are met:

(i) (semi-group property).

The operator A + B generates a C 0 semigroup on X. Besides, for any z ∈ C such that Re(z) > 0, the linear operators zA and zB, generate C 0 semigroups on X. Lastly, for any given z ∈ C with Re(z) > 0, there exists a real number ω such that ∀t ≥ 0, e t(A+B) ≤ e ωt , and e tzA + e tzB ≤ e ω|z|t .

(3.3)

(ii) (smoothness assumption).

For any operator E p+1 that is obtained as the product of exactly p + 1 factors chosen amongst A and B, there is a constant C > 0 such that Assume this splitting method is a p-th order approximation of e h(A+B) , meaning that whenever A and B are replaced by finite dimensional matrices M and N we have ∏ s j=1 e b j hM e a j hN = e h(M+N) + O(h p+1 ), in the sense of matrix norms. Under all these assumptions, the following holds. There exists a constant C > 0 such that for any integer n ≥ 0 and any time step h > 0 satisfying nh ≤ T , we have, S ne nh(A+B) u 0 ≤ Ch p .

∀0 ≤ t ≤ T, E p+1 e t(A+B) u 0 ≤ C. ( 3
(3.5)

In particular, under the present assumptions on A, B, and u 0 , the methods Φ

(p) h , Φ h (p)
, and Ψ (p) h discussed in the previous paragraphs satisfy estimate (3.5) Before ending this paragraph, we propose a slightly different formulation, adopting a Hilbert space setting: so-called mα-dissipative operators are well-adapted in the present context, where we eventually wish to derive splitting methods that are adapted to parabolic equations.

Let α belong to [0, π/2] and define the sector S α in the complex plane by

S α = {z ∈ C, z = 0 or | arg z| ≤ α}. 0 α S α
Let H be a complex Hilbert space with scalar product denoted by (•, •). Take a linear, unbounded operator A on H, with domain D(A), a dense subspace of H. We recall that A is said mα-dissipative whenever for all u in D(A), the quantity (-Au, u) belongs to S α , and if for all complex z / ∈ S α , the operator zId + A is an isomorphism from D(A) to H. A nice introduction to mα-accretive operators 1 can be found in [START_REF] Crouzeix | Approximation of parabolic equations[END_REF]. It is known that an mα-dissipative operator generates a C 0 semigroup on H, denoted by e tA (t ≥ 0), and e tA is a contraction operator from H to H. Besides, if an operator A is such that there exists a real number c for which A + cId is mα-dissipative for some α ∈ [0, π/2], then A generates a C 0 semi-group as well, and we have the estimate e tA ≤ e +ct whenever t ≥ 0. Corollary 3.2 (Hilbert space formulation) Let H be a complex Hilbert space with scalar product (•, •) and associated norm . . Take an initial datum u 0 ∈ H, a time T ≥ 0, and an integer p. Assume the following:

(i) Let A resp. B

Application to linear diffusion equations

In this parapraph, we apply the above results in the case of the linear heat equation with potential

∂ t u(t, x) = ∆ u(t, x) + V (x)u(t, x), u(0, x) = u 0 (x).
To fix ideas, consider the case when x belongs to the whole space R d . There are many settings adapted to this equation, and one may either seek solutions u having Hölder smoothness, or L p (R d ) smoothness, or Sobolev regularity (in turn based either on L p -p = 2 -or on L 2 ). In order to keep a simple presentation, we choose to work in an L 2 -based Sobolev space setting. Therefore, we introduce the Hilbert space H = L 2 (R d ), and the two operators A : u → ∆ u and B : u → Vu. The operator A with domain D

(A) = H 2 (R d ) is m0-dissipative on H. Whenever V ∈ L ∞ (R d ), the operator B is bounded on H. Lastly, the operator A + B is such that A + B -V L ∞ Id is m0-dissipative.
Hence assumption (i) in corollary 3.2 is met.

To ensure assumption (ii), namely the smoothness assumption, we take an integer p, and assume that u 0 ∈ H 2(p+1) (R d ) and V ∈ W p,∞ (R d ). This ensures that assumption (ii) in corollary 3.2 is met.

We are in position to state the 1 An operator B is said mα-accretive whenever A = -B is mα-dissipative. Needless to say, the similar result holds when the heat equation is considered on the Torus as well, or on a bounded domain with appropriate boundary conditions and smooth enough boundary, etc. , p = 2,4,6,8. Right picture: "quadruple jump" composition methods Ψ h (p) , p = 2,4,6,8. For all these pictures, solid lines: basic method is the Strang splitting with exponential maps (2.3) -dashed lines: basic method is Peaceman-Rachford formula (2.4) .

The nonlinear case

At least formally, all above results immediately extend to the nonlinear situation, provided all exponentials e hA etc. are replaced by the appropriate nonlinear flows (see introduction).

In that perspective, we consider the one-dimensional, non-linear heat equation on the Torus T (identified with

[0, 1]) ∂ t u(t, x) = ∆ u(t, x) + F(u(t, x)), u(0, x) = u 0 (x),
where F(u) is a non-linear reaction term and, for the purpose of testing our methods, we have retained Fisher's potential

F(u) = u(1 -u). The differential equation ∂ u ∂ z = u(1 -u), u(0) = u 0 can be solved analytically as u(z) = u 0 + u 0 (1 -u 0 ) (e z -1) 1 + u 0 (e z -1)
, which is well defined for small complex time z. We discretize the equation in space as in the linear case. The original nonlinear heat equation becomes u = Au + F(u), (4.1)

where the vector u(t) is as in the linear case, the Laplacian ∆ is approximated by the above N × N matrix A, and the vector F(u) stands for F(u) = u 1 (1u 1 ), . . . , u N (1u N ) .

The experimental conditions are the ones we used to illustrate the linear case. In Figure 4.2, we compare the accuracy of the composition methods introduced in this article ("triple" (2.10) and "quadruple" (2.13) jump compositions). We plot for many stepsizes the solution error at time T as a function of the number of evaluations of the basic method. As a basic method, we consider (in solid lines) alternatively the Strang splitting Φ h involving exact flows (where the terms e hB/2 are replaced by half a time step of the exact flow of the nonlinear differential equation ẏ = F(y)), and (in dashed lines) the Peaceman-Rachford formula Φ P h (where the implicit Euler and explicit Euler operators are replaced by the appropriate nonlinear operators obtained by applying the explicit and/or implicit Euler algorithm to the nonlinear equation ẏ = F(y)).

Conclusion

We have constructed new high-order compositions methods and splitting methods using complex coefficients for parabolic linear and non-linear parabolic partial differential equations. Based on the results of Hansen & Osterman [START_REF] Hansen | Exponential splitting for unbounded operators[END_REF], a convergence analysis is provided in the linear case. The methods we have derived are all based on a composition procedure with complex time steps; they actually allow to build up high order methods which are based on the composition of possibly cheap low order algorithms.

Going a bit further, we may stress here that it is also possible to construct high-order splitting methods for which only one operator carries complex time steps. Such methods are however not based on the composition technique we have presently developed. For instance, the following splitting method is symmetric and of order 4, e b 1 hV e a 1 hA e b 2 hV e a 2 hA e b 3 hV e a 2 hA e b 2 hV e a 1 hA e b 1 hV (5.1)

where b 1 = 1/10i/30, b 2 = 4/15 + 2i/15, b 3 = 4/15i/5 are complex, and a 1 = a 2 = a 3 = a 4 = 1/4 are all reals. Such a decomposition is interesting when the evolution along A carries the most computational cost and the evolution with V is cheap to compute (e.g. when V is a diagonal matrix): in that case the extra cost due to the complex numbers is marginal. This type of splitting method is also of great interest in the case where one operator has its eigenvalues close to the imaginary axis, like e.g. the Ginzburg-Landau equation. Note however that in this very case, and because this algorithm is not based on a composition procedure, the true exponentials cannot by replaced, in general, by low order approximations: the exponential needs to be approximated to within the appropriate order in any circumstance here. Independently, note that a systematic study of optimal composition methods (i.e. methods with optimal error constants) is out of the scope of this paper and will be the subject of a future article by the same authors. It requires the resolution in C of the polynomial systems of order conditions for composition methods and splitting methods. Also, a theoretical analysis in the case of a non-linear source is in preparation.

. 4 )

 4 (iii) (splitting method). Next, consider s complex numbers a 1 , . . . a s , b 1 , . . . , b s , and take a splitting method S of the form S = s ∏ j=1 e b j hA e a j hB .

Fig. 4 . 1

 41 Fig. 4.1 Plot: Linear potential V (x) = 2 + sin(2πx). Error (L 2 norm at time T ) of composition methods versus number of evaluations of the basic method Φ h . Left picture: "triple jump" composition methods Φ h (p)

  Values of max i=1...s |argγ i | for various composition methods.

	π/2	max i=1...s	|arg(γ i )|
	π/4					triple jump Φ triple jump Φ h (p) h (2.9) (p) (2.10)
						quadruple jump Ψ h (p)	(2.13)
							order p
	2 0	4	6	8	10	12	14	16
	Fig. 2.2

  for 0 ≤ nh ≤ T , with S replaced by Φ , whenever p = 2, 4, 6, 8, resp. p = 2, 4, 6, 8, 10, 12, 14, resp. p = 2, 4, 6, 8, 10, 12.

	Φ h	(p)	, resp. Ψ	(p)	(p) h , resp.

h

  be such that there exist two real numbers a and b for which A + aId resp. B + bId are mα-dissipative resp. mβ -dissipative operator for some α ∈ [0, π/2] resp. β ∈ [0, π/2]. Lastly, assume that there is a real number c for which A + B + cId is mγ-dissipative for some γ ∈ [0, π/2].

	(ii) Assume that for any operator E p+1 that is obtained as the product of exactly p + 1 factors chosen
	amongst A and B, that there is a constant C > 0 such that
	∀0 ≤ t ≤ T,	E p+1 e t(A+B) u 0 ≤ C.	(3.6)
	Then, the methods Φ h , Φ h (p) 0 ≤ nh ≤ T , with S replaced by Φ (p) , and Ψ h (p) (p) h , resp. Φ h discussed in the previous paragraphs satisfy estimate (3.5) for (p) , resp. Ψ (p) h , whenever p = 2, 4, 6, 8, resp. p = 2, 4, 6, 8, 10, 12, 14, resp. p = 2, 4, 6, 8, 10, 12.

  Theorem 3.2 Under all these assumptions, considering either the method Φ ≤ p ≤ 12), estimate(3.5) with S replaced by one of the above methods holds true.

	p ≤ 14) or Ψ h (p)	(2	(p) h (2 ≤ p ≤ 8) or Φ h	(p)	(2 ≤

This number of discretization points is naturally arbitrary. The fact that, with this value of N, the extrapolation method (4.2) diverges, as expected -see below -, while it converges for smaller values of N, indicates that the value N = 100 is enough for our purposes.

Acknowledgement. The authors are grateful to Ernst Hairer for useful discussions and comments. The fourth author was partially supported by the Fonds National Suisse, project No. 200020-109158.

Numerical results

In this section, we numerically illustrate the above convergence results.

The linear case

We consider the one-dimensional linear heat equation with potential on the Torus T (identified with [0, 1])

where the potential V is taken as V (x) = 2 + sin(2πx).

In order to discretize the equation in space, we take a (large) integer N and choose a finite differences procedure on the regular grid 0, 1/N, 2/N, . . ., N/N. The original heat equation then becomes

where the vector u(t) belongs to R N , and has the form

and u j (t) is an approximation of u t, j N , while the Laplacian ∆ is approximated by the N × N matrix A given by

, and the vector Bu stands for

We take the C ∞ initial condition u 0 (x) = sin(2πx), and consider a spatial discretization with N = 100 points 2 . In Figure 4.1, we compare the accuracy of the composition methods introduced in this article ("triple" (2.10) and "quadruple" (2.13) jump compositions) on the time interval [0, T ], where T = 0.2. We plot for many stepsizes the solution error at time T as a function of the number of evaluations of the basic method. As a basic method, we consider (in solide lines) alternatively the Strang splitting Φ h involving exact flows (where the terms e hB/2 are replaced by half a time step of the exact flow of the nonlinear differential equation ẏ = F(y)), and (in dashed lines) the Peaceman-Rachford method Φ P h . The 'exact' solution is computed with a very small time step. We observe the expected orders (lines of slopes 2, 4, 6, 8). Surprisingly, composition methods using the Peaceman-Rachford formula are slightly more accurate than the one using exponentials. h , p = 2,4,6,8. For all these pictures, solid lines: basic method is the Strang splitting with exponential maps (2.3) -dashed lines: basic method is the Peaceman-Rachford formula (2.4).

In Figure 4.3, we compare the "quadruple jump" composition method of order 4 with two extrapolation methods. We also give the results for the Strang splitting of order 2. We use the same initial data and parameters as before. The first extrapolation formula we consider is

where for the basic method Φ h , we take alternatively the Strang splitting Φ h with the exact flows, see left picture in Figure 4.3, or the -conveniently adapted, see above-Peaceman-Rachford formula, see right picture. However, as pointed out in [START_REF] Schatzman | Numerical integration of reaction-diffusion systems, Numerical methods for ordinary differential equations (Auckland[END_REF]Sect. 6] this scheme is not stable and does not converge in the second case (see dashed-dotted line in right picture). Another extrapolation method (dashed lines) is considered in [START_REF] Schatzman | Numerical integration of reaction-diffusion systems, Numerical methods for ordinary differential equations (Auckland[END_REF] and taken from [START_REF] Dia | Méthodes de directions alternées d'ordre élevé en temps[END_REF], 45 64

Although the formal order of this method is 4, it is said in [START_REF] Schatzman | Numerical integration of reaction-diffusion systems, Numerical methods for ordinary differential equations (Auckland[END_REF] that the true order of convergence of this method is not clearly understood, and in the numerical experiments for linear problems in [START_REF] Dia | Méthodes de directions alternées d'ordre élevé en temps[END_REF], "the formal order in not reached ; the experimental precision is smaller than the theoretical precision, and the difference is smaller that 1". Finally, and for a fair comparison in Figure 4.3, it should be mentioned that computations using complex numbers are actually about four times more expensive than computations with reals numbers (because of the cost of a multiplication). h (solid line), extrapolation method (4.2) (dashed-dotted line), extrapolation method (4.3) (dashed lines). Left picture: Basic method is the Strang splitting with exponential maps (2.3). Right picture: Peaceman-Rachford formula (2.4).


Denis Mounier
email: denis.mounier@univ-lemans.fr

Christophe Poilane

Haithem Khelfa

Pascal Picart

Christophe Poil Âne

Gigahertz laser resonant ultrasound spectroscopy for the evaluation of transverse elastic properties of micrometric bers

Keywords: fiber, laser ultrasonics, resonant ultrasound spectroscopy (RUS), finite element modeling, laser interferometry, cylinder, eigenmodes, Rayleigh modes, wave gallery modes, design of experiment (DOE)

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Resonant ultrasound spectroscopy (RUS) has been extensively used during the past two decades to evaluate the elastic constants of various anisotropic solid materials with an ulti-mate precision [1][2][3][4]. From the measurement of the mechanical resonance frequencies of a sample, the elastic constants of the material can be determined.

The PZT-RUS method, i.e. the RUS technique using piezoelectric transducers for both the excitation and the measurement of vibrations, is the most widely used technique. The optimum size of the sample for PZT-RUS is about 1 cm, but the technique can be used for small samples whose dimensions are of the order of 0.5 mm with specially designed transducers [3]. For such small samples, the RUS technique may require the measurement of resonance frequencies up to 50 MHz [5].

However, even in the case of samples of several millimeters in size, contact forces between the ultrasonic transducers and the sample may lead to spurious frequency shifts for the resonance frequencies in addition to the widening of resonances [6][START_REF] Sedlák | 1st International Symposium on Laser Ultrasonics: Science, Technology and Applications[END_REF][START_REF] Park | [END_REF]. The problem of contacts comes up more and more crucial as the size of samples decreases.

The suppression of contact forces by using lasers for both the excitation and the detection of vibrations results in getting sharp resonances, which is beneficial to the accuracy of measurement [START_REF] Sedlák | 1st International Symposium on Laser Ultrasonics: Science, Technology and Applications[END_REF][START_REF] Park | [END_REF]. The excitation of vibrations can be carried out by using a periodic laser modulation, whose frequency is swept over the required range [9]. Alternatively, short laser pulses can be used to excite the vibrations of the sample. The advantage of the pulsed method lies in the possibility of exciting simultaneously many eigenmodes over a very large frequency range. The broad band excitation provided by short optical pulses is particularly appropriate to the study of very small samples. The vibrations are measured in the time domain and the eigenfrequencies are obtained from the Fourier transform of the measured signal. The pulsed Laser-RUS method is preferable to determine the elastic constants of small samples [START_REF] Sedlák | 1st International Symposium on Laser Ultrasonics: Science, Technology and Applications[END_REF][START_REF] Park | [END_REF]10,11] or to detect flaws in a small composite structure [12][START_REF] Amziane | Proc. of SPIE[END_REF][START_REF] Petit | [END_REF][15]. In this paper, we present a pulsed LRUS technique applied to the evaluation of elastic constants of micrometric fibers, whose diameters are in the range 5-50 µm. For such small dimensions, the eigenfrequencies must be measured in the 0.05-1 GHz range.

The LRUS allows one to evaluate the elastic constants on a single fiber, in contrast with indirect methods of evaluation [16]. From the evaluation of the elastic constants of a unidirectional fiber composite, the elastic constants of the fibers are calculated. This indirect method is based on specific assumptions about the nature of the reinforcing effect of fibers in the polymer matrix. These assumptions may not be valid for all kinds of reinforced polymers, so the application of this method may lead to erroneous results. Therefore, it is is much more reliable to evaluate the elastic constants directly on a single fiber.

A direct evaluation of elastic constants of carbon fibers by a non-contact technique was first achieved using the technique of laser picosecond ultrasonics [17]. The transverse elastic constant C 11 of carbon fibers was evaluated through the measurement of the time of flight of acoustic pulses that propagate back and forth across a fiber diameter, which determines the longitudinal acoustic velocity and then the C 11 elastic constant. However, the evaluation of the C 44 elastic constant, which requires the measurement of the transverse acoustic velocity, could not be achieved because of the difficulty to generate and to detect picosecond transverse acoustic waves.

Fibers are commonly used in reinforced polymers, such as glass or carbon fibers. The interest of the method is great for anisotropic fibers (carbon and vegetal fibers) because of the lack of reliable and accurate methods for the evaluation of the elastic constants of fibers, in particular the elastic constants characterizing the transverse elastic properties.

The evaluation of the longitudinal Young modulus of fibers is generally carried out using the tensile test method [18,19]. But this method is impracticable to measure the transverse elastic properties.

We illustrate the application of the LRUS method to the evaluation of the transverse Young modulus and Poisson ratio of a micrometric aluminum wire. The evaluation of the transverse elastic properties of fibers requires the excitation of the cross section eigenmodes.

In our approach, the propagation of guided acoustic waves along the fiber axis is neglected, so the mechanical displacements are restricted to a plane perpendicular to the fiber axis.

The restriction in 2D of the cylinder vibrations results in a discrete vibration spectrum, which is required for the application of the RUS method. In section 2, we will describe the experimental method used to excite and detect cross-section modes of a fiber. In section 3, we will present the method used to solve the inverse problem, which leads to the evaluation of the transverse elastic properties of an aluminum fiber. We will also discuss the conditions of validity RUS applied to a cylinder.

Measurement of ultrasonic vibrations

Laser ultrasonic setup

The irradiation of a cylinder with a laser pulse can excite many acoustic guided modes.

The propagation of guided modes along the cylinder axis is determined by the dispersion relation ω(k), where k is the propagation wave vector along the cylinder axis and ω/(2π) is the frequency of the guided mode. The dispersion relation determines both the phase velocity v φ = ω/k and the group velocity v g = dω/dk of the guided mode.

Both the geometry of the laser spot focused on the cylinder surface and the duration of the laser pulse determine the spectral range of guided modes that can be excited. If the pump beam is focused to form a linear spot parallel to one of the generating line of the cylinder, then a linear acoustic source will be formed and a cylindrical acoustic wave propagating perpendicularly to the cylinder axis will be generated. Consequently, a linear acoustic source parallel to the cylinder axis, which is very long compared to the fiber diameter will favor the excitation of guided modes which very small wave vectors k ≈ 0, i.e. with very large wavelengths λ = 2π/k compared to the fiber diameter. Moreover, such a linear acoustic source will force the displacements in the directions perpendicular to the cylinder axis and the displacements will be almost independent of the z-coordinate around the center of the acoustic source.

Figure 1 shows the LRUS setup, which includes two parts: the first part is devoted to the excitation of vibrations and the second part, the interferometer, serves to measure ultrasonic vibrations. The vibrations of fibers are excited using a Q-switched Nd:YAG microchip laser at 1064 nm from TEEM photonics whose pulses have a maximum energy of 10 µJ and a duration of 0.6 ns at a maximum repetition rate of 5 kHz. With a 0.6 ns pump pulse, it is possible to generate vibrations up to the 1-2 GHz frequency range, which is sufficient to excite many eigenmodes of a fiber whose diameter is in the 5-50 µm range.

The fibers are mounted on a 1 cm diameter plastic ring (Inset of Fig. 1, at the bottom right). In order to excite acoustic waves which propagate mainly perpendicularly to the fiber axis, the length of the pump spot must be much larger than the fiber diameter. The pump laser is therefore shaped to form a narrow line, parallel to the fiber axis. By using a cylindrical lens placed in the optical path of the pump beam, we obtain an elliptical spot with dimensions 6 µm × 100 µm at the fiber surface, with the main axis of the ellipse aligned with the fiber axis (see left inset of Fig. 1). By a proper adjustment of the pulse energy with an variable attenuator (not represented on Fig. 1), the vibrations of the fiber can be excited in the non destructive thermoelastic regime [START_REF] Scruby | Laser ultrasonics: techniques and applications[END_REF].

For the measurement of radial ultrasonic vibrations, a stabilized homodyne Michelson interferometer is used [START_REF] Scruby | Laser ultrasonics: techniques and applications[END_REF]. The probe beam at 532 nm is focused onto the fiber surface using an aspherical lens (x20) having 8 mm focal length and 0.5 numerical aperture. The resulting focused spot centered on the pump line has a diameter of about 2 µm (inset of Fig. 1).

An incident probe power of 1-2 mW is sufficient to achieve a measurement sensitivity of 1 mV nm -1 on a metallic surface.

The optical interferometric signal is received by a high-speed photodetector with a bandpass of 1 GHz. Using a digital oscilloscope with a band-pass of 3 GHz, the signals are recorded in synchronization with the pump pulses during about 2 µs after the pump laser pulse. Thanks to the averaging of 10 000 acquisitions, which necessitates a only two seconds at the pulse repetition rate of 5 kHz, we obtain a detection threshold of about 5 pm for the measurement of out-of-plane displacements, which is sufficient to detect the vibrations with a satisfactorily signal to noise ratio.

The high frequency ultrasonic vibrations are superimposed to a slow and huge displacement of the fiber; a displacements of about 5-10 nm which is probably induced by the sudden expansion of the air in contact with the heated fiber surface, thus inducing the recoil of the fiber. This slow drift of the signal is filtered out in order to magnify the high frequency vibrations corresponding to the fiber cross-section eigenmodes. The vibration spectrum is obtained with 0.5 MHz resolution by calculating the Fast Fourier Transform (FFT) of the digitized signal (Fig. 2).

Vibration spectrum of an aluminum wire

The test sample is an aluminum bonding wire with 1% silicon, manufactured by HER-AEUS [START_REF] Heraeus | AlSi 1%, The Aluminum Fine Wire Solution[END_REF]. With a scanning electron microscope (SEM), we measured the wire diameter at several points of the wire. A mean diameter of (32.7 ± 0.1) µm was determined over a length of 1 cm. The ultrasonic measurements were made at several points each millimeter over a wire length of 5 mm.

Evaluation of the transverse elastic properties of an aluminum fiber

Two-dimensional finite element modeling

In the following, we will focus our attention on k = 0 guided modes, whose displacements are restricted within the xy-plane. From the experimental conditions depicted before, we expect that the excited eigenmodes are mainly the k = 0 guided modes. We will refer to such guided modes as cross-section modes. The eigenfrequencies of cross-section modes depend only on the transverse Young modulus and Poisson ratio of the cylinder. Consequently, from the measured eigenfrequencies, the transverse Young modulus and Poisson ratio can be determined.

A 2D mechanical model of the aluminum wire is carried out. It is assumed that the elastic properties are homogeneous and invariant by rotation around the fiber axis, i.e. the mechanical properties of the material are transverse isotropic. Thus, the eigenmodes are determined by: the cylinder diameter d, the Young modulus E, the Poisson ratio ν and density ρ. The vibration eigenmodes of an aluminum wire are determined with the software COMSOL Multiphysics. The mode shapes of the first 20 eigenmodes are illustrated in Fig.

3.

Nomenclature of modes

The guided modes of a cylinder are named according to the nomenclature of Silk and Bainton [START_REF] Silk | [END_REF]. The modes are labeled X(m, n) (or simply Xm, n as in Fig. 4 or Table 1) where X = L, T, or F and m, n are positive integers. The letters L, T or F are the initials of longitudinal, torsional and flexural modes, respectively. The number m of a X(m, n) mode characterizes the rotational symmetry of the mode shape and n the ordering of the Xm mode in the frequency scale. In cylindrical coordinates (r, θ, z), the u r and u θ components of the displacement field vary according to the functions cos mθ or sin mθ, which means that the mode shape is invariant after a rotation of 2π/m rad around the cylinder z-axis. The longitudinal and torsional modes, which have both a rotational symmetry for an arbitrary rotation angle around the cylinder axis have m = 0. The longitudinal and the torsional modes are characterized by the displacements u θ = 0 and u r = 0, respectively. The modes with m > 0 are flexural modes.

According to Victorov [START_REF] Viktorov | Rayleigh and Lamb waves[END_REF][START_REF] Clorennec | [END_REF], the cross-section modes belong to one of the two categories: the Rayleigh modes (R) or the Wave Gallery modes (WG). The Rayleigh modes are flexural modes of the n = 0 series, which we denote more explicitly R(m, 0) instead of F (m, 0). The Rayleigh modes series starts with the R(2, 0) mode, as the R(1, 0) mode is the trivial Rayleigh mode corresponding to the translation of the cylinder along a radial direction, having thus the frequency zero. The WG modes are denoted W G(m, n) instead of F (m, n), with n > 0 and m > 0. The first longitudinal mode, denoted L(0, 1), is often called the breathing mode. The first torsional mode is denoted T (0, 1). The L(0, 0) and T (0, 0) modes are trivial modes with their eigenfrequency equal to zero.

The distinction between Rayleigh and some series of Wave Gallery modes will be useful to consider in the perspective of mode identification, as all Rayleigh modes have significant radial components u r and weak orthoradial components u θ at the surface of the cylinder.

On the contrary, the n = 1 series of Wave Gallery modes is characterized by weak radial components compared to the orthoradial ones and the displacements are mainly concentrated in the vicinity of the surface. The displacements for WG modes of the series n > 1 may have very weak components at the surface. In contrast, for some particular WG modes (for example see mode 12 in Fig. 3 and Table 1), the magnitude of the radial displacement may be significant. Considering the correlation between the relative magnitude of radial displacements predicted by the modeling and the amplitudes of the spectral lines may help to identify the modes (cf. section 3.2).

The longitudinal modes L(0, n) and the torsional modes T (0, n) are not degenerate.

On the contrary, the m > 0 modes are two-fold degenerate. This property arises from the existence of the two angular functions cos mθ and sin mθ, expressing the angular dependence of a mode shape. The two degenerate mode shapes can be deduced from each other by a rotation of π/(2 m) rad. If a fiber deviates a little from the circular symmetry, a degeneracy lifting may occur. If the spectral resolution is sufficient, it could be possible to appreciate an elliptical shape of the fiber cross-section. In the spectrum of Fig. 2, the absence of doublets in the spectrum demonstrates the almost perfect circular cross-section of the aluminum wire.

Modes identification

The RUS technique requires the correct mode identification of each measured eigenfrequency. This can be done by determining experimentally the mode shapes, which can be then compared to the calculated ones [10,25]. A raster-scanning of the probe laser on the surface can image the mode shape.

However, in the case of a non-planar surface, like a fiber, it is not possible to scan the probe laser and thus to determine experimentally the mode shapes. Nevertheless, the digital holography technique could be an alternative to determine experimentally the mode shapes of fibers [26,27], but this technique is not yet implemented for small object as fibers. therefore, we will try to guess which eigenmode corresponds to each measured eigenfrequency.

The Rayleigh eigenmodes are characterized by strong radial components at the surface and would thus be more likely attributed to the strongest lines of the spectrum. Indeed, Raleigh modes can be observed in the spectrum up to the order m = 15 (eigenfrequency at 464 MHz). Conversely, WG modes which have predominant orthoradial components at the surface would correspond to spectral lines of smaller amplitude. The torsional modes, with zero radial components, are missing modes.

Eigenfrequencies up to about 200 MHz are well separated in the spectrum and can be identified without ambiguity, except for the two close lines at 135 MHz and 137 MHz, which may be attributed either to the first longitudinal mode L(0, 1) (the breathing mode) or to the W G(2, 1) mode. The strongest line at 137 MHz would correspond more likely to the L(0, 1) mode (mode 5), as this mode has exclusively radial components. The W G(2, 1) mode (mode 4 in Fig. 3) displays very weak radial components, so it would thus more likely corresponds to the weakest spectral line.

The mode identification can be carried out without ambiguity for 9 eigenfrequencies among the first 11 eigenfrequencies shown in Table 1. The values in Table 1 corresponds to the mean experimental frequencies measured on several spectra. Most of the frequencies displayed in Fig. 2 are close to the frequencies of Table 1 within 0.5 MHz.

Solving the inverse problem

In order to determine the mechanical parameters, we have to search the set of mechanical parameters which minimizes the root mean square of the distance between the experimental frequencies f exp i and the calculated ones f calc i , which is:

σ res = (1/N ) N i=1 (f exp i -f calc i ) 2 (1)
A non linear regression algorithm (Gauss-Newton) where N is the number of the experimental frequencies that are considered for the fit. Only frequencies below 210 MHz are considered, as these lines can be identified without ambiguity. Eigenfrequencies beyond 210 MHz have been omitted, as the overlap of some unresolved spectral lines may lead to an ambiguous identification. Only 9 eigenfrequencies among the 11 below 210 MHz can be identified unambiguously. We do not consider the weak line around 135 MHz, probably the W G(2, 1) mode, which is visible only a few times among all the measured spectra, and thus this line cannot be measured accurately. The torsional T (0, 1) mode, which is missing in the spectra, is not considered.

Solving an inverse problem requires in general many iterations. For each iteration, the eigenfrequencies are calculated for a set of input parameters. The use of finite element calculation is not very practicable for iterative processes. For a cylindrical fiber, an analytic calculation of eigenfrequencies would be probably the best method to calculate the eigenfrequencies. However, we preferred to use a more general method which can be applied to a fiber of arbitrary cross-section shape. In order to calculate quickly the eigenfrequencies for a given set of the parameters E, ν and ρ, we determined first the coefficients of polynomial interpolation functions, which can be used to calculate the eigenfrequencies. For a predetermined diameter d 0 , each interpolation function can be expressed in the form of a Taylor series of the three variables E, ν and ρ, which can be express in the form:

f i (E, ν, ρ) = α,β,γ C α β γ i (E -E 0 ) α (ν -ν 0 ) β (ρ -ρ 0 ) γ (2)
where E 0 = 70 GPa and ν 0 = 0.33 and ρ 0 = 2700 kg/m 3 .

A second-order Taylor series is sufficient to get an accurate prediction of the frequencies, so that only 10 coefficients are determined, those verifying the inequality α + β + γ ≤ 2. For each frequency i, the 10 coefficients C α β γ i were determined from a set of 15 finite element models with the parameters (E, ν, ρ) suitably chosen inside the 3-dimensional domain D defined by the intervals: E = (70 ± 3) GPa, ν = 0.33 ± 0.03 and ρ = (2700 ± 75) kg/m 3 .

The method of design of experiment (DOE) for quadratic models [START_REF]Nist/sematech e-handbook of statistical methods[END_REF] was used for this purpose.

The interpolation functions of Eq. 2 are used to predict the eigenfrequencies with an error inferior to 0.05 MHz in the D-domain, which is sufficient compared to the experimental uncertainties.

For a different diameter d of the cylinder, the eigenfrequencies were calculated by first using Eq. 2 and by multiplying the eigenfrequencies vector by the scale factor d 0 /d.

Results

With the 9 interpolation functions, we searched the set of parameters (E, ν, ρ) which minimizes the standard residual σ res . For the first run, we fixed the diameter, the Young modulus and the density: 32.7 µm, 70.0 GPa and 2700 kg/m 3 . We found that the best fit is obtained with σ res = 0.47 MHz when the Poisson ratio ν = 0.3496. If the diameter and the density are fixed to the preceding values, the best fit is E = 69.6 GPa and ν = 0.350, which is obtained with σ res = 0.17 MHz. No significant improvement of the fit can be obtained by varying the density. In order to evaluate the measurement uncertainty for the Young modulus and the Poisson ratio, due to the measurement uncertainty of the fiber diameter, we supposed that the fiber diameter is now 32.8 µm. The best evaluation of the parameters is E = 70.0 GPa and ν = 0.350 with σ res = 0.17 MHz. Therefore, an error in the evaluation of the fiber diameter of 0.1 µm propagates an error of 0.4 GPa in the Young modulus. But it is remarkable that the Poisson ratio is not sensitive to the diameter uncertainty. The evaluation of the Young modulus is also sensitive to the density, with a sensitivity ∂E/∂ρ = 0.026 GPa/(kg/m 3 ).

If the uncertainty of the fiber density is neglected, the result is E = (69.6 ± 0.4) GPa and ν = 0.350 ± 0.001.

The Young modulus and the Poisson ratio were determined by considering only 9 experimental frequencies below 210 MHz. We must verify now if determined parameters give the correct prediction of the other eigenfrequencies. Table 1 shows the comparison between the experimental and calculated frequencies up to 300 MHz with the parameters d = 32.7 µm, E = 69.6 GPa, ν = 0.350 and ρ = 2700 kg/m 3 . Each experimental frequency f exp is an average of several measured eigenfrequencies. The relative errors ∆ = (f exp -f calc )/f calc remain under 0.5% for all the calculated frequencies below 300 MHz. The W G(2, 1) mode is located before the "breathing" L(0, 1) mode, in agreement with our initial guess.

Moreover, the correlation between the calculated and the experimental frequencies is quite good for 31 modes up to the frequency of the Rayleigh mode R(15, 0), which is calculated at 466 MHz and measured at 464 MHz. The relative residuals ∆ = (f exp -f calc )/f calc are plotted in Fig. 4. For all the measured eigenfrequencies the relative error is below 0.5%.

The limits of the 2-dimensional approach

Though the 2D model is capable to fit rather well most of the measured eigenfrequencies, it cannot explain some particular features observed in the experimental spectrum, in particular the unexpected weakness of the "breathing mode" line at 137 MHz.

In the 2D approach, the propagation of acoustic waves along the fiber axis is neglected.

We implicitly assume that the group velocities of the excited guided modes are zero. In order to evaluate an upper limit for the group velocity, we must take into account of the finite length of the pump line, which determines a cutoff in the wave vector spectrum of the excited guided modes. The cutoff wave number is roughly 1/L p = 1/0.1 = 10 mm -1 , where L p is the length of the pump spot. The propagation of guided modes can be ignored during the acquisition time T a = 2 µs only if their group velocity is much lower than v max g = L p /T a = 50 m s -1 . In consequence, the 2D approach is valid only if the group velocity is inferior to 50 m s -1 for all wave vectors k 2π • 10 mm -1 .

The group velocity is zero for pure cross-section modes, i.e. with the k = 0 guided modes. But for some modes, the variation of the group velocity with k may be too rapid between the zero and the cutoff wave vector k max = 2π • 10 mm -1 and this would result in a quick vanishing signal and thus a strong apparent damping of the measured vibrations. In the spectral domain, this will result in weak and broad spectral line.

The fact that most of the spectral lines displayed in Fig. 2 have relatively high amplitudes proves a posteriori that the group velocities of most guided modes are not much larger than the above estimated value 50 m s -1 . However, the condition v g < 50 m s -1 may not be fulfilled for some modes, in particular for the breathing mode.

The 2D model does not take into account the possible Zero Group Velocity modes (ZGV) that may be excited in a cylinder. Such modes are known in structures capable of guiding acoustic waves, such as plates [START_REF] Prada | [END_REF]. The existence of ZGV modes in a cylinder is predicted theoretically [START_REF] Prada-Julia | Institut Langevin[END_REF]. The ZGV modes can exhibit very high quality factors in the range of Q = 1000 -10000. Moreover, it is known that ZGV eigenfrequencies are different from the cross-section eigenfrequencies. The differences between cross-section and ZGV eigenfrequencies could explain some discrepancies between the experimental and the calculated eigenfrequencies in the 2D modeling.

Conclusion

This paper presents a method to evaluate the transverse Young modulus and Poisson ratio of a micrometric metallic wire by using laser resonant ultrasound spectroscopy. A precision of 0.7% for the Young modulus and 0.3% for the Poisson ratio of an aluminum fiber was achieved.

The pump laser was shaped as a line spot parallel to the fiber axis in order to excite cross section eigenmodes. Using a stabilized homodyne Michelson interferometer, the cross section eigenfrequencies were measured up to 500 MHz. The 2D vibration model of the wire proved to be in good agreement for about 30 measured eigenfrequencies within an accuracy of 0.5% in the 50-500 MHz range.

Nevertheless, the spectral resolution can be improved by increasing the acquisition time windows, as some modes are not completely damped 2 µs after the excitation. The improved resolution could be useful to detect minute geometrical defects of the cylinder through the observation of the degeneracy lifting of some modes. Furthermore, the measurement of high frequency Rayleigh modes could be useful to probe the mechanical properties of the material in the vicinity of the surface, as the displacements of high frequency Rayleigh modes are localized mainly near the surface.

In order to achieve ultimate precision for the evaluation of the elastic parameters of fibers, it is essential to take into account the propagation of guided modes along the fiber and therefore to carry out the 3-dimensional modeling of the mechanical system and then to determine the eigenfrequencies of the Zero Group Velocity modes (ZGV). Such modes could be useful to improve the precision in the determination of the elastic parameters and could give access to the determination of the intrinsic damping of the material, i.e. the imaginary part of the Young modulus.

The application of the LRUS method is not restricted to metallic fibers. We intend to apply more extensively our LRUS method to characterize fibers used in reinforced composite materials, such as carbon, glass, kevlar and vegetal fibers. The study of vegetal fibers, in particular flax fiber, which is a hollow structure approximately cylindrical, requires the association of LRUS with a method to determine the actual geometry of the fiber, such as digital holographic tomography.

modes in cylinders (personal communication).

Figure 1 :

 1 Figure 1: Experimental setup of laser ultrasonic resonance spectroscopy (L-RUS) and support of the fiber. BS: beam splitter, PBS : polarizing beam splitter, VA: variable attenuator, λ/4: quarter-wave plate, λ/2: half-wave plate, M1, M2, M3, M4: mirrors, PZT: piezoelectric actuator, GD: ground potential.

Figure 2

 2 Figure2shows one of the measured vibration spectra of the aluminum wire. The inset of Fig.2represents the high frequency ultrasonic component used to calculate the vibration spectrum. Each eigenfrequency can be determined with an accuracy of about 0.5 MHz. For each spectral line, the mean frequency and standard deviation were calculated from the vibration spectra measured at different points of the wire. The standard deviations was found of the order of the spectral resolution 0.5 MHz. The repeatability of the vibration

Figure 2 :

 2 Figure 2: Eigenfrequency spectrum of an aluminum wire of diameter 32.7 µm calculated from the signal shown in the inset, which shows the high frequency component of the recorded signal, corresponding to the vibrations of the wire cross-section.

Figure 3 :

 3 Figure 3: Mode shapes of the first 20 calculated eigenmodes of an aluminum cylinder. Numbers correspond to the ordering numbers of Table 1. Category of each eigenmode: Rayleigh (R), Wave Gallery (WG), longitudinal (L) or torsional (T) is indicated. Rayleigh modes are characterized by predominant radial components. Colors represent the magnitude of the displacements: the deep blue color means a zero displacement, i.e. the nodal areas of the mode shape and the red color is for the maximum magnitude. The arrows indicate both the direction and the magnitude of displacements.

Figure 4 :

 4 Figure 4: Percentage of relative error between experimental and calculated eigenfrequencies with the parameters: d = 32.7 µm, E = 69.6 GPa, ν = 0.350 and ρ = 2700 kg/m 3 .

Table 1 :

 1 Measured and calculated eigenfrequencies in MHz, corresponding to the mode shapes of Figure3.The right column shows the differences: ∆f = f exp -f calc .

		mode	f exp	f calc	∆f		mode	f exp	f calc	∆f
	1	R2,0 70.54	70.67 -0.13 11	R6,0 204.71 204.44 0.27
	2 WG1,1 86.15	86.54 -0.39 12 WG1,3 233.02 233.37 -0.35
	3	R3,0 109.07 108.99 0.08 13 WG4,1 234.28 233.55 0.73
	4 WG2,1 134.72 135.10 -0.38 14	R7,0 235.11 234.33 0.79
	5	L0,1 136.97 136.92 0.06 15 WG2,2 241.98 240.70 1.28
	6	R4,0 142.47 142.49 -0.01 16	T0,2	NA	253.17 NA
	7	T0,1	NA	154.47 NA	17	R8,0 264.02 263.83 0.19
	8	R5,0 173.87 173.97 -0.10 18 WG5,1 280.83 280.62 0.22
	9 WG3,1 184.64 184.69 -0.05 19 WG3,2 281.62 282.05 -0.43
	10 WG1,2 198.33 198.27 0.05 20	R9,0 292.70 293.09 -0.39