Erwann Aubry
email: eaubry@math.unice.fr

APPROXIMATION OF THE SPECTRUM OF A MANIFOLD BY DISCRETIZATION

Keywords: Spectral theory, graphs, finite elements, Riemannian geometry, finite elements, discretization

We approximate the spectral data (eigenvalues and eigenfunctions) of compact Riemannian manifold by the spectral data of a sequence of (computable) discrete Laplace operators associated to some graphs immersed in the manifold. We give an upper bound on the error that depends on upper bounds on the diameter and the sectional curvature and on a lower bound on the injectivity radius.

Introduction

We prove that the spectral data (eigenvalues and eigenfunctions) of any closed Riemannian manifold can be approximated by the corresponding spectral data of the Laplace operator of some graphs geodesically immersed in the manifold. It is an extension of the finite elements method to the Riemannian setting. The two main points of our method are the following.

(1) The error made on the spectral data are bounded above by universal functions of some bounds on the geometry of the manifold (i.e. bounds on the injectivity radius, the sectional curvature and the diameter) and of the graph (i.e. bounds on the thinness and mesh of the graph). This errors tend to 0 as the mesh of the graph tends to 0. (2) The Laplace operator of a metric graph is a universal and explicitly computable function of its edge-lengths.

Before stating our main results, we need a few definitions and notations.

1.1. Definitions and notations. We will work with a special kind of immersed graphs, that we will call geodesic triangulations (see the definition in Section 2). They are not necessarily actual triangulations of M (for instance the simplices of dimension greater than 1 are not necessarily realized as subset of M ) but are more easier to construct. A geodesic triangulation T of an n-dimensional Riemannian manifold (M n , g) is a set of points (x i ) 1 i N of M endowed with a structure of abstract simplicial complex K which satisfies the two properties 2.1 and 2.2 of section 2. We denote by S p the set of closed p-simplices of K. We identify the edges of K with some minimizing, geodesic segment between their vertexes. For any closed simplex σ of K, we set St(σ) (resp. St p (σ)) the set of the closed simplices (resp. of dimension p) of K that contains σ. The vertices of any σ ∈ S p are naturally ordered by their indices (σ = {x iσ (0) , . . . , x iσ (p) } with i σ (0) < . . . < i σ (p)). We set X σ = x i σ(0) and for any distinct 0 ≤ k ≤ n, we set v σ k a vector of T Xσ M such that x iσ (k) = exp Xσ (v σ k ). We set also A σ the associated Gramm matrix g

(v σ k , v σ l ) 1 ≤ l ≤ n 1 ≤ k ≤ n
. Given a geodesic triangulation of (M n , g),

we note m T its mesh (the maximal length of its edges) and Θ T its thinness, i.e the quantity Θ T = max max

σ ∈ Sn 0 ≤ k ≤ n m T (det A σ ) -1 2n , max (e1,e2)∈S1
length(e 1 ) length(e 2 ) .

Eventually, on the set R N of functions y : T → R, (where N is the number of vertices of T and we identify T with S 0 ), we define two quadratic forms by the formulae (1.1)

|y| 2 T = 2 (n + 2)! 1≤i≤j≤N y i y j σ∈Stn([xi,xj ]) √ det A σ , (1.2) q T (y) = 1 n! σ∈Sn √ det A σ n k,l=1
(A σ ) kl (y iσ (k) -y iσ (0) )(y iσ (l) -y iσ (0) ).

Note that if K is a sub-complex of R n then |y| T and q T (y) give respectively the L 2 -norm and Dirichlet energy of the affine-by-parts expansion of y.

1.2. Main results. For any closed, Riemannian n-manifold (M, g), we denote by δ M its diameter, by R an upper bound of all its sectional curvatures and by i M its injectivity radius. We denote also by 0 = λ 0 (T ) ≤ • • • ≤ λ N -1 (T ) the eigenvalues of q T with respect to | • | 2 T and 0 = λ 0 (M ) < λ 1 (M ) ≤ • • • ≤ λ i (M ) ≤ • • • the eigenvalues of (M n , g).

Theorem 1.1. Let n ≥ 2 be an integer, and ǫ ∈]0, 1[ be a real number. There exists a constant C(n) such that if i) (M n , g) is a closed, Riemannian n-manifold which satisfies δ 2 M |R| ≤ Λ 2 , ii) T is a geodesic triangulation of M which satisfies mT δM ≤ C(n) iM δM ΘT e e Λ p 3n 3 ǫ, then we have

(1 -ǫ)λ p (T ) ≤ λ p (M ) ≤ (1 + ǫ)λ p (T ).
Remark 1.2. The constant C(n) is computable. The second condition says that any finite number of eigenvalues can be approximated provided the mesh of the graph is small enough and the thinness is controlled.

Remark 1.3. The matrices A σ depend on the angle between some edges of T issued from a same vertex, but the same result is valid if we replace the coefficient g(v σ k , v σ l ) by 1 2

d 2 (X σ , x iσ (k) ) + d 2 (X σ , x iσ (l) ) -d 2 (x iσ (k) , x i σ(l) )
in the definition of the matrix A σ . This gives approximation of the eigenvalues of M by the eigenvalues of a discrete Laplace operator whose coefficients are universal functions of the lengths of a geodesically immersed graph of M . Note that in [START_REF] Burago | A graph discretization of the Laplace-Beltrami operator[END_REF], the authors get the same result for another geometric quadratic form q T , whose coefficients depend on the volume of the Voronoï cells associated to a lattice (x i ) i∈I which need not to be the vertices of a geodesic triangulation.

We denote by (f T i ) the eigenvectors of q T with respect to | • | 2 T and let (f i ) i∈N be a L 2 orthonormal family of eigenfunctions of (M n , g) such that ∆f i = λ i f i for all i ∈ N. For some integers p < q, we set E p,q (resp. F p,q ) the sum of the eigenspaces of ∆(M ) (resp. q T ) associated to the eigenvalues (λ i (M )) p+1 i q (resp. (λ i (T )) p+1 i q and P p,q (resp. Q p,q ) the normal projection on E p,q (resp. F p,q ). Theorem 1.4. Under the assumptions of Theorem 1.1, if there exist some integers p < q and η > 0 such that λ p + η λ p+1 and λ q + η λ q+1 , then for any f ∈ E p , we have

R(f )-P p,q •R(f ) 2 T C(q,n,Λ, δ M i 0 ) √ η ( mT δM ) 1 6n 2 R(f ) 2
T , and for any

(y i ) ∈ F p , we have W (y i )- Q p,q • W (y i ) 2 T C(q, n, Λ, δM i0 , η)( mT δM ) 1 6n 2 W (y i ) 2 T .
To get an approximation of the spectral data of (M n , g) by those of q T in Theorems 1.1 and 1.4, we need some geodesic triangulations with arbitrary small mesh but bounded thinness. The existence of such fat triangulation is often admitted or conjectured but we do not know complete published proof of this fact. For sake of completeness, we give a constructuve proof of the following result (based on the previous work of J.Cheeger, S.Müller and R.Schräder).

Theorem 1.5. Let n ≥ 2 be an integer and D, i 0 and Λ be some positive real numbers. There exist some constants β(n) and C(i 0 /D, Λ, n) such that for any Riemannian manifold (M n , g) with diameter δM D, sectional curvature δ 2 M |R| Λ 2 and injectivity radius i M i 0 and for any ǫ ∈]0, C(i 0 /D, Λ, n)[, there exists a geodesic triangulation T of M with mesh m T ≤ ǫ and thinness Θ T ≤ 1/C(i 0 /D, Λ, n).

Remark 1.6. Once again the constants of Theorem 1.5 are explicitly computable. Combining Theorems 1.5 and 1.1, for any compact manifold (M n , g), any N ∈ N and any ε > 0 we get a method to construct a geodesic triangulation T of M such that we have (1ǫ)λ p (T ) ≤ λ p (M ) ≤ (1 + ǫ)λ p (T ) for any p N .

1.3.

Main steps of the proof. Let (E, • , • ) be a Euclidean space endowed with a bilinear symmetric form q, and λ 0 ≤ • • • ≤ λ dimE-1 be the eigenvalues of q with respect to • , • . Using the min-max principle we readily infer the following spectral comparison principle.

Proposition 1.7 (small eigenvalue principle). Let E 1 , •, • 1 and E 2 , •, • 2 be two Euclidean spaces endowed respectively with quadratic forms q 1 , q 2 . If there exists a linear map Φ : E 1 → E 2 and two positive real numbers α, β such that

Φ(x), Φ(x) 2 ≥ α x, x 1 and q 2 Φ(x), Φ(x) ≤ β q 1 (x, x),
then we have λ k (q 2 ) ≤ β α λ k (q 1 ) for any k.

Proposition 1.7 is usually used to compare spectra under small perturbations on the metric or on the manifold. It is the key tool of our eigenvalues approximation method.

Given a manifold M and a geodesic triangulation T of M , we denote by (x i ) 1≤i≤N the vertices of a T , by E p the subspace of H 1,2 (M ) spanned by the p + 1 first eigenfunctions

(f i ) 0≤i≤p of M , by R : E p → R N the natural restriction map R(f ) = (f (x i )) 1≤i≤N , by •,
• the scalar product on E p induced by the L 2 -norm on M , and we set q(f ) = M |df | 2 . The spectrum of q with respect to •, • is given by the p + 1 first eigenvalues (counted with multiplicities) of M . We then proceed in two steps.

(1) A Moser's iteration scheme gives bounds of the quotients ∇df ∞ f 2 on E p \ {0} by a universal function of λ p , δ and Λ (see proposition 3.1). This Hessian bounds imply the following estimates (see Propositions 4.1 and 4.2)

R(f ), R(f ) T - M f 2 ≤ Cm T M f 2 , (1.3) q T R(f ) ≤ (1 + Cm T ) M |df | 2 , (1.4)
for any f ∈ E p , where | • | T and q T are the discrete quadratic forms on R N given by the formulae (1.1) and (1.2), and where C is a constant which depends on bounds on λ p and on the geometries of (M n , g) and T . Proposition 1.7 gives then some lower bounds on the spectrum of (M n , g) of the form (see Theorem 4.7) λ k (M ) ≥ 1 -Cm T λ k (T ) for any k ≤ p.

(2) In Section 4.3, we construct an expansion (Withney) map W : R N → C ∞ (M ) such that R • W = Id R N and which satifies the following estimates (1.6) for any (y i ) ∈ R N . From Proposition 1.7 again we infer that for any k ≤ p we have 1 -Cm T λ k (M ) ≤ λ k (T ) (see Theorem 4.8). The construction of the withney map is the main technical difficulty of the proof. It is done by local mean of the affine expansions obtained by identifying the simplex of the geodesic triangulation with Euclidean simplicex through the Riemannian exponential maps at the vertices of the simplex. Note that J. Dodziuk [START_REF] Dodziuk | Riemannien structures and triangulations of manifolds[END_REF] developed another generalization of the finite element method to compact Riemannian manifolds in which, to any smooth triangulation of (M, g) is associated the subspace of H 1,2 (M n , g) of the continuous functions on M which are affine on each simplex (this subspace has finite dimension). The authors consider on it the quadratic forms induced by • 2 and the ambient Dirichlet form q(f ) = M |df | 2 . They prove that the spectrum of the discrete Dirichlet form with respect to the the discrete L 2 norm converges to the spectrum of (M n , g) when the mesh of the triangulations tend to 0 with controlled thinness. However, they do not prove that the error is bounded by geometrical bounds on the manifolds, and moreover, the discrete quadratic forms cannot be explicitly computed as function of the geometric data (edge's lengths, edge's angle) of the triangulations.

(y i ) 2 T - M W (y i ) 2 ≤ Cm T (y i ) 2 T + q T (y i ) , (1.5) M |dW (y i )| 2 ≤ 1 + Cm T q T (y i ),
The proof of Theorem 1.4 is done in section 5 using the above estimates and the technique developed by Y. Colin de Verdière in [START_REF] De Verdière | Sur la multiplicité de la première valeur propre non nulle du Laplacien[END_REF].

Aknowledgement We thank S.Gallot for fruitful discussions and C.Vernicos for bringing our attention to the paper [START_REF] Burago | A graph discretization of the Laplace-Beltrami operator[END_REF].

Geodesic triangulations

2.1. Definition. A geodesic triangulation T of an n-dimensional Riemannian manifold (M n , g) is a set of points (x i ) 1 i N of M endowed with a structure of abstract n-dimensional simplicial complex K whose simplices are all contained in a n-dimensional simplex of K and which satisfies two more properties for which we need to complete the notations of the introduction.

For σ ∈ S p and any distinct 0 ≤ k, l ≤ p, v σ kl is a vector of

T x iσ (k) M such that x iσ (l) = exp x iσ (k) (v σ kl ) and |v σ kl | = L σ kl = d(x iσ (k) , x iσ (l) ). We set A σ k the matrix g(v σ kl , v σ km ) l = k m = k .
For any vertex x of T and any σ ∈ St(x), C σ is the cone of T k M spanned by the vectors (v σ (i σ ) -1 (x)l ) 0≤l≤p . Note that for any σ ∈ S n , we have

X σ = x iσ (0) , L σ k = L σ 0k and A σ = A σ 0
(according to the definitions given in the introduction section). For any simplex σ of K, we denote by N p (σ) the set of all the simplices of dimension p that intersect σ.

We set

F = {(θ i ) ∈ R n+1 / i θ i = 1}. For any λ ∈ R, ∆ n λ is the closed n-simplex F ∩ [1 -λ, +∞[ n+1 (we will denote ∆ n = ∆ n 1 )
. Given σ ∈ S n and 0 ≤ k ≤ n, we get some local barycentric coordinates on M by the formula

B σ k : (θ i ) ∈ F → exp x iσ (k) l =k θ l v σ kl .
We set λ • T σ = B σ 0 (∆ n λ ). Eventually, a geodesic triangulation T has to satisfy the following two properties.

(2.1) For any vertex x of T , (C σ ) σ∈St(x) induces a triangulation of the unit sphere of T x M .

(2.2) For any disjoint σ, σ ′ ∈ S n and any 0 ≤ k, k ′ ≤ n, we have

B σ k (∆ n ) ∩ B σ ′ k ′ (∆ n ) = ∅.
Eventually, a geodesic triangulation with boundary T of an n-dimensional Riemannian manifold (M n , g) is a set of points (x i ) 1 i N of M endowed with a structure of abstract ndimensional simplicial complex K whose simplices are all contained in a n-dimensional simplex of K and which satisfies condition (2.2) but condition (2.1) only for vertexes not on the boundary of K, where we call boundary of K the complex of the simplices of K that are contained in a n -1 simplex of K itself contained in only one n-dimensional simplex of K.

Metric estimates.

We now study some metric properties of the geodesic triangulations in bounded geometry. We first recall some estimates on the Riemannian exponential map whose proofs can be found in [START_REF] Buser | Gromov's almost flat manifolds[END_REF].

Theorem 2.1. Let (M n , g) be a compact, Riemannian manifold with δ 2 M |σ| ≤ Λ 2 . Let v ∈ T x M be fixed, y = exp x (v) and for any w ∈ T x M , let w(t) be the parallel translation of w along t → exp x (tv). If we define two maps from T x M to M by F (w) = exp x (v + w) and G(w) = exp y w(1) , then they satisfy the following estimates

d F (w), G(w) ≤ 1 3 d(x, y)|w| gx Λ δ M sinh Λ δ M (d(x, y) + |w| gx ) , |d v exp x (w) -w(1)| ≤ |w| sinh(Λ |v| δM ) Λ |v| δM -1 . Theorem 2.2. Let (M n , g) be a manifold with δ 2 M |σ| ≤ Λ 2 and ǫ < inf i M , δM 2Λ be a positive real. Then for any x ∈ M , the map exp x is a diffeomorphism from B(0 x , ǫ) ⊂ T x M to B(x, ǫ)
and for any u, v ∈ B(0 x , ǫ) we have that

1 -Λ 2 ( ǫ δ M ) 2 |u -v| gx ≤ d exp x (u), exp x (v) ≤ 1 + Λ 2 ( ǫ δ M ) 2 |u -v| gx .
If T is a geodesic triangulation of (M n , g) with mesh smaller than i M /10, then for any σ ∈ S n and any 0 ≤ k ≤ n, the map B σ k gives some coordinates on a neighbourhood of ∆ n 10 . We can compare these coordinates for fixed σ but different values of k.

Lemma 2.3. Let (M n , g) be a manifold with δ 2 M |σ| ≤ Λ 2 and T be a geodesic triangulation with mesh 10 m T ≤ inf i M , δM 2Λ . For any σ ∈ S n , any 0 ≤ k 1 , k 2 ≤ n and any (θ i ) ∈ ∆ n 10 , we have

d B σ k1 (θ i ), B σ k2 (θ i ) ≤ 10 Λ 2 ( m T δ M ) 2 m T i |θ i |. Proof. Let (θ i ) ∈ ∆ n 10 . We set v = v σ k1k2 and w θ = l θ l (v σ k1l -v σ k1k2 ). Then we have B σ k1 (θ i ) = exp x iσ (k 1 ) (v + w θ ). If w θ (t) is the parallel transport of w along s → exp x iσ (k 1 ) (sv σ k1k2 ) then Theorem 2.1 implies that d B σ k1 (θ i ), exp x iσ (k 2 ) (w θ (1)) = d exp x iσ (k 1 ) (v + w θ ), exp x iσ (k 2 ) (w θ (1)) ≤ ( Λ δ M ) 2 m 3 T i |θ i |.
For

θ i = δ il we get d B σ k1 (δ il ), exp x iσ (k 2 ) (w δ il (1)) = d exp x iσ (k 2 ) (v σ k2l ), exp x iσ (k 2 ) (w δ il (1)) ≤ ( Λ δ M ) 2 m 3 T , so, by Theorem 2.2, we have |u σ k2l -v σ k2l | gx iσ (k 2 ) ≤ 2( Λ δM ) 2 m 3 T , where u σ k2l is the parallel transport from x iσ (k1) to x iσ (k2) of the vector v σ k1l -v σ k1k2 . Hence we get d B σ k1 (θ i ), B σ k2 (θ i ) ≤ d B σ k1 (θ i ), exp x iσ (k 2 ) (w θ (1)) + d exp x iσ (k 2 ) (w θ (1)), B σ k2 (θ i ) ≤ ( Λ δ M ) 2 m 3 T i |θ i | + d exp x iσ (k 2 ) ( l θ l u σ k2l ), exp x iσ (k 2 ) ( l θ l v σ k2l ) ≤ ( Λ δ M ) 2 m 3 T + 3( Λ δ M ) 2 m 3 T l |θ l |.
For any τ ∈ S p , with 1 ≤ p ≤ n -1, we set

T ′ τ = σ∈Stn(τ ) ∪ p k=0 exp x iτ (k) { p j=0 θ j v σ ki -1 σ (x iτ (j) ) / θ j ≥ 0, θ j ≤ 1} , T τ = σ∈Stn(τ ) ∪ n k=0 exp x iσ (k) { p j=0 θ j v σ ki -1 σ (x iτ (j) ) / θ j ≥ 0, θ j ≤ 1} .
For any subset A ⊂ M , we set B(A, r) the tubular neighbourhood of A and radius r.

If A is empty we set B(A, r) = ∅.
Lemma 2.3 implies the following result.

Corollary 2.4. There exists a constant

C(n) such that if (M n , g) and T satisfy δ 2 M |σ| ≤ Λ 2 and 10m T ≤ inf i M , δM C(n)Θ 2n
T Λ then we have the following properties

(1) δ Tσ ≤ m T 1 + 10( Λ δM ) 2 m 2 T for any σ ∈ S n , (2) the (1 -η) • T σ σ∈Sn are disjoints and the (1 + η) • T σ σ∈Sn cover M , (3) for any σ, τ ∈ K such that σ ∩ τ = ∅, we have d(T σ , T τ ) ≥ mT C(n)Θ 2n T , (4
) for any σ ∈ S p and τ ∈ K, the tubular neighbourhoods B T τ , α p+1 mT δM m T and B T σ \ B T σ∩τ , α p mT δM m T , α p+1 mT δM m T are disjoint, where

η = C(n) Θ 2n T ( Λ δM ) 2 m 2 T and α = 1 C(n)Θ 2n T . Proof. We set m = m T and Θ = Θ T . Since the points x iσ (k) are in B X σ , m , the n-simplex ∆ σ = (0 Xσ , v σ 01 , . . . , v σ 0n ) of T Xσ M has a diameter less than m 1-( Λ δ M ) 2 m 2 ≤ m(1 + 2( Λ δM ) 2 m 2 ). Moreover, if σ ′ is a (n -1)-face and H is the iso-barycentre of ∆ σ , then the distance from H to σ ′ is equal to n Vol ∆σ (n+1) Vol σ ′ ≥ m C(n)Θ 2n .
Hence, Theorems 2.3 and 2.2 imply that we can choose

C(n) large enough so that T σ have diameter less than m 1 + 4( Λ δM ) 2 m 2 , the (1 + η)T σ contain all the B σ k (∆ n ) (0 ≤ k ≤ n) and (1 -η)T σ ∩ B σ ′ k (∆ n ) = ∅ for all σ ′ ∈ S n \ {σ} and all 0 ≤ k ≤ n. In Particular the (1 -η)T σ σ∈Sn are disjoint.
Let τ ∈ S p be a simplex of T . We now show by recurrence on p that

N τ = σ ∈ N n (τ ) 0 ≤ k ≤ p B σ i -1 σ (x iτ (k) ) (∆ n ) is a neighbourhood of T τ in M and that d(T τ , ∂N τ ) ≥ m C(n)Θ 2n . Note that for any σ ∈ S n , ∆ σ has heights greater than n Vol ∆σ Vol σ ′ ≥ m C(n)Θ 2n
, where σ ′ is the face of σ with dimension n -1 and smallest volume. So the case p = 0 derives from the first axiom of geodesic triangulations and from Lemma 2.3.

If τ is of dimension p ≥ 1 then exp x iτ (k) { p j=0 θ j v σ ki -1 σ (x iτ (j) ) / θ j > 0, θ j < 1} is interior to σ∈Stn(x iτ (k) ) B σ i -1 σ (iτ (k)) (∆ n ) for any 0 ≤ k ≤ p,
by the first axiom of the geodesic triangulation and the smallness of m. Its boundary is also included in N τ by the recurrence assumption and Lemma 2.

3. So T ′ τ is a subset of N τ . It remains to show that d(T ′ τ , ∂N τ ) ≥ m C(n)Θ 2n , which combined with Lemma 2.3 will imply that T τ is included in N τ and that d(T τ , ∂N τ ) ≥ m C(n)Θ 2n
for C(n) large enough. By pulling back the vertices of N τ to T xi τ (0) under the map exp x iτ (0) , we can assume that (M n , g) = (R n , eucl) (by Lemma 2.3 this operation does not change m and Θ too much for C(n) large enough). In that case T τ is a real simplex and by convexity argument on the distance function, d(T τ , ∂N τ ) is bounded from below by the infimum of the distances between disjoint faces of a n-simplex of T τ . For a n-simplex with mesh m and thinness Θ, an easy computation, based on multi-linearity of the determinant, gives that this distance is bounded from below by m C(n)Θ 2n . We easily infer ( 3) from what precedes and from the second axiom of the definition of the geodesic triangulations. We also have

σ∈Sn ∪ n k=0 B σ k (∆ n ) both closed
and open in M , and so equal to M . This implies that the sets (1 + η)T σ σ∈Sn cover M . The property ( 4) is obvious when σ ⊂ τ or τ ⊂ σ and follows from (3) when σ ∩ τ = ∅. In particular, ( 4) is true when σ or τ is a vertex. We now suppose that τ and σ intersect and no one is a subset of the other. As for Point (3), we pull back σ and τ in T z M , where z is a vertex of σ ∩ τ . By Lemma 2.3 it remains to show that ( 4) is satisfied in the Euclidean case.

Let α(n) > 0 such that for any k-face σ and any face τ of ∆ n we have

B σ \ B(σ ∩ τ, α k ), α k+1 ∩ B(τ, α k+1 ) = ∅
By dilation based on a vertex of τ ∩ σ and rate r ≤ 1 we get

B σ \ B(σ ∩ τ, α k r), α k+1 r ∩ B(τ, α k+1 r) = ∅
Since the linear map which maps ∆ n to any T s for s ∈ S n is auto-adjoint with eigenvalues in [ C(n)m Θ , C(n)m] and by Lemma 2.3, we get point (4). Given a geodesic triangulation T of M , we set, for any x ∈ M , dim T (x) = inf{p ≥ 0/ ∃σ ∈ S p , d(x, T σ ) ≤ mT α p+1 δM m T }. This is well defined by Point (2) of Corollary 2.4 as soon as m T ≤ δM C(n)Θ 4n 4 Λ 2 . For any simplex σ ∈ S p , we set

S σ = B T σ , m T α p+1 δM m T \ τ ⊂∂σ B T τ , m T α dim τ +1 δM m T S σ = τ ⊂σ S τ .
The following properties follow readily from Lemma 2.3 and Corollary 2.4. It fundamental for our application to spectral approximations. It says that, even if a geodesic triangulation is not an actual triangulation of the manifold, you can decompose the manifold into some thickening of the generalized faces T σ .

Proposition 2.5. Let (M n , g) be a manifold with

δ 2 M |σ| ≤ Λ 2 and T be a geodesic triangulation of M with 10m T ≤ inf i M , δM C(n)Θ 4n 2 T Λ(1+Λ)
. We have the following properties (1) for any σ ∈ S p , we have

S σ = {x ∈ T σ / dim T (x) = p},
(2) M is the disjoint union of the (S σ ) σ∈K and for any (σ, τ ) ∈ K 2 , we have S σ ∩ S τ = S σ∩τ . (3) for any σ ∈ S p and τ ∈ St n (σ) we have that Vol S σ ≤ C(n) mT δM n-p Vol S τ .

Construction of good geodesic triangulations.

Given a compact manifold, we can use the Riemannian exponential maps to construct some geodesic triangulations at the neighbourhood of any point with bounded thinness and arbitrarily small mesh (image of some Euclidean triangulations of the tangent space) and then adapt the Cheeger-Müller-Schräder ( [START_REF] Cheeger | On the curvature of piecewise flat space[END_REF]) procedure to interpolate these local triangulations in a global, controlled triangulation of the manifold.

To make easier the control of the thinness in our construction, we will work with an alternative (fortunately equivalent in bounded geometry) thinness ΘT of triangulations. In that purpose we replace the Gramm matrix A σ by the matrix

Ãσ = 1 2 d 2 (X σ , x iσ (k) ) + d 2 (X σ , x iσ (l) ) -d 2 (x iσ (k) , x i σ(l) ) k,l
in the definition of the thinness given in the introduction. By Theorem 2.2, if (M n , g) satisfies

δ 2 M |σ| ≤ Λ 2 and m T < inf i M , δM 2Λ then we have 1 -Λ 2 ( m T δ M ) 2 2 |v σ k | 2 gX ≤ d(X, x iσ (k) ) 2 ≤ 1 + Λ 2 ( m T δ M ) 2 2 |v σ k | 2 gX , 1 -Λ 2 ( m T δ M ) 2 2 |v σ l | 2 gX ≤ d(X, x iσ (l) ) 2 ≤ 1 + Λ 2 ( m T δ M ) 2 2 |v σ l | 2 gX , 1 -Λ 2 ( m T δ M ) 2 2 |v σ l -v σ k | 2 gX ≤ d(x iσ (l) , x iσ (k) ) 2 ≤ 1 + Λ 2 ( m T δ M ) 2 2 |v σ l -v σ k | 2 gX , which easily gives A σ kl -Ãσ kl 6m 2 T (Λ mT δM ) 2 .
This easily implies the existence of some functions

C 1 , C 2 such that ΘT C 1 (Θ T , n) and Θ T C 2 ( ΘT , n) as soon we have m T < inf i M , δM 2Λ
. Note that this two thinnesses coincide for Euclidean simplicial complexes and the our thinness is essentially the inverse of the fatness used in [START_REF] Cheeger | On the curvature of piecewise flat space[END_REF].

Let (x i ) i∈I be a maximal family of points of M such that the balls B xi (10

√ ε) are disjoint. Let I 1 , • • • , I k be
k Vol B xi k (120 √ ε) j Vol B xi j (160 √ ε) j Vol B xi j (10 √ ε) max j Vol B xi j (160 √ ε) Vol B xi j (10 √ ε) C(n) Vol B xi k (120 √ ε)
and so we have k C(n) for any ε c(n) Λ . By iteration, we will construct a family of geodesic triangulations, possibly with boundary,

C 1 , • • • , C k in M with mesh less than εC(i, n), thinness ΘCi C(i, n
) and whose vertices of the boundary are outside the set

∪ i j=1 ∪ l∈Ij B x l (30 √ ε -c(i, n)ε) for any ε C(i M , n, δ Λ ). For i = k and ε C(i M , n, δ Λ ), C k will be a geodesic triangulation of M (without boundary since ∪B xi (20 √ ε) = M ) with mesh less than C(k, n)ε and thinness less than Θ(n) = C(k, n).
For any ε > 0, there exists a constant C(n) > 0 such that R n admits a triangulation with mesh less than ε and thinness less than C(n). For any i ∈ I, let T i be such a triangulation of T xi M and T ′ i the subcomplex whose simplices are those of T i contained in B 0x i (30 √ ε) ⊂ T xi M . For any j ∈ {1, • • • , k}, we set K j the simplicial complex ∪ i∈Ij T ′ i . K j is naturally identified with an abstract simplicial complex of M with vertices {exp xi (y), y ∈ T ′ i , i ∈ I j }. By Theorem 2.2, we have m Kj ε/2 and ΘKj 2C(n) for any ε C(n, i M , δM Λ ). We set C 1 = K 1 . Assume that C i is constructed. We now construct C i+1 by interpolation of C i with K i+1 . For any l ∈ I i+1 , we consider in T x l M the complex T ′ l and the complex S l whose vertices are the pull back by exp x l of the vertices of C i that are contained in B 0x l (40 √ ε), and whose simplices have the same combinatorial than in C i . Using Theorem 2.2 as above we get that S l is an Euclidean complex with mesh less than 2C(i, n)ε and thinness less than 2C(i, n) for ε

1 η(i,n) inf(i M , δ 2Λ ) for η(i, n) large enough.
In what follows, for any Euclidean complex T , we denote by T its support, i.e. the union of its simplices. Let A l (respectively A ′ l ) be the complex formed by the simplices of T ′ l (respectively not) contained in S l , and B l (respectively B ′ l ) be the complex formed by the simplices of S l (respectively not) contained in T ′ l . We set D l the set of simplices of T x l M that are intersection of a simplex of A l and of a simplex of B l . Let E l be the set of the simplices σ of A ′ l disjoint from B ′ l and whose intersection with D l is either empty or a union of faces of σ. Similarly, we set F l the set of the simplices σ of B ′ l disjoint from A ′ l and whose intersection with D l is either empty or a union of faces of σ. D l is a polyhedral complex. We obtain a triangulation D ′ l of D l by first barycentric subdivision as follows. Let C be a cell of D l . It is a closed, convex polyhedral cell, as an intersection σ ∩ σ ′ with σ in A l and σ ′ in B l . For each face σ α ⊂ ∂C, we set p α the isobarycentre σ α . The simplices of D ′ l are those spanned by all sets p α1 , ..., p αt , with σ αi σ αi+1 for any 1 i t -1. By the proof of Lemma 6.3, p.439-440, and by Lemma 7.1 3) of [START_REF] Cheeger | On the curvature of piecewise flat space[END_REF], there exists some constants f (Θ, m, n), g(Θ, n) and h(Θ, n) such that, up to a move of the vertices of T ′ l by at most ε2C(i, n)f (2C(i, n), n), the thinness of D ′ l is less than g(2C(i, n), n) and its mesh less than εh(2C(i, n), n). Actually, we perform this deformation of the complex T ′ l before the definition of the complexes A l , A ′ l , B l , B ′ l , D l , E l and F l . Since [START_REF] Cheeger | On the curvature of piecewise flat space[END_REF] allows a control of the thinness of the first barycentric subdivision of any intersection cell of a simplex of (the deformation of) T ′ l and of a simplex of S l , there is no circular definition.

To extend this triangulation to E l ∪ F l , we keep unchanged the simplices of E l ∪ F l that do not intercept D ′ l and we subdivides all the n-simplices of E l (or F l ) with non-empty intersection with D l . For such a simplex σ of E l , we have σ ∩ D l = ∂σ ∩ ∂D l , since any cell of D l is covered by some simplex of A ′ l , and so different from σ. The triangulation D ′ l induces a partition of σ∩D l , which by definition of E l is a triangulation (∂σ) ′ of ∂σ. We then subdivides σ by forming all the simplices spanned its barycentre and by a face of (∂σ) ′ . Thus we get a triangulation G l of E l ∪ F l . Once again, by Lemma 7.1 3) of [START_REF] Cheeger | On the curvature of piecewise flat space[END_REF], the thinness of G l is less than g(2C(i, n), n) and its mesh less than h(2C(i, n), n).

We set

C ′ i+1,l = G l ∪ D ′ l and C(i + 1, n) = max 1, g(2C(i, n), n), h(2C(i, n), n) . Then C ′ i+1,l
is a simplicial complex of T x l M with mesh less than εC(i+1, n) and thinness less than C(i+1, n). Indeed, since G l and D ′ l are simplicial complexes, we just need to remark that a simplex of G l and a simplex of D ′ l intersect along a simplex of ∂D ′ l by what precedes. Moreover, any point of

T ′ l ∪ S l is either in D l or contained in a n-simplex of A ′ l or of B ′ l . So (T ′ l ∪ S l ) \ C ′ i+1,l is covered by the simplices of A ′ l (respectively of B ′ l ) that intersect B ′ l (respectively A ′ l )
or whose intersection with D l is not a union of their faces. In the former case, the simplex obviously intersects both ∂S l and ∂T ′ l and so is at distance from ∂S l and from ∂T ′ l less than 2m C ′ i+1,l . It is the same in the latter case, since if σ is a simplex of (for instance) A ′ l whose non empty intersection with D l is not a union of face of σ, then a face σ ′ of σ must satisfy σ ′ ∩ D l = ∅ and σ ′ \ D l = ∅. It first gives that σ intersect S l (but is not contained in S l ) and so is at distance from ∂S l less than m C ′ i+1,l . By definition of D l , σ ′ ∩ D l is covered by some simplices of T ′ l contained in S l , and since T ′ l is a simplicial complex, we get that σ ′ ⊂ S l . So σ ′ is covered by some simplices of S l . Since σ ′ is not in D l , we infer that σ ′ intersect some simplex of B ′ l . Hence, σ ′ is at distance from ∂T ′ l less than m C ′ i+1,l and σ is at distance from ∂T ′ l less than 2m C ′ i+1,l . We infer that (T ′ l ∪ S l ) \ C ′ i+1,l is covered by the intersection of the 2m C ′ i+1,l -tubular neighbourhoods of ∂T ′ l and ∂S l . We now set C i+1,l the union of the vertices of C i \B x l (35 √ ε) and of the image by exp x l of the vertices of G l ∪ D ′ l . We endow it with the abstract structure of complex obtained by gathering that of C i \ B x l (30 √ ε) and that of G l ∪ D ′ l . Since the complex S l is not deformed during the previous interpolation and since the only vertices of S l that disappears during the interpolation are in B x l (35 √ ε), we really get an abstract structure of simplicial complex on C i+1,l such that any simplex is contained in a n-dimensional simplex. Moreover, by Theorem 2.2, C i+1,l is a triangulation with boundary of (M n , g) with thinness less than 2C(i + 1, n) and mesh less than 2C(i + 1, n)ε for any ε C(n, i M , δ Λ ). Finally, by what precedes the vertices of the boundary of C i+1,l are the same as the vertices of the boundary of C i outside B x l (30 √ ε+C(i, n)ε) and are at distance less than C(i, n)ε from ∂B x l (30 √ ε) inside B x l (30 √ ε+C(i, n)ε) for any ε C(n, i M , δ Λ ) (once again by Theorem 2.2). From this, we get that the vertices of the boundary of C i+1,l are outside

∪ i j=1 ∪ p∈Ij B xp (30 √ ε -c(i + 1, n)ε) ∪ B x l (30 √ ε -c(i + 1, n)ε). C i+1,l
is just the interpolation between C i and T ′ l . But since the family of balls B x l (40 √ ε), l ∈ I i+1 are disjoint, we can interpolate C i with all the T ′ l (l ∈ I i+1 ) simultaneously to get C i+1 . Note that the constant c(i + 1, n) and C(i + 1, n) will be the same whatever the cardinal of I i+1 is since the operations done during the interpolation of two different T ′ l do not interact. So we get the geodesic triangulation C i+1 with all the needed estimates.

Note that the image of any simplex of our geodesic triangulation by the barycentric coordinates map associated to its vertices gives an embedded simplex of M and thus a true triangulation of M whose edges are minimizing geodesic segments.

Estimates on the eigenfunctions

The following proposition gives bounds on the gradient and Hessian of the eigenfunctions. Proposition 3.1. Let (M, g) be a compact Riemannian manifold with delta 2 M |σ| ≤ Λ 2 . For any f ∈ E p , we have that Let ∆ = D * D be the rough Laplace operator on covariant tensors of (M n , g) (where D * is the L 2 adjoint of D). The rough Laplace operator coincides on the functions with the usual Laplace operator and, on the 1-forms it is related to the Hodge Laplace operator by the Bochner formula ∆ = ∆ + Ric. For any tensor T on M it satisfies the equality ∆T,

f ∞ ≤ C(n)e n 2 2 Λ (1 + δ 2 M λ p ) n 4 f 2 , (3.1) δM df ∞ ≤ C(n)e n 2 2 Λ (1 + δ 2 M λ p ) n 4 δM df 2 ≤ C(n)e n 2 2 Λ (1 + δ 2 M λ p ) n+2 4 f 2 , (3.2) δM Ddf ∞ ≤ C(n)e n 3 Λ (1 + δ 2 M λ p ) n 2 2 df ∞ ≤ C(n)e 2n 3 Λ (1 + δ 2 M λ p )
T = |DT | 2 + 1 2 ∆ |T | 2 .
Let T be any tensor on M and V be a field of symmetric endomorphisms on the covariant tensor bundle of M . Fix a real V ≥ 0 such that <V (T ), T > ≥ -V |T | 2 for all tensors T . We set u = |T | 2 + ǫ 2 , then we have that

u∆u = 1 2 ∆(u 2 ) + |du| 2 ≤ 1 2 ∆(|T | 2 ) + |DT | 2 = ∆T, T ≤ |(∆ + V )T |u + V u 2 .
This inequality and the Green formula gives, for any real k > 1/2

d(u k ) 2 2 = k 2 2k -1 M du, d(u 2k-1 ) Vol M = k 2 2k -1 M (u∆u)u 2k-2 Vol M ≤ k 2 2k -1 (∆ + V )T 2k u 2k-1 2k + V u 2k 2k
We apply the above Sobolev inequality to the function u k and make then ǫ tends to 0. This gives us the following inequality

(3.4) T k 2kn n-2 ≤ T k 2k + C(n)k δM e n 2 Λ √ 2k -1 (∆ + V )T 2k T 2k-1 2k + V T 2k 2k
We have E p = Vect{f i , i ≤ p}. For any k ≥ 1, we set

A k = sup f ∈Ep\{0} f k f 2 (resp. B k = sup f ∈Ep\{0} df k df 2 )
. Since E p is stable by ∆ we have

∆f 2k ≤ A 2k ∆f 2 ≤ A 2k λ p f 2 , ∆df 2k ≤ B 2k ∆df 2 ≤ B 2k λ p df 2 .
Hence, by applying inequality 3.4 to T = f and V = 0 (resp. to T = df and V = Ric), for any f ∈ E p , we get

A 2nk n-2 ≤ 1 + C(n)e n 2 Λ k δM λ p √ 2k -1 1/k A 2k B 2kn n-2 ≤ 1 + C(n)e n 2 Λ k √ 2k -1 (n -1)Λ 2 + δ 2 M λ p 1/k B 2k
We multiply the inequalities obtained by setting successively k = ν j with ν = n n-2 > 1 and j ∈ N. Since A m tends to A ∞ (resp. B m tends to B ∞ ) when m tends to ∞, we get

max(A ∞ , B ∞ ) ≤ ∞ j=0 1 + C(n)e n 2 Λ ν j √ 2ν j -1 (n -1)Λ 2 + δ 2 M λ p 1 ν j
.

To get a more convenient upper bound, note that 1 + a

√ b ≤ √ 1 + b (1 + a) and that the infinite product ∞ j=0 1 + ν j √ 2ν j -1 1 ν j converges, hence max(A ∞ , B ∞ ) ≤ C ′ (n)e n 2 4 Λ 1 + (n -1)Λ 2 + δ 2 M λ p n 4 .
This gives us the first three inequalities of proposition 3.1. For any f ∈ E p , we have (see for instance [START_REF] Aubry | Variétés de courbure de Ricci presque minorée, inégalités géométriques optimales et stabilité des métriques extrémales[END_REF])

(3.5) ∆Ddf, Ddf ≤ (D * R)df, Ddf + C(n)Λ 2 δ 2 M |Ddf | 2 + D∆df, Ddf . Now, if we set u = |Ddf | 2 + ǫ 2 , we have u∆u ≤ ∆Ddf, Ddf . From Lemma 3.5, we infer M |d(u k )| 2 ≤ k 2 2k -1 M (u∆u)u 2(k-1) ≤ k 2 2k -1 C(n)Λ 2 δ 2 M u 2k 2k + M D∆df, Ddf u 2(k-1) + M D * Rdf, Ddf u 2(k-1) (3.6)
The divergence theorem applied to the field u 1) .

2(k-1) ∆df, D • df # , gives (∀k ≥ 1) M D∆df, Ddf u 2(k-1) = M |∆df | 2 u 2(k-1) -2(k -1) i ∆df, Ddf (e i ) du(e i )u 2k-3 ≤ M |∆df | 2 u 2(k-1) + 2(k -1) M |∆df ||du|u 2(k-1) ≤ k -1 2 M |du| 2 u 2(k-1) + (2k -1) M |∆df | 2 u 2(k-1) ≤ k -1 2 M |du| 2 u 2(k-1) + 2(2k -1) M (|∆df | 2 + | Ric(df )| 2 )u 2(k-1) ≤ k -1 2 M |du| 2 u 2(k-1) + 2(2k -1)( ∆df 2 ∞ + C(n)Λ 4 δ 4 M df 2 ∞ ) M u 2(k-
The same process applied to u 2(k-1) (tr

1,3 ( R (•,•) df, D • Ddf )) # , combined to the equality i,j
Rdf (e i , e j ),

D 2 df (e i , e j ) = 1 2 |Rdf | 2 ,
gives:

M D * Rdf, Ddf u 2(k-1) = M 1 2 |Rdf | 2 u 2(k-1) + 2(k -1) i,j M R E (e i , e j )df, D ej df du(e i )u 2k-3 ≤ k -1 2 M |du| 2 u 2(k-1) + (2k -1) M |Rdf | 2 u 2(k-1) ≤ k -1 2 M |du| 2 u 2(k-1) + (2k -1) C(n)Λ 4 δ 4 M df 2 ∞ M u 2(k-1) Since M |du| 2 u 2(k-1) = 1 k 2 M |d(u k )| 2 , inequality (3.6) implies: d(u k ) 2 2 ≤ 4k 2 C(n)Λ 2 δ 2 M u 2k 2k + C(n)Λ 4 δ 4 M df 2 ∞ u 2k-2 2k-2 + ∆df 2 ∞ u 2k-2 2k-2 Now, since ∆df ∞ ≤ B ∞ ∆df 2 ≤ B ∞ λ p df ∞ , we get d(u k ) 2 2 ≤ k 2 C(n)e Λn 2 δ 2 M u 2k-2 2k-2 u 2 ∞ + df 2 ∞ δ 2 M (1 + δ 2 M λ p ) n+4 2
.

We can now apply the Sobolev inequality to get:

Ddf k 2kn n-2 ≤ Ddf k 2k + kC(n)e Λn 2 Ddf k-1 2k-2 Ddf ∞ + df ∞ δM (1 + δ 2 M λ p ) 1+ n 4 ≤ Ddf k-1 2k-2 Ddf ∞ 1 + C(n)e Λn 2 k 1 + (1 + δ 2 M λ p ) n+4 4 df ∞ δM Ddf ∞ ≤ Ddf k-1 2k-2 Ddf ∞ 1 + C(n)e Λn 2 k(1 + δ 2 M λ p ) n+4 4
since we can suppose that δM Ddf ∞ ≥ df ∞ . Hence (set k = a l /2 + 1)

F l+1 ≤ 1 + kC(n)e Λn 2 (1 + δ 2 M λ p ) n+4 4 2 ν F l , where ν = n n-2 > 1, F l = Ddf a l Ddf ∞ a l
ν l and (a l ) is the sequence defined by a 0 = 2 and a l+1 = n n-2 (a l + 2). Since the sequence a l /ν l tends to n, we get

Ddf ∞ Ddf 2 ≤ ∞ i=0 1 + C(n)e Λn 2 a i (1 + δ 2 M λ p ) n+4 4
1/ν i .

The previous inequality gives

Ddf ∞ ≤ C(n)e n 3 2 Λ (1 + δ 2 M λ p ) (n+4)n 8
Ddf 2

But if we integrate the Bochner formula ∆df, df

= 1 2 ∆|df | 2 + |Ddf | 2 + Ric(df, df ) we easily get δM Ddf 2 ≤ δ 2 M λ p + (n -1)Λ df 2 So we have δM Ddf ∞ ≤ C(n)e n 3 Λ (1 + δ 2 M λ p ) (n+2) 2 8
df ∞ .

Proposition 3.1 implies that at small scale, the eigenfunctions are almost affine.

Lemma 3.3. Let (M, g) be a compact n-manifold which satisfies δ 2 |σ| ≤ Λ 2 , and T be a geodesic triangulation of M such that 10m T ≤ inf(i M , δM 2Λ ). For any σ ∈ S n , we define a function L f σ on 10 • T σ by

L f σ exp Xσ ( n j=1 θ j v σ j ) = f (X σ ) + n j=1 θ j f (x iσ (j) ) -f (X σ ) .
Then we have the following estimates on 10

• T σ f -L f σ ∞ ≤ C(n)e 2n 3 Λ (1 + δ 2 M λ p ) (n+1) 2 2 m T δM 2 f 2 df -dL f σ ∞ ≤ C(n)Θ 2 e 2n 3 Λ (1 + δ 2 λ p ) 3n 2 4 m T δ M df 2 .
Proof. We set v = n j=1 θ j v σ j and γ(t) = exp Xσ (tv), then we have that

f γ(1) -f (X σ ) -g x Df (x), γ(0) = 1 0 s 0 d 2 dt 2 f • γ(t) dt ds ≤ 1 0 s 0 |Ddf | • γ(t) dt ds | γ(0)| 2 ≤ Ddf ∞ d X σ , γ (1) 2 2 .
This implies that

f exp Xσ ( n j=1 θ j v σ j ) -f (X σ ) - n j=1 θ j g x (Df (x), v σ j ) ≤ C(n)e 2n 3 Λ (1 + δ 2 M λ p ) (n+1) 2 2 m T δM 2 f 2 .
If we set θ j = δ jk , then this inequality gives that

(3.7) f (x iσ (k) ) -f (X σ ) -g x (Df (x), v σ k ) ≤ C(n)e 2n 3 Λ (1 + δ 2 M λ p ) (n+1) 2 2 m T δM 2 f 2 ,
which combined with the previous inequality gives the first result.

Set L the linear form on

T Xσ M such that L f σ = L • exp -1
Xσ . For any w ∈ T Xσ M we set w(t) the parallel transport of w along γ. Then by Theorem 2.1, we have that

(3.8) dL f σ w(1) -L(w) ≤ L d exp -1 Xσ w(1) -w ≤ 4 L |w|( Λm δ ) 2 .
On the other hand, we have that

df w(1) -df (w) ≤ 1 0 d dt df w(t) dt ≤ 1 0 Ddf γ, w(t) dt ≤ Ddf ∞ |v||w|.
Hence we get (3.9)

dL f σ -df w(1) ≤ 2|df -L| + C(n)e n 3 Λ (1 + δ 2 λ p ) n 2 2 df ∞ m δ |w|.
Now, by Inequality (3.7) we have that

df (v σ j ) -L(v σ j ) = df (v σ j ) -f (x iσ (j) ) + f (X σ ) ≤ C(n)Θe n 3 Λ (1 + δ 2 M λ p ) n 2 2 m 2 δ df ∞ ,
and so

|df -L| 2 = jk (A -1 σ ) jk (df -L)(v σ j )(df -L)(v σ k ) ≤ A -1 σ j [(df -L)(v σ j )] 2 ≤ Θ 4 C(n)e 2n 3 Λ (1 + δ 2 M λ p ) n 2 ( m δ ) 2 df 2 ∞ .
If we combine this inequality with Inequality (3.9), we get

|dL f σ -df | ≤ C(n)Θ 2 e n 3 Λ (1 + δ 2 λ p ) n 2 2 m δ df ∞ ,
which gives the result by Proposition 3.1.

4.

Proof of Theorem 1.1

4.1. A discrete L 2 norm. We prove that |R(f )| 2 T gives an approximation of f 2 2 on E p . Proposition 4.1. Let (M n , g) be a compact Riemannian manifold with δ 2 M |σ| ≤ Λ 2 and T a geodesic triangulation with 10m T ≤ inf i M , δM C(n)ΛΘ n T . Then, for all f ∈ E p , we have M f 2 -|R(f )| 2 T ≤ C(n)Θ 2n e 5n 3 Λ 1 + δ 2 M λ p 3n 2 m T δM 2 M f 2 .
Proof. Corollary 2.4 (2) implies that

σ∈Sn (1-η)Tσ f 2 ≤ M f 2 ≤ σ∈Sn (1+η)Tσ f 2 .
We have

(1+η)Tσ f 2 ≤ (1+η)Tσ (L f σ ) 2 + (2 f ∞ + L f σ -f ∞ ) L f σ -f ∞ .
f 2 ≤ (1+η)Tσ (L f σ ) 2 + C(n)e 5n 3 Λ (1 + δ 2 M λ p ) 3n 2 ( m δ ) 2 f 2 2 Vol (1 + η) • T σ
By Theorem 2.2 we have that

(1+η)Tσ (L f σ ) 2 = (1+η)∆σ f (X σ ) + n j=1 θ j f (x iσ (j) ) -f (X σ ) 2 dv exp * Xσ g ≤ (1 + 4Λ 2 ( m δ ) 2 ) n (1+η)∆σ f (X σ ) + n j=1 θ j f (x iσ (j) ) -f (X σ ) 2 dv gX σ ≤ ∆σ f (X σ ) + n j=1 θ j f (x iσ (j) ) -f (X σ ) 2 dv gX σ + 2 n 4 f 2 ∞ Λ 2 ( m δ ) 2 + η Vol ∆ σ ≤ 2 √ det A σ (n + 2)! 0≤j≤j ′ ≤n f (x iσ (j) )f (x iσ (j ′ ) ) + C(n)Θ 2n ( Λm δ ) 2 f 2 ∞ Vol(1 -η)T σ
So, we have

(1+η)Tσ f 2 ≤ 2 √ det A σ (n + 2)! 0≤j≤j ′ ≤n f (x iσ (j) )f (x iσ (j ′ ) ) +C(n)e 5n 3 Λ (1 + δ 2 M λ p ) 3n 2 Θ 2n ( m δ ) 2 f 2 2 Vol (1 -η) • T σ By summing on σ ∈ S n , we get M f 2 -|R(f )| 2 ≤ C(n)e 5n 3 Λ (1 + δ 2 M λ p ) 3n 2 ( m δ ) 2 Θ 2n M f 2
. By the same way, we get the reverse inequality.

4.2.

A discrete Dirichlet energy. We prove that q T R(f ) approximates df 2 2 on E p .

Proposition 4.2. Let (M n , g) be a compact Riemannian manifold with δ 2 M |σ| ≤ Λ 2 and T a geodesic triangulation of M with 10m T ≤ inf i M , δM C(n)Θ n Λ . Then for any f ∈ E p , we have that

q T R(f ) - M |df | 2 ≤ C(n)Θ 2n+4 e 3n 3 Λ (1 + δ 2 λ p ) 2n 2 m T δM M |df | 2 .
Proof. As in the proof of Proposition 4.1, we have

(1+η)Tσ |df | 2 ≤ (1+η)Tσ |dL f σ | 2 + C(n)Θ 4 e 2n 3 Λ (1 + δ 2 M λ p ) n 2 m δM df 2 ∞ Vol (1 + η) • T σ .
Using Lemma 3.3 and the fact that

∆σ |L| 2 = √ det Aσ n! n k,l=1 A kl σ L(v k σ )L(v l σ )
, where L is the linear map defined in the proof of Proposition 4.1 (i.e. L(v k σ ) = f (x iσ (k) )f (X σ )) we get the result.

4.3.

A Withney map. We construct an extending (i.e. Withney type) map which to each (y i ) ∈ R N associates a function f : M → R such that f (x i ) = y i for all 1 ≤ i ≤ N . This function f has to be such that f 2 and |df | 2 be close to |y i | 2

T and q T (y i ). In that purpose we take f almost linear by part.

We need first some controlled partitions of the unity on M associated to the geodesic triangulations. We have obviously ϕ σ = 0 outside (1 + 2ζ)T σ By a volume argument, the number of non zero term in the sum τ d x ψ τ is bounded from above by 2 n Θ 2n , and so we have that

|dϕ σ | = dψ σ τ ψ τ - ψ σ τ dψ τ ( τ ψ τ ) 2 ≤ 5n2 n Θ 4n ζm . Definition 4.4. Let f T = (y i ) ∈ R N . For any σ ∈ K, we set L σ : (θ k ) ∈ F → k θ k y iσ(k) . For m T small enough, the function L fT σ = L σ • (B σ 0 ) -1 is defined on S σ and W (f T ) = σ∈Sn ϕ σ L fT σ is well defined on M .
This extending map W satisfies the following properties.

Proposition 4.5.

If 10m T ≤ inf i M , δM C(n)Λ(1+Λ)Θ 4n 2 T then (1) R • W = Id R N , (2) for any f T ∈ R N , we have that M W (f T ) 2 -|(f T )| 2 T ≤ C(n)(1 + Λ)Θ 2n ( m δ ) 2 |f T | 2 T + δ 2 M q T (f T ) , (3)
for any (f T ) ∈ R N , we have that

M |dW (f T )| 2 -q T (f T ) ≤ C(n)Θ 20n 2 m δ q T (f T ).
Proof. Point (1) is obvious. We now prove point [START_REF] Aubry | Diameter pinching in almost positive Ricci curvature[END_REF]. By Inequality (3.8), we have

|dL fT σ | 2 - |L σ | 2 ≤ 16|L σ | 2 Λm δ 2 . Since |L σ | 2 = n k,l=1 (A σ ) kl (y iσ (k) -y iσ (0) )(y iσ (l) -y iσ (0) ), we have M |dW (f T )| 2 ≥ σ∈Sn Sσ |dW (f T )| 2 = σ∈Sn Sσ |dL fT σ | 2 ≥ σ∈Sn Vol T σ 1 -16( Λm δ ) 2 n+1 1 - mα n δM n |L σ | 2 ≥ 1 -C(n) m δ q T (f T ).
On the other hand, we have that

|dW (f T )| 2 = σ∈Sn Sσ |dW (f T )| 2 + n-1 p=0 τ ∈Sp Sτ |dW (f T )| 2 ≤ 1 + C(n) m δ q T (f T ) + n-1 p=0 τ ∈Sp Sτ |dW (f T )| 2 ,
and for any τ ∈ K

Sτ |dW (f T )| 2 = Sτ | σ∈Stn(τ ) d(L fT σ ϕ σ )| 2 ≤ 2 Sτ σ∈Stn(τ ) ϕ σ |dL fT σ | 2 + 2 Sτ σ∈Stn(τ ) L fT σ dϕ σ 2 . Since |dL fT σ | ≤ (1 + 8( Λm δ ) 2 )|L σ |,
ϕ σ |dL fT σ | 2 ≤ n-1 p=0 τ ∈Sp (1 + 16( Λm δ ) 2 ) σ∈Stn(τ ) |L σ | 2 Vol S τ ≤ C(n) m δM σ∈Sn |L σ | 2 Vol S σ ≤ C(n) m T δM q T (f T ).
To bound the remaining term Sτ σ∈Stn(τ ) L fT σ dϕ σ 2 , we set σ 0 ∈ St n (τ ) and τ 0 the realisation of τ associated to exp Xσ 0 . Then there for any x ∈ S τ , there exists x ′ ∈ τ 0 such that

d(x, x ′ ) ≤ m 2 δ α p+1 . Hence we have |L fT σ0 (x) -L fT σ0 (x ′ )| ≤ 2|L σ0 | m 2 δ α p+1
. Note that the barycentric coordinates (θ k ) of x ′ in the simplex σ satisfy θ k = 0 if x iσ 0 (k) / ∈ τ . For any σ ∈ St n (τ ), we set x σ the point whose barycentric coordinates in σ are θ

′ i -1 σ (x) = θ i -1 σ 0 (x) if
x is a vertex of τ and θ ′ k = 0 otherwise. By applying Lemma 2.3 at most twice, we get d

(x ′ , x σ ) ≤ 200Λ 2 m 3 δ 2 . Since L fT σ (x σ ) = L fT σ0 (x ′ ), we have that |L fT σ (x) -L fT σ0 (x)| ≤ |L fT σ (x) -L fT σ (x σ )| + |L fT σ0 (x ′ ) -L fT σ0 (x)| ≤ C(n)(|L σ | + |L σ0 |) m 2 δ . Since σ dϕ σ = 0, Lemma 4.3 gives us n-1 p=1 τ ∈Sp Sτ σ∈Stn(τ ) L fT σ dϕ σ 2 = n-1 p=1 τ ∈Sp Sτ σ∈Stn(τ ) (L fT σ -L fT σ0 )dϕ σ 2 ≤ n-1 p=1 τ ∈Sp C(n)Θ 10n 2 #St n (τ )( σ∈Stn (τ ) |L σ | 2 + |L σ0 | 2 ) Vol S τ ≤ C(n)Θ 20n 2 m T δM σ∈Sn |L σ | 2 Vol S σ ≤ C(n)Θ 20n 2 m T δM q T (f T ).
We now prove point [START_REF] Aubry | Pincement sur le spectre et le volume en courbure de Ricci positive[END_REF]. As in the proof of point (3), we have that

σ∈Sn Sσ (L fT σ ) 2 ≤ M |W (f T )| 2 ≤ σ∈Sn Sσ (L fT σ ) 2 + n-1 p=0 τ ∈Sp Sτ |W (f T )| 2 and n-1 p=0 τ ∈Sp Sτ |W (f T )| 2 = n-1 p=0 τ ∈Sp Sτ | σ∈Stn(τ ) L fT σ ϕ σ | 2 ≤ n-1 p=0 τ ∈Sp σ∈Stn(τ ) C(n)Θ 2n Sτ |L fT σ | 2 = C(n)Θ 2n σ∈Sn S σ \Sσ |L fT σ | 2 .
Hence we have that

M |W (f T )| 2 - σ Tσ (L fT σ ) 2 ≤ C(n)Θ 2n σ∈Sn S σ \Tσ |L fT σ | 2 ,
and since by Theorem 2.2, we have

1 -( Λm δ ) 2 n |f T | 2 T = 1 -( Λm δ ) 2 n σ ∆σ (L σ ) 2 ≤ σ Tσ (L fT σ ) 2 ≤ 1 + ( Λm δ ) 2 n σ ∆σ (L σ ) 2 ≤ 1 + ( Λm δ ) 2 n |f T | 2 T , it only remains to bound from above Sσ\Sσ |L fT σ | 2 . If we set ζ = nΘ 2n α n+1 m δM , then S σ \ T σ ⊂ (1 + ζ)∆ σ \ (1 -ζ)∆ σ and so S σ \Tσ |L fT σ | 2 ≤ 1 + ( Λm δ ) 2 n (1+ζ)∆σ\(1-ζ)∆σ |L σ | 2
Let H λ be the dilation in F of factor λ and centred at the iso-barycentre of Λ n . By the fundamental theorem of the calculus we have

1 (1 -(n + 1)ζ) n (1-ζ)∆σ (L σ ) 2 - ∆σ (L σ ) 2 ≤ ∆σ (L σ ) 2 H 1-(n+1)ζ (y) -(L σ ) 2 (y) . We set h = (n + 1) m δM ζ. Since (L σ ) 2 H 1-(n+1)ζ (y) -(L σ ) 2 (y) ≤ 2|L σ (y)| dL σ h δM + dL σ 2 δ 2 M h 2 ≤ C(n)( m δM ) 2 |L σ (y)| 2 + δ 2 M dL σ 2 then we have 1 (1 -(n + 1)ζ) n (1-ζ)∆σ (L σ ) 2 - ∆σ (L σ ) 2 ≤ C(n)( m δM ) 2 ∆σ (L σ ) 2 + δ 2 M L σ 2 ,
By the same way we get

1 (1 + (n + 1)ζ) n (1+ζ)∆σ (L σ ) 2 - ∆σ (L σ ) 2 ≤ C(n)( m δM ) 2 ∆σ (L σ ) 2 + δ 2 M L σ 2 Hence σSn (1+ζ)∆σ\(1-ζ)∆σ |L σ | 2 ≤ C(n)( m δM ) 2 (|f T | 2 T +δ 2 M q T (f T ))
, which combined with the previous inequalities gives the result. 4.4. Conclusion. By Propositions 4.1, 4.2 and 1.7 we can almost bound from below the first eigenvalues of (M n , g) by the eigenvalues of q T with respect to •, • T . To have an error bound that depends on p and not on λ p , we use the following rough version of a well known result due to S. Cheng. Lemma 4.6. Let (M n , g) be a compact Riemannian manifold such that δ 2 M |σ| ≤ Λ 2 . For any k ∈ N, we have δ

2 M λ k ≤ C(n)( δM iM ) 2 e ne Λ 2 k 2 .
We infer the following theorem.

Theorem 4.7. Let ǫ ∈]0, 1[ be a real number, (M n , g) be a compact Riemannian manifold such that δ 2 M |σ| ≤ Λ 2 and T be a geodesic triangulation of (M n , g) such that mT δM ≤ ǫ C(n) iM δM ΘT e e Λ p 3n 3

, then we have λ p (q T ) ≤ λ p (M n , g)(1 + ǫ).

Once again, by Propositions 4.5 and 1.7 we can bound from above the eigenvalues of (M n , g) by the eigenvalues of q T with respect to •, • T . Note that to bound λ p (q T ), we just have to use Theorem 4.7.

Theorem 4.8. Let ǫ ∈]0, 1[ be a real number, (M n , g) be a compact Riemannian manifold such that δ 2 M |σ| ≤ Λ 2 and T be a geodesic triangulation of (M n , g) such that mT δM ≤ ǫ C(n) iM δM ΘT e e Λ p 3n 3

. Then, we have that λ p (M n , g) ≤ λ p (q T )(1 + ǫ).

Proof. When applying Proposition 1.7 with E 1 the space spanned by the p+1 first eigenfunction of q T , Proposition 4.5 gives α ≥ 1-C(n)e Λ Θ 2n ( mT δM ) 2 1+δ 2 M λ p (q T ) and β ≤ 1+C(n)Θ 20n 2 m δM . Now, by Theorem 4.7 and Lemma 4.6, we have δ 2 M λ p (q T ) ≤ Proof. We use the same idea as in [START_REF] De Verdière | Sur la multiplicité de la première valeur propre non nulle du Laplacien[END_REF]. We consider in Λ p+1 R N the operator

A(v 0 ∧ • • • ∧ v p ) = p i=0 v 0 ∧ • • • ∧ ∆ T (v i ) ∧ • • • ∧ v p ,
where ∆ T is the symmetric operator such that q T (x) = ∆ T (x), x T . The eigenvalues of A are exactly the sum λ i1 (T ) + • • • + λ ip+1 (T ) with 0 i 1 < • • • < i p+1 N , and so its first eigenvalue is λ 0 (T )+• • •+λ p (T ). We set R :

Λ p+1 E p → Λ p+1 R N defined by R(f 0 ∧ • • • ∧ f p ) = R(f 0 ) ∧ • • • ∧ R(f p ). Since |R(f 0 ∧ • • • ∧ f p )| 2 T = (det G) 2 ,
|R(f 0 ∧ • • • ∧ f p )| 2 T -1 C(p, n, Θ, Λ, δM i M )( m T δM ) 2
By Propositions 4.1 and 4.2 and by Theorem 1.1, we get

AR(f 0 ∧ • • • ∧ f p ), R(f 0 ∧ • • • ∧ f p ) - p i=0 λ i (T )|R(f 0 ∧ • • • ∧ f p ) 2 = p i=0 R(f 0 ) ∧ • • • ∧ (A(R(f i )) -λ i (T )R(f i )) ∧ • • • ∧ R(f p ), R(f 0 ) ∧ • • • ∧ R(f p ) = p i,j=0 R(f 0 ) ∧ • • • ∧ A(R(f i )) -λ i (T )R(f i )), R(f j ) G ij R(f j ) ∧ • • • ∧ R(f p ), R(f 0 ) ∧ • • • ∧ R(f p ) = p i=0 G ii (det G) 2 (q T (R(f i )) -λ i (T )|R(f i )| 2 T ) p i=0 G ii (det G) 2 (q T (R(f i )) -df i 2 ) + p i=0 G ii (det G) 2 (λ i (M ) -λ i (T )) + p i=0 G ii (det G) 2 λ i (T )(1 -|R(f i )| 2 T ) C(p, n, Θ, Λ, δM i M )( m T δM ) 1 3n 2
Let (y i ) an orthonormal family of eigenfunctions of q T associated to the eigenvalue (λ i (T )). We

√ η

The proof of the other estimate is exactly the same, but we first have to bound from below the gap λ p+1 (T )λ p (T ) using the bound on the gap λ p+1 (M )λ p (M ) and Theorem 1.1.

We easily infer Theorem 1.4 from the previous Lemma. Indeed, by Proposition 4.1, the quadratic form |P p • R| 2 on E q takes only values less than 1 + C( mT δM ) 1 6n 2 on the unit sphere of E q and so its trace with respect to •, • T is less than p + C( mT δM ) 1 6n 2 (complete an orthonormal basis of the kernel of P p • R). But the previous lemma, applied for the spectral gap at p, implies that the trace of |P p • R| 2 on E p is close to p and so p -C( mT δM )

1 6n 2 + p+1 i q |P p • R(f i )| 2 = tr|P p • R| 2 p + C( mT δM )
1 6n 2 , and so p+1 i q |P p • R(f i )| 2 C( mT δM ) 1 6n 2 . This gives the result when combined with the previous lemma applied to the spectral gap at q.

  a partition of I into (non empty) parts such that each (B xi (40 √ ε)) i∈Ij is a maximal family of disjoint balls among the (B xi (40 √ ε)) i∈I\∪ k<j I k . Since for i k ∈ I k and for each j < k, B xi k (40 √ ε) has to intersect at least one ball B xi j (40 √ ε) with x ij ∈ I j , the Bishop-Gromov inequality gives us

Remark 3 . 2 . 2 M 2 Λ

 3222 The three first inequalities of Proposition 3.1 are valid under the weaker assumption δ Ric ≥ -Λ 2 g. Proof. The proof of Proposition 3.1 is based on a Moser iteration schema. For any f ∈ H 1 (M ), we have the Sobolev inequality f 2n n-2 ≤ C(n)e n δM df 2 + f 2 (if M is a surface, this inequality, and what follows, is valid with n = 4).

Lemma 4 . 3 . 2 T. 2 T δM m 2 T

 43222 Let (M n , g) be a compact, Riemannian n-manifold with δ(M ) 2 |σ| ≤ Λ 2 and T be a geodesic triangulation of M with 10m T ≤ i M , δM C(n)Λ 2 Θ 4n There exists some smooth functions (ϕ σ ) σ∈Sn such that(1) ϕ σ : M → [0, 1], ϕ σ = 1 on S σ and Supp ϕ σ ⊂ S σ for any σ ∈ S n , (2) σ∈Sn ϕ σ = 1 and |dϕ σ | ≤ C(n)Θ 5nfor any σ ∈ S n .Proof. We set ζ = α n+1 mT 2 δM . Since the distance in (T Xσ M, g Xσ ) between (1 + ζ)∆ σ and T Xσ M \ (1 + 2ζ)∆ σ is larger than ζm nΘ 2n (see the proof of Corollary 2.4) and since exp Xσ is a 1 + Λ 2 ( m δ ) 2 -Lipschitzian map, there exists a functionψ σ : M → [0, 1] such that ψ σ = 1 on (1 + ζ)T σ , ψ σ = 0 outside (1 + 2ζ)T σ and |dψ σ | ≤ 4nΘ 2nζm . We set ϕ σ = ψσ τ ∈Sn ψτ , which is well defined since by Corollary 2.4 we have ∪ τ ∈Sn (1 + ζ)T τ = M . By the same kind of arguments as in the proof of Corollary 2.4, we have (1 -4ζ)T σ ∩ (1 + 2ζ)T τ = ∅ for τ = σ, and so ϕ σ = ψ σ = 1 on (1 -4ζ)T σ .

C(n) δ 2 M e ne Λ 2 p 2 iM 2 . 5 . 4 5. 1 .Lemma 5 . 1 . 2 T 1 6n 2 R(f ) 2 T 2 T 1 6n 2 W

 22541512122212 Proof of Theorem 1.Approximation of the eigenfunctions. To get the relations between the eigenfunctions of (M n , g) and the discrete eigenfunctions, we first prove the following result, where the notations are the same as in the introduction. Let δ > 0 and assume that λ p (M ) + η λ p+1 (M ).For any f ∈ E p , we have R(f ) -P p • R(f ) where P p is the orthogonal projection from R N to the space F p spanned by the first p eigenfunctions of q T .For any(y i ) ∈ F p , we have W (y i ) -Q p • W (y i ) C(p, n, Λ, δM i0 , η)( mT δM ) (y i ) 2T , where Q p is the orthogonal projection from L 2 (M ) to E p .

  where G is the Gramm matrix of the family R(f 0 ), • • • , R(f p ), and since there exists C p such that |(det G) 2 -1| C p |G -I p+1 | for G near I p+1 , Proposition 4.1 and Lemma 4.6 give us

2 T

 2 set R(f 0 ∧ • • • ∧ f p ) = αy 0 ∧ • • • ∧ y p + n with n orthogonal to y 0 ∧ • • • ∧ y p . The above estimates give us |α 2 + |n| 2 -1| C(p, n, Θ, Λ, δM iM )( mT δM ) 2 and δ|n| 2 C(p, n, Θ, Λ, δM iM )( mT δM ) i ) -P p (R(f i )) |R(f 0 ∧ • • • ∧ f p )y 0 ∧ • • • ∧ y p | 2 C(p, n, Θ, Λ, δM iM )( mT δM )

  and by Proposition 2.5, we have

	n-1
	2
	p=0 τ ∈Sp	Sτ σ∈Stn(τ )


Modelling the long-term evolution of permeability in full-scale municipal MBRs: a multivariate statistical modelling approach

Introduction

Over the last decade, filtration in membrane bioreactors (MBR) has been extensively studied since fouling is one of the main bottlenecks for the development of this technology [START_REF] Santos | Membrane bioreactors: Two decades of research and implementation[END_REF]. To help filtration management and design, different models have been proposed and calibrated, essentially on lab-scale and pilot-scale experiments. On the one hand mechanistic models based on physical filtration laws were proposed, resistances in series coupled with mass balances being the most common approach (e.g. [START_REF] Jiang | Characterization and Modelling of Soluble Microbial Products in Membrane Bioreactors[END_REF]Busch et al. 2006;[START_REF] Li | Modelling of membrane fouling in a submerged membrane bioreactor[END_REF][START_REF] Sarioglu | Dynamic in-series resistance modeling and analysis of a submerged membrane bioreactor using a novel filtration mode[END_REF]. These models contain a large number of parameters and variables, such as sludge supernatant composition and especially soluble microbial products (SMP) concentration, shear stress, or size distribution of particles and membrane pores. On the other hand, empirical models based on statistical analysis of data can integrate unknown mechanisms in multivariate linear correlations [START_REF] Jang | Steady-state modeling of bio-fouling potentials with respect to the biological kinetics in the submerged membrane bioreactor (SMBR)[END_REF][START_REF] Tian | Correlating membrane fouling with sludge characteristics in membrane bioreactors: An especial interest in EPS and sludge morphology analysis[END_REF]. In this study, a multivariate statistical approach was used to establish links between long-term fouling and variables usually measured or calculated on full-scale plants. The purpose is to identify relevant variables or groups of variables explaining long-term permeability evolution.

Materials and Methods

A 60,000 PE MBR plant (fig. 1) in the Paris area was extensively monitored over nearly one year. Mixed liquor from the two biological treatment trains (anoxic, anaerobic and sequenced aerobic tanks) is mixed and then filtered by hollow fibre modules (Zeeweed 500d, Zenon) installed in four membrane tanks (MT). For 8 month (January-august 2011, fig. 2), carbon (COD), nitrogen (NH 4 , NO 3 ), and phosphorus concentrations in influent and permeate were measured by online sensors and laboratory analysis and the MLSS concentration was analysed in the reactors. Flows, TMP and other SCADA data were collected and analysed during the whole year 2011. A time series analysis program was implemented in R to isolate stable flux periods and improve the permeability calculations so as to integrate a correction of permeability for (i) temperature, (ii) head loss in permeate tubes, and (iii) pressure sensor drift. A multivariate statistical study was performed to find links between daily operational variables (flux, iron dose, coarse bubble aeration time, SRT, MLSS, MLVSS, temperature, organic load) and factors indicating permeability variations (daily average permeability, impact of backwashes, instantaneous permeability evolution, impact of maintenance cleanings). First, a Principal Component Analysis (PCA) was performed on all explanatory variables to reveal correlations between them and identify independent variables. Then, principal components were used as independent variables to model the long-term permeability evolution. Differences between predictions by the different principal components separately and by combinations of principal components were compared to emphasize which variables are more predictive for the evolution of the long-term permeability.

Results and discussion

The four membrane tanks had significantly different fouling behaviour. Similar trends can be observed during some periods but never for the four MT simultaneously (fig. 3). These differences cannot be explained by differences in concentrations of the fouling components in the sludge because the sludge of the two aeration tanks is mixed prior to filtration. The main trend is a cyclic evolution of permeability: it increases from January to June, and reduces from July on, with an acceleration of the fouling between October and November before intensive chemical cleanings.

The PCA results show the fouling correlation between explanatory variables: the first principal component represents the solids contents of the sludge (SRT, MLSS) and temperature. The second component essentially represents membrane flux and aeration. The third component represents the iron injection (fig. 4). When a multivariate model is calibrated on these three variables, the first principal component seems better correlated with fouling than the eight others (the results with the first three components are shown in fig. 5), while a combination of the first three principal components gives a better result. Unfortunately, the measured variables inside this first component are too correlated to allow discerning effects of each on the membrane fouling. Furthermore, principal components weights in multivariate models differ from one MT to another. Some authors found a correlation between temperature and membrane fouling (Rosenberger et al. 2006, van den Brink et al. 2011) to explain seasonal variations of permeability. In our case, the membrane permeability variations seem not to be correlated with it (for example a permeability decrease occurs before the temperature decrease in MT1, 3 and 4). For other authors, a high SRT which is supposed to influence the SMP concentration (Rosenberger et al. 2006) can have a positive impact on permeability. In the monitored plant the variations in SRT were limited to the range of 50-70d and out of the range 10-50d that are supposed to impact quantities of SMP [START_REF] Kimura | Changes in characteristics of soluble microbial products in membrane bioreactors associated with different solid retention times: Relation to membrane fouling[END_REF][START_REF] Miyoshi | Seasonal variation in membrane fouling in membrane bioreactors (MBRs) treating municipal wastewater[END_REF]Liang et al. 2009[START_REF] Al-Halbouni | Correlation of EPS content in activated sludge at different sludge retention times with membrane fouling phenomena[END_REF][START_REF] Malamis | Fractionation of proteins and carbohydrates of extracellular polymeric substances in a membrane bioreactor system[END_REF]. Even if severe fouling periods seem to occur in the four membrane tanks, they don't occur at the same time. It emphasizes the limits of statistical analysis for determining the causes of fouling. To deal with these limits, a more mechanistic approach could be chosen: for instance one could explain the permeability decreases by the spatial heterogeneity of membrane fouling with complete pore blocking models. Some authors (Li&Wang 2006 or Jiang et al. 2007) emphasize the fact that when parts of a membrane are severely fouled, the local flux on another area of the membrane increases until reaching a critical flux that leads to severe fouling of membrane. The history of each filtration unit should probably be included in the analysis for a better prediction. Further experimental work is required to confirm this hypothesis.

Conclusion

An improved permeability calculation procedure was used to better understand full-scale filtration behaviour in a full-scale MBR. Long-term monitoring of four membrane filtration modules in parallel allows emphasizing phenomena that are not visible on the short-term and are not reproducible on the four membrane tanks. It appears that among the explanatory variables used, a combination of SRT, MLSS, F:M ratio, pH and temperature is best to represent the long-term drifts of permeability. But this approach fails at quantifying Author-produced version of the article published in Procedia Engineeering, 2012, 44, 574-580. The original publication is available at http://www. sciencedirect.com/ doi : 10.1016/j.proeng.2012.08.490 homogenously permeability of the four MT and at predicting rapid fouling events. A complementary mechanistic approach based on historical data could help to explain phenomena like pore blocking or fiber network clogging.

Figures

Figure

  Figure 1: Schematic plant layout

Figure 2 :Figure 3 :Figure 4 :Figure 5 :

 2345 Figure 2: Evolution of main operational variables during the monitored period (for MLSS, the solid line is the on-line measurement and the dashed line is the estimation from dry matter analysis)

  Author-produced version of the article published in Procedia Engineeering, 2012, 44, 574-580. The original publication is available at http://www.sciencedirect.com/ doi : 10.1016/j.proeng.2012.08.490

		rate)	38 h
		Membrane type	Hollow fibre ZEEWEED 500d (Zenon)
	Membranes	Mean pore size Material Membrane area Filtration flux	0.04µm PVDF 4 x 4550 m 2 10 -50 LMH
		Instantaneous SAD m	0.62 Nm 3 /m 2 /h (50% time)
	Acivated sludge	Mean organic loading Mean F:M MLSS (aeration tanks) MLSS (membrane tanks) Temperature Iron dose	4200 kgCOD/d 0.052 kgCOD/kgMLVSS/d 5 -9 g/L 6 -10 g/L 10 -21 °C 400 kgFeCl 3 /d (winter) 1400 kgFeCl 3 /d (summer)


G Gandanegara

X Roboam

B Sareni

Storage device sizing for a hybrid railway traction system by means of bicausal bond graphs

Keywords: bond graph, bicausality, model inversion, dimensioning, electrical engineering, hybrid systems

In this paper, the application of bicausal bond graphs for system design in electrical engineering is emphasized. In particular, it is shown how this approach is very useful for model inversion and parameter dimensioning. To illustrate these issues, a hybrid railway traction device is considered as a case study. The synthesis of a storage device (a supercapacitor) included in this system is then discussed.

INTRODUCTION

complex energetic systems and is complementary to the analysis process. Two main issues can be addressed from this approach: model inversion and Electrical engineering systems are more and more complex and heterogeneous, consisting of elements of parameter sizing. On bond graph formalism, the different natures, strongly coupled between different properties of bicausality can help to solve these physical fields. Within this framework, the system issues [13][14][15][16][17][18][19]. analysis becomes complicated so that a unified

This paper aims at showing how bond graph formalism such as the bond graph (BG) [1][2][3][4] is formalism associated with the bicausal approach can particularly useful. This graphical method illustrates be useful in constructing inverse models and solving the energetic transfers in the system. By analysing the system design problems. A hybrid railway traction causal bond graph, it is possible to know whether device is considered as a case study. Coupling 'energy' the association is physically (and energetically) and 'power' sources is firstly investigated when convenient and then 'numerically consistent'. Homothe consumption cycles vary significantly. Indeed, the geneous modelling and system analysis methods size of the energy source can be adjusted with directly applicable for bond graphs are the major the average power consumption, the power peaks interest of bond graph formalism [2][3][4][5][6][7][8][9][10][11][12].

being provided by the storage device. Furthermore, Having chosen the system architecture and paraif breaking modes appear during the driving mission, meter sizing, the system analysis process consists of it can be possible to store the breaking energy. These verifying whether the device fulfils the requirements: issues are well known for hybrid vehicles, but it has this is usually done from a system's model and also recently appeared in railway traction. its simulation. On the other hand, the synthesis In the paper, the example of a supercapacitor process consists of choosing the system structure dimensioning by means of the bicausal approach is and its sizing, directly starting from the requirements. considered. To present this application, bicausality This 'inverse process' is essential for the design of and model inversion are briefly discussed in the next section. Then, the model of a hybrid traction device is presented and the storage element dimensioning is test this property on the bond graph before applying the bicausality inversion process. For this purpose, several definitions are employed [17]. In a causal bond graph, each bond has one causal stroke. In fact, it can be divided into two causal half-Definition 1 strokes: one for the flow variable and one for the effort variable. Thus, the assignment on the model Two single-input/single-output (SISO) causal paths can be examined by applying this type of causality, are disjoint if they do not pass through any common which is called bicausality.

v(S)= ∑ m i=1 v p i
(2) initial states are given and if a set of parameters are fixed, the other parameters values can be By using these definitions on the direct bond graph, the invertibility condition for an MIMO synthesized. The number of parameters that can (multi-input/multi-output) model with m inputs be determined depends on the degrees of freedom and m outputs is: in the system (number of provided inputs/outputs).

1. If there is no choice for the set of m disjoint SISO In a bicausal bond graph, the causality of each causal paths, the model is structurally not variable is separately examined. With causal halfinvertible. strokes, there are four possibilities. The conventional 2. If there is only one choice for the set of m disjoint or physical causality is a particular case of the SISO causal paths, the model is structurally bicausality where both causal half-stokes are placed invertible. on the same side of the bond. Note that, if the con-3. If there are several choices, the modified sequential ventional causality gives considerable information causality procedure for inversion (MSCAPI) should about the physical meaning of associations, the be applied [16, 17]: bicausality is purely conceptual and is consequently (a) Determine a set S 0 whose order is the smallest only useful for the design process.

among them. In a bicausal bond graph, it is preferred to replace (b) Replace all sources and detectors associated the notation Sf, Se, De, and Df by SS (source-sensor with the control variables or inputs and outelements) [20], even if their causality is not changed puts by SS elements. (both causal half strokes are on the same side of (c) Assign the effort source-flow source causality the bond). In this case, these SS elements are called to the SS output elements and propagate the the effort source-flow sensor (for Se and Df) or the causal information to the SS input elements. flow source-effort sensor (for Sf and De).

This propagation has to arrive and impose An inverse model can only be obtained if the direct the effort detector-flow detector causality on the input SS elements. Other elements take model is invertible. Therefore, it is necessary to firstly the causality due to the bicausality propagation Generally, some derivative causalities appear when the inversion process is applied. An implicit differential and junction conventions. These conventions are:

algebraic equation (DAE) model is obtained and the inverse formulation with integral (int) and derivative $ For 1 junctions:

-effort side: only one bond without a (der) causality is as follows half-stroke close to the junction; -flow side: only one bond without a half-d dt x int = f inv (x int , y spec , … , yn spec ) stroke close to the junction. $ For 0 junctions: u inv =g inv (x int , y spec , … , yn spec ) -effort side: only one bond with a halfwhere x int is the state variable corresponding to the stroke near the junction; element that keeps an integral causality in the inverse -flow side: only one bond with a halfbond graph. The specified inputs ( y spec ) of the stroke near the junction.

inverse model have to be sufficiently differentiable. (d) If there is at least one causality conflict, the

In particular, if the inversion is processed along an model is not invertible. In the opposite case, SISO causal path with a length of n, the correspondit is invertible.

ing input of the inverse model must be n times When the model is invertible, the following prodifferentiable if no symbolic manipulations are made. cedure can be applied to construct the inverse model Some propositions have been made in reference [21] (or the synthesis model) using bicausality [16].

to facilitate the numerical derivation. This issue can also be solved by filtering the specified inputs of the 1. Replace all source and detector elements by SS inverse model. From a specified data file, several elements.

methods of filtering can be used with different cut-2. Apply the bicausality effort source-flow source on off frequencies. In particular, causal (lowpass filter) the output elements and the effort sensor-flow or acausal (centred weighted moving average) filtersensor on the input elements or on the elements ing methods can be applied. In the first case (causal) whose parameters have to be synthesized, in the signal phase is shifted following the cut-off relation to the degrees of freedom.

frequency. In order to avoid this delay an acausal 3. Propagate the bicausality from outputs to inputs.

filtering can be used. Other elements take the causality due to the bicausality propagation with respect to junction conventions.

APPLICATION ON A RAILWAY TRACTION

The obtained bond graph is then called a bicausal DEVICE bond graph.

In order to apply these methods, a railway traction Model derivation system is considered. Basically, the model of this device is a 'simplified vision' of the traction part of A direct formulation with integral causality can be the Alstom BB36000 locomotive [22] (see Fig. 1). The expressed as considered structure is composed of a direct current (d.c.) voltage source, an RLC input filter feeding d dt x= f (x, u); y=g(x, u) an induction machine that drives the mechanical f 12 f 13 f 14 f 16 -C acc -e 16 -e 14 -e 17 In this paper, the induction motor with a d.c.

-e 18 differentiation formula (BDF) [18] is used as an integration method. Note that by considering the driving mission (CBD cycles) as requirements, the electrical con-As a validation of the inversion process, Fig. 5 shows the duty cycle obtained from the inverse model straints can 'directly' be synthesized by means of this bicausal approach from the model inversion.

-I 5 -f 18 -f 19 -f 20 -f 21 -K 4 -e
choice: v 1 (a, f 42 )=12 -f 40 -f 41 -f 42 L 1 =a-e 7 -e 9 -L DC -f 9 -f 10 -e 11 -e 12 -I 6 Given that several choices of the SISO causal path -f 12 -f 13 -f 14 -f 15 -K acc -e
superimposed with the one obtained from the direct model. In this particular case, the validation process This example emphasizes the design capacity of this methodology in electrical engineering in the context consists in using the output of the direct model (i.e. the train velocity) as the input of the inverse of a 'top-down' systemic approach.

Simulations are carried out by means of the 20 Sim model. Therefore the velocity is differentiable at least 10 times and is obtained in the direct model simulation software [24, 25] by using a modified library developed in the authors' laboratory, in order through 10 dynamic elements with integral causality.

In the general case, the specified input signal must be to take bicausality into account. In the case study, having an SISO causal path length of 10, a set of DAE 'sufficiently filtered' to inverse and differentiate the model. Some symbolic manipulations of the derived is obtained for which the output (i. rather simplified, it is energetically significant and 3.2.1 Hybrid systems sufficiently accurate to demonstrate the synthesis In numerous devices, such as the one considered process. This storage element is connected to the here as a case study, the loading power is very timeprimary bus with a d.c./d.c. boost converter including variable (see here the CBD cycle of Fig. 3) so that the 'boost inductance'. The storage current (I store ) is the peak power is far from the average consumed then controllable (see Fig. 6). power. For such systems, the main source must be

The proposed energy management strategy is the 'oversized' to take into account the maximum peak following. power demand, which presents a great drawback from 1. For a given and nearly constant d.c. bus voltage an economic point of view. The general idea in hybrid (V bus ), the loading power variations are equivalent systems is well known for electrical vehicles but it to the loading current variations (I load ). also becomes applicable in other fields. For railway 2. A low pass filtering (LPF) of the loading current is traction systems, hybridization can offer advantages considered, the idea being to control the system such as a reduction of energy consumption and of in such a way that the primary energy source pollution (carbon emissions) when diesel-electric only provides the lowpass filtered loading currents devices are used.

(ILPF load ). As a result, the primary source only pro-In hybrid systems, it is necessary to associate the vides the sum of the average loading power, the main energy source with a storage device, such as average value of the system losses, and the lowsupercapacitors, batteries of accumulators, or inertia frequency power harmonics. On the contrary, the wheels. With such components, the primary energy storage device must provide the high-frequency source will only have to furnish the average loading loading currents (IHF load ). The bandwidth of the LPF power and the average system losses. The variations is determined according to the primary source of the consumed power can be provided by the specifications. storage element(s). In this type of system, the con-3. The system is controlled to ensure that the storage cept of energy sources is used for the sources supplyelement provides the high-frequency (HF) part of ing the average loading power and power sources for the loading current. Note in this case that the those furnishing the loading power variations in average power produced by the storage element short periods according to the driving mission time.

is null on a time interval superior to the inverse of the filtering cut-off frequency. In this way, the 3.2.2 Energy management and storage control loading state of the storage device remains almost constant: the loading management of the storage A supercapacitor will be considered here as the storage element. This device is modelled by a low element is therefore really simplified. T i pB (4) power). This reduction value can be obtained by tuning the filter bandwidth. To accomplish the with objective, a second-order filter is applied with a bandwidth (cut-off frequency) of 0.01 Hz. This T i = L store R store and K= L store T CL frequency is also convenient for filtering the power variations of the considered system with the CBD where T CL is the desired closed-loop time constant. cycles.

The control of the storage element is defined in Fig. 7. The storage element should provide the HF

Storage synthesis part of the consumed current (I load

). Respecting the In the particular case of a supercapacitor, an analytical power transfer orientation gives approach is possible to synthesize the C store value. This analytical process is used as a comparison with

I source =I cf +I in +I load (3)
the bicausal approach which is the main aim of this paper. where

I source #ILPF load Analytical approach.
The energy that can be stored or supplied by the storage element is defined as

I cf #0
The HF current I in is indirectly set by controlling the E store inverse BG model developed in section 3.1. For this 1. The case where V store is maximal (100 per cent): purpose, the train velocity has been filtered with a causal filtering so that the inverse system can be max(V store )=V bus when b max =1 derived with the BDF algorithm of 20 Sim. From the The associated energy is inverse model presented in Fig. 4, the current f 6 and the associated voltage e 6 can be obtained at the

(t)= P t 0 [V bus (t)I in (t)] dt+E store (0) (
E max = 1 2 C(V bus )2 (9)
output of the LRC filter (see Fig. 9). It can be seen 2. The case where V store is minimal (50 per cent): that the voltage of the input filter is almost constant. Therefore, the power can be tuned if the current is min(V store )=0.5V bus when b min =0.5 itself controlled. The associated energy is

E min = 1 2 C A V bus 2 B 2 = 1 8 C(V bus )2 (10)
To generalize, if the discharge depth is x%, then the minimal energy will be

E min = 1 2 C(x{%}V bus )2 (11)
From these assumptions, The bicausal BG of Fig. 4 allows the direct calcu-appropriate and it is interesting to integrate a storage element. lation of the stored energy to be made at any given time. It is sufficient to take the maximal and minimal Numerical applications. Applying both methods values of the stored energy during the CBD cycles.

E max -E min V bus = 3 8 CV bus #max CP
previously presented with a discharge depth of 50 For a 50 per cent discharge depth, the supercapacitor per cent gives: value is obtained using the following formula With the analytical approach

C min = 8 3(V bus )2 (E max -E min ) (14) max AP t 0 I in (t) dt B -min AP t 0 I in (t) dt B =2081.514
The general expression for a discharge depth of x% is

Then C min = 2 (1-x{%})2(V bus )2 (E max -E min ) (15)
C min = 8 3 2081.514 2750 =2.018F Note that the obtained minimal value C min is only valid for a certain driving mission (here the CBD With the bicausality approach cycle). If the power profile is changed, this value C min =2.022F should be recalculated.

With the CBD cycles (cf. Fig. 3), the loading power Note that these methods give quite similar results. The difference is related to the accuracy of both is extremely time variable (from -400 kW to 1 MW). For this kind of application, a hybrid structure is integration methods. analysis of a railway traction system. In Proceed-Second, the bicausality has been used to synings of the 9th European Conference on Power thesize a storage element in the wide application electronics and applications (EPE'01), Graz, Austria, case of hybrid systems. For such systems, the idea is September 2001. to size the primary energy source following average 11 Gandanegara, G., Roboam, X., Sareni, B., and power consumption. For this purpose, a storage Dauphin-Tanguy, G. Modeling and multi-time scale analysis of railway traction systems using bond element, here a supercapacitor, is added in order to graphs. In Proceedings of the International Conprovide the intermittent power. By using bicausality, ference on Bond graph modeling and simulation the profile required for the storage power can be

Storage synthesis validation

Fig. 1

 1 Fig. 1 The BB36000 simplified block diagram

Fig. 2

 2 Fig. 2 Direct causal bond graph of a railway traction device with an equivalent d.c. motor

  e. the duty cycle a) is expressed versus the input (i.e the train velocity). model can also be made, as given in reference [21]. Note that this includes the control strategy, contrary Without any symbolic manipulation [21] of the equations, this input is differentiated 10 times in the to the inverse model. This result proves that the inverse model is validated. inverse model (10 I and C elements are in derivative

Fig. 4 Fig. 5

 45 Fig. 4 Inverse model to obtain the duty cycle a

Fig. 6

 6 Fig. 6 Synoptic of a hybrid system with a storage element

  ) should be provided by the storage element and can easily be deduced from the If the discharge state of the supercapacitor is 50 per bicausal bond graph of Fig.8. In this BG, the loading cent of the nominal voltage, there are two extreme current (SS element) has been synthesized from the cases.

Fig. 8 Fig. 9

 89 Fig.8Inverse model with a submodel to calculate the

Fig. 9

 9 Fig. 9(a). The storage element is represented by a supercapacitor of C store = 3 F with internal resistive In order to validate the bicausal approach based losses R store = 5 mV. The boost inductance is fixed at on storage synthesis, the direct model of Fig. 6 is L store = 1 mH in order to smooth the current ripple correctly. The responses of this direct BG model are simulated with the loading current illustrated in

  21 motor is replaced in order to simplify the bicausal approach of the model inversion. This 'equivalent

	motor' is sized to present the same power balance (identical efficiency) as the one actually obtained		-e 22 -K 5 -e 27 -e 23 -e 24 -e 28 -I 4 -f -e 29 -e 30 24 -f 25 -I 3 -f -f 26 30 -f 27 -f 31
	with the induction motor. The previous voltage source inverter is then also replaced by a d.c.-d.c. chopper. The direct causal bond graph correspond-		-f 32 -f 36	-f 33 -f 37	-K jac -e 33 -f 38 -K ess -e 38 -e 32	-e 35 -e 39 -e 36 -e 40 -I 9 -I eq
	ing to this system is presented in Fig. 2. To test the device behaviour with real conditions, the Central Business District (CBD) cycle has been retained as	-f 40 Third choice: v 3 -f 41 (a, f 42 -f 42	)=11
	the system driving mission [23]. The CBD cycle is considered as a reference for the design of traction systems. Each cycle includes a velocity acceleration	L 3	=a-e 7 -f 12 -f 13 -e 9	-L DC -f -f 14 -f 15 -K acc -e 15 9 -f 10 -e 11 -e 12 -e 14 -I 6 -e 17
	phase, a constant velocity phase at 20 mile/h, a velocity deceleration (or braking) phase, and a phase with zero velocity. The curves of velocity and		-e 18 -e 22	-I 5 -f -e 23 -e 24 18 -I 4 -f -f 19 -f 24 20	-f 21 -f 25	-K 4 -e 21 -f 26 -f 27
	power source applied to the present case study are illustrated in Fig. 3. It can be seen that negative powers (i.e. the regenerative phase) are obtained.		-K 5 -e 27 -f 32 -f 33 -C jac -e 33 -e 28 -e 29 -e 30 -e 32 -I 3 -f 30 -e 35 -e 36 -f 31 -I 9
			-f 36	-f 37	-f 38	-K ess -e 38	-e 39	-e 40	-I eq
	3.1 The inverse model		-f 40	-f 41	-f 42
	In the direct causal bond graph, the d.c. voltage source (U cont ), the duty cycle a of the d.c./d.c. con-) are considered as verter and the resistive force (F res the inputs and the obtained source current and velocity as the outputs. In particular, consider the SISO causal path L from the control variable a to	Fourth choice: v 4 L 4 =a-e 7 -e 9 -L DC -f )=10 (a, f 42 9 -f 12 -f 13 -f 14 -f 15 -C acc -e 15 -f 10 -e 11 -e 12 -e 14 -I 6 -e 17 -e 18 -I 5 -f 18 -f 19 -f 20 -f 21 -K 4 -e 21
	the flow variable (train velocity) at the bond 42 f 42 (see Fig. 2): there are four choices of SISO causal paths, where X indicates that the element X is		-e 22 -K 5 -e 27 -e 23 -e 24 -e 28 -I 4 -f -e 29 -e 30 24 -f 25 -I 3 -f -f 26 30 -f 27 -f 31
	crossed by the causal path:		-f 32	-f 33	-C jac -e 33	-e 32	-e 35	-e 36	-I 9
	First		-f 36	-f 37	-f 38	-K ess -e 38	-e 39	-e 40	-I eq

25 Broenink, J. F. and Kleijn, C. Computer

		APPENDIX
	dynamics using bicausal bond graphs. J. Franklin
	Inst., 2000, 337, 743-769 (Elsevier Science Ltd.). 15 Mechin, O., Marquis-Favre, W., Scavarda, S., and	Notation
	Ferbach, P. Sizing of an electric power steering	C	bond graph generalized capacitance-
	system on dynamic and energetic criteria. In Pro-ceedings of the International Conference on Bond graph modeling and simulation (ICBGM'03), vol. 35, no. 2, Orlando, Florida, January 2003, pp. 137-147 (Society for Computer Simulation).	C( p) De (or Df)	type storage element controller transfer function with p as the Laplace parameter bond graph effort (or flow) detector
	16 Ngwompo, R. F., Scavarda, S., and Thomasset, D.		element
	Inversion of linear time-invariant SISO systems	e	bond graph generalized effort
	modelled by bond graph. J. Franklin Inst., 1996, 333,		variable
	157-174 (Elsevier Science Ltd). 17 Ngwompo, R. F. Contribution au dimensionne-ment des syste `mes sur des crite `res dynamiques et e ´nerge ´tiques. PhD Thesis, INSA de Lyon, 1997. 18 Ngwompo, R. F. and Gawthrop, P. J. Bond graph-	f GY I	bond graph generalized flow variable bond graph gyrator element bond graph generalized inertia-type
	based simulation of non-linear inverse systems using		storage element
	physical performance specifications. J. Franklin Inst., 1999, 336, 1225-1247 (Elsevier Science Ltd). 19 Ngwompo, R. F. and Scavarda, S. Dimensioning problems in system design using bicausal bond graphs. Simulation Practice and Theory, 1999, 7, 577-587 (Elsevier Science BV).	I cf I in I load IHF load ILPF load	filter (bus) capacitor current input current of the storage converter element loading current high-frequency loading current lowpass filtered loading current
	High speed railway traction system modelling for	SS	generic bond graph source-sensor
	simulating electromechanical interactions. In Pro-		element
	ceedings of the World Conference on Railway research (WCRR'97), Firenze, Italy, November 1997. 23 Emission test cycles: Central Business District (CBD), Ecopoint Inc., URL: http://www.dieselnet.com/ standards/cycles/cbd.html.	TF V bus V store x	bond graph transformer element bus voltage supercapacitor voltage average value of x(t)
	24 20-Sim, Controllab Products BV, URL: http://
	www.20sim.com.	a or b	converter duty cycle
	-aided	0	bond graph common effort junction
	design of mechatronic systems using 20-Sim 3.0. In Proceedings of the 2nd Workshop on European	1	bond graph common flow junction
	scientific and industrial collaboration (WESIC'99),
	Newport, UK, 1999, pp. 27-34.
	directly deduced and a value of the supercapacitor can be efficiently synthesized following the choice of the discharge depth. Note that the minimal value	(ICBGM'01), vol. 33, no. 1, Phoenix, Arizona, January 2001, pp. 62-67 (Society for Computer Simulation). 12 Gandanegara, G., Roboam, X., Sareni, B., and Dauphin-Tanguy, G. One model for one frequency
	of the supercapacitor is only valid for the considered	range: comparison of bond graph based simplifi-
	mission cycle ensuring that the finality character	cation methods. In Proceedings of the International
	(i.e. the actual functioning mode of the system) is	Conference on Bond graph modeling and simulation
	directly linked with the design process, as for any systemic approach. More generally, this application example has emphasized the capacity of the bicausal methodology	(ICBGM'03), vol. 35, no. 2, Orlando, Florida, January 2003, pp. 53-58 (Society for Computer Simulation). 13 Gawthrop, P. J. Bicausal bond graphs. In Proceed-ings of the International Conference on Bond graph modeling and simulation 1995 (ICBGM'95), Las
	with a view to a systemic approach using electrical	Vegas, USA, vol. 27, January 1995, pp. 83-88 (Society
	engineering design.	for Computer Simulation).

20 Gawthrop, P. J. and Smith, L. Causal augmentation of bond graphs with algebraic loops. J. Franklin Inst., I source main source current 1992, 329(2), 291-303 (Elsevier Science Ltd). I store supercapacitor current 21 van Dijk, J. and Breedveld, P. C. Simulation of R bond graph generalized dissipated system models containing zero-order causal path -I. element Classification of zero-order causal paths. J. Franklin Se (or Sf) bond graph effort (or flow) source Inst., 1991, 328(5/6), 959-979 (Elsevier Science Ltd). element 22 Lochot, C., Roboam, X., de Fornel, B., and Moll, F.