Thibault Bourgeron

Vincent Calvez

Jimmy Garnier

Thomas Lepoutre

Existence of recombination-selection equilibria for sexual populations

Keywords: recombination, selection, infinitesimal model, non-homogeneous environment, adaptation, structured populations, eigenproblem, quantitative genetics 2

We study a birth and death model for the adapatation of a sexual population to an environment. The population is structured by a phenotypical trait, and, possibly, an age variable. Recombination is modeled by Fisher's infinitesimal operator. We prove the existence of principal eigenelements for the corresponding eigenproblem. As the infinitesimal operator is 1-homogeneous but nor linear nor monotone, the general Kreȋn-Rutman theory cannot be applied to this problem.

1 Introduction

The infinitesimal model

We first recall the definition and some properties of the infinitesimal model, see [START_REF] Bürger | The Mathematical Theory of Selection, Recombination, and Mutation[END_REF][START_REF] Doebeli | Multimodal pattern formation in phenotype distributions of sexual populations[END_REF][START_REF] Metz | How should we define fitness for general ecological scenarios?[END_REF] for biological motivation, and [START_REF] Mirrahimi | Dynamics of sexual populations structured by a space variable and a phenotypical trait[END_REF][START_REF] Diekmann | The dynamics of adaptation: an illuminating example and a Hamilton-Jacobi approach[END_REF][START_REF] Desvillettes | On selection dynamics for continuous structured populations[END_REF][START_REF] Lorz | Dirac mass dynamics in multidimensional nonlocal parabolic equations[END_REF] for mathematical properties of this kind of models. The infinitesimal model provides a simple and robust model describing inheritance of quantitative traits. It can be rigorously derived, see [START_REF] Barton | The infinitesimal model[END_REF], from a discrete model with Mendelian inheritance, mutation and environmental noise, when the genetic component of the trait is purely additive (non-epistasis case) or when epistasis is not consistent in direction, as the number of underlying loci tends to infinity.

In this model males and females are not distinguished. This makes senses, for instance, if they have the same distribution (e.g. for hermaphrodites). Mating is assumed to be random and uniform among the population. The number of offsprings is assumed to be proportional to the density (of females). Let us consider a density of individuals f (z) structured by a continuous phenotypical trait z ∈ R. Sexual reprodution involves two parents, having the traits z and z who give birth to an offspring having the trait z. In the infinitesimal model the genetic component of offspring traits follows a normal distribution centered the average of the parental traits z +z 2 , while this distribution has a variance that is independent of the parental values and which remains constant. The density of descendants is normalized to be proportional to the density of individuals f . This leads to consider the following mixing operator:

Gf (z) := 1 R f (x) dx R 2 f (z )f (z ) G z - z + z 2 dz dz , (1) 1
where G is the gaussian distribution with the given variance V LE /2, where V LE is the variance at linkage equilibrium, see [START_REF] Bulmer | The mathematical theory of quantitative genetics[END_REF][START_REF] Turelli | Genetic and statistical analyses of strong selection on polygenic traits: what, me normal?[END_REF][START_REF] Tufto | Quantitative genetic models for the balance between migration and stabilizing selection[END_REF][START_REF] Barfield | Evolution in stage-structured populations[END_REF][START_REF] Huisman | Comparison of non-gaussian quantitative genetic models for migration and stabilizing selection[END_REF].

As noticed in [START_REF] Turelli | Genetic and statistical analyses of strong selection on polygenic traits: what, me normal?[END_REF], and used in [START_REF] Mirrahimi | Dynamics of sexual populations structured by a space variable and a phenotypical trait[END_REF] for numerical simulations, the operator G can be defined as a double convolution product:

Gf (z) = 4 R f (x) dx (G * f (2•) * f (2•)) (z), (2)
or equivalently through its Fourier transform, with the convention f (ξ) = R e -ixξ f (x) dx:

Gf (ξ) = 1 f (0) G(ξ) f (ξ/2) 2 . ( 3
)
Beyond the biological background, we do not assume that the probability density function G is a gaussian function but we consider slightly more general cases.

The recombination-selection model without age structure

To model the adaptation of a population to an environment, let us consider a population density f (t, z) where t 0 is the time variable and z the trait variable. We assume the time variable t to be continuous. This implies that overlapping generations exist that is more than one breeding generation is present at any time. For the sake of simplicity we assume that the birth rate is constant and that selection is taken into account through the heterogeneous mortality rate µ(z).

These assumptions lead us to consider the following recombination-selection equation:

∂ t f (t, z) = -µ(z)f (t, z) + β G(f (t, •))(z), (4)
where G is defined by [START_REF] Abrams | The fitness costs of senescence: the evolutionary importance of events in early adult life[END_REF]. The model (4) will be referred to as the homogeneous model.

The death rate µ(z) is assumed to be a symmetric function in the trait variable z, so that the mortality is minimal at z = 0. We further assume that µ(0) = 0, without loss of generality. The term µf models selection: individuals having the optimal trait z = 0 have the minimal mortality rate, so are more likely to be maintained in the population. The second term β G(f ) is the birth term and models inheritance of the traits by the sexual reproduction, at birth. For µ = β constant, one can check that the size of the population is constant and the mean phenotypical trait is conserved. It has been shown, see [START_REF] Bulmer | The mathematical theory of quantitative genetics[END_REF][START_REF] Turelli | Genetic and statistical analyses of strong selection on polygenic traits: what, me normal?[END_REF], that (4) has a unique steadystate, that is the Gaussian distribution of variance σ 2 = V LE . Furthermore, it can be shown that the distribution f (t, •) converges exponentially fast to the Gaussian distribution in the sense of Wasserstein distances, see [START_REF] Magal | Dynamics of a kinetic model describing protein exchanges in a cell population[END_REF].

The age-structured recombination-selection model

The antagonistic pleiotropy theory explains senescence by the selection of genes with positive effects at young ages even though they can be deleterious later in life, see [START_REF] Williams | Pleiotropy, natural selection, and the evolution of senescence[END_REF]. This idea was formalized by Hamilton, [START_REF] Hamilton | The moulding of senescence by natural selection[END_REF], and developped by Charlesworth, cf. [START_REF] Charlesworth | Evolution in age-structured populations[END_REF][START_REF] Charlesworth | Evolution in age-structured populations[END_REF][START_REF] Charlesworth | Patterns of age-specific means and genetic variances of mortality rates predicted by the mutation-accumulation theory of ageing[END_REF], and recently numerically studied in [START_REF] Cotto | Maladaptation as a source of senescence in habitats variable in space and time[END_REF]. To model ageing of individuals we thus introduce an age variable a 0 in Equation (4), see [START_REF] Perthame | Transport equations in biology[END_REF] for classical age-structured PDEs. This leads to consider the age-structured recombination-selection equation:

     ∂ t f (t, a, z) + ∂ a f (t, a, z) + µ(a, z)f (t, a, z) = 0, f (t, 0, z) = G ∞ 0 β(a)f (t, a, •)da (z), (5)
where G is defined by [START_REF] Abrams | The fitness costs of senescence: the evolutionary importance of events in early adult life[END_REF]. The death rate µ(a, z) and the birth rate β(a) can depend on the age variable a. The first equation models that the age is increasing at speed one while the mortality is acting as for the homogeneous model [START_REF] Baudisch | Hamilton's indicators of the force of selection[END_REF]. The second equation models the inheritance of the trait at birth (age 0).

Examples. For different models, the age-specific effect of mutations are considered to affect a single age class or a range of age class. Different assumptions about the age-specific effect of mutations translate into different patterns of age-specific mortality, see [START_REF] Abrams | The fitness costs of senescence: the evolutionary importance of events in early adult life[END_REF][START_REF] Charlesworth | Patterns of age-specific means and genetic variances of mortality rates predicted by the mutation-accumulation theory of ageing[END_REF][START_REF] Baudisch | Hamilton's indicators of the force of selection[END_REF][START_REF] Moorad | A theory of age-dependent mutation and senescence[END_REF][START_REF] Wachter | The age-specific force of natural selection and biodemographic walls of death[END_REF]. In the previous age-continuous setting, some of these assumptions can be translated into considering death rates of the form:

µ(a, z) = |z| α δ a=a * , or |z| α 1 a a * , (6)
where α > 0 and a * > 0.

The principal eigenelements

As the operators f → βf and f → µG(f ) are 1-homogeneous, special solutions to (4) of the form f (t, z) = e λt F (z), resp. solutions to (5) of the form f (t, a, z) = e λt F (a, z), are characterized by pairs (λ, F ) solving the following spectral problems. The real number λ is the rate at which the number of individuals grows (λ 0) or decays (λ 0).

For the homogeneous model ( 4) the spectral problem is to find a pair (λ, F ) solution of:

     λF (z) + µ(z)F (z) = G (F ) (z), R F (z) dz = 1. (7)
For the age-structured model ( 5) the spectral problem is to find a pair (λ, F ) solution of:

               λF (a, z) + ∂ a F (a, z) + µ(a, z)F (a, z) = 0, F (0, z) = G ∞ 0 β(a)F (a, •) da (z), R F (0, z) dz = 1. (8)
As Problem ( 8) is 1-homogeneous with respect to F , we chose to normalize the density of newborns to be 1.

The operator G defined by (1) is not linear, so that the classical Kreȋn-Rutman theory, see [START_REF] Kreȋn | Linear operators leaving invariant a cone in a Banach space[END_REF][START_REF] Dautray | Analyse mathématique et calcul numérique pour les sciences et les techniques[END_REF], cannot be applied to this problem. The operator G is positively 1-homogeneous and compact on subsets of L 1 (R) with bounded variance, cf. proof of Theorem 2.1. But the operator G is not monotone, because of the division by R f dx, so the generalization of Kreȋn-Rutman theory to 1-homogeneous operators cannot be applied neither, see for instance [START_REF] Mahadevan | A note on a non-linear Krein-Rutman theorem[END_REF].

1. The death rate µ(z) is non-negative, symmetric, non-decreasing for z 0 and satisfies µ(0) = 0.

2. We further assume integrability conditions at 0 and ∞:

1 µ ∈ L 1 (-1, 1), (9)
1 µ ∈ L 1 (1, ∞). (10)
A typical example is µ(z) = |z| α with α > 1. Assumption (9) will be referred to as weak selection, in contrast to strong selection, cf. §3.3. It can be rephrased as: H µ (-1,1) is zero, where H is the harmonic mean defined by:

H(f ) = R dz f (z) -1 .
A similar weak selection assumption in the case of asexual reproduction is assumed in [START_REF] Martin | The Nonstationary Dynamics of Fitness Distributions: Asexual Model with Epistasis and Standing Variation[END_REF].

Theorem 2.1 (weak selection). Let G be a function satisfying the Assumptions on G. We consider the resulting operator G defined by (1). Let β be a positive number and µ a function satisfying the Assumptions on µ. Then, there exists a solution (λ, F ) to Problem [START_REF] Charlesworth | Evolution in age-structured populations[END_REF] with λ > 0 and F a nonnegative L 1 (R) function.

Assumption ( 9) can be relaxed. In this case it may happen that the limiting value λ = 0 is reached and F has a Dirac mass at the origin, see Figure 1. To include more singular behaviour at the origin, we consider the normed vector space:

E = Rδ + L 1 (R) = xδ + f, x ∈ R, f ∈ L 1 (R) , ( 11
)
where δ is the Dirac measure at {z = 0}. It is endowed with a product norm on R × L 1 (R).

Theorem 2.2 (strong selection).

Let G be a function satisfying the Assumptions on G. We consider the resulting operator G defined by (1). Let β be a positive number and µ a function satisfying the Assumption [START_REF] Cotto | Maladaptation as a source of senescence in habitats variable in space and time[END_REF]. Then, there exists a solution (λ, F ) to Problem (7) with λ 0 and F ∈ E, defined by [START_REF] Coville | Singular measure as principal eigenfunction of some nonlocal operators[END_REF], the singular and the regular parts of which are non-negative.

For a convolution operator a similar result is obtained in [START_REF] Coville | Singular measure as principal eigenfunction of some nonlocal operators[END_REF].

These recombination-selection equilibria are built as fixed point of some operator linked to G. Assuming that G and µ are symmetric allows us to build symmetric equilibria. This operator is shown to stabilize the variance, thus G is needed to have a finite second moment, cf. Lemma 3.1. Also, the key Lemma 3.2 shows that the multiplication by the mortality rate µ does not affect this property if µ is a symmetric and non-decreasing (for z 0) function. Lastly, the uniform continuity of G and the integrability of 1/µ at infinity are needed to establish the compacity of this operator, cf. Proofs of the Theorems.

The age-structured model, assumptions and statements

Assumptions on G. The probability density function G is symmetric and it has a finite second moment.

Assumptions on β and µ.

1. The reproduction rate β and the death rate µ(a, z) are non-negative.

2. The death rate z → µ(a, z) is symmetric, non-decreasing for z 0 and satisfies µ(a, 0) = 0.

3. Our last and main assumption links the reproduction rate β and the death rate µ. There exists λ 0 ∈ R such that:

∀z ∈ R 1 ∞ 0 β(a)e -λ0a e -a 0 µ(a ,z) da da ∞ 0 β(a)e -λ0a da < ∞. (12)
Assumption ( 12) is a balance assumption between birth and death rates which is uniform in the trait variable. As the weak selection assumption (9) for the homogeneous case, Assumption ( 12) avoids concentration at z = 0. For the death rates given by ( 6) and a constant birth rate β, Assumption ( 12) is satisfied as soon as a * > 0. From a technical viewpoint, Assumption [START_REF] Dautray | Analyse mathématique et calcul numérique pour les sciences et les techniques[END_REF] ensures the existence of eigenvalues bounded from below by λ 0 , cf. Lemma 3.3. Theorem 2.3. Let G be a function satisfying the Assumptions on G. We consider the resulting operator G defined by [START_REF] Abrams | The fitness costs of senescence: the evolutionary importance of events in early adult life[END_REF]. Let β(a) and µ(a, z) be functions satisfying the Assumptions on β and µ. Then, there exists a solution (λ, F ) to system (8) with λ λ 0 .

Proofs of the existence theorems

We first prove Theorem 2.1 because some estimates are easier and its proof contains all the steps necessary for the proof of Theorem 2.3. Namely, the dispersion relation gives implicitly the principal eigenvalue corresponding to a principal eigenvector. This allows us to define a mixing operator, the fixed points of which are principal eigenelements. Proving the theorem comes down to build a convex set conserved by this operator in order to apply the Schauder fixed point theorem.

Confinement and technical lemmas

The proofs of these theorems rely on the variance stabilizing property of the operator G, i.e. G has a confinement effect in the trait variable. The proofs do not use the confinement effect of the selection operator f → µf . Lemma 3.1 describes how the infinitesimal operator acts on the first moments. Lemma 3.2, which is the key lemma, implies that the multiplication by a symmetric, non-increasing function reduces the variance. The third and last lemma is a short technical lemma concerning the Laplace transform. It is only used in the proof of Theorem 2.3. Lemma 3.1. (first moments of the operator G)

1. Let f ∈ L 1 (R) with R f (z) dz = 0. Then: R Gf (z) dz = R f (z) dz, and R zGf (z) dz = R zf (z) dz.
2. Let f be an L 1 (R) function satisfying: R f (z) dz = 1. The following equality holds:

R (Gf (z) -f (z)) z 2 dz = Var G - 1 2 Var f,
where

Var ϕ = R z 2 ϕ(z) dz - R zϕ(z) dz 2 .
In addition if R zf (z) dz = 0, one obtains:

Var G(f ) = Var G + 1 2 Var f. As a consequence, if Var f 2Var G then Var G(f ) 2Var G.
Proof. For any function ϕ, from (1), using R G(z) dz = 1, one obtains:

R (Gf (z) -f (z)) ϕ(z) dz = 1 R f (x) dx R 3 f (z ) f (z ) G(z) ϕ z + z + z 2 - 1 2 (ϕ(z ) + ϕ(z )) dz dz dz. (13)
1. The first point is obtained for ϕ(z) = 1, and ϕ(z) = z using R zG(z) dz = 0. 2. Again, for ϕ(z) = z 2 , equality (13) leads to the result. After integrating the equality:

z + z + z 2 2 - 1 2 z 2 + z 2 = z 2 - 1 4 z 2 - 1 4 z 2 + zz + zz + 1 2 z z ,
one gets:

R (Gf (z) -f (z)) z 2 dz = R z 2 G dz R f dz -2 1 4 R z 2 f dz R G dz + 2 R zG dz R zf dz + 1 2 R zf dz 2 R f dz R G dz = Var G - 1 2 Var f. Lemma 3.2.
Let ν be a measurable, non-negative, symmetric on R, non-increasing on [0, +∞) function. Let W be a W 1,1 loc , non-negative, symmetric on R, non-decreasing on [0, +∞) function. The quadratic form

Q(ϕ) = R νϕW dx R ϕ dx - R ϕW dx R νϕ dx
is non-positive on the set:

D = ϕ 0 : ϕ, νϕ ∈ L 1 (R, dx) ∩ L 1 (R, W dx) .
Proof. The form Q is linear with respect to W and vanishes if W is a constant function. Thus considering W -W (0) we can assume that W (0) = 0.

Let us define the bilinear symmetric form:

B(ϕ, ψ) = 1 2 R νψW dx R ϕ dx + R νϕW dx R ψ dx - R ϕW dx R νψ dx - R ψW dx R νϕ dx ,
such that: Q(ϕ) = B(ϕ, ϕ). As: B(ϕ, ψ) = 0, if ϕ is odd, considering the decomposition of a function ϕ = ϕ e + ϕ o as the sum of its even and its odd components, leads to:

Q(ϕ) = Q(ϕ e ) + B(ϕ o , ϕ o ) + 2B(ϕ e , ϕ o ) = Q(ϕ e ).
Thus the functions in D can be assumed to be even. If ϕ is even, then:

1 4 Q(ϕ) = ∞ 0 νϕW dx ∞ 0 ϕ dx - ∞ 0 ϕW dx ∞ 0 νϕ dx. (14)
For a non-negative function ϕ : [0, ∞) → R which belongs to L 1 (R) let us denote:

F ϕ (x) = ∞ x ϕ(z) dz.
For ϕ ∈ L 1 (W dx), as W is non-decreasing, in the limit x → ∞, we have:

F ϕ (x)W (x) ∞
x ϕW dz → 0. By integration by parts, using W (0) = 0, we get:

∞ 0 ϕW dx = ∞ 0 F ϕ dW dx dx.
Thus, ( 14) can be changed to:

1 4 Q(ϕ) = ∞ 0 (AF νϕ (x) -BF ϕ (x)) dW dx dx, with A = ∞ 0 ϕ dx, B = ∞ 0 ϕν dx.
As dW dx 0, it amounts to show that I = AF νϕ -BF ϕ is non-positive. The derivative of -I is:

J := - dI dx = (Aν -B)ϕ.
The function J satisfies ∞ 0 J(z) dz = I(0) -I(∞) = 0 so that the function Aν -B changes of sign. As Aν -B is non-increasing there exists x 0 ∈ [0, ∞) such that: J(x) 0 for x x 0 and J(x) 0 for x x 0 . Hence, I is non-increasing on [0, x 0 ] and non-decreasing on [x 0 , ∞), hence non-positive because I(0) = I(∞) = 0. Lemma 3.3. Let ϕ : R + → R be a measurable function, which is non-negative and non-zero. For p ∈ R let us define:

L(p) = ∞ 0 e -pz ϕ(z) dz ∈ [0, ∞]. If there exists λ 0 ∈ R such that 1 L(λ 0 ) < ∞, then there exists a unique p ∈ R such that L(p) = 1.
Proof. The function L(p) is finite for p λ 0 . Let us define the abscissa of convergence

p 0 = inf{p ∈ R | L(p) < ∞} ∈ [-∞, λ 0 ].
Clearly, the function L is decreasing, analytic on (p 0 , ∞) and: lim p→+∞ L(p) = 0. Therefore there is at most one p ∈ R such that L(p) = 1. The existence of such an element is equivalent to L(p 0 ) 1.

The recombination-selection equation without age structure, weak selection case

Proof of theorem 2.1.

1. For λ > 0 and z ∈ R, or λ 0 and z = 0, let us denote:

ν(λ, z) = β λ + µ(z) . (15)
A solution (λ, F ) to [START_REF] Charlesworth | Evolution in age-structured populations[END_REF], with λ > 0, can be expressed as a fixed point using ν:

F = ν(λ, •) G(F ).
Let us consider the following convex subset of L 1 (R):

C = F ∈ L 1 (R) F 0, R F (z) dz = 1, F (-z) = F (z), R |z| 2 F (z) dz 2Var G . (16)
The set C is a nonempty, bounded, closed, convex subset of L 1 (R). Regularization properties of the operator G, defined as a double convolution (2), makes the application of the Schauder fixed point theorem possible to the set C and the operator M defined by [START_REF] Magal | Dynamics of a kinetic model describing protein exchanges in a cell population[END_REF]. The construction of M is based on the following Claim.

2. Claim. For any function F in C there exists a unique real number λ such that:

R ν(λ, z) G(F )(z) dz = 1. ( 17
)
This defines a function: F → λ(F ). The function λ is bounded from below by a constant ε > 0 which depends only on G and µ.

The function (0, ∞) λ → R ν(λ, z) G(F )(z) dz ∈ (0, ∞) is decreasing and tends to 0 at infinity. So, it is enough to prove that its limit at 0 is infinite. Uniformly for F in C:

|z|>r F (z) dz 1 r 2 |z|>r |z| 2 F (z) dz 1 r 2 2Var G → 0, (18)
so that there exists r > 0 satisfying: |z| r F dz

1 2
. By definition of G, (1), for all F in C and z such that |z| r, we obtain:

G F (z) (z ,z )∈[-r,r] 2 F (z )F (z ) G z - z + z 2 dz dz inf |z| 2r G |z| r F (z) dz 2 inf |z| 2r G × 1 4 .
So, we get:

R ν(λ, z) G(F )(z) dz |z| r ν(λ, z) G(F )(z) dz 1 4 inf |z| 2r G |z| r β λ + µ(z) dz.
Thanks to Fatou's lemma and Assumption ( 9) the last quantity tends to +∞ as λ goes to 0 + . In particular there exists ε > 0 such that:

1 4 inf |z| 2r G |z| r β ε + µ(z) dz 1, ( 19
)
and we get the uniform lower bound for λ: λ(F ) ε.

3. We define the mixing operator:

M(F ) = ν(λ(F ), •) G(F ). ( 20
)
Proving Theorem 2.1 is equivalent to show that the operator M has a fixed point.

•

(z) = z 2 : Var M(F ) Var G(F ) = Var G + 1 2
Var F 2Var G.

• The set M(C) is relatively compact thanks to the Fréchet-Kolmogorov theorem. As M(C) ⊂ C, tails are uniformly bounded thanks to [START_REF] Kreȋn | Linear operators leaving invariant a cone in a Banach space[END_REF]. Let τ h be the translation by h defined by:

τ h F = F (• -h).
From the definition of M, (20), we get:

τ h (MF ) -MF L 1 (R) = (τ h (ν(λ(F ), •)) -ν(λ(F ), •)) τ h (GF ) + ν(λ(F ), •) (τ h (GF ) -GF ) L 1 (R) , τ h (ν(λ(F ), •)) -ν(λ(F ), •) L 1 (R) GF L ∞ (R) + ν(λ(F ), •) L 1 (R) τ h (GF ) -GF L ∞ (R) (21)
As λ(F ) ε, we get:

|τ h ν(λ(F ), •) -ν(λ(F ), •)| = |µ -τ h µ| (λ(F ) + τ h µ)(λ(F ) + µ) |µ -τ h µ| (ε + τ h µ)(ε + µ) = |τ h ν(ε, •) -ν(ε, •)| . (22)
In addition, by definition (2), we have:

τ h (GF ) -GF = 4(τ h G -G) * F (2•) * F (2•). ( 23
)
Using Young's inequality we get the uniform upper bound:

τ h (MF ) -MF L 1 (R) τ h ν(ε, •) -ν(ε, •) L 1 (R) G L ∞ (R) + ν(ε, •) L 1 (R) τ h G -G L ∞ (R) ,
which goes to 0 as h tends to 0 because: ν(ε, •) is in L 1 (R), see Assumption [START_REF] Cotto | Maladaptation as a source of senescence in habitats variable in space and time[END_REF], and G is uniformly continuous.

• The operator M : C → C is continuous. Let F 1 , F 2 ∈ C, and let us denote:

λ k = λ(F k ), ν k = ν(λ k , •), for k = 1, 2.
Using [START_REF] Lorz | Dirac mass dynamics in multidimensional nonlocal parabolic equations[END_REF], we have λ k ε, and:

ν k L 1 (R) ν(ε, •) L 1 (R) < ∞.
Writing:

M(F 1 ) -M(F 2 ) L 1 (R) = 4 (ν 1 -ν 2 ) G(F 1 ) + ν 2 (G(F 1 ) -G(F 2 )) L 1 (R) 4 ν 1 -ν 2 L 1 (R) G L ∞ (R) + 4 ν(ε, •) L 1 (R) G(F 1 ) -G(F 2 ) L ∞ (R) . ( 24
)
it remains to show that the maps

C F → G(F ) ∈ L ∞ (R) and C F → ν ∈ L 1 (R) are continuous.
Thanks to definition (2), we have:

G(F 1 ) -G(F 2 ) = 4 G * (F 1 -F 2 )(2•) * (F 1 + F 2 )(2•). ( 25
)
This implies that the first map is continuous thanks to Young's inequality.

For the continuity of the second map we prove, as an intermediate stage, that C F → λ ∈ R is continuous. To do so let us fix F 1 , define λ = λ 2 -λ 1 , and:

I(λ) := R (ν 1 -ν 2 ) G(F 1 ) dz . (26)
By definitions of the

ν k , F k : R ν 1 G(F 1 ) dz = 1 = R ν 2 G(F 2
) dz, and:

I(λ) = R (ν 1 -ν 2 ) G(F 1 ) dz = R ν 2 (G(F 2 ) -G(F 1 )) dz ν(ε, •) L 1 (R) G(F 2 ) -G(F 1 ) L ∞ (R) (27)
tends to 0 as F 1 -F 2 L 1 (R) goes to 0.

We have:

I(λ) = R (ν 1 -ν 2 ) G(F 1 ) dz = R |λ| (λ 1 + µ)(λ 1 + µ + λ) G(F 1 ) dz, (28)
where the last equality holds because the integrand has a constant sign.

It is easily checked from (28) that the function I satisfies: I(0) = 0, it is continuous, strictly decreasing for λ 0 and strictly increasing for λ 0. In particular: as I(λ) goes to 0, |λ| tends to 0.

From [START_REF] Tufto | Quantitative genetic models for the balance between migration and stabilizing selection[END_REF], we have: |λ| = |λ 2 -λ 1 | goes to 0. This allows us to conclude thanks to Lebesgue's dominated convergence theorem and Assumption [START_REF] Cotto | Maladaptation as a source of senescence in habitats variable in space and time[END_REF], because:

ν 1 -ν 2 = (λ 1 + µ) -1 -(λ 2 + µ) -1 . (29)

The recombination-selection equation without age structure, strong selection case

Proof of theorem 2.2. 1. A solution (λ, F ) to ( 7) can be expressed as a fixed point using ν defined by ( 15): F = ν(λ, •) G(F ). Indeed the pair (0, xδ + f ) ∈ E, where E is defined by [START_REF] Coville | Singular measure as principal eigenfunction of some nonlocal operators[END_REF], is a solution to [START_REF] Charlesworth | Evolution in age-structured populations[END_REF] if and only if µf = G(xδ + f ) because µ(0) = 0.

2. Under the weak selection assumption [START_REF] Charlesworth | Patterns of age-specific means and genetic variances of mortality rates predicted by the mutation-accumulation theory of ageing[END_REF], equation ( 17) for λ 0, has a solution for every

F ∈ C ⊂ L 1 (R) because: R ν(0, z) G(F )(z) dz = +∞.
If assumption (9) does not hold true, it may hapen that R ν(0, z) G(F )(z) dz < 1 for some F ∈ L 1 (R). Thus equation ( 17) may have no solution. In such a case λ = 0 and F may have a Dirac part at the origin. This Dirac mass compensates the lack of mass of the regular part.

First, this leads to consider the following convex subset of E:

C = mδ + F ∈ E m 0, F 0, m + R F (z) dz = 1, F (-z) = F (z), R |z| 2 F (z) dz 2Var G . (30)
For xδ + f in E such that x + R f dz = 0, we have:

G(xδ + f ) = 1 x + R f dz x 2 G + 4x G * f (2•) + 4 G * f (2•) * f (2•) . (31)
Thus, G(xδ

+ f ) is an L 1 (R) ∩ L ∞ (R) function.
Second, for F in C, we define the mixing operator M by:

MF = ν(λ, •) G(F ) if m = 1 - R ν(0, •) G(F ) dz < 0, ν(0, •) G(F ) + mδ if m 0, (32)
where, in the first case, λ > 0 is still defined by equation [START_REF] Huisman | Comparison of non-gaussian quantitative genetic models for migration and stabilizing selection[END_REF], while in the second case, λ = 0.

As a by-product of this definition, we have m × λ = 0, so that m acts as a degree of freedom when λ = 0.

The Schauder fixed point theorem is applied to the operator M on the convex set C, to prove Theorem 2.2. Let us denote p 1 , p 2 the canonical projections:

p 1 : E → R, p 2 : E → L 1 (R).
• The set C is a nonempty, bounded, closed, convex subset of E. It is preserved by M thanks to definition (32).

• The set p 1 ( M( C)) ⊂ [0, 1] is bounded and p 2 ( M( C)) ⊂ p 2 ( C) is relatively compact thanks to the Fréchet-Kolmogorov theorem. Thanks to [START_REF] Kreȋn | Linear operators leaving invariant a cone in a Banach space[END_REF], the tails of the elements of p 2 ( M( C)) are uniformly bounded. As in the previous case, from the definition of M, (32), we get:

τ h p 2 ( MF ) -p 2 ( MF ) L 1 (R) = τ h ν(λ, •) GF -ν(λ, •) GF L 1 (R) τ h ν(0, •) -ν(0, •) L 1 (R) G L ∞ (R) + ν(0, •) L 1 (R) τ h G -G L ∞ (R) ,
because λ 0. This expression goes to 0 as h tends to 0 because: ν(0, •) is in L 1 (R), see Assumption [START_REF] Cotto | Maladaptation as a source of senescence in habitats variable in space and time[END_REF], and G is uniformly continuous.

• The operator M : C → C is continuous. Let F 1 , F 2 ∈ C, and let us denote:

λ k = λ(F k ), ν k = ν(λ k , •), for k = 1, 2. We have: λ k 0 and ν k L 1 (R) ν(0, •) L 1 (R) < ∞.
First, the map C F → G(F ) ∈ L ∞ (R) is continuous as a consequence of (31), ( 25) and Young's inequality. Then, definition (32) yields:

p 1 ( MF 1 ) -p 1 ( MF 2 ) = 1 -ν(0, •) GF 1 L 1 (R) + -1 -ν(0, •) GF 2 L 1 (R) + ν(0, •) GF 1 L 1 (R) -ν(0, •) GF 2 L 1 (R) ν(0, •) (GF 2 -GF 1 ) L 1 (R) ν(0, •) L 1 (R) GF 2 -GF 1 L ∞ (R) , thus, p 1 • M : C → R is continuous.
Thanks to [START_REF] Mirrahimi | Dynamics of sexual populations structured by a space variable and a phenotypical trait[END_REF], proving the continuity of p 2 • M comes down to check the continuity of the maps

C F → G(F ) ∈ L ∞ (R) and C F → ν ∈ L 1 (R).
The continuity of the first map has just been proved. Using [START_REF] Wachter | The age-specific force of natural selection and biodemographic walls of death[END_REF], the continuity of the second map is a consequence of the continuity of C F → λ ∈ R, which is obtained using [START_REF] Turelli | Genetic and statistical analyses of strong selection on polygenic traits: what, me normal?[END_REF].

The age-structured recombination-selection equation

The proof of Theorem 2.3 follows the same steps. The differences are that a first reduction step is needed but no lower bound on the eigenvalue has to be proved because it is provided by Assumption [START_REF] Dautray | Analyse mathématique et calcul numérique pour les sciences et les techniques[END_REF].

Proof of theorem 2.3.

1. A solution to equation ( 8) is determined by the phenotypical distribution at age 0:

P = F (0, •), (33)
through: F (a, z) = F (0, z)s(a, z)e -λa , where:

s(a, z) = e -a 0 µ(a ,z) da , (34)
is the survival function.

Furthermore the function P can be expressed as a fixed point: P = G (ν(λ, •) P ) , where:

ν(λ, z) = ∞ 0 β(a)s(a, z)e -λa da. (35)
2. Claim. For any non-negative function P satisfying R P dz = 1, there exists a unique λ ∈ R such that: R ν(λ, z)P (z) dz = 1.

(36)

This defines a function: P → λ(P ) ∈ R.

This claim is a consequence of Lemma 3.3 applied to the function a → β(a) R s(a, z)P (z) dz, because Assumption (12) implies:

1 = R P (z) dz R ∞ 0 e -λ0a β(a)s(a, z)P (z) dz da ∞ 0 e -λ0a β(a) da < ∞,
using s(a, z) 1. We get the uniform lower bound λ(P ) λ 0 .

3. We define the mixing operator:

M(P ) = G (ν(λ(P ), •) P ) . (37)
As previously, the Schauder theorem on C, defined by [START_REF] Williams | Pleiotropy, natural selection, and the evolution of senescence[END_REF], is applied to conclude.

• The set C is preserved by M. The last estimate, is a consequence of Lemmas 3.1, 3.2 with ν = ν(λ, •) and W (z) = z 2 :

Var

M(P ) = Var G + 1 2 Var ν(λ(P ), •) P Var G + 1 2
Var P 2Var G.

• The set M(C) is relatively compact thanks to the Fréchet-Kolmogorov theorem. The estimate on τ h (MP ) -MP L 1 (R) is easier to obtain:

τ h (MP ) -MP L 1 (R) τ h G -G L 1 (R) ν(λ(P ), •)P L 1 (R) ν(λ(P ), •)P L 1 (R) = τ h G -G L 1 (R) → 0.
• The operator M : C → C is continuous. Let P 1 , P 2 ∈ C, and let us denote:

λ k = λ(P k ), ν k = ν(λ k , •), for k = 1, 2. Assumption (12) implies λ k λ 0 and ν k L ∞ (R) ν(λ 0 , •) L ∞ (R) ∞ 0 e -λ0a β(a) da < ∞.
Thanks to (2) and ( 25) we have:

M(P 1 ) -M(P 2 ) L 1 (R) G L 1 (R) (ν 1 P 1 ) * (ν 1 P 1 ) -(ν 2 P 2 ) * (ν 2 P 2 ) L 1 (R) = [(ν 1 -ν 2 )P 1 + ν 2 (P 1 -P 2 )] * (ν 1 P 1 + ν 2 P 2 ) L 1 (R) ν 1 -ν 2 L ∞ (R) P 1 L 1 (R) + ν 2 L ∞ (R) P 1 -P 2 L 1 (R) 2 ∞ 0 e -λ0a β(a) da It remains to prove that L 1 (R) P → ν ∈ L ∞ (R) is continuous. Let us fix P 1 , define λ = λ 2 -λ 1 , and: i(x) = |1 -e -x |.
As previously, we have:

R (ν 1 -ν 2 )P 1 dz = R ν 2 (P 1 -P 2 ) dz P 1 -P 2 L 1 (R) ∞ 0 e -λ0a β(a) da. (38)
Let us define the function I(λ) by:

I(λ) := R ∞ 0 i(λa) e -λ1a P 1 (z)s(a, z)β(a) da dz = R ∞ 0 1 -e -(λ2-λ1)a e -λ1a P 1 (z)s(a, z)β(a) da dz = R (ν 1 -ν 2 )P 1 dz . (39)
The equalities hold because a → e -λ1a -e -λ2a has a constant sign.

It is easily checked that the function I enjoys the following properties of i: I(0) = 0, it is continuous, decreasing for λ 0 and increasing for λ 0. In particular: as I(λ) goes to 0, |λ| tends to 0.

As P 1 -P 2 L 1 (R) goes to 0 we get that |λ| = |λ 2 -λ 1 | tends to 0 thanks to (38), (39). Lastly:

ν 1 -ν 2 L ∞ (R)
∞ 0 e -λ1a -e -λ2a β(a) da tends to 0 as |λ 1 -λ 2 | goes to 0.

Conclusions and perspectives

As the Kreȋn-Rutman theory could not be applied to Problems [START_REF] Charlesworth | Evolution in age-structured populations[END_REF], [START_REF] Charlesworth | Evolution in age-structured populations[END_REF], because the operator G is nor linear nor monotone, we used the Schauder fixed point theorem to prove the existence of principal eigenelements. Another difficulty was that the dispersion relation was not explicit. Using regularizing effects of the operator G and properties of moments we were able to build some convenient convex subset of L 1 (R) on which the Schauder theorem was applied, to prove the joint existence of a principal eigenvector and of the corresponding principal eigenvalue.

The previous approach can be straightforwardly adapted to the case of a genetic component of offspring traits z being a convex combination of the parental traits z , z : z = αz + βz with α + β = 1 and 0 < α, β < 1. The resulting operator G is defined as:

Gf (z) = 1 R f (x) dx R 2 f (z )f (z ) G (z -(αz + βz )) dz dz , = 1 αβ 1 R f (x) dx G * f • α * f • β (z).
This operator conserves the 0th order and the 1st order moments, and the main confinement property remains true: for ϕ ∈ L 1 (R) satisfying R ϕ(z) dz = 1 and R zϕ(z) dz = 0, we have:

Var G(ϕ) = Var G + (α 2 + β 2 )Var ϕ, with α 2 + β 2 < 1.

It would be interesting to include more general birth rates than the constant ones in our analysis. Weither the principal eigenelements are unique or not is also an interesting question. As concerns the asymptotic behaviour of Equations ( 4), [START_REF] Bulmer | The mathematical theory of quantitative genetics[END_REF], convergence to natural attractors could be investigated: e -λt f (t, z) → F (z), resp. e -λt f (t, a, z) → F (a, z), in some weighted norm, where F (z), F (a, z) are defined by [START_REF] Charlesworth | Evolution in age-structured populations[END_REF], [START_REF] Charlesworth | Evolution in age-structured populations[END_REF].

In the context of a changing environment, see [START_REF] Cotto | Maladaptation as a source of senescence in habitats variable in space and time[END_REF] for instance, we can assume that the optimal trait moves at a speed c, by changing the mortality rate µ(z) by µ(z -ct) in Equations ( 4), [START_REF] Bulmer | The mathematical theory of quantitative genetics[END_REF]. For instance, for the homogeneous model, in the moving frame z = z -ct, Problem (7) writes:

     λF (z) -c ∂ z F (z) + µ(z)F (z) = G (F ) (z), R F (z) dz = 1.
For this problem, the function F cannot be supposed to be symmetric. Thus our confinement proof cannot be applied to this case. Other models could also be studied: cyclic changes of environment, space structured models with migration, see [START_REF] Cotto | Maladaptation as a source of senescence in habitats variable in space and time[END_REF].

Figure 1 :

 1 Figure 1: The selection-mutation equilibrium F (z) obtained from (4), for µ(z) = |z| α + |z| 2 . It may have a Dirac mass at z = 0 when the limiting eigenvalue λ = 0 is reached.

  The set C is preserved by M. Let F ∈ C. Clearly from (1), MF 0 and thanks to the definition of λ(F ), R MF (z) dz = 1. The equality M(F )(-z) = M(F )(z) is a consequence of the symmetry of the death rate µ and the probability density function G. The last point is a consequence of Lemmas 3.1, 3.2 with ν = ν(λ, •) and W

Acknowledgments. The authors thank Dr. É. Bouin for fruitful discussions. T.B's, V.C.'s, J.G.'s research is supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme, ERC Starting Grant MESOPROBIO No. 639638.


Gourdon
email: gourdon@insp.jussieu.fr

S Shihab

L Thevenard

A Lemaître

C Gourdon

Counter-rotating standing spin-waves: a magneto-optical illusion

Keywords: numbers: 75.78.Jp, 75.30.Ds, 75.50.Pp

Since pioneering work on nickel [1], laser-induced magnetization dynamics has been widely used to investigate ultrafast magnetic processes not only in magnetic metals [2,3], but also in magnetic semiconductors [4][5][6][7][8][9] and insulators [10,11], exploring the fundamentals of light-spin interaction in view of a full and ultrafast optical control of magnetic order.

Ultrashort pulses can trigger a wide variety of processes, including ultrafast demagnetization [1], full magnetization reversal [12] as well as coherent precession [2,4,10]. In magnetically ordered materials, ferro-and ferrimagnets, as well as in antiferromagnets, the coherent magnetization dynamics arises from collective spin excitations, spin waves (SW) (or magnons, their quanta), which attract a considerable interest motivated by their possible use as information carriers in magnonics applications [13,14]. Magnons are versatile excitations since their dispersion curves, comprising magnetostatic and exchange modes [15], can be tuned by a magnetic field or by micro-or nanostructuring the material in any of its dimensions [14]. For instance, perpendicular standing spin-wave (PSSW) modes in a single nanometric layer of thickness L with free boundary conditions (no surface anisotropy) will have their wave-vector quantized by an integer p (k = pπ/L), and their energy by p 2 , proportionally to the spin stiffness D.

SWs can be studied in the frequency domain by, ments using two pump beams can be performed [16,17].

55

Furthermore the possibility to fully reconstruct the mag-56 netization trajectory using two different magneto-optical 57 effects [18,19] brings a deep insight into magnetization 58 dynamics.

59

Whereas the excitation mechanism of SWs by optical, 60 acoustical or magnetic field pulses [14,[20][21][22], has been 61 thoroughly discussed, their optical detection has been 62 much less addressed [4,23]. In particular, in contrast to 63 the cavity FMR detection of SWs which was modeled a 64 long time ago [24], the respective amplitudes of optically 65 detected SWs remained unexplained [2,9,14,25]. In this 66 Letter, we provide a comprehensive model to explain the 67 large amplitude of these non-uniform SWs that should 68 not be detectable in the framework of simple models [4].

69

We moreover present intriguing experimental results of The existence of circular and linear magnetic birefringence/dichroism [31] makes the TRMO signal sensitive to both the out-of-plane δθ and the in-plane δφ components of the transient magnetization (Fig. 1(a)). This allows for the reconstruction of the magnetization trajectory using the expected dependency of the rotation angle δβ r on the probe polarization angle β (Ref. 19 and Suppl-info):

δβ exp r (t) =K r δθ exp (t) + 2V r δφ exp (t) cos 2(β -φ 0 ) -2V r δM (t) M sin 2(β -φ 0 ) , (1)
where K r and V r are the static Kerr and Voigt rotation coefficients, M is the magnetization vector modulus, and φ 0 is the in-plane equilibrium angle of the magnetization.

146

The SW space and time profiles are obtained by solving the LLG equation within the small precession angle approximation for in-plane static magnetization:

δ φ = γ F θθ δθ -D ∂ 2 δθ ∂z 2 + α G δ θ δθ = -γ F φφ δφ -D ∂ 2 δφ ∂z 2 + δB exc (z, t) -α G δ φ . (2)
F is equal to E/M where E is the magnetic anisotropy 147 energy [28]. α G is the Gilbert damping. profiles at any t for this mode (see Suppl-info). Under symmetric conditions the p-PSSW eigenmodes (p=0, 1, 2...) are even (odd) with respect to the layer mid-plane for even (odd) p. δB exc (z, t) is the optically induced effective field that launches the magnetization precession.

It arises from the in-plane rotation of the magnetic easy axis induced by transient thermal effects or by the optical spin-orbit torque [32]. δB exc (z, t) is taken as a product of time and space functions f (t)g(z). f (t) is chosen so that the calculated δθ(t) and δφ(t) match the experimental ones. g(z) is taken as a Fourier series over the PSSW eigenmodes (see Suppl-info). The depth dependence of the magnetization trajectory for the p=0, 1, 2 SW modes (with time dependence cos(2πf p t+φ p )) is shown in Fig. 2 (a-c), respectively. It is seen that inside the layer the direction of rotation is the same for the three modes. For modes 1 and 2 the magnetization vector experiences a π-shift at each node, but its direction of rotation does not change with z (Fig. 2(b,c)).

The TRMO signal is then calculated using a multilayer and transfer matrix model to obtain the Kerr and Voigt rotation angle and ellipticity. The magnetic layer of thickness L is divided in N sub-layers with magnetization components M (m x,y + δm x,y,z (z i , t)) (Fig. S1 of Supplinfo). The calculation is performed for normal incidence of light along the z-direction and linear polarization making an angle β with the x-axis. The theoretical dynamical polarization rotation δβ th r (t) is obtained by taking the limit of an infinite number of sub-layers (N → ∞). δβ th r (t) is the sum of the Kerr and Voigt rotation angles, δβ th r (t) = δβ th Kr (t) + δβ th Vr (t) cos 2(β -φ 0 ) with:

δβ th Kr (t) = - 4π λ Re n 2 Q n 2 -1 L 0 e i 4πnz λ δθ(z, t)dz δβ th Vr (t) = 4π λ Im n 2 B (n 2 -1) L 0 e i 4πnz λ δφ(z, t)dz , (3)
where pointed out [33,34] and evidenced experimentally in the 198 90s [35] and recently in GaMnAs layers [36]. Similar 199 ideas were at play when conceiving magneto-optical sen-200 sors using magnetic quantum wells in optical cavities [37] 201 or magneto-photonic crystals with enhanced Faraday ro-202 tation [38,39]. However, these ideas had so far not been

B = B 1 + Q 2 . Q (∝ M ) and B 1 (∝ M 2 )

203

applied to the time-resolved optical detection of PSSWs 204 in ferromagnetic layers.

205

The dynamical rotation angles δβ th Kr (t) and δβ th Vr (t) are calculated according to Eq. 3 with the real and imaginary parts of Q and B extracted from the static rotation and ellipticity using Eqs.

S12 and S13 of Suppl-info. In order to compare the optically detected magnetization trajectory and the theoretical one we define δθ opt (z, t) and δφ opt (z, t) so that

δβ th r (t) = K r δθ opt z (t) + 2V r δφ opt z (t) cos 2(β -φ 0 ) where • • • z = (1/L) L 0 • • • (z)
dz denotes the average value over the layer thickness, K r and V r are given by Eqs. S11 and S18 of Suppl-info, respectively and

δθ opt (z, t) = - 1 K r Re φ opt nQ (n 2 -1) e i 4πnz λ δθ(z, t) , (4)
δφ opt (z, t) = 1 2V r Im φ opt nB (n 2 -1) e i 4πnz λ δφ(z, t) , (5)
with φ opt = 4πnL/λ the complex optical phase.

224

In order to explain why odd modes may exhibit an apparently inverted direction of rotation when detected optically, we take a simplified model and calculate the sign of the amplitude ratio δθ opt (t)/δφ opt (t). The damped p-PSSW mode at frequency f p is expressed as δθ p (z, t) = a θ p exp (-χ p t) cos (2πf p t) cos (pπz/L), δφ p (z, t) = a φ p exp (-χ p t) sin (2πf p t) cos (pπz/L). Ne- glecting absorption, the optical precession amplitudes normalized to the excitation amplitudes a θ p and a φ p are:

δθ opt p = - 4π K r η η 2 -1 0 Re Q e i4πu cos (pπu/ ) du δφ opt p = 4π 2V r η η 2 -1 0
Im B e i4πu cos (pπu/ ) du , (6) where = L/(λ/η). It is straightforward to show that for even p (even modes) the ratio r p = δθ opt p /δφ opt p is positive and equal to 1. For odd modes, r p can on the contrary be positive or negative depending on the layer thickness and the ratios B i /B r and Q i /Q r of the imaginary and real parts of B and Q, respectively. It is moreover independent of p and is given by r

odd p = -( Bi Br C l + S l )( Qi Qr C l + S l )(C l -Bi Br S l ) -1 (C l -Qi Qr S l ) -1
with C l = cos(2π ) and S l = sin(2π ). Therefore the possibility to change the sign of only one of the δθ and δφ components (r p <0) and hence to observe a change of the direction of rotation is achieved exclusively for the odd SW modes. This result is not fully conserved when taking into account absorption as can be seen in Fig. 3 where the direction of rotation (CW, CCW) given from the sign of r p is plotted in (dark/light) gray scale. However, given our parameters, r 1 is always negative for L in the range 26-72 nm encompassing the layer thickness of our sample (50 nm) while r 0 and r 2 are positive. This is an important result of this paper as it provides a tool to identify PSSW modes.

This model also accounts very well for the large amplitude ratio of the high/low frequency modes observed

  e.g., ferromagnetic resonance (FMR) experiments or Brillouin light scattering as well as in the time-domain by time-resolved laser pump-probe (PP) experiments using magneto-optical effects. The latter distinguish themselves in various ways from the former. The coherent excitation of several PSSWs with different frequencies, from the sub-GHz to the THz range, is made possible by the wide frequency spectrum of the pulsed excitation induced by the femtosecond pump laser pulses. In a single time-scan PP experiments can detect several coherent SW modes, and provide their time period, their respec-tive phases and their damping. Coherent control experi-54

70

  apparent different chiralities for SWs of different pari-71 ties. We show that they can only be explained by this 72 theoretical description of the magneto-optical detection 73 through the Kerr and Voigt effects that takes into ac-74 count the absorption depth and the optical phase shift 75 inside the layer. Important consequences are expected 76 when the SW wavelength (determined by the layer thick-77 ness) becomes a few tenth of the wavelength of light in 78 the material, λ/η, where η is the refractive index. In par-79 ticular, ignoring this effect can lead to an erroneous de-80 termination of the SW stiffness, and of the relative mode 81 amplitude, a signature of the up-to-now elusive magnon 82 excitation mechanisms. We show that the optical phase

83

  shift can lead to a striking and non-intuitive optical ef-84 fect in the detection of SWs: an apparent reversal of the 85 magnetization rotation direction for SWs of odd parity 86 with respect to the layer mid-plane. A key result of this 87 paper is that the optical phase shift provides a unique 88 tool for the determination of the SW mode number, or 89 in other words its parity. 90 We study thin layers of the ferromagnetic semiconduc-91 tor alloy GaMnAsP. Most samples show only one PSSW mode in the FMR spectra while one or two modes are optically detected in the TRMO signal [8, 29]. The results presented here are obtained in an in-plane magnetized GaMnAsP layer with thickness L=50 nm and Phosphorus concentration 4.3 % grown on a (001) GaAs substrate by molecular beam epitaxy and annealed 1 hour at 250 • C. The effective Mn concentration is 4 % and the Curie temperature is 85 K. The anisotropy constants were determined by FMR. PP experiments are carried out at T =12 K in zero external magnetic field after a 60 mT in-plane initialization of the magnetization direction. The laser source is a 76 MHz Ti:Sapphire laser at a wavelength λ=700 nm. To limit thermal effects, low pump and probe fluence are used (1 µJ cm -2 and 0.4 µJ cm -2 , respectively) [30]. The pump beam is modulated at 50 kHz. The pump-induced magnetization dynamics is detected as a function of the pump-probe delay through the rotation of the probe beam linear polarization detected by a balanced optical diode bridge and a lock-in amplifier. The static rotation and ellipticity signals are obtained with the probe beam only.

Figure 1 (

 1 Figure 1(b) shows the dependence of the TRMO signal on the incident probe polarization. The signal is fitted with u(t) + v(t) sin 2β + w(t) cos 2β from which we obtain the δθ exp (t) and δφ exp (t) functions taking into account the magnetization equilibrium angle φ 0 (M close to [100]) as shown in Fig. 1(c) (see Suppl-info). The plot of the trajectory, δθ exp (t) versus δφ exp (t) shown in Fig. 1(d) reveals a very complex dynamics that actually results from the contributions of two SWs that clearly appear in δθ exp (t) and δφ exp (t) and their Fourier transform (Fig. 1(c) and inset). Let us note the large ratio (≈ 0.6) of the amplitude of the second SW with respect to the first one.δθ exp (t) and δφ exp (t) are fitted with two damped oscillating signals and a sum of two exponentials that reflects the shape of the laser-induced pulsed excitation (τ 1 =0.03 ns, τ 2 =1 ns). Plotting separately the trajectory for each SW (f 0 =2.36 GHz and f 1 =3.90 GHz) in Fig.1(e, f) results in a much clearer picture of the magnetization precession. Surprisingly, the magnetization seems to rotate in opposite directions for the two SWs. As we shall see

FIG. 1 .

 1 FIG. 1. (a) Reference frame. (b) Dependence of the TRMO signal δβ on the probe beam linear polarization with respect to the sample crystallographic axes. (c) Time dependence of δθexp and δφexp. Inset: Fourier transform amplitude of δθexp. (d) Optically detected magnetization trajectory. (e), (f), (g) Decomposition of the experimental magnetization trajectory (d) into two oscillating signals at frequencies f0=2.36 GHz and f1=3.90 GHz and a non-oscillating signal, respectively.

  F i,j = ∂ 2 F ∂i∂j are 148 the second derivatives of F with respect to the angles 149 (i, j). Since it was shown that in thin (Ga,Mn)(As,P) lay-150 ers the lowest frequency PSSW is a nearly uniform mode, 151 independent of the layer thickness[8,9], the boundary 152 conditions at z=0 and z = L were chosen to ensure very 153 weak surface pinning, giving nearly flat δθ(z) and δφ(z) 154

  are the elements of the dielectric permittivity tensor describing the Kerr and Voigt effects, respectively (see Suppl-info). n = η + iκ is the layer mean complex refractive index and e i 4πnz λ = e -αz e i 4πηz λ with α = 4πκ λ the absorption coefficient. The important result is the modulation of the spatial dependence of the δθ(z, t) and δφ(z, t) magnetization components by the optical phase factor e i 4πηz λ that reflects the propagation of light from the surface to the depth z and back. The phase factor is damped by the e -αz absorption factor. Therefore, in the case of strong absorption as in metallic layers, the TRMO signal will be sensitive only to the SW amplitude very close the surface within the absorption depth. In the case of weak absorption and layer thickness L comparable to a fraction of the light wavelength inside the material λ/η, the optical phase shift plays a crucial role. It is precisely the case of the sample studied here where L ≈ λ/4η. 193 Actually, for static magnetization, the importance of a 194 phase shift factor that makes the magneto-optical effects 195 sensitive to the magnetization at a specific depth inside 196 single or multiple ferromagnetic layers was theoretically 197

206Figures 2

 2 Figures2(d-f) show the depth dependence of the 207

FIG. 2 .

 2 FIG. 2. (left)Magnetization trajectory in the depth of the ferromagnetic layer for the p=0 (a), p=1 (b), p=2 (c) PSSW modes and corresponding detected trajectory in (g), (h), (i) assuming that the optical signal would result from a depth-averaged amplitude. (right) (d), (e), (f) Theoretical effective magnetization trajectory in the depth of the ferromagnetic layer taking into account the optical phase factor and corresponding optically detected trajectory in (j), (k), (l). L=50nm, λ=700 nm, the other parameters are given in Suppl-info.

FIG. 3 .

 3 FIG. 3. Theoretical optically detected direction of rotation of the magnetization vector for p=0, 1, and 2 PSSWs modes from the amplitude ratio rp of δθ opt p and δφ opt p (see text). The dashed line indicates the layer thickness. λ=700 nm, η=3.67, κ=0.1.

Phys. D: Appl. Phys.36, R277 (2003).


Silvia Bertoluzza
email: silvia.bertoluzza@imati.cnr.it

Valérie Perrier
email: lerie.perrier@imag.fr

A new construction of boundary interpolating wavelets for fourth order problems

In this article we introduce a new mixed Lagrange-Hermite interpolating wavelet family on the interval, to deal with two types (Dirichlet and Neumann) of boundary conditions. As this construction is a slightly modification of the interpolating wavelets on the interval of Donoho, it leads to fast decomposition, error estimates and norm equivalences. This new basis is then used in adaptive wavelet collocation schemes for the solution of one dimensional fourth order problems. Numerical tests conducted on the 1D Euler-Bernoulli beam problem, show the efficiency of the method.

Introduction

The properties of wavelet bases are nowadays well understood, and they allow the design of methods of arbitrary high order for the solution of elliptic partial differential equations, that display an optimal lowest complexity coupled with the good local behaviour necessary for local refinement [START_REF] Cohen | Wavelet methods in numerical analysis[END_REF]. Wavelet collocation methods were introduced in the mid nineties [START_REF] Bertoluzza | A Wavelet Collocation Method for the Numerical Solution of Partial Differential Equations[END_REF][START_REF] Vasilyev | A Dynamically Adaptive Multilevel Wavelet Collocation Method for Solving Partial Differential Equations in a Finite Domain[END_REF] and have since proven their strength on second order problems arising in fluid-dynamics [START_REF] Bertoluzza | Adaptive Wavelet Collocation Method for the Solution of Burgers Equation[END_REF][START_REF] Vasilyev | An Adaptive Multilevel Wavelet Collocation Method for Elliptic Problems[END_REF][START_REF] Schneider | Wavelet Methods in Computational Fluid Dynamics[END_REF] and in structural mechanics. In particular, recent works show their effectiveness in the study of buckling [START_REF] Ferreira | A wavelet collocation approach for the analysis of laminated shells[END_REF] or vibration [START_REF] Ferreira | A high order collocation method for the static and vibration analysis of composite plates using a first order theory[END_REF] phenomena.

The wavelet collocation method is based on two (independent) ingredients: expanding the unknown solution of the equation considered in a Deslaurier-Dubuc interpolating basis [START_REF] Deslauriers | Symmetric Iterative Interpolation Processes[END_REF], and determining the unknown coefficients by requiring that the equation is exactly verified at the corresponding "collocation points". The construction of interpolating wavelets on the interval by Donoho [START_REF] Donoho | Interpolating Wavelet Transform[END_REF], allows to efficiently handle Dirichlet boundary conditions [START_REF] Bertoluzza | Multiresolution Analysis and Wavelets for the Numerical Solution of Partial Differential Equations[END_REF][START_REF] Bertoluzza | Adaptive Wavelet Collocation Method for the Solution of Burgers Equation[END_REF]. Adaptive versions of the method on non uniform grids have been also proposed, based on a quite simple adaptive strategy where the wavelet coefficients themselves are used to decide where to refine or coarsen the underlying grid.

With respect to a Galerkin approach, collocation has some advantage, in particular when non linear equations are considered [START_REF] Bertoluzza | Adaptive Wavelet Collocation for nonlinear BVPs[END_REF]. Thanks to the particular structure of the Interpolating Wavelet Transform, the evaluation of local nonlinearities turns out to be quite easy, also in an adaptive framework, and collocation allows to avoid the need of computing integrals involving non linear expressions, which, in the wavelet context, can be quite cumbersome.

In this article we aim at extending such a method to the numerical solution of fourth order problems. Such problems occur in a variety of applications including structural mechanics (beams, plates, ...), fluid dynamics (ice formation, fluids in the lungs), nanotechnologies, or image denoising. We should point out that, in this context, Galerkin methods using orthogonal wavelet bases have been already successfully employed in previous works [START_REF] Joly | Wavelet preconditioning of the Stokes Problem in ψ-ω formulation[END_REF][START_REF] Ma | A study of the construction and application of a Daubechies wavelet-based beam element[END_REF][START_REF] Diaz | Daubechies wavelet beam and plate finite elements, Daubechies wavelet beam and plate finite elements[END_REF]. The first problem that we need to face is the handling of the boundary conditions. Focusing, for the moment, on the one dimensional case, we will have to impose two boundary conditions at each of the two extrema of the interval. In order to maintain the balance between the number of equations and the number of degrees of freedom, we will then need to enrich the space with one more function per extremum. To this aim, we introduce a new family of interpolating wavelets on the interval, that we will call Mixed Lagrange-Hermite interpolatory wavelet bases. It consists in adding a boundary interpolating function at each extremum, satisfying an Hermite interpolation property (vanishing at all nodes and with first derivative equals to one at one node) and in modifying the other (Lagrange) interpolating boundary functions, in such a way that their derivative vanishes at the extrema. This new interpolating wavelet family will be used to develop our adaptive collocation schemes for the solution of one-dimensional fourth order problems.

The rest of the paper is organized as follows: Section 2 introduces the construction of the boundary mixed Lagrange-Hermite interpolating wavelets, and states the main properties satisfied by this new wavelet basis. Section 3 describes their use in the framework of an adaptive collocation wavelet schemes, and presents first numerical tests in 1D on two Euler-Bernulli beam problems.

Interpolating boundary Hermite wavelets on the Interval

To construct our new boundary modified interpolating wavelet basis on the interval, we start with an Interpolating Multiresolution Analysis (MRA) on the line [START_REF] Donoho | Interpolating Wavelet Transform[END_REF][START_REF] Mallat | A wavelet tour of signal processing[END_REF]. Let φ be a Deslaurier-Dubuc scaling function on R, which, we recall, satisfies an interpolation property: φ(0) = 1 and φ(k) = 0, ∀k = 0. Such function φ can be constructed by autocorrelation from a Daubechies orthonormal compactly scaling function Φ with N + 1 vanishing moments [START_REF] Daubechies | Ten Lectures on Wavelets, CBMS-NSF Regional Conference Series in Applied Mathematics[END_REF]:

φ(x) = Φ * Φ(-•)(x). (1)
We let L = 2N +1 denote the polynomial reproduction allowed by φ, i.e., for each polynomial p of degree lower or equal than L we have the identity p(x) = k p(k)φ(x -k). Moreover, we let R > 1 be the regularity of φ, such that φ is Hölder continuous of order R, which, for simplicity, we will assume to be integer. We observe that supp φ = [-L, L]. Figure 1 shows different functions φ obtained from different Daubechies scaling functions.

We will construct in the following a MRA on the interval (0, 1), with interpolating functions at interior points and allowing an Hermite interpolation at boundaries 0 and 1. This MRA will be a slight modification of interpolating MRA on the interval introduced by [START_REF] Donoho | Interpolating Wavelet Transform[END_REF], obtained by adding at each scale one function at each boundary, allowing the interpolation of the derivative at 0 and 1, while modifying the other boundary functions in such a way that their derivative vanishes at 0 and 1. We remark that this MRA does not fall in the framework of Hermite Multiwavelets constructed in [START_REF] Dahmen | Biorthogonal Multiwavelets on the Interval: Cubic Hermite Spines[END_REF].

We aim at constructing our MRA on (0, 1) in such a way that any continuously differentiable function can be expanded as a linear combination of interpolating functions as follows:

f = 2 j 0 +1 k=-1 β j 0 ,k ϕ j 0 ,k + j≥j 0 2 j k=0 α j,k ψ j,k (2)
where j 0 will be specified in the following, and where the functions ϕ j 0 ,k and ψ j,k satisfy a mixed Lagrange-Hermite interpolating condition:

ϕ j 0 ,k (k 2 -j ) = 2 j 0 /2 δ k,k , ∀k = 0, • • • , 2 j 0 and ϕ j 0 ,k (0) = ϕ j 0 ,k (1) = 0, φ j 0 ,-1 (k 2 -j ) = 0, ∀k = 0, • • • , 2 j 0 and (φ j 0 ,-1 ) (0) = 1 , (φ j 0 ,-1 ) (0) = 0, φ j 0 ,2 j +1 (k 2 -j ) = 0, ∀k = 0, • • • , 2 j 0 and (φ j 0 ,2 j 0 +1 ) (1) = 0 , (φ j 0 ,2 j 0 +1 ) (1) = 0, ψ j,k (k 2 -j-1 ) = 2 3j/2 δ 2k+1,k , ∀j ≥ j 0 and ∀k = 0, • • • , 2 j+1 .
The convergence of the series [START_REF] Bertoluzza | A Wavelet Collocation Method for the Numerical Solution of Partial Differential Equations[END_REF] in Sobolev spaces will be a consequence of Theorem 2.2.

Construction of the basis

The construction follows the same approach as in [START_REF] Donoho | Interpolating Wavelet Transform[END_REF]. Let j 0 = log 2 (L -1) , and let j ≥ j 0 . Let f ∈ C 0 ([0, 1]) be continuously differentiable at 0 and 1. We construct a sequence { βj,k ; k ∈ Z} in such a way that

k∈Z βj,k φ j,k (n/2 j ) = f (n/2 j ), n = 0, • • • , 2 j , k∈Z βj,k φ j,k (0) = f (0 + ), k∈Z βj,k φ j,k (1) = f (1 -),
where, throughout this section we set, as usual, φ j,k (x) = 2 j/2 φ(2 j x -k). To do so, we let {β j,k = 2 -j/2 f (k2 -j ) ; k = 0, . . . 2 j } be normalized samples of f in [0, 1], and

{β h j = 2 -3j/2 f (0 + ), β h j = 2 -3j/2 f (1 -
)} its normalized derivative values (resp. right derivative and left derivative) at boundaries 0 and 1. The extrapolated samples { βj,k ; k < 0} are defined by introducing the Lagrange-Hermite polynomial H j,L of degree L interpolating f at the L leftmost samples, as well as its derivative at 0:

H j,L (k2 -j ) = f (k2 -j ) for 0 ≤ k ≤ L -1 and (H j,L ) (0) = f (0 + )
and using this to extrapolate to k < 0:

βj,k = 2 -j/2 H j,L (k2 -j ) , k < 0. (3)
Similarly the samples { βj,k ; k > 2 j } are obtained by introducing the Lagrange-Hermite polynomial H j,L of degree L interpolating f at the L rightmost samples, as well as its derivative at 1:

H j,L (k2 -j ) = f (k2 -j ) for 2 j -L + 1 ≤ k ≤ 2 j and (H j,L ) (1) = f (1 -).
Using this to extrapolate to k > 2 j :

βj,k = 2 -j/2 H j,L (k2 -j ) , k > 2 j , (4)
the extrapolation operator E j,L of degree L at level j is then defined by:

E j,L (β j,k ) 2 j k=0 , β h j , β h j = ( βj,k ) j∈Z , (5)
where βj,k coincides with β j,k for k = 0, • • • , 2 j , and is defined by (3) for k < 0 and (4) for k > 2 j .

Definition 2.1. Let V j be the collection of the restrictions to [0, 1] of all functions f defined by: f

= k∈Z βj,k φ j,k (6)
where the coefficients βj,k are the extrapolated samples, defined by relation [START_REF] Bertoluzza | Adaptive Wavelet Collocation for nonlinear BVPs[END_REF], of the normalized node values of some function f .

The following properties hold.

Property 2.1. The space V j is a mixed Lagrange-Hermite interpolating space in the sense that for each function f ∈ V j , defined by [START_REF] Chiavassa | On the Effective Construction of Compactly Supported Wavelets Satisfying Homogeneous Boundary Conditions on the Interval[END_REF] with coefficients obtained by extrapolation of the normalized samples of a function f , it holds that f (k2 -j ) = 2 j/2 β j,k = f (k2 -j ) ∀k = 0, • • • 2 j , and f (0

) = f (0 + ), f (1) = f (1 -).
Proof. f (k2 -j ) = 2 j/2 βj,k = 2 j/2 β j,k for k = 0, • • • 2 j due to the interpolation property of scaling functions: φ j,k (n2 -j ) = 2 j/2 δ k,n . The property on the derivatives follows from the fact that f is a C 1 -function (φ is assumed to be C R with R > 1), and polynomial outside [0, 1], and has for right-derivative f (0 + ) at 0 and left-derivative f (1 -) at 1. By C 1 -continuity, one gets:

f (0) = f (0 + ), f (1) = f (1 -).
Together with the definition of V j , property 2.1 gives us a condition which is necessary and sufficient for a function f to belong to

V j : f ∈ C 1 ([0, 1]) belongs to V j if and only if f = k βj,k φ j,k | [0,1] ,
where βj,k are the extrapolated samples, defined by relation [START_REF] Bertoluzza | Adaptive Wavelet Collocation for nonlinear BVPs[END_REF], of the normalized node values of f itself. Property 2.2. V j reproduces the polynomials of degree ≤ L: if Π L denotes the space of polynomials of degree less than L on [0, 1], one has Π L ⊂ V j .

Proof. Property 2.2 follows from the unicity of interpolating polynomial, and the reproduction of polynomials allowed by the φ j,k on R [START_REF] Donoho | Interpolating Wavelet Transform[END_REF].

We now exhibit the basis functions of the space V j . We first define the left extrapolation weights: ∀n ∈ {1, • • • , L -1} and ∀k ∈ Z,

a k,n = L j,n (k2 -j ), a k,0 = H j,0 (k2 -j ), b k = 2 j H j,1 (k2 -j ),
where the Lagrange-Hermite polynomials L j,n , H j,0 , H j,1 are the interpolation polynomials needed to define H j,L , that are defined respectively as follows:

• for n = 1, • • • , L -1, L j,
n is the polynomial assuming the value 1 in n/2 j , vanishing at k/2 j for k = 0, • • • , L -1, k = n and whose derivative vanishes in 0,

• H j,0 is the polynomial vanishing at k/2 j , k = 1, • • • , L -1, assuming value 1 at 0 and whose derivative vanishes in 0,

• H j,1 is the polynomial vanishing at k/2 j , k = 0, • • • , L-1, and whose derivative assumes the value 1 in 0. The exact analytic expression for such polynomials are given in Annexe 4 (see equations (34,35,36)). Remark that the weights are indeed independent of j. With such weights we have

βj,k = β h j b k + L-1 n=0 a k,n β j,n ∀k < 0.
For j ≥ j 1 = log 2 (L -1) + 2, the left boundary scaling functions are then defined by:

ϕ j,n = φ j,n + -1 k=-L+1 a k,n φ j,k , 0 ≤ n ≤ L -1 (7)
ϕ h j = φ j,0 + -1 k=-L+1 b k φ j,k (8)
Remark that equations (7,8) also make sense for j = 0, so that they allow us to define functions φ 0,n , n = 0, • • • , L -1, and φ h 0 , such that using the the scaling invariance of the φ j,n , we obtain ϕ j,n (x) = 2 j/2 ϕ 0,n (2 j x) and ϕ h j,n (x) = 2 j/2 ϕ h 0,n (2 j x).

Similarly, the right extrapolation weights are defined by: ∀n ∈ {1, • • • , L-1} and ∀k ∈ Z,

a k,n = L j,n ((2 j -k)2 -j ), a k,0 = H j,0 ((2 j -k)2 -j ), b k = 2 j H j,1 ((2 j -k)2 -j )
where the Lagrange-Hermite polynomials L j,n , H j,0 , H j,1 are defined analogously to L j,n , H j,0 and L j,1 and are used to build H j,L (see once again Annexe 4). Then

βj,2 j -k = β h j b k + L-1 n=0 a k,n β j,2 j -n ∀k < 0
The right boundary scaling functions are then defined by:

ϕ j,2 j -n = φ j,2 j -n + -1 k=-L+1 a k,n φ j,2 j -k , 0 ≤ n ≤ L -1 (9) ϕ h j = φ j,2 j + -1 k=-L+1 b k φ j,2 j -k (10)
These boundary scaling functions are completed by the "interior" scaling functions

ϕ j,k = 2 j/2 φ(2 j x -k), k = L, • • • , 2 j -L
to span the space V j , and the following proposition holds.

Proposition 2.1. For j ≥ j 1 , V j is a finite dimensional space of dimension 2 j + 3 and it satisfies:

V j = span{ϕ j,k } L-1 k=0 ⊕ span{ϕ h j } ⊕ span{ϕ j,k } 2 j -L k=L ⊕ span{ϕ j,k } 2 j k=2 j -L+1 ⊕ span{ϕ h j }
Moreover, every f ∈ V j has the representation:

f = L-1 k=0 β j,k ϕ j,k + β h j ϕ h j + 2 j -L k=L β j,k ϕ j,k + 2 j k=2 j -L+1 β j,k ϕ j,k + β h ϕ h j with the coefficients: β j,k = 2 -j/2 f (k2 -j ), β h j = 2 -3j/2 f (0 + ), β h j = 2 -3j/2 f (1 -).
Remark 2.1. For j ≥ j 1 we have that 2 -j (L -1) < 1 -2 -j (L -1) and then the samples used in the left hand side extrapolation (k < 0) are disjoint from the samples used in the right hand side extrapolation (k > 2 j ). For j ∈ [j 0 , j 1 -1] the left and the right hand side extrapolations share some of the nodes. As a result we obtain some basis functions which have contributions both from the left and the right side of the interval. More precisely, for 2 j -L + 1 ≤ n ≤ L -1 we can introduce the function

ϕ j,n = φ j,n + -1 k=-L+1 a k,n φ j,k + -1 k=-L+1 a k,2 j -n φ j,2 j -k .
We then have the representation, for j 0 ≤ j < j 1

V j = span{ϕ j,k } 2 j -L k=0 ⊕ span{ϕ h j } ⊕ span{ϕ j,k } L-1 k=2 j -L+1 ⊕ span{ϕ j,k } 2 j k=L ⊕ span{ϕ h j }.
We have the following Lemma.

Lemma 2.1. The sequence {V j } j≥j 0 is nested:

V j ⊂ V j+1 (11)
Proof. We start by recalling a fundamental property of the Deslaurier-Dubuc MRA [START_REF] Donoho | Interpolating Wavelet Transform[END_REF].

Let U j = span{φ j,k ; k ∈ Z} ∈ L 2 loc (R). Let f = k∈Z 2 -j/2 f (k2 -j )φ jk ∈ U j .
Then the values of f at dyadic points at level j + 1 can be obtained from the values at the dyadic points at level j by centered polynomial interpolation of order L = 2N + 1. More precisely we have that

f ((2k + 1)2 -(j+1) ) = π j k ((2k + 1)2 -(j+1) ) (12)
with π j k ∈ Π L defined by

π j k (n2 -j ) = f (n2 -j ), n = k -N, • • • , k + N + 1. ( 13
)
We next recall that V j can be characterized as the restriction to [0, 1] of the subspace of U j (which by abuse of notation we will also denote by V j ) whose coefficients are obtained through the extrapolation procedure E j,L defined in [START_REF] Bertoluzza | Adaptive Wavelet Collocation for nonlinear BVPs[END_REF]:

f ∈ V j if and only if f = k β jk φ jk with (2 j/2 β jk ) k = E j,L (f (0), f (2 -j ), • • • , f (1 -2 -j ), f (1), 2 -j f (0), 2 -j f (1)). ( 14
)
Let then f ∈ V j be the restriction to [0, 1] of f = k β jk φ jk ∈ U j satisfying [START_REF] Daubechies | Ten Lectures on Wavelets, CBMS-NSF Regional Conference Series in Applied Mathematics[END_REF]. Since

U j ⊆ U j+1 we have that f ∈ U j+1 and f = k β j+1,k φ j+1,k with coefficients β j+1,k = 2 -(j+1)/2 f (k2 -(j+1)
), given by [START_REF] Dahmen | Biorthogonal Multiwavelets on the Interval: Cubic Hermite Spines[END_REF][START_REF] Dahmen | Biorthogonal Spline-wavelets on the interval. Stability and moment conditions[END_REF]. In order to prove that f ∈ V j+1 it is then sufficient to prove that (β j+1,k

) k = ( βj+1,k ) k with (2 j+1 2 βj+1,k ) k := E j+1,L (f (0), f (2 -(j+1) ), • • • , f (1 -2 -(j+1) ), f (1), 2 -(j+1) f (0), 2 -(j+1) f (1)).
By the definition of E j+1,L , for k = 0, • • • , 2 j+1 we have that βj+1,k = 2 -(j+1)/2 f (k/2 j+1 ) = β j+1,k . We then need to prove the identity for k < 0 and k > 2 j+1 . To fix the ideas let us consider the k < 0 case. Let p = H j,L denote the polynomial of degree

L = 2N + 1 that verifies p(k2 -j ) = f (k2 -j ), k = 0, • • • , L -1, p (0) = f (0), (used in the extrapolation E j,L ). Since f ∈ V j we have that f (k/2 j ) = p(k/2 j ) for k ≤ L -1.
Using [START_REF] Dahmen | Biorthogonal Multiwavelets on the Interval: Cubic Hermite Spines[END_REF] it is then immediate to verify that f (k/2 j+1 ) = p(k/2 j+1 ), for k ≤ L -1: for even values of k, this is trivial, while for odd values of k = 2n + 1 we have

f (k/2 j+1 ) = π j k (k/2 j+1 ) = p(k/2 j+1 ),
since π j k coincides with f (and therefore with p) at the

L + 1 nodes (k -N )/2 j , • • • , (k + N + 1)/2 j . Letting now q = H j+1,L denote the polynomial of degree L that verifies q(k2 -(j+1) ) = f (k2 -(j+1) ), k = 0, • • • , L -1, q (0) = f (0), it is immediate to check that q = p (since q(k/2 j+1 ) = p(k/2 j+1 ), k = 0, • • • , L -1 and q (0) = p (0)). It is now easy to conclude that for k < 0 βj+1,k = 2 -(j+1) q(2 -(j+1) k) = 2 -(j+1) p(2 -(j+1) k) = 2 -(j+1) f (2 -(j+1) k) = β j+1,k .
An analogous argument yields βj+1,k = β j+1,k for k > 2 j+1 .

Figure 3: The basis functions for the spaces V j 0 (top-left) and V j 1 (top-right), and the multiscale basis for V j 1 obtained from the splitting of the form (20) (bottom).

Interpolation and wavelets

From now on, in order to simplify the notation, we will denote the boundary scaling functions of V j by ϕ j,k := ϕ j,k for k = 0, . . . , L -1 and ϕ j,k := ϕ j,k for k = 2 j -L + 1, . . . , 2 j . The "Hermite" scaling functions ϕ h j and ϕ h j will be respectively denoted by ϕ j,-1 and ϕ j,2 j +1 . This leads to the uniform representation:

V j = span{ϕ j,k ; k = -1, . . . , 2 j + 1}
where, of course, the scaling functions ϕ j,k for k < L and k > 2 j -L do not coincide those on R:

ϕ j,k (x) = φ j,k (x) = 2 j/2 φ(2 j x -k).
We may then define the interpolation operator

I j : C 1 ([0, 1]) → V j I j f = 2 -3j/2 f (0) ϕ j,-1 + 2 -j/2 2 j k=0 f (k2 -j ) ϕ jk + 2 -3j/2 f (1) ϕ j,2 j +1 . ( 15
)
Thanks to the polynomial reproduction property the following error estimate holds.

Theorem 2.2. If f ∈ H s (0, 1) with 3/2 < s < L + 1, 0 ≤ r ≤ R -1, r ≤ s then f -I j f H r (0,1) 2 -j(s-r) |f | H s (0,1) . (16)
Proof. We prove the result for j ≥ j 1 (the proof for j < j 1 being analogous). Let r be an integer. We have

|f -I j f 2 H r (0,1) = 2 j -1 k=0 |f -I j f | 2 H r (K j,k ) , with K j,k = (k/2 j , (k + 1)/2 j ).
Now, for p ∈ Π L arbitrary we can write:

|f -I j f | H r (K j,k ) = |(f + p) -I j (f + p)| H r (K j,k ) ≤ |f + p| H r (K j,k ) + |I j (f + p)| H r (K j,k ) .
We distinguish three cases:

k ≤ L -1, L ≤ k ≤ 2 j -L, and k ≥ 2 j -L + 1.
Let us consider the first case. For k ≤ L -1 the modified right boundary functions ϕ h j and ϕ j,2 j -n , 0 ≤ n ≤ L -1, vanish on K j,k and then, for x ∈ K j,k , we have

I j g(x) = 2 j+1 n=-1 β j,n ϕ j,n (x) = k+L-1 n=-1 β j,n ϕ j,n (x) = β j,-1 2 j/2 ϕ h 0 (2 j x) + L-1 n=0 β j,n 2 j/2 ϕ 0,n (2 j x) + k+L-1 n=L β j,n 2 j/2 φ(2 j x -n), (17)
with β j,-1 = 2 -3j/2 g (0), β j,n = 2 -j/2 g(n/2 j ), n ≥ 0.

For g = f + p we let ĝ be such that

g(x) = 2 j/2 ĝ(2 j x)
and we set ĥ

(x) = β j,-1 ϕ h 0 (x) + L-1 n=0 β j,n ϕ 0,n (x) + k+L-1 n=L β j,n φ(x -n),
so that I j g(x) = 2 j/2 ĥ(2 j x). We have

|g -I j g| H r (K j,k ) = 2 jr |ĝ -ĥ| H r (K 0,k ) ≤ 2 jr (|ĝ| H r (K 0,k ) + | ĥ| H r (K 0,k ) ).
Let us consider the contribution of ĥ. We easily see that

| ĥ| H r (K 0,k ) max{|β j,n |, n = -1, • • • , k + L -1} ĝ W 1,∞ ([0,2L-2]) ĝ H s ([0,2L-2]) , yielding |g -I j g| H r (K j,k ) 2 jr ĝ H s ([0,2L-2]) .
Letting f and p be such that f (x) = 2 j/2 f (2x) and p(x) = 2 j/2 p(2 j x), so that ĝ = f + p, and given the arbitrariness of p we immediately obtain

|f -I j f | H r (K j,k ) 2 jr inf p∈Π L f + p H s ((0,2L-2)) 2 jr | f | H s (0,2L-2) ,
where the last inequality is obtained by applying Theorem 3.1.1 in [START_REF] Ciarlet | The finite element method for elliptic problem[END_REF]. A change of variable finally yields

|f -I j f | H r (K j,k ) 2 j(r-s) |f | H s (0,2 -j (2L-2)) .
Analogously, for L ≤ k ≤ 2 j -L we prove

|f -I j f | H r (K j,k ) 2 j(r-s) |f | H s (2 -j (k-L+1),2 -j (k+L-1)) ,
and for k ≥ 2 j -L + 1 we prove

|f -I j f | H r (K j,k ) 2 j(r-s) |f | H s (2 -j (2 j -2L+2),1)
.

Squaring and adding up we obtain

|f -I j f | H r (]0,1[) 2 j(r-s) |f | H s (]0,1[) .
which easily yields the thesis.

In addition to the direct type inequality given by Theorem 2.2, the space V j inherits from U j an inverse type inequality.

Proposition 2.3. For all 0 ≤ r ≤ s < R, for all f ∈ V j it holds that

f H s (]0,1[) 2 j(r-s) f H r (]0,1[) .
Let us now introduce the wavelet space W j ⊂ V j+1 defined as

W j = (I j+1 -I j ) V j = span{ϕ j+1,2k+1 ; k = 0, . . . , 2 j -1}. ( 18
)
For analogy with orthonormal wavelets we introduce the notation

ψ jk = ϕ j+1,2k+1 . (19)
We may remark that the function ψ jk verifies

ψ jk ((2n + 1)/2 j+1 ) = 2 j+1 2 δ kn , and ψ jk (n/2 j ) = 0, ∀j ≤ j, n = 0, • • • , 2 j
It is easy to check that

V j+1 = V j ⊕ W j . ( 20
)
Such kind of multiresolution analysis has been extensively investigated in the Lagrangeinterpolating case [START_REF] Donoho | Interpolating Wavelet Transform[END_REF]. In the present case, it is possible to define a Lagrange-Hermite interpolating wavelet transform mapping any continuously differentiable function into the sequence of its coefficients

f → {β j 0 ,k , k = 0, • • • , 2 j 0 }, {α j 0 ,k , k = 1, • • • , 2 j 0 }, {α j 0 +1,k , k = 1, • • • , 2 j 0 +1 }, . . . .
The following proposition, which can be proven analogously to the corresponding proposition in [START_REF] Donoho | Interpolating Wavelet Transform[END_REF], holds.

β j 0 ,k ϕ j 0 ,k + j≥j 0 2 j k=0 α j,k ψ j,k = I j 0 f + j≥j 0 Q j f,
the operator I j f being defined in [START_REF] Diaz | Daubechies wavelet beam and plate finite elements, Daubechies wavelet beam and plate finite elements[END_REF], and with Q j = I j+1 -I j , with uniform convergence of partial sums (with respect to j). The scaling wavelet coefficients β j 0 ,k correspond to:

β j 0 ,k = 2 -j 0 /2 f (k2 -j 0 ) , ∀k = 0, . . . 2 j β j 0 ,-1 = 2 -3j 0 /2 f (0) and β j,2 j 0 +1 = 2 -3j 0 /2 f (1)
whereas the wavelet coefficients α j,k are given by:

α j,k = 2 -j/2 f ((k + 1/2) 2 -j ) -I j f ((k + 1/2) 2 -j ) .

Norme equivalences

One of the relevant properties of wavelet bases is the possibility of expressing equivalent norms for a wide range of functional spaces in terms of wavelet coefficients. This also holds true for the new wavelet basis we just introduced, provided, of course, the functional space considered is continuously embedded in C 1 . In view of the application to the solution of PDEs, we concentrate here on Sobolev spaces. The following theorem holds.

Theorem 2.4. Let f ∈ C 1 ([0, 1]). Then f ∈ H s (]0, 1[), 3/2 < s < R if and only if f 2 s := 2 j 0 +1 k=-1 |β j 0 k | 2 + j≥0 2 j k=1 2 2js |α jk | 2 < +∞. ( 21
)
Moreover • s is an equivalent norm for H s (]0, 1[).

Proof. For j ≥ j 0 let Q j = I j+1 -I j . We start by proving that if f s < +∞ then f ∈ H s (]0, 1[). In fact we have

f 2 H s (]0,1[) = I j 0 f + j≥j 0 Q j f 2 H s (]0,1[) I j 0 f 2 H s (]0,1[) + j≥j 0 Q j f 2 H s (]0,1[) Now we have j≥j 0 Q j f 2 H s (]0,1[) = j≥j 0 Q j f 2 H s (]0,1[) + 2 j≥j 0 k>j Q j f, Q k f H s (]0,1[) . (22)
We recall that the following bound holds for all f ∈ H s (]0, 1[), g ∈ H s+2ε (]0, 1[), provided ε > 0 is sufficiently small:

f, g H s (]0,1[) f H s-2ε (]0,1[) g H s+2ε (]0,1[) . (23)
Then we can write

j k>j Q j f, Q k f H s (]0,1[) j k>j Q j f H s+2ε (]0,1[) Q k f H s-2ε (]0,1[) j 2 js Q j f L 2 (]0,1[) k>j 2 -2ε|j-k| 2 ks Q k f L 2 (]0,1[) j 2 2js Q j f 2 L 2 (]0,1[) j |b j | 2 (24) with b j = k>j 2 -2ε|j-k| 2 ks Q k f L 2 (]0,1[) .
Using the 2 boundedness of the convolution product with 1 convolution kernel we have

(b j ) j 2 2 j≥j 0 2 2js Q j f 2 L 2 (]0,1[) . (25)
Using ( 22) and ( 24) we obtain

j≥j 0 Q j f 2 H s (]0,1[) j≥j 0 2 2js Q j f 2 L 2 (]0,1[) .
In order to prove the converse inequality we consider a decomposition f = fm 0 + m≥m 0 dm of the function f ∈ H s (]0, 1[) in a Daubechies orthonormal wavelet basis with a sufficiently high number of zero moments, in such a way that for all r, 0 ≤ r ≤ L + 1 the following norm equivalence holds:

f 2 H r (]0,1[) fm 0 2 L 2 (]0,1[) + m≥m 0 2 2mr dm 2 L 2 (]0,1[) .
We can then bound (1

-I j )f L 2 (]0,1[) as follows (1 -I j )f L 2 (]0,1[) ≤ (1 -I j ) fm 0 L 2 (]0,1[) + m≥m 0 (1 -I j ) dm L 2 (]0,1[) .
Let us bound the different terms on the right hand side. Let σ < L + 1 -s. Using Theorem 2.2, we have

(1 -I j ) fm 0 L 2 (]0,1[) 2 -j(s+σ) fm 0 H s+σ (]0,1[) 2 -j(s+σ) fm 0 L 2 (]0,1[) .
Let us now bound the term (1 -I j ) dm L 2 (]0,1[) . Let us at first consider the case m ≤ j. We have

(1 -I j ) dm L 2 (]0,1[) 2 -j(s+σ) dm H s+σ (]0,1[) 2 (m-j)s 2 -σ|m-j| dm L 2 (]0,1[) .
For m ≥ j we have

(1 -I j ) dm L 2 (]0,1[) 2 -j(s-σ) dm H s-σ (]0,1[) 2 (m-j)s 2 -σ|m-j| dm L 2 (]0,1[) . Since Q j f = (1 -I j )f -(1 -I j+1 )f , it is not difficult to verify that 2 sj Q j f L 2 (]0,1[) 2 -jσ fm 0 L 2 (]0,1[) + m 2 -σ|m-j| 2 ms dm L 2 (]0,1[)
Using once again the boundedness of the convolution operator we can show that

j 2 2js Q j f 2 L 2 (]0,1[) fm 0 2 L 2 (]0,1[) + m 2 2ms dm 2 L 2 (]0,1[) f 2 H s (]0,1[) .
With a similar argument we obtain

I j 0 f 2 L 2 (]0,1[) f 2 L 2 (]0,1[) + (1 -I j 0 )f 2 L 2 (]0,1[) f 2 H s (]0,1[) .
Finally we have proved that:

f 2 H s (]0,1[) I j 0 f 2 L 2 (]0,1[) + j≥j 0 Q j f 2 L 2 (]0,1[) .
We conclude by observing that

I j 0 f 2 L 2 (]0,1[) 2 j 0 +1 k=-1 |β j 0 k | 2
(since the space V j 0 is finite dimensional and all norms are equivalent), and that

Q j f 2 L 2 (]0,1[) 2 j k=0 |α jk | 2 ,
which can be obtained by using Lemma 7.6 and Lemma 7.7 of [START_REF] Donoho | Interpolating Wavelet Transform[END_REF].

Numerical solution of Fourth order problems

We want now to demonstrate the potential of the new basis for the solution of first order problems. In order to do so, in the following, we will show on classical test-cases, how our construction of mixed Lagrange-Hermite basis functions used within a collocation scheme allows to define an accurate numerical scheme, provided a multiple precision is used in selected computations.

Collocation wavelet method for fourth order problem

Let us consider a one-dimensional fourth-order boundary value problem of the form:

Find u : [0, 1] → R such that Au = f in ]0, 1[, B u(0) = g , B u(1) = g ( 26
)
where A is a fourth order operator, and

B = (B 1 , B 2 ) : C 1 ([0, 1]) → R 2 , B = (B 1 , B 2 ) : C 1 ([0, 1]) → R 2
denotes the boundary operators on the left and right boundary respectively.

In this section we adapt to the case of fourth order problems the classical wavelet collocation method introduced in [START_REF] Bertoluzza | A Wavelet Collocation Method for the Numerical Solution of Partial Differential Equations[END_REF] for second order problems, by considering the modified interpolated wavelets introduced in section 2, and the corresponding multiresolution analysis spaces (V J ). For a fixed J > 0, we introduce the grid G J of dyadic points at resolution 2 J in [0, 1]:

G J = {k/2 J : k = 0, • • • , 2 J }
We observe that each interior point k/2 J ∈ G J , k = 1, • • • , 2 J -1 corresponds to the basis function ϕ Jk , while each of the two extrema 0 = 0/2 J and 1 = 2 J /2 J corresponds to two basis functions, namely for x = 0 we have the two basis functions ϕ J,-1 and ϕ J,0 while for x = 1 we have the two functions ϕ J,2 J and ϕ J,2 J +1 . The collocation formulation of Problem [START_REF] Monasse | Orthogonal Wavelet Bases Adapted For Partial Differential Equations With Boundary Conditions[END_REF] writes: find u J ∈ V J such that

Au J (k2 -J ) = f (k2 -J ) k = 1, • • • , 2 J -1, (27)
B u J (0) = g , B u J (1) = g .
The unknown u J = 2 J +1 k=-1 u Jk ϕ Jk is sought by solving the linear system A J ūJ = fJ [START_REF]Multiprecision Computing Toolbox for MATLAB, vX.X.X.XXXX[END_REF] where the matrix A J is defined as follows

A n+2,k+2 = Aϕ Jk (n2 -J ), n = 1, • • • , 2 J -1, k = -1, • • • , 2 J + 1 A i,k+2 = B i ϕ Jk (0), A 2 J +i+1,k+2 = B i ϕ Jk (1), i = 1, 2, k = -1, • • • , 2 J + 1,
and where the vectors ūJ and fJ are, respectively, the vector of the coefficients of the unknown function u J and the vector of the values of the function f at the interior dyadic points and of the functions g at 0 and g at 1.

Collocation on non-uniform grids

It has been shown, in the case of second order problems, that the wavelet collocation method is well suited to be used in the framework of adaptivity. By using the wavelet coefficients as a guide to decide whether to refine or coarsen the grid, it is possible to design non uniform grids and corresponding non uniform wavelet spaces, well suited to the approximation of a given problem. This holds true also in the case of fourth order problem. In order to define the solution of a fourth order problem on a given non uniform grid (to be selected by a suitable adaptive refinement procedure), it will be convenient to introduce the following notation for dyadic points corresponding to the scaling functions ϕ jk and ψ jk

x jk = k2 -j , k = 0, • • • , 2 j , y jk = (2k -1)2 -(j+1) , k = 1, • • • , 2 j .
We observe that the grid G J can be obtained as

G J = G j 0 ∪ J-1 j=j 0 {y jk , k = 1, • • • , 2 j } = G j 0 ∪ {y λ , λ ∈ Λ}, with Λ = {(j, k), j ≥ j 0 , k = 1, • • • , 2 j }.
If we let Λ h ⊂ Λ be any finite set of Λ we can consider the grid

G h = G j 0 ∪ {y λ , λ ∈ Λ h }
and the corresponding finite dimensional space

V h = V j 0 ⊕ span {ψ λ , λ ∈ Λ h } .
We then consider the problem: find

u h ∈ V h such that Au h (λ) = f (λ) ∀λ ∈ G h \ {0, 1}, B u h (0) = g , B u h (1) = g . (29)

A simple adaptive strategy

The simplest way to construct a non-uniform approximation space well suited to the solution of a given problem is to use the simple adaptive strategy already in use for second order problem [START_REF] Bertoluzza | Multiresolution Analysis and Wavelets for the Numerical Solution of Partial Differential Equations[END_REF][START_REF] Bertoluzza | Adaptive Wavelet Collocation Method for the Solution of Burgers Equation[END_REF], which is essentially the one initially proposed in [START_REF] Maday | Adaptativité dynamique sur bases d'ondelettes pour l'approximation d'équations aux dérivées partielles[END_REF]. The idea, heuristically justified by the norm equivalence [START_REF] Grivet-Talocia | Wavelets on the interval with optimal localization[END_REF], is to use the absolute value of the coefficients themselves as indication of the possible lack of regularity (and therefore need of employing a finer grid). A very small coefficient tells us that the corresponding point is not relevant, that is the solution is already approximated at a coarser grid. The point can then be removed without deteriorating the quality of the approximation. A big coefficient, on the other hand, indicates that either the function or one of its derivatives has a sharp variation in a neighborhood of the corresponding point, and that a locally finer grid might be necessary to achieve a good approximation. The adaptive strategy is then defined as follows. We choose two (small) tolerance parameters ε R and ε C with 0 < ε C < ε R , and we let G h and V h and u h = u 0 h + λ∈Λ h u λ ψ λ denote the actual grid and corresponding space and the actual solution (satisfying ( 29)). The next grid is defined according to the following strategy. For each λ ∈ Λ h three things can occur:

1. |u λ | ≤ ε C ; 2. ε C < |u λ | < ε R ; 3. |u λ | ≥ ε R .
In the case 1. we remove the point y λ from the grid. In the case 2. we leave y λ in the grid as it is. In the third case we leave y λ in the new grid, and we add two neighbouring points at the next level. More precisely, if for λ = (j, k) we set N λ = {(j + 1, 2k -1), (j + 1, 2k + 1)}, the new grid G * h is defined as

G * h = G j 0 ∪ {λ : |u λ | ≥ ε C } ∪ λ:|u λ |≥ε R N λ .
Using this refinement strategy we can think of a basic adaptive method where at each iteration the solution is computed and the grid is updated according to the strategy we just described. Though this is not an optimal method in terms of computational complexity, it gives us a good idea of the performance that we can expect from a more refined method in this context.

Implementation of the method

The implementation of the method is carried out in two stages. In the first stage, the values of the function φ as well as those of its derivatives up to order 4 are computed according at integers to the cascade algorithm described in [START_REF] Daubechies | Ten Lectures on Wavelets, CBMS-NSF Regional Conference Series in Applied Mathematics[END_REF]. This implies solving an eigenvector problem for a matrix of order 2L -1, L being the polynomial reproduction allowed by φ. By applying a refinement scheme the value of φ and of its derivatives at dyadic points in the support of φ are computed. Next we use extrapolation to compute the values of the modified functions ϕ 0,n and ϕ h 0 at positive dyadic points, as well as the values of the modified functions ϕ 0,n and ϕ h 0 at negative dyadic point of a grid at meshsize 2 -Jmax . This stage can be performed once and for all and the values stocked.

In the second stage the actual stiffness matrix is assembled by using the stocked values. In this procedure one has to take care of the precision with which the different step are computed. The high order extrapolation procedure results in an amplification of the roundoff error, which is further amplified by the multiplication by the factor 2 4j deriving from computing the fourth derivative of the rescaled functions. This can have devastating effects on the results, as can be seen by considering the trivial problem u = 0, u(0) = u(1) = 1, u (0) = u (1) = 0 whose solution is the constant function u = 1. In order to demonstrate the issue we solved this equation on different uniform grids by performing both stages of the procedure in double precision. The results are reported in table 1.

A better (but still unsatisfactory, considering that the basis used exactly reproduces polynomials of order 11) result is obtained when using a wavelet basis (see Table 2).

The solution that we propose to this problem is to carry out the first stage of the procedure (which, we recall, is carried out once and for all), and in particular the solution of the initial small eigenvector problem and the high order extrapolation, in multiple precision. The solution of the linear system is instead carried out in double precision. In order to do it we used the Matlab Multiprecision Toolbox [START_REF]Multiprecision Computing Toolbox for MATLAB, vX.X.X.XXXX[END_REF]. The results are displayed in table [START_REF] Bertoluzza | Adaptive Wavelet Collocation Method for the Solution of Burgers Equation[END_REF]. In table 3 we present the result obtained by performing the whole computation in multiple precision, with 30 exact decimal digits.

Numerical tests

The Euler-Bernoulli equation describes the relationship between the beam's deflection ω and the applied load q. The equation reads where E is the elastic modulus and I is the second moment of area of the beam cross section. The curve ω(x) describes the deflection of the beam in the z direction at some position x.

d 2 dx 2 E(x)I(x) d 2 dx 2 ω(x) = q(x), (30)

Test 1.

The first test we carried out is proposed in [START_REF] Twizell | A sixth order method for two beam problems[END_REF]. The equation we consider is the following:

w (iv) + 4w = 1, in (-1, 1) (31)
w(-1) = w(1) = 0, w (-1) = -w (1) = sinh(2) -sin(2) 4(cosh(2) + cos [START_REF] Bertoluzza | A Wavelet Collocation Method for the Numerical Solution of Partial Differential Equations[END_REF])

.

The exact solution takes the following form w = 1 -2(sin(1) sinh( 1) sin(x) sinh(x) + cos(1) cosh( 1) cos(x) cosh(x) 4(cosh(2) + cos( 2)

)
The problem is solved on a uniform grid with different mesh sizes. In table 5 we report the error, both in the L 2 and in the L ∞ norm. In the first two columns, we report the error obtained by using multiple precision also for the solution of the linear system, while the error reported on the two columns on the right is obtained by solving the linear system in double precision. The degradation of the error as the mesh size decreases is a consequence of the ill conditioning of the matrix. This issue will of course need to be faced. However we observe on the one hand that, when solving in multiple precision, our method outperforms the sixth order method of [START_REF] Twizell | A sixth order method for two beam problems[END_REF], and, on the other hand, that, also when solving the linear system in double precision the results on relatively coarse grids (h = .0625 and h = .0313), on which the roundoff error appears to be of the same order of magnitude as the approximation error, are quite good.

Test 2. Clamped-Clamped beam, half-loaded

The second example, taken from [START_REF] Diaz | Daubechies wavelet beam and plate finite elements, Daubechies wavelet beam and plate finite elements[END_REF][START_REF] Ma | A study of the construction and application of a Daubechies wavelet-based beam element[END_REF], is a non-uniformly loaded clamped beam of length 2L with equal cross section (see Figure 4), and with constant bending rigidity E 0 I 0 . The load, limited to the right half of the beam, takes the form

q(x) = 0 if 0 ≤ x ≤ L k 0 q 0 ( x L -1) if L ≤ x ≤ 2L
Then the exact solution of In the experiment, we took k 0 = 480, and q 0 E 0 I 0 = 1 and L = 0.5. The solution for such values of the parameters is shown on the right hand side of Figure 4. We apply to this example the adaptive strategy described in section 3.2. In Table 6 we display, on the right hand side, the error in the L 2 norm for the result obtained with a uniform discretization for different values of the meshsize h = 2 -J . On the right hand side we display the error (once again in the L 2 norm) between the true solution and the solution obtained by the adaptive strategy for different values of the tolerance ε R and for ε C = 10 -2 • ε R . Since the solution ω is not C ∞ but displays a jump in the 5-th derivative, the errors for the uniform discretization are definitely bigger than the errors displayed in the previous tests (whose solution is, instead, C ∞ ). The adaptive scheme manages to capture this lack of regularity and refines the grid at the mid-point of the segment, as shown in Figure 5. In figure 6 we show, for both uniform discretization and adaptive method, the error as a function of the number of degrees of freedom.

E 0 I 0 d 4 dx 4 ω(x) = q(x), ω(0) = ω (0) = ω(2L) = ω (2L) = 0 (33) is given by ω(x) =    q 0 k 0 E 0 I 0 -L 120 x 3 + 7L 2 480 x 2 0 ≤ x ≤ L q 0 k 0 E 0 I 0 1 120L x 5 -1 24 x 4 + 3L 40 x 3 -11L 2 160 x 2 + L 3 24 x -L 4 120 L ≤ x ≤ 2L 20

Annexe: on Lagrange-Hermite interpolation

The Lagrange-Hermite polynomial H L of degree L, interpolating f at samples (x 0 , x 1 , • • • , x L-1 ) and f at x 0 :

H L (x i ) = f (x i ) ∀i = 0, L -1, H L (x 0 ) = f (x 0 ), is given by H L (x) = f (x 0 )H 0 (x) + f (x 0 )H 1 (x) + L-1 i=1 f (x i )L i (x),
where: l i is the Lagrange interpolating polynomial of degree L-1 associated to samples (x 0 , x 1 , • • • , x L-1 ):

L i (x) = x -x 0 x i -x 0 l i (x) ∀i ≥ 1, (34)
H 0 (x) = (1 -2(x -x 0 )l 0 (x 0 )) l 0 (x), (35) H 1 (x) = (x -x 0 ) l 0 (x). ( 36
l i (x) = L-1 q=0,q =i
x -x q x i -x q , l 0 (x 0 ) = L-1 q=1 1 x 0 -x q .

Example 4.1. At a given level j > j 1 , for nodes of the form i2 -j , the interpolating polynomials L i , H 0 and H 1 will be denoted by L j,i , H j,0 and H j,1 . We separate the case x 0 = 0 (left), from the case x 0 = 1 (right) with different notations for the interpolating polynomials.

-Left extrapolation weights: at level j, the nodes are given by: x i = i2 -j , ∀i = 0, • • • , L -1. In this case:

l i (k2 -j ) = L-1 q=0,q =i k -q i -q , l 0 (0) = - L-1 q=1 1 q2 -j
Then for n ∈ {1, • • • , L -1} and ∀k ∈ Z,

a k,n = L j,n (k2 -j ) = k n 2 L-1
q=1,q =i k -q n -q a k,0 = H j,0 (k2 -j ) = 1 + 2k

L-1 q=1 1 q L-1 q=0,q =i k -q -q b k = 2 j H j,1 (k2 -j ) = k L-1 q=0,q =i k -q -q -Right extrapolation weights: at level j, the nodes are now:

x i = (2 j -i)2 -j , ∀i = 0, • • • , L -1. and ∀k ∈ Z,

l i ((2 j -k)2 -j ) = L-1 q=0,q =i k -q i -q , l 0 (1) = L-1 q=1 1 q2 -j Then ∀n ∈ {1, • • • , L -1}, a k,n = L j,n ((2 j -k)2 -j ) = k n 2 L-1
q=1,q =i k -q n -q a k,0 = H j,0 ((2 j -k)2 -j ) = 1 + 2k

L-1 q=1 1 q L-1 q=0,q =i k -q -q b k = 2 j H j,1 ((2 j -k)2 -j ) = -k L-1 q=0,q =i k -q -q
Remark that, for k = 0, • • • , L -1, the right extrapolation weights satisfy:

a k,n = a k,n , a k,0 = a k,0 , b k = -b k
and do not depend on the level j.

Proof of bound [START_REF] Ma | A study of the construction and application of a Daubechies wavelet-based beam element[END_REF] We distinguish between s integer and not. For integer s we have f, g H s (0,1) = f, g H s-1 (0,1) + 1 0 f (s) (x)g (s) (x) dx f H s-1 (0,1) g H s-1 (0,1)

+ 1 0 f (s) (x)g (s) (x) dx
We can write, for ε < 1/2.

1 0 f (s) (x)g (s) (x) dx ≤ f (s) H -ε (0,1) g (s)

H ε (0,1)

f H s-ε g H s+ε ,
which, combined with the previous bound, since s -1 ≤ s -ε ≤ s + ε, yields [START_REF] Ma | A study of the construction and application of a Daubechies wavelet-based beam element[END_REF].

Let us now consider s ∈ (0, 1). This time we have, for s = [s] and ω = s -[s],

f, g H s (0,1) = f, g H s(0,1) + which, combined with the previous bound, yields the the thesis.

Figure 1 :

 1 Figure 1: Deslaurier-Dubuc functions φ for N = 0, • • • , 5

Figure 2 :

 2 Figure 2: The Hermite functions ϕ h j (in blue) and ϕ h j (in red), for j = 2 (left) and j = 3 (right)

Proposition 2 . 1 .

 21 (IWT) Any continuously differentiable function f ∈ C 1 ([0, 1]) may be reconstructed from its transform by means of f = 2 j 0 +1 k=-1

Figure 4 :

 4 Figure 4: Clamped-clamped beam and the exact solution

12 Figure 5 :Figure 6 :

 1256 Figure 5: Clamped-clamped beam and the exact solution (in blue). Non uniform grid red) for values of ε = 1.e -4, 1.e -6, 1.e -8, 1.e -9 (from left to right and top to bottom)

  s) (x) -f (s) (y))(g (s) (x) -g (s) (y) |x -y| 2ω+1 f H s(0,1) g H s(0,1) s) (x) -f (s) (y))(g (s) (x) -g (s) (y) |x -y| 2ω+1Now, for ε < ω, we can writes) (x) -f (s) (y))(g (s) (x) -g (s) (y) |x -y| 2ω+1 s) (x) -f (s) (y))(g (s) (x) -g (s) (y) |x -y| ω+ε+1/2 |x -y| ω-ε+1/2 |f (s) | H ω-ε (0,1) |g (s) | H ω+ε (0,1)

Table 1 :

 1 Solution of trivial problem in double precision -scaling basis

	h	L 2 error	L ∞ error
	.0625 9.696e-06 1.5677e-05
	.0313 .00015748 .00025088
	.0156	.0025464	.0040258
	.0078	.042886	.067578
	.0039	3.0851	4.8962

Table 2 :

 2 Solution of trivial problem in double precision -wavelet basis

	h	L 2 error	L ∞ error
	.0625 9.696e-06 1.5677e-05
	.0313 9.8408e-06 1.5677e-05
	.0156 9.9162e-06 1.5677e-05
	.0078 9.9545e-06 1.5677e-05
	.0039 9.9734e-06 1.5677e-05

Table 3 :

 3 Solution of trivial problem in multiple precision -30 exact decimal digits for both stage 1 and stage 2 of the computation

	h	L 2 error	L ∞ error
	.0625 2.3871e-25 3.8769e-25
	.0313 3.8769e-24 6.1761e-24
	.0156 6.2506e-23 9.8818e-23
	.0078 1.004e-21 1.5811e-21
	.0039 1.6095e-20 2.5297e-20

Table 4 :

 4 Solution of trivial problem -30 exact decimal digits for stage 1 and double precision for stage 2

	h	L 2 error	L ∞ error
	.0625 3.3959e-13 9.952e-13
	.0313 3.1068e-13 4.6829e-13
	.0156 6.1219e-13 8.5487e-13
	.0078 1.6684e-13 5.0271e-13
	.0039 1.0915e-13 2.4181e-13

Table 5 :

 5 Results for Test 1: multiple precision (left) and double precision (right)

	h	L 2 error	L ∞ error	L 2 error	L ∞ error
	.

Table 6 :

 6 Results for the clamped-clamped beam, half loaded -On the left results with uniform grids, on the left adaptive scheme

		(a)		(b)
	N	L 2 error	ε R	N	L 2 error
	35	0.0039	1.e-4 35	.0039
	67 9.8449 e-04	1.e-6 51 9.8896 e-04
	131 2.4704 e-04	1.e-8 83 2.4859 e-04
	259 6.1876 e-05	1.e-9 120 6.1896 e-05