François Boulier
email: francois.boulier@univ-lille1.fr

François Lemaire
email: francois.lemaire@univ-lille1.fr

Finding First Integrals Using Normal Forms Modulo Differential Regular Chains

Keywords: First integral, linear algebra, differential algebra, nonlinear system

This paper introduces a definition of polynomial first integrals in the differential algebra context and an algorithm for computing them. The method has been coded in the Maple computer algebra system and is illustrated on the pendulum and the Lotka-Volterra equations. Our algorithm amounts to finding linear dependences of rational fractions, which is solved by evaluation techniques.

Introduction

This paper deals with the computation of polynomial first integrals of systems of ODEs, where the independent variable is t (for time). A first integral is a function whose value is constant over time along every solution of a system of ODEs. First integrals are useful to understand the structure of systems of ODEs. A well known example of first integral is the energy of a mechanical conservative system, as shown by Example 1.

Example 1. Using a Lagrange multiplier λ(t), a pendulum of fixed length l, with a centered mass m submitted to gravity g can be coded by:

x y l m g Σ      m ẍ(t) = λ(t) x(t) m ÿ(t)
= λ(t) y(t) + mg x(t) 2 + y(t) 2 = l 2 .

(1)

A trivial first integral is x(t) 2 + y(t) 2 since x(t) 2 + y(t) 2 = l 2 on any solution. A less trivial first integral is m 2 ẋ(t) 2 + ẏ(t) 2 -mg y(t) which corresponds to the total energy of the system (kinetic energy + potential energy).

When no physical considerations can be used, one needs alternative methods. Assume we want to compute a polynomial first integral of a system of ODEs. Our algorithm findAllFirstIntegrals, which has been coded in Maple, relies on the following strategy. We choose a certain number of monomials µ 1 , µ 2 , . . . , µ e built over t, the unknown functions and their derivatives (namely t, x(t), y(t), ẋ(t), ẏ(t), . . . on Example 1), and look for a candidate of the form q = α 1 µ 1 +• • •+α e µ e satisfying dq (t) dt = 0 for all solutions, where the α i are in some field K. If the α i are constant (i.e. dαi dt = 0 for each α i ), then our problem amounts to finding α i such that dq (t) dt = α 1 dµ1 dt + • • • + α e dµe dt is zero for all solutions. On Example 1, µ 1 , µ 2 and µ 3 could be the monomials ẋ(t) 2 , ẏ(t) 2 and y(t) and α 1 , α 2 and α 3 could be m/2, m/2 and -mg.

Anticipating Section 4, differential algebra techniques combined with the Nullstellensatz Theorem permit to check that a polynomial p vanishes on all solutions of a radical ideal A by checking that p belongs to A. Moreover, in the case where the ideal A is represented 1 by a list of differential regular chains M = [C 1 , . . . , C f ], checking that p ∈ A can be done by checking that the normal form (see Section 4) of p modulo the list of differential regular chains M , denoted NF(p, M ), is zero. Since the normal form is K-linear (see Section 4), finding first integrals can be solved by finding a linear dependence between NF( dµ1 dt , M ), NF( dµ2 dt , M ), . . . , NF( dµe dt , M ). This process, which is in the spirit of [START_REF] Boulier | Efficient computation of regular differential systems by change of rankings using Kähler differentials[END_REF][START_REF] Faugère | Efficient computation of Gröbner bases by change of orderings[END_REF], encounters a difficulty: the normal forms usually involve fractions. As a consequence, we need a method for finding linear dependences between rational functions. A possible approch consists in reducing all fractions to the same denominator, and solving a linear system based on the monomial structure of the numerators. However, this has disadvantages. First, the sizes of the numerators might grow if the denominators are large. Second, the linear system built above is completely modified if one considers a new fraction. Finding linear dependences between rational functions is addressed in Sections 2 and 3 by evaluating the variables occurring on the denominators, thus transforming the fractions into polynomials. The main difficulty with this process is to choose adequate evaluations to obtain deterministic (i.e. non probabilistic) and terminating algorithms.

This paper is structured as follows. Section 2 presents some lemmas around evaluation and Algorithm findKernelBasis (and its variant findKernelBasisMatrix) which computes a basis of the linear dependences of some given rational functions. Section 3 presents Algorithm incrementalFindDependence (and its variant incrementalFindDependenceLU) which sequentially treats the fractions and stops when a linear dependence is detected. Section 4 recalls some differential algebra techniques and presents Algorithm findAllFirstIntegrals which computes a basis of first integrals which are K-linear combinations of a predefined set of monomials.

In Sections 2 and 3, k is a commutative field of characteristic zero containing R, A is an integral domain with unit, a commutative k-algebra and also a k-vector space equipped with a norm • A . Moreover, K is the total field of fractions of A. In Section 4, we will require, moreover, that K is a field of constants. In applications, k is usually R, A is usually a commutative ring of multivariate polynomials over k (i.e. A = k[z 1 , . . . , z s ]) and K is the field of fractions of A (i.e. the field of rational functions k(z 1 , . . . , z s )). This section presents Algorithms findKernelBasis and findKernelBasisMatrix which compute a basis of the linear dependences over K of e multivariate rational functions denoted by q 1 , . . . , q e . More precisely, they compute a K-basis of the vectors (α 1 , . . . , α e ) ∈ K e satisfying e i=1 α i q i = 0. Those algorithms are mainly given for pedagogical reasons, in order to focus on the difficulties related to the evaluations (especially Lemma 2). The rational functions are taken in the ring K(Y )[X] where Y = {y 1 , . . . , y m } is a finite set of indeterminates and X is another (possibly infinite) set of indeterminates such that Y ∩ X = ∅. The idea consists in evaluating the variables Y , thus reducing our problem to finding the linear dependences over K of polynomials in K[X], which can be easily solved for instance with linear algebra techniques. If enough evaluations are made, the linear dependences of the evaluated fractions coincide with the linear dependences of the non evaluated fractions.

Even if it is based on evaluations, our algorithm is not probabilistic. A possible alternative is to use [START_REF] Stoutemyer | Multivariate partial fraction expansion[END_REF] by writing each rational function into a (infinite) basis of multivariate rational functions. However, we chose not to use that technique because it relies on multivariate polynomial factorization into irreducible factors, and because of a possible expression swell.

Preliminary Results

This section introduces two basic definitions as well as three lemmas needed for proving Algorithm findKernelBasis.

Definition 1 (Evaluation). Let D = {g 1 , . . . , g e } be a set of polynomials of K[Y ] where Y = {y 1 , . . . , y m }. Let y 0 = (y 0 1 , y 0 2 , . . . , y 0 m ) be an element of K m , such that none of the polynomials of D vanish on y 0 . One denotes σ y 0 the ring homomorphism from K[X, Y, g -1 1 , . . . , g -1 e ] to K[X] defined by σ y 0 (y j ) = y 0 j for 1 ≤ j ≤ m and σ y 0 (x) = x for x ∈ X. Roughly speaking, the ring homomorphism σ y 0 evaluates at y = y 0 the rational functions whose denominator divides a product of powers of g i .

Definition 2 (Linear combination application). Let E be a K-vector space. For any vector v = (v 1 , v 2 , . . . , v e ) of E e , Φ v denotes the linear application from

K e to E defined by α = (α 1 , α 2 , . . . , α e ) → Φ v (α) = e i=1 α i v i .
The notation Φ v defined above is handy: if q = (q 1 , . . . , q e ) is (for example) a vector of rational functions, then the set of the vectors α = (α 1 , . . . , α n ) in K e satisfying e i=1 α i q i = 0 is simply ker Φ q . The following Lemma 1 is basically a generalization to several variables of the classical following result: a polynomial of degree d over an integral domain admitting d + 1 roots is necessarily the zero polynomial.

Lemma 1. Let p ∈ K[Y ] where Y = {y 1 , . . . , y m }. Assume deg yi p ≤ d for all 1 ≤ i ≤ m. Let S 1 , S 2 , . . . , S m be m sets of d + 1 points in Q. If p(y 0 ) = 0 for all y 0 ∈ S 1 × S 2 × . . . × S m , then p is the zero polynomial.
Proof. By induction on m. When m = 1, p has d + 1 distinct roots and a degree at most d. Since p is a polynomial over a field K, p is the zero polynomial.

Suppose the lemma is true for m. One shows it is true for m + 1. Seeing p as an element of K(y 1 , . . . , y m )[y m+1 ] and using the Lagrange polynomial interpolation formula [9, where S m+1 = {s 1 , . . . , s d+1 }.

For each 1 ≤ i ≤ d+1, p(y 1 , . . . , y m , s i ) vanishes on all points of S 1 ×S 2 ×. . .×S m . By induction, all p(y 1 , . . . , y m , s i ) are thus zero, and p is the zero polynomial.

The following lemma is quite intuitive. If one takes a nonzero polynomial g in K[Y ], then it is possible to find an arbitrary large "grid" of points in N m where g does not vanish. This lemma will be applied later when g is the product of some denominators that we do not want to cancel. Lemma 2. Let g be a nonzero polynomial in K[Y ] and d be a positive integer. There exist m sets S 1 , S 2 , . . . , S m , each one containing d + 1 consecutive nonnegative integers, such that g(y 0 ) = 0 for all

y 0 in S = S 1 × S 2 × • • • × S m .
The proof of Lemma 2 is the only one explicitly using the norm • A . The proof is a bit technical but the underlying idea is simple and deserves a rough explanation. If g is nonzero, then its homogeneous part of highest degree, denoted by h, must be nonzero at some point ȳ, hence on a neighborhood of ȳ. Since g "behaves like" h at infinity (the use of • A will make precise this statement), one can scale the neighborhood to prove the lemma.

Proof. Since K is the field of fraction of A, one can assume with no loss of generality that g is in A[Y ].

Denote g = h + p where h is the homogeneous part of g of (total) degree e = deg g. Since, g is nonzero, so is h. By Lemma 1, there exists a point ȳ ∈ R m >0 , such that h(ȳ) = 0, where R >0 denotes the positive reals. Without loss of generality, one can assume that ȳ ∞ = 1 since h is homogeneous. There exists an open neighborhood V of ȳ such that V ⊂ R m >0 and 0 / ∈ h(V ). Moreover, this open neighborhood V also contains a closed ball B centered at ȳ for some positive real i.e. B = {y ∈ R m >0 | y -ȳ ∞ ≤ } ⊂ V . Since B is a compact set and the functions h and t are continuous, the two following numbers are well-defined and finite: m = min y∈B h(y) A and M = max y∈B ( i∈I a i A m i (y)) (where p = i∈I a i m i with a i ∈ A and m i is a monomial in Y ). Moreover, m > 0 (since by a compactness argument there exists y ∈ B such that h(y) = m).

Take y ∈ B and s > 1 ∈ R. Then g(sy) = h(sy) + p(sy) = s e h(y) + p(sy). By the reverse triangular inequality and the homogeneity of h, one has g(sy) A ≥ s e h(y) A -p(sy) A . Moreover

p(sy) A ≤ i∈I a i A m i (sy) ≤ i∈I a i A s e-1 m i (y) ≤ s e-1 M .
Consequently, g(sy) A ≥ s e m-s e-1 M . If one takes s sufficiently large to ensure s e m-s e-1 M > 0 and s ≥ (d+1)/2, the ball B obtained by uniformly scaling B by a factor s contains no root of g. Since the width of the box B is at least d + 1, the existence of the expected S i is guaranteed.

Roughly speaking, the following lemma ensures that if a fraction q in K(Y )[X] evaluates to zero for well chosen evaluations of Y , then q is necessarily zero. Lemma 3. Take an element q in K(Y )[X]. Then, there exist an integer d, and m sets S 1 , S 2 , . . . , S m , each one containing d+1 nonnegative consecutive integers, such that σ y 0 (q) is well-defined for all y 0 ∈ S = S 1 × S 2 × • • • × S m . Moreover, if σ y 0 (q) = 0 for all y 0 in S, then q = 0.

Proof. Let us write q = p/g where p ∈ K[X, Y ] and g ∈ K[Y ]. Consider d = max 1≤i≤m deg yi (p). By Lemma 2, there exist m sets S 1 , . . . , S m of d+1 consecutive nonnegative integers such that g(y 0 ) = 0 for all y 0 in S = S 1 ×S 2 ו • •×S m . Consequently, σ y 0 (q) is well-defined for all y 0 in S. Let us assume that σ y 0 (q) = 0 for all y 0 in S. As a consequence, one has σ y 0 (p) = 0 for any y 0 in S. By Lemma 1, one has p = 0, hence q = 0.

Algorithm findKernelBasis

We rely on the notion of iterator, which is a tool in computer programming that enables to enumerate all values of a structure (such as a list, a vector, . . . ). An iterator can be finite (if it becomes empty after enumerating a finite number of values) or infinite (if it never gets empty).

In order to enumerate evaluation points (that we take in N m for simplicity), we use two basic functions to generate one by one all the tuples, which do not cancel any element of D, where D is a list of nonzero polynomials of K[Y ]. The first one is called newTupleIteratorNotCancelling(Y, D). It builds an iterator I to be used with the second function getNewTuple(I). Each time one calls the function getNewTuple(I), one retrieves a new tuple which does not cancel any element of D. The only constraint we require for getNewTuple is that any tuple of N m (which does not cancel any element of D) should be output after a finite number of calls to getNewTuple. To ensure that, one can for example enumerate all integer tuples by increasing order (where the order is the sum of the entries of the tuple) refined by the lexicographic order for tuples of the same order. Without this constraint on getNewTuple, Algorithm findKernelBasis may not terminate.

Example 2. Take K = Q(z), Y = {a} and X = {x}. Consider q = (q 1 , q 2 , q 3 ) = ( ax+2 a+1 , z(a-1-ax) 1+a , 1
). One can show that the only (up to a constant) K-linear dependence between the q i is -q 1 -(1/z)q 2 + q 3 = 0.

Apply Algorithm findKernelBasis. One has D = [a + 1] at Line 3, so one can evaluate the indeterminate a on any nonnegative integer. At Line 7, S contains the tuple s 1 = (0), corresponding to the evaluation a = 0. One builds the vectors v 1 , v 2 and v 3 at Line 7 by evaluating q 1 , q 2 and q 3 on a = 0. One obtains polynomials in K Algorithm 1: findKernelBasis algebra over K, which contains z. The test at Line 12 fails because the vector b 1 = (-1/2, 0, 1) of B does not yield a linear dependence over the q i . Indeed, -(1/2)q 1 + 0q 2 + q 3 = 0. Consequently, S is enriched with the new tuple s 1 = (1) at Line 14 and one performs a second iteration. One builds the vectors v 1 , v 2 and v 3 at Line 7 by evaluating q 1 , q 2 and q 3 on a = 0 and a = 1. Thus v 1 = (2, 1 + x/2), v 2 = (-z, -xz/2) and v 3 = (1, 1). A basis B computed at Line 10 could be (-1, -1/z, 1). This time, the test at Line 12 succeeds since -q 1 -(1/z)q 2 +q 3 = 0. The algorithm stops by returning the basis (-1, -1/z, 1).

[x]. Thus v 1 = (2), v 2 = (-z) and v 3 = (1). A basis of ker Φ v computed at Line 10 could be B = {(-1/2, 0, 1), (z/2, 1, 0)}. Note that it is normal that both v = (v 1 , v 2 , v 3 )
Proof (Algorithm findKernelBasis). Correctness. If the algorithm returns, then ker Φ q ⊃ ker Φ v . Moreover, at each loop, one has ker Φ q ⊂ ker Φ v since α i q i = 0 implies α i σ(q i ) = 0 for any evaluation which does not cancel any denominator. Consequently, if the algorithm stops, ker Φ q = ker Φ v , and B is a basis of ker Φ q . Termination. Assume the algorithm does not terminate. The vector space ker Φ v is smaller at each step since the set S grows. By a classical dimension argument, the vector space ker Φ v admits a limit denoted by E, and reaches it after a finite number of iterations. From ker Φ q ⊂ ker Φ v , one has ker Φ q ⊂ E. Since the test at Line 12 always fails (the algorithm does not terminate), ker Φ q E. Take α ∈ E \ ker Φ q and consider the set S obtained by applying Lemma 3 with q = e i=1 α i q i . Since the algorithm does not terminate, S will eventually contain S. By construction of v, one has e i=1 α i σ y 0 (q i ) = 0 = σ y 0 (q) for all y 0 in S. By Lemma 3, one has q = e i=1 α i q i = 0 which proves that α ∈ ker Φ q . Contradiction since α / ∈ ker Φ q .

Input: Two lists of variables Y and X Input: A vector q = (q1, . . . , qe) of e rational fractions in K(Y )[X] Output: A basis of ker Φq 1 begin 2 For each 1 ≤ i ≤ e, denote the reduced fraction qi as fi/gi Algorithm 2: findKernelBasisMatrix

with fi ∈ K[Y, X], gi ∈ K[Y ];

A Variant of findKernelBasis

A variant of algorithm findKernelBasis consists in incrementally building a matrix M with values in K, encoding all the evaluations. Each column of M corresponds to a certain q j , and each line corresponds to a couple (s, m), where s is an evaluation point, and m is a monomial of K[X]. This yields Algorithm findKernelBasisMatrix.

Example 3. Apply Algorithm findKernelBasisMatrix on Example 2. The only difference with Algorithm findKernelBasis is that the vectors v i are stored vertically in a matrix M that grows at each iteration. At the end of the first iteration of the while loop, one has y 0 = (0), q = (2, -z, 1), L = [START_REF] Boulier | The BLAD libraries[END_REF], and

M = N = 2 -z 1 , and B = {(-1/2, 0, 1), (z/2, 1, 0)}. Another iteration is needed since -(1/2)q 1 + 0q 2 + q 3 = 0.
At the end of the second iteration, one has y 0 = (1), q = (1 + x/2, -xz/2, 1),

L = [1, x], N = 1 0 1 1/2 -z/2 0 , and M =   2 -z 1 1 0 1 1/2 -z/2 0   . A basis B is (-1, -1/z, 1
) and the algorithm stops since -q 1 -(1/z)q 2 + q 3 = 0.

Proof (Algorithm findKernelBasisMatrix). The proof is identical to the proof of findKernelBasis. Indeed, the vector v i in findKernelBasis is encoded vertically in the i th column of the M matrix in findKernelBasisMatrix. Thus, the kernel of Φ v in findKernelBasis and the kernel of M in findKernelBasisMatrix coincide.

In terms of efficiency, Algorithm findKernelBasisMatrix is far from optimal for many reasons. First, the number of lines of the matrix M grows excessively. Second, a basis B has to be computed at each step of the loop. Third, the algorithm needs to know all the rational functions in advance, which forbids an incremental approach. Next section addresses those issues.

Incremental Computation of Linear Dependences

The main idea for building incrementally the linear dependences is presented in Algorithm incrementalFindDependence. It is a sub-algorithm of Algorithms findFirstDependence and findAllFirstIntegrals. Assume we have e linearly independent rational functions q 1 , . . . , q e and a new rational function q e+1 : either the q 1 , . . . , q e+1 are also linearly independent, or there exists (α 1 , . . . , α e+1 ) ∈ K e+1 such that e+1 i=1 α i q i = 0 with α e+1 = 0. It suffices to iterate this idea by increasing the index e until a linear dependence is found. When such a dependence has been found, one can either stop, or continue to get further dependences.

Algorithm incrementalFindDependence

The main point is to detect and store the property that the q 1 , . . . , q e are linearly independent, hence the following definition. Moreover, for efficiency reasons, one should be able to update this property at least cost when a new polynomial q e+1 is considered.

Definition 3 (Evaluation matrix)

. Consider e rational functions q 1 , . . . , q e in K(Y )[X], and e couples (s 1 , ω 1 ), . . . , (s e , ω e ) where each s i is taken in N m and each ω i is a monomial in X. Assume that none of the s i cancels any denominator of the q i . Consider the e × e matrix M with coefficients in K, defined by M ij = coeff(σ si (q j ), ω i ) where coeff(p, m) is the coefficient of the monomial m in the polynomial p. The matrix M is called the evaluation matrix of q 1 , . . . , q e w.r.t (s 1 , ω 1 ), . . . , (s e , ω e ).

Example 4. Recall q 1 = ax+2 a+1 and q 2 = z(a-1-ax) 1+a from Example 2. Consider (s 1 , ω 1 ) = (0, 1) and (s 2 , ω 2 ) = (1, 1). Thus, σ s1 (q 1 ) = 2, σ s1 (q 2 ) = -z, σ s2 (q 1 ) = x/2 + 1 and σ s2 (q 2 ) = -zx/2. Then, the evaluation matrix for (s

1 , ω 1 ), (s 2 , ω 2 ) is M = 2 -z 1 0
• If w 2 were the monomial x instead of 1, the evaluation matrix

would be M = 2 -z 1/2 -z/2
• In both cases, the matrix M is invertible.

The evaluation matrices computed in Algorithms incrementalFindDependence, findFirstDependence and findAllFirstIntegrals will be kept invertible, to ensure that some fractions are linearly independent (see Proposition 1 below).

Proposition 1. Keeping the same notations as in Definition 3, if the matrix M is invertible, then the rational functions q 1 , . . . , q e are linearly independent. Proof. Consider α = (α 1 , . . . , α e ) such that e i=1 α i q i = 0. One proves that α = 0. For each 1 ≤ j ≤ e, one has e i=1 α i σ sj (q i ) = 0, and consequently e i=1 α i coeff(σ sj (q i ), ω j ) = 0. This can be rewritten as M α = 0, which implies α = 0 since M is invertible.

By Proposition 1, each evaluation matrix of Example 4 proves that q 1 and q 2 are linearly independent. In some sense, an invertible evaluation matrix can be viewed as a certificate (in the the computational complexity theory terminology) that some fractions are linearly independent.

Proposition 2. Take the same notations as in Definition 3 and assume the evaluation matrix M is invertible. Consider a new rational function q e+1 . If the rational functions q 1 , . . . , q e+1 are linearly independent then one necessarily has q e+1 + e i=1 α i q i = 0 where α is the unique solution of M α = -b, with b = (. . . , coeff(σ sj (q e+1 ), ω j ), . . .) 1≤j≤e .

Proof. Since M is invertible and by Proposition 1, any linear dependence involves q e+1 with a nonzero coefficient, assumed to be 1. Assume q e+1 + e i=1 α i q i = 0.

Then for each j, one has e i=1 α i coeff(σ sj (q i ), ω j ) = -coeff(σ sj (q e+1 ), ω j ) which can be rewritten as M α = -b.

Example 5. Consider the q i of Example 2. Take D = [a+1], and (s 1 , ω 1 ) = (0, 1), and the 1 × 1 matrix M = (2) which is the evaluation matrix of q 1 w.r.t. (s 1 , ω 1 ). Apply algorithm incrementalFindDependence on Y , X, D, [q 1 ], [(s 1 , ω 1 )], M and q 2 . The vector b at Line 2 equals (-z) since q 2 = z(a -ax -1)/(a + 1) evaluates to -z when a = 0. Solving M α = -b yields α = (z/2). Then h = q 2 + (z/2)q 1 = az(2-x) 2(a+1) = 0. One then iterates the repeat until loop until h evaluates to a non zero polynomial. The value a = 0 is skipped, and the repeat until loop stops with s 2 = 1. Choosing the monomial w 2 = 1 yields the matrix M = 2 -z 1 0 • Proof (Algorithm incrementalFindDependence). Correctness. Assume the fractions q 1 , . . . , q e+1 are not linearly independent. Then Proposition 2 ensures that h must be zero and Case 1 is correct. It the q 1 , . . . , q e+1 are linearly independent, then necessarily h is non zero. Assume the repeat until loop terminates (see proof below). Since σ se+1 (h) = 0, a monomial ω e+1 such that coeff(σ se+1 (h), ω e+1 ) = 0 can be chosen. By construction, the matrix M is the evaluation matrix of q 1 , . . . , q e+1 w.r.t (s 1 , ω 1 ), . . . , (s e+1 , ω e+1 ). One just needs to prove that M is invertible to end the correctness proof. Assume M v = 0 with a non zero vector v = (α 1 , . . . , α e , β). If β = 0, then M α = 0 where α = (α 1 , . . . , α e ) and α = 0 since v = 0 and β = 0. Since M is invertible, α = 0 hence a contradiction. If β = 0, then one can assume β = 1. The e first lines of M v = 0 imply M α = -b where α is the vector computed in the algorithm. The last line of M v = 0 implies that l(α 1 , . . . , α e , 1) = 0, which implies coeff(σ se+1 (h, ω e+1 )) = 0 and contradicts the choice of ω e+1 .

Termination. One only needs to show that the repeat until loop terminates. This follows from Lemma 3 in the case of the single fraction q e+1 .

Improvement Using a LU-decomposition

In order to optimize the solving of M α = -b in incrementalFindDependence, one can require a LU-decomposition of the evaluation matrix M . The specification of Algorithm incrementalFindDependence can be slightly adapted by passing a LU-decomposition of M = LU (with L lower triangular with a diagonal of 1, and U upper triangular), and by returning a LU-decomposition of M = L U in Case 2. Note that a PLU-decomposition (where P is a permutation matrix) is not needed in our case as shown by Proposition 4 below. // E is a list of (s, ω), where s is an evaluation and w is a monomial

Finding the First Linear Dependence

The main advantage of Algorithm incrementalFindDependence is to limit the number of computations if one does not know in advance the number of rational fractions needed for having a linear dependence. Imagine the rational functions are provided by a finite iterator (i.e. a iterator that outputs a finite number of fractions), one can build the algorithm findFirstDependence which terminates on the first linear dependence it finds, or fails if no linear dependence exists. Example 6. Let us apply Algorithm findFirstDependence on Example 2. The first fraction q to be considered is ax+2 a+1 . The call to incrementalFindDependence returns the 1 × 1 matrix (2) and the couple (s 1 , ω 1 ) = (0, 1). One can check that q evaluated at a = 0 yields 2 and that (2) is indeed the evaluation matrix for the couple (s 1 , ω 1 ) = (0, 1).

Continue with the second iteration. One now considers q = z(a-ax-1) a+1 . Example 5 shows that the call to incrementalFindDependence returns the 2×2 invertible matrix 2 -z 1 0 and the couple (s 2 , ω 2 ) = (1, 1).

Finally, the third iteration makes a call to incrementalFindDependence. Line 2 builds the vector b = (1, 1). Line 3 solves M α = -b, yielding α = (-1, -1/z). Line 4 builds h = q 3 -q 1 -(1/z)q 2 which is zero, so incrementalFindDependence returns at Line 6. As a consequence, findFirstDependence detects that a linear dependence has been found at Line 11 and returns ((-1, -1/z), [q 1 , q 2 , q 3 ]) at Line 15.

Proof (Algorithm findFirstDependence). Termination. The algorithm obviously terminates since the number of loops is bounded by the number of elements output by the iterator J.

Correctness. One first proves the following loop invariant: M is invertible, and M is the evaluation matrix of Q w.r.t. E. The invariant is true when first entering the loop (even if the case is a bit degenerated since M , Q and E are all empty). Assume the invariant is true at some point. The call to incrementalFindDependence either detects a linear dependence α, or returns an evaluation matrix M and a couple (s, ω). In the first case, M , Q and E are left unmodified so the invariant remains true. In the second case, M , Q and E are modified to incorporate the new fraction q, and the invariant is still true thanks to the specification of incrementalFindDependence. The invariant is thus proven.

When the algorithm terminates, it returns either (α, Q) or FAIL. If it returns (α, Q), this means that the variable bool has changed to true at some point. Consequently a linear dependence α has been found, and the algorithm returns (α, Q). The dependence is necessarily the one with the smallest e because of the invariant (ensuring M is invertible) and Proposition 1. This proves the Case 1 of the specification.

If the algorithms returns FAIL, the iterator has been emptied, and the variable bool is still false. Consequently, all elements of the iterator have been stored in Q, and because of the invariant and Proposition 1, the elements of Q are linearly independent. This proves the Case 2 of the specification.

Remark 1. Please note that Algorithm findFirstDependence can be used with an infinite iterator (i.e. an iterator that never gets empty). However, Algorithm findFirstDependence becomes a semi-algorithm since it will find the first linear dependence if it exists, but will never terminates if no such dependence exists.

Complexity of the Linear Algebra

When findFirstDependence terminates, it has solved at most e square systems with increasing sizes from 1 to e. Assuming the solving of each system M α = b has a complexity of O(n ω ) [START_REF] Bunch | Triangular factorization and inversion by fast matrix multiplication[END_REF] arithmetic operations, where n is the size of the matrix and ω is the exponent of linear algebra, the total number of arithmetic operations for the linear algebra is O(e ω+1 ). If the LU-decomposition variant is used in Algorithm incrementalFindDependence, then the complexity of drops to O(e 3 ), since solving a triangular system of size n can be made in O(n2 ) arithmetic operations. As for the space complexity of algorithm findFirstDependence, it is O(e 2 ), whether using the LU-decomposition variant or not.

Application to Finding First Integrals

In this section, we look for first integrals for ODE systems. Roughly speaking, a first integral is an expression which is constant over time along any solution of the ODE system. This is a difficult problem and we will make several simplifying hypotheses. First, we work in the context of differential algebra, and assume that the ODE system is given by polynomial differential equations. Second, we will only look for polynomial first integrals.

Basic Differential Algebra

This section is mostly borrowed from [START_REF] Boulier | On the Regularity Property of Differential Polynomials Modulo Regular Differential Chains[END_REF] and [START_REF] Boulier | A normal form algorithm for regular differential chains[END_REF]. It has been simplified in the case of a single derivative. The reference books are [START_REF] Ritt | Differential Algebra[END_REF] and [START_REF] Kolchin | Differential Algebra and Algebraic Groups[END_REF]. A differential ring R is a ring endowed with an 2 abstract derivation δ i.e. a unary operation which satisfies the axioms δ(a

+ b) = δ(a) + δ(b) and δ(a b) = δ(a) b + aδ(b)
for all a, b ∈ R. This paper considers a differential polynomial ring R in n differential indeterminates u 1 , . . . , u n with coefficients in the field K. Moreover, we assume that K is a field of constants (i.e. δk = 0 for any k ∈ K). Letting U = {u 1 , . . . , u n }, one denotes R = K{U }, following Ritt and Kolchin. The derivation δ generates a monoid w.r.t. the composition operation. It is denoted: Θ = {δ i , i ∈ N} where N stands for the set of the nonnegative integers. The elements of Θ are the derivation operators. The monoid Θ acts multiplicatively on U , giving the infinite set ΘU of the derivatives.

If A is a finite subset of R, one denotes (A) the smallest ideal containing A w.r.t. the inclusion relation and [A] the smallest differential ideal containing A. Let A be an ideal and S = {s 1 , . . . , s t } be a finite subset of R, not containing zero. Then A : A differential polynomial q is said to be partially reduced w.r.t. p if it does not depend on any proper derivative of the leading derivative v of p. It is said to be reduced w.r.t. p if it is partially reduced w.r.t. p and deg(q, v) < deg(p, v). A set of differential polynomials of R \ K is said to be autoreduced if its elements are pairwise reduced. Autoreduced sets are necessarily finite [10, chapter I, section 9]. To each autoreduced set C, one may associate the set L = ld C of the leading derivatives of C and the set N = ΘU \ ΘL of the derivatives which are not derivatives of any element of L (the derivatives "under the stairs" defined by C).

S ∞ = {p ∈ R | ∃ a 1 , . . . , a t ∈ N, s a1 1 • • • s at t p ∈ A} is
In this paper, we need not recall the (rather technical) definition of differential regular chains (see [START_REF] Boulier | A normal form algorithm for regular differential chains[END_REF]Definition 3.1]). We only need to know that a differential regular chain C is a particular case of an autoreduced set and that membership to the ideal [C] : H ∞ C can be decided by means of normal form computations, as explained below.

Normal Form Modulo a Differential Regular Chain

All the results of this section are borrowed from [START_REF] Boulier | A normal form algorithm for regular differential chains[END_REF] and [START_REF] Boulier | On the Regularity Property of Differential Polynomials Modulo Regular Differential Chains[END_REF] The interest of [START_REF] Boulier | On the Regularity Property of Differential Polynomials Modulo Regular Differential Chains[END_REF] is that it provides a normal form algorithm which always succeeds (while [START_REF] Boulier | A normal form algorithm for regular differential chains[END_REF] provides an algorithm which fails when splittings occur).

Recall that the normal form algorithm relies on the computation of inverses of differential polynomials, defined below. Definition 5. Let f be a nonzero differential polynomial of R. An inverse of f modulo C is any fraction p/q of nonzero differential polynomials such that p ∈ K[N ∪ L] and q ∈ K[N ] and f p ≡ q mod A.

Normal Form Modulo a Decomposition

This subsection introduces a new definition which in practice is useful for performing computations modulo a radical ideal [Σ] expressed as an intersection of differential regular chains (i.e.

[Σ] = [C 1 ] : H ∞ C1 ∩ • • • ∩ [C f ] : H ∞ C f ).
Such a decomposition can be computed with the RosenfeldGroebner algorithm [START_REF] Boulier | The BLAD libraries[END_REF][START_REF] Boulier | Computing representations for radicals of finitely generated differential ideals[END_REF]. Definition 6 (Normal form modulo a decomposition). Let Σ be a set of differential polynomials, such that [Σ] is a proper ideal. Consider a decomposition of [Σ] into differential regular chains C 1 , . . . , C f for some ranking, that is differential regular chains satisfying

[Σ] = [C 1 ] : H ∞ C1 ∩ . . . ∩ [C f ] : H ∞ C f .

First Integrals in Differential Algebra

Definition 7 (First integral modulo an ideal). Let p be a differential polynomial and A be an ideal. One says p is a first integral modulo A if δp ∈ A.

For any ideal A, the set of the first integrals modulo A contains the ideal A. If A is a proper ideal, the inclusion is strict since any element of K is a first integral. In practice, the first integrals taken in K are obviously useless.

Example 7 (Pendulum). Take K = Q(m, l, g). Consider the ranking

• • • > l > ẍ > ÿ > l > ẋ > ẏ > l > x > y.
Recall the pendulum equations Σ in Equations [START_REF] Boulier | The BLAD libraries[END_REF]. Algorithm RosenfeldGroebner [START_REF] Boulier | Computing representations for radicals of finitely generated differential ideals[END_REF] shows that

[Σ] = [C 1 ] : H ∞ C1 ∩ [C2] : H ∞ C2
where C 1 and C 2 are given by:

-C 1 = [ λ = -3 ẏgm l 2 , ẏ2 = --λy 2 l 2 +λl 4 -y 3 gm+ygml 2 ml 2 , x 2 = -y 2 + l 2 ] ; -C 2 = [λ = -ygm l 2 , x = 0, y 2 = l 2 ]
. Remark that the differential regular chain C 2 corresponds to a degenerate case, where the pendulum is vertical since x = 0. Further computations show that

NF(δ( m 2 ( ẋ2 + ẏ2 ) -mg y), [C 1 , C 2 ]) = [0, 0]), proving that p = m 2 ( ẋ2 + ẏ2 ) -mg y is a first integral modulo [Σ]. Remark that x 2 + y 2 is also a first integral. This is immediate since δ(x 2 + y 2 ) = δ(x 2 + y 2 -1) ∈ [Σ].
Definition 7 is new to our knowledge. It is expressed in a differential algebra context. The end of Section 4.4 makes the link with the definition of first integral in a analysis context, through analytic solutions using results from [START_REF] Boulier | A computer scientist point of view on Hilbert's differential theorem of zeros[END_REF]. Definition 8 (Formal power solution of an ideal). Consider a n-uple ū = (ū 1 (t),. . . ,ū n (t)) of formal power series in t over K. For any differential polynomial, one defines p(ū) as the formal power series in t obtained by replacing each u i by ūi and interpreting the derivation δ as the usual derivation on formal power series. The n-uple ū = (ū 1 , . . . , ūn ) is called a solution of an ideal A if p(ū) = 0 for all p ∈ A.

Lemma 4. Take a differential polynomial p and n-uple ū = (ū 1 (t),. . . ,ū n (t)) of formal power series. Then (δp)(ū) = dp(ū) dt . If p is a first integral modulo an ideal A and ū is a solution of A, then the formal power series p(ū) satisfies dp(ū) dt = 0.

Proof. Since δ is a derivation, (δp)(ū) = dp(ū) dt is proved if one proves it when p equals any u i . Assume that p = u i . Then (δp)

(ū i ) = (δu i )(ū i ) = dp(ū)
dt . Assume p is a first integral modulo A and ū is a solution of A. Then δp ∈ A and (δp)(ū) = 0. Using (δp)(ū) = dp (ū) dt , one has dp(ū) dt = 0 Take a system of differential polynomials Σ. By [5, Theorem and definition 3], a differential polynomial p in R vanishes on all analytic solutions (over some open set with coordinates in the algebraic closure of the base field) of Σ if and only if p ∈ [Σ]. Applying this statement to p = δq for some first integral q w.r.t. [Σ], then δq vanishes on all analytic solutions of [Σ] so p is a first integral in the context of analysis, if one excludes the non analytic solutions.

Algorithm findAllFirstIntegrals

In this section, one considers a proper ideal A given as A

= [C 1 ]:H ∞ C1 ∩• • •∩[C f ]: H ∞ C f where the C i are differential regular chains. Denote M = [C 1 , . . . , C f ]
. Take a first integral modulo A of the form p = α i µ i , where the α i 's are in K and the µ i 's are monomials in the derivatives. Computing on lists componentwise, we have 0 = NF(δp, M ) = α i NF(δµ i , M ). Consequently, a candidate α i µ i is a first integral modulo A if and only if the NF(δµ i , M ) are linearly dependent over K.

Since the µ i have no denominators, every irreducible factor of the denominator of any NF(w i , C j ) necessarily divides (by Proposition 5) the denominator of the inverse of a separant or initial of C j . As a consequence, the algorithms presented in the previous sections can be applied since we can precompute factors which should not be cancelled.

Algorithm findAllFirstIntegrals is very close to Algorithm findFirstDependence and only requires a few remarks. Instead of stopping at the first found dependence, it continues until the iterator has been emptied and stores all first dependences encountered. It starts by precomputing a safe set D for avoiding cancelling the denominators of any NF(δw i , C j ). The algorithm introduces some dummy When increasing the degree bound, one finds more and more spurious first integrals like ẏ y 2 + ẋ x y (which is in the ideal) or some powers of the first integral -2g y + ẋ2 + ẏ2 .

Example 9 (Lotka-Volterra equations).

C ẋ(t) = a x(t) -b x(t) y(t) ẏ(t) = -c y(t) + d x(t) y(t)

x(t) u(t) = ẋ(t)

y(t) v(t) = ẏ(t) (2)
Take K = Q(a, b, c, d) and the ranking

• • • > u > v > ẋ > ẏ > u > v > x > y.
One can show that C is a differential regular chain for the chosen ranking.

The two leftmost equations of C corresponds to the classical Lotka-Volterra equations, and the two rightmost ones encode the logarithms of x(t) and y(t) in a polynomial way. A call to findAllFirstIntegrals with the monomials of degree at most 1 built over x, y, u, v yields [1, -av d -cu d + by d + x] which corresponds to the usual first integral -a ln(y(t)) -c ln(x(t) + by(t) + dx(t).

Remark 2. The choice of the degree bounds and the candidate monomials in the first integrals is left to the user, through the use of the iterator J. This makes our algorithm very flexible especially if the user has some extra knowledge on the first integrals or is looking for specific ones. Finally, this avoids the difficult problem of estimating degree bounds, which can be quite high. Indeed, the simple system ẋ = x, ẏ = -ny, where n is any positive integer, admits x n y as a first integral, which is minimal in terms of degrees.

Complexity

The complexity for the linear algebra part of Algorithm findAllFirstIntegrals is the same as for Algorithm findFirstDependence: it is O(e 3 ), where e is the cardinal of the iterator J, if one uses the LU-decomposition variant. However, e can be quite large in practice. For example, if one considers the monomials of total degree at most d involving s derivatives, then e is equal to s+d s .

  page 101] over the ring K(y 1 , . . . , y m )[y m+1 ], one has p = d+1 i=1 p(y 1 , . . . , y m , s i )

2 For each 1 3 D 4 I 5 S 6 while true do 7 For each 1 ≤ 8 / 9

 2134567189 and B involve z since one performs linear Input: Two lists of variables Y and X Input: A vector q = (q1, . . . , qe) of e rational fractions in K(Y )[X] Output: A basis of ker Φq 1 begin ≤ i ≤ e, denote the reduced fraction qi as fi/gi with fi ∈ K[Y, X], gi ∈ K[Y ]; ← [g1, . . . , ge]; ← newTupleIteratorNotCancelling(Y, D); ← [getNewTuple(I)]; i ≤ e, denote vi the vector (σs 1 (qi), σs 2 (qi), . . . , σs r (qi)) where S = [s1, s2, . . . , sr]; / each vi is a vector of r = | S| elements of K[X], obtained by evaluating qi on all points of S Denote v = (v1, . . . , ve); 10 Compute a basis B of the kernel of Φv using linear algebra; 11 // if ker Φq ⊃ ker Φv, one returns B 12 if e i=1 biqi = 0 for all b = (b1, . . . , be) ∈ B then return B; 13 ; 14 Append to S the new evaluation getNewTuple(I);

3 D 4 I 5 M 6 while true do 7 y 0 8 // evaluate the vector q at y 0 9 q

 34567089 ← [g1, . . . , ge]; ← newTupleIteratorNotCancelling(Y, D); ← the 0 × e matrix; ← getNewTuple(I); ← σ y 0 (q); 10 build L = [ω1, . . . , ω l ] the list of monomials involved in q; 11 build the l × e matrix N = (Nij) where Nij is the coefficient of qj in the monomial ωi ; 12 M ← M N ; 13 Compute a basis B of the right kernel of M ; 14 if e i=1 biqi = 0 for all b = (b1, . . . , be) ∈ B then return B; 15 ;

Input: 2 b 4 h 5 if h = 0 then 6 return 8 I 10 ; 11 choose a monomial

 245681011 Two lists of variables Y and X, and a list D of elements of K[Y ] Input: A list Q = [q1, . . . , qe] of e rational functions in K(Y )[X] Input: A list E = [(s1, ω1), . . . , (se, ωe)], with si ∈ N m and ωi a monomial in X Input: M an invertible eval. matrix of q1, . . . , qe w.r.t. (s1, ω1), . . . , (se, ωe) Input: qe+1 a rational function Assumption: denote qi = fi/gi with fi ∈ K[X, Y ] and gi ∈ K[Y ] for 1 ≤ i ≤ e + 1. Each gi divides a power of elements of D. Moreover, σs i (dj) = 0 for any si and dj ∈ D. Output: Case 1 : (α1, . . . , αe, 1) s.t. qe+1 + e i=1 αiqi = 0 Output: Case 2 : M , (se+1, ωe+1) such that M is the evaluation matrix of q1, . . . , qe+1 w.r.t. (s1, ω1), . . . , (se+1, ωe+1), with M invertible 1 begin ← (. . . , coeff(σs j (qe+1), ωj), . . .) 1≤j≤e ; 3 solve M α = -b in the unknown α = (α1, . . . , αe); ← qe+1 + e i=1 αiqi; (α1, . . . , αe, 1) ; // Case 1: a linear dependence has been found 7 else ← newTupleIteratorNotCancelling(Y, D); 9 repeat se+1 ← getNewTuple(I) until σs e+1 (h) = 0; ωe+1 such that coeff(σs e+1 (h), ωe+1) = 0 ; 12 l ← (coeff(σs e+1 (q1), ωe+1) • • • coeff(σs e+1 (qe+1), ωe+1)); (se+1, ωe+1) // Case 2: q1, . . . , qe+1 are linearly independent Algorithm 3: incrementalFindDependence

Proposition 3 ( 2 M 3 Q

 323 Solving α). Solving M α = -b is equivalent to solving the two triangular systems Ly = -b, then U α = y. Proposition 4 (The LU-decomposition of M ). The LU-decomposition of M can be obtained by: solve γU = l 1:e (in the variable γ), where l 1:e denotes the e first components of the line vector l computed in findFirstDependence Input: Two lists of variables Y and X Input: A list D of elements of K[Y ] Input: A finite iterator J which outputs rational functions q1, q2, . . . in K(Y )[X] Assumption: Denote the reduced fraction qi = fi/gi with fi ∈ K[X, Y ] and gi ∈ K[Y ]. Each gi divides a power of elements of D. Output: Case 1: a shortest linear dependence i.e. a vector (α1, . . . , αe) and a list [q1, . . . , qe] with e i=1 αiqi = 0 and e the smallest possible. Output: Case 2: FAIL if no linear dependence exists 1 begin ← the 0 × 0 matrix ; // M is an evaluation matrix ← the empty list [] ; // Q is the list of the qi 4

-

  solve Ly = -b (in the variable y) l e+1 + γy

  called the saturation of A by the multiplicative family generated by S. The saturation of a (differential) ideal is a (differential) ideal [10, chapter I, Corollary to Lemma 1]. Fix a ranking, i.e. a total ordering over ΘU satisfying some properties [10, chapter I, section 8]. Consider some differential polynomial p / ∈ K. The highest derivative v w.r.t. the ranking such that deg(p, v) > 0 is called the leading derivative of p. It is denoted ld p. The leading coefficient of p w.r.t. v is called the initial of p. The differential polynomial ∂p/∂v is called the separant of p. If C is a finite subset of R \ K then I C denotes its set of initials, S C denotes its set of separants and H C = I C ∪ S C .

Definition 4 .Proposition 5 .

 45 . Let C be a regular differential chain of R, defining a differential ideal A = [C]:H ∞ C . Let L = ld C and N = ΘU \ ΘL. The normal form of a rational differential fraction is introduced in [2, Definition 5.1 and Proposition 5.2], recalled below. Let a/b be a rational differential fraction with b regular modulo A.A normal form of a/b modulo C is any rational differential fraction f /g such that 1 f is reduced with respect to C, 2 g belongs to K[N ] (and is thus regular modulo A), 3 a/b and f /g are equivalent modulo A (in the sense that a g -b f ∈ A). Let a/b be a rational differential fraction, with b regular modulo A. The normal form f /g of a/b exists and is unique. The normal form is a K-linear operation. Moreover 4 a belongs to A if and only if its normal form is zero, 5 f /g is a canonical representative of the residue class of a/b in the total fraction ring of R/A, 6 each irreducible factor of g divides the denominator of an inverse of b, or of some initial or separant of C .

Proposition 6 .

 6 For any differential fraction a/b with b regular modulo each [C i ] : H ∞ Ci , one defines the normal form of a/b w.r.t. to the list [C 1 , . . . , C f ] by the list [NF(a/b, C 1 ), . . . , NF(a/b, C f )] . It is simply denoted by NF(a/b, [C 1 , . . . , C f ]). Since it is an open problem to compute a canonical (e.g. minimal) decomposition of the radical of a differential ideal, the normal form of Definition 6 depends on the decomposition and not on the ideal. With the same notations as in Definition 6, for any polynomial p ∈ R, one has p ∈ [Σ] ⇐⇒ NF(p, [C 1 , . . . , C f ]) = [0, . . . , 0]. Moreover, the normal form modulo a decomposition is K-linear.

f 1 begin 2 result 3 for i = 1 to f do 4 I 6 Y 7 X 8 M 9 E 10 while the iterator J is not empty do 11 µ

 23467891011 Input: A list of differential regular chains M = [C1, . . . , C f ] Input: A finite iterator J which outputs differential monomials Output: A K-basis of the linear combinations of momomials of J which arefirst integrals w.r.t. [C1] : H ∞ C 1 ∩ • • • ∩ [C f ] : H ∞ C ←the empty list [] ; D ← the empty list [] ; ← the inverses of the initials and separants of Ci modulo Ci ; 5 append to D the denominators of the elements of I ; ← the list of derivatives occurring in D ; ← the list of dummy variables [d1, . . . , d f ] ; ← the 0 × 0 matrix ; Q ← the empty list [] ; ← the empty list [] ; F ← the empty list []; ← getNewDerivativeMonomial(J); 12 q ← f i=1 diNF(δµ, Ci); 13 append to X the variables of the numerator of q, which are neither in X nor Y ; 14 r ← incrementalFindDependence(Y, X, D, Q, E, M, q); 15 if r is a linear dependence then 16 // A new first integral has been found 17 (α1, . . . , αe, 1) ← r ; 18 append α1µ1 + • • • + αeµe + w to result, where F = (µ1, . . . , µe); 19 else 20 append µ to the end of F ; append q to the end of Q ; 21 M , (s, ω) ← r ; M ← M ; append (s, ω) to the end of E; 22 return result ; Algorithm 5: findAllFirstIntegrals variables d 1 , . . . , d f for storing the normal form NF(δµ, M ), which is by definition a list, as the polynomial d 1 NF(δµ, C 1 )+• • •+d f NF(δµ, C f ). This alternative storage allows us to directly reuse Algorithm incrementalFindDependence which expects a polynomial. Example 8 (Pendulum). Take the same notations as in Example 7. Take an iterator J enumerating the monomials 1, y, x, ẏ, ẋ, y 2 , xy, y 2 , ẏy, ẏx, ẏ2 , ẋy, ẋx, ẋ ẏ and ẋ2 . Then Algorithm findAllFirstIntegrals returns the list [1, x 2 + y 2 , ẏy + ẋx, -2g y + ẋ2 + ẏ2 ]. Note the presence of ẏy + ẋx which is in the ideal since it is the derivative of (x 2 + y 2 -1)/2. The intermediate computations are too big be displayed here. As an illustration, one gives the normal forms of δy, δ( ẏ2 ), δ( ẋ ẏ) and δ( ẋ2 ) modulo [C 1 , C 2 ] which are respectively [ ẏ, 0] , 2(λy ẏ + mg ẏ) m , 0 , x ẏ(λ(2y 2 -l 2 ) + mgy) m(y 2 -l 2 ) , 0 and -2λy ẏ m , 0 .

In the general setting, differential ring are endowed with finitely many derivations


Joseph Gril
email: joseph.gril@univ-montp2.fr

Bertrand Marcon

Fabrice Brémand

Linda Cocchi

Paolo Dionisi-Vici

Jean-Christophe C Dupré

Cécilia Gauvin

Giacomo Goli

Franck Hesser

Delphine Jullien

P Mazzanti

E Ravaud

M Togni

V Valle

L Uzielli

The Mona Lisa Project: an update on the progress of measurement, monitoring, modelisation and simulation

Introduction

Since 2004 an international research group of wood technologists has been given by the Louvre museum the task of analysing the mechanical situation of the wooden panel on which Leonardo da Vinci painted his "Mona Lisa" (Fig. 1a), possibly between 1503 and 1506. The general purpose of such study was to evaluate the influences that could possibly derive from environmental fluctuations in the showcase where the painting is exhibited as well as outside the showcase for occasional checks, and develop measurements and models to improve its conservation conditions. To acquire data on the mechanical behaviour of the panel, and to feed and calibrate appropriate simulation models, the team has set up a continuous monitoring by means of automatic equipment (Fig. 1c). The "Mona Lisa" is painted on a poplar panel (Populus alba L.) ~790 x 530 mm, ~13 mm thick, which is inserted in an oak frame (called châssis-cadre), and is slightly forced against it by means of four crossbars, holding it flatter than it would be if unconstrained. In turn the châssis-cadre is inserted in a wooden gilded frame (Fig. 1c). In 2006 a book was published [START_REF] Mohen | Mona Lisa: Inside the Painting[END_REF] offering a wealth of scientific studies and researches concerning the "Mona Lisa". Also in 2006 a report about finite elements model was presented [START_REF] Gril | Mona Lisa saved by the Griffith theory: assessing the crack propagation risk in the wooden support of a panel painting[END_REF] and in 2011 a scientific article on the modelisation was published [START_REF] Dureisseix | A partitioning strategy for the coupled numerical simulation of hygromechanical wood structures of Cultural Heritage[END_REF].

Development of measurement and monitoring techniques

Initially the panel shape has been determined through manual shape measurement by means of a mechanical comparator and a reference bar. This technique is slow, and allows surveying only a limited number of points. Then optical techniques shadow moiré and fringe pattern profilometry -FPP-have been used to measure [START_REF] Brémand | Mechanical structure analysis by Digital Image Correlation and Fringe Pattern Profilometry -Proceedings of COST Actions IE0601 Non-destructive techniques to study Wooden Cultural Heritage Objects[END_REF], on both front and rear faces, relief (Fig. 1b) and out of plane deformation field. To obtain the 3D surface displacements in some points of the panel stereo-optical tracking [START_REF] Brémand | Mechanical structure analysis by Digital Image Correlation and Fringe Pattern Profilometry -Proceedings of COST Actions IE0601 Non-destructive techniques to study Wooden Cultural Heritage Objects[END_REF] has been used (following image contrast and craquelure pattern). Accurate and reliable data about the forces exerted between panel and crossbars are obtained thanks to self-designed equipment including four sub-miniature load cells incorporated into the crossbar thickness (next to the panel's four corners). The deflection variations along time are measured by means of three deformation transducers, located inside a thin reference aluminium profiled crossbar (carrying data-loggers, transducers and instrumentation), fixed on the châssis-cadre and providing records of both (a) transversal deflection at the panel centre with reference to the lateral edges, and (b) longitudinal deflection with reference to the châssis-cadre, see Fig. 1c. The contact forces between panel and châssis-cadre have been localized and estimated on the basis of local contact pressures through a pressure-sensitive foil as described in [START_REF] Goli | Locating contact areas and estimating contact forces between the "Mona Lisa" wooden panel and its frame[END_REF].

Modelisation and simulation strategies

The results from the above mentioned measurements are being processed to be included into a 3D numerical finite elements model to simulate the panel behaviour under environmental fluctuations, see Fig. 1d. We focus herein on a numerical strategy, taking into account the panel specificities, including its shape and thickness, its sawing pattern (the elastic behaviour is orthotropic, roughly with a cylindrical symmetry), the boundary conditions imposed by the châssis-cadre (contact area and forces [START_REF] Goli | Locating contact areas and estimating contact forces between the "Mona Lisa" wooden panel and its frame[END_REF]), the crack at the upper edge of the panel and the remedial wooden butterfly, the unpainted back face of the panel responsible for the out-ofplane movements caused by moisture gradients across the panel's thickness and resulting moisture-induced expansion during relative humidity fluctuations. Such elements -mostly based on measurements-are essential to propose a model that can simulate with accuracy the Mona Lisa behaviour and are exposed in detail in [START_REF] Dureisseix | A partitioning strategy for the coupled numerical simulation of hygromechanical wood structures of Cultural Heritage[END_REF] and [START_REF] Marcon | Hygromécanique des panneaux en bois et conservation du patrimoine culturel[END_REF]. The aim of this strategy is to improve Mona Lisa's preventive conservation by virtually testing and predicting, via numerical simulation, if the masterpiece will be safe under various scenarios.

Fig. 1 (

 1 Fig. 1 (a) The painted face and the back face; (b) Shape of the panel obtained by FPP; (c) Completed assembly of the Mona Lisa: panel, crossbars, châssis-cadre, gilded frame and the monitoring systems; (d) Mona Lisa numerical finite elements model.

Acknowledgments

The authors wish to thank the Louvre museum and V. Delieuvin, conservator at the Department of Paintings, in charge of the Mona Lisa; they also acknowledge P. Mandron and D. Jaunard, independent restorers, for their essential contribution and fruitful discussions.


Elodie Mattio

Fabien Robert-Peillard

Catherine Branger

Kinga Puzio

André Margaillan

Christophe Brach-Papa

Joël Knoery

Jean Luc Boudenne

Bruno Coulomb
email: bruno.coulomb@univ-amu.fr

Jean-Luc Boudenne

3D-printed flow system for determination of lead in natural waters

Keywords: Lead determination, 3D-printed MPFS system, stereolithography, natural waters

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Lead is considered as one of the most toxic heavy metals [START_REF] Garnier | Toxicité du plomb et de ses dérivés[END_REF] in the light of its environmental [START_REF] Dikilitas | Chapter 3 -Effect of Lead on Plant and Human DNA Damages and Its Impact on the Environment[END_REF] and health [START_REF] Juberg | Lead and human health: An update[END_REF][START_REF] Johnson | The genetic effects of environmental lead[END_REF] impacts; it may cause irreversible neurological effects and digestive and kidney malfunctions. Its presence in the environment and more particularly in water can be mainly explained by anthropogenic sources like paints, arms and electronic [START_REF] Pecht | The impact of lead-free legislation exemptions on the electronics industry[END_REF] industries products. For these reasons, the World Health Organization has identified lead as one of the 10 chemicals of major public health concern [START_REF]World Health Organization | Lead poisoning and health[END_REF] and has recommended a guideline value of 10 μg.L -1 of lead in drinking water. It is therefore necessary to quantify lead by a rapid and efficient method to avoid toxic consumption.

Many methods are already available for lead determination such as spectrophotometry [START_REF] Rahman | Determination of lead in solution by solid phase extraction, elution, and spectrophotometric detection using 4-(2-pyridylazo)-resorcinol[END_REF],

voltammetry [START_REF] Fang | Determination of trace lead and cadmium using stripping voltammetry in fluidic microchip integrated with screen-printed carbon electrodes[END_REF][START_REF] Abate | Complexation of Cd (II) and Pb (II) with humic acids studied by anodic stripping voltammetry using differential equilibrium functions and discrete site models[END_REF], graphic furnace atomic absorption spectrometry [START_REF] Bruland | Analysis of seawater for dissolved cadmium, copper and lead: An intercomparison of voltammetric and atomic absorption methods[END_REF], or inductively coupled plasma spectroscopy [START_REF] Ho | Determination of trace metals in seawater by an automated flow injection ion chromatograph pretreatment system with ICPMS[END_REF][START_REF] Zougagh | Automatic on line preconcentration and determination of lead in water by ICP-AES using a TSmicrocolumn[END_REF], but they require costly and sophisticated devices and do not allow real-time and on-site measurements. In this context, mesofluidic and microfluidic systems can help to meet these needs thanks to their intrinsic advantages: miniaturization and low energy consumption, decreased reagents consumption and waste generation. Flow analysis offers many opportunities [START_REF] Cerdà | Flow techniques in water analysis[END_REF][START_REF] Trojanowicz | Advances in Flow Analysis[END_REF] and Flow Injection Analysis (FIA), Sequential Injection Analysis (SIA)

or MultiSyringe Flow Injection Analysis (MSFIA) based systems have been successfully applied for the determination of a wide range of substances in environmental matrices. Among flow analysis system, pulse flow systems using solenoid diaphragm micropumps (MPFS: Multi Pumping Flow System) present high flexibility, easy configuration and low cost. Another advantage is that the signal peaks are higher compared to other flow techniques, due to turbulences created by pump diaphragm strokes, which improve mixing between reagents and sample [START_REF] Santos | Multipumping flow systems: The potential of simplicity[END_REF]. Recently, such MPFS systems were developed for phosphorous [START_REF] González | Low cost analyzer for the determination of phosphorus based on open-source hardware and pulsed flows[END_REF] or boron [START_REF] González | Multi-pumping flow system for the determination of boron in eye drops, drinking water and ocean water[END_REF] determination in aqueous samples.

Several units can be combined to create a flow system in full compliance with the analytical needs: solid phase extraction for pre-concentration, photo-oxidation or digestion, membrane or membrane-less separation, detection and many others. To increase versatility and ease of fabrication of flow systems, 3D printing is increasingly used in flow analysis [START_REF] Bhattacharjee | The upcoming 3D-printing revolution in microfluidics[END_REF][START_REF] Au | 3D-Printed Microfluidics[END_REF]. This printing technology can be divided into three main categories: Fused Deposition Modeling (FDM) where a thermoplastic material is heated and extruded from a XYZ positionable nozzle, Multi Jet Modeling (MJM) which is based on an inkjet head to deposit liquid photopolymer (plastic resin or casting wax) layer by layer and finally stereolithography (SL) [START_REF] Au | Mail-order microfluidics: evaluation of stereolithography for the production of microfluidic devices[END_REF] invented in the 1980s and which is a layer manufacturing process with liquid materials [START_REF] Bártolo | Stereolithography: Materials, Processes and Applications[END_REF] as MJM. SL is based on a mobile platform which dives in a resin tray, where a laser polymerizes the resin layer by layer. The most frequently used material is poly(methyl methacrylate) (PMMA) which allows to decrease fabrication cost and improve resolution [START_REF] Polzin | Characterization and evaluation of a PMMA-based 3D printing process[END_REF]. 3D-printing technologies emergence enables creation of new type of units [START_REF] Frizzarin | A portable multi-syringe flow system for spectrofluorimetric determination of iodide in seawater[END_REF][START_REF] Cerdà | Chips: How to build and implement fluidic devices in flow based systems[END_REF] for microfluidic systems which can be imagined and combined according to the analytical needs.

In this paper, the development of a novel 3D-printed multi-pumping flow system for the determination of lead in natural waters is presented. The system is composed of three 3Dprinted units: a resin column for lead solid phase extraction, a mixing coil and a classical flow cell for UV-Spectroscopy. The 3 modules are connected to one another by a screw system to limit tubings between the three units. Lead solid phase extraction is based on a commercial crown-ethers resin (TrisKem Pb resin) and the chromogenic reagent used for the spectrophotometric detection of lead is 4-(2-pyridylazo)resorcinol (PAR). The system was applied to the determination of lead in real river water samples.

Materials and methods

2.1.Reagents and samples

All chemicals used were of analytical grade and used without further purification. Freshwater samples were collected at purposely chosen points in a coastal river "The Arc" in the south-east of France, near industrial or urban effluents discharge points. All samples were UV-photo-oxidized with a 254 nm low-pressure mercury lamp (UVP PenRay, USA) during 30 minutes before analysis allowing liberation of lead linked to natural organic ligands, inorganic ligands or even to anthropogenic organic ligands.

2.2.Apparatus

Flow system

The three units of the MPFS system (Fig. 1) were designed with Rhinoceros® 5.0 3D software (Robert McNeel & Associates Europe, Spain), then printed with a poly(methyl methacrylate) resin on the 3D printer Form1+ (Formlabs, USA). The first unit consisted of a resin column, tightly closed by a screw piece with three entry channels for injection of sample and reagents.

The resin was packed in this unit between two layers of glass wool. To optimize mixing of eluate and chromogenic reagent used for lead detection (PAR), a second unit composed of a serpentine mixing coil (1.5 mm internal diameter, 50 cm length) and of a connection for PAR inlet was added to the system. Finally, the detection step was performed with a classical spectroscopic flow-cell with a 5-cm optical pathlength. Sample and reagents were introduced inside the system by means of four solenoid micro-pumps (Bio-ChemValve Inc., USA) that had a stroke volume of 20 µL and a high frequency of 250 cycles/min. These pumps were computercontrolled by a MCFIA/MPFS system (Sciware, Spain) with eight digital 12V output channels.

For the detection, two FC-UV600 optical fibres (Ocean Optics, USA) were connected at the ends of the optical pathlength, and isolated from the reaction mixture with two tailor-made quartz discs, to guide the light from the source to the spectrophotometric detector. The radiation of the halogen bulb of a DH-2000 UV-Vis light source (Ocean Optics, USA) was transmitted to a USB2000 miniature spectrometer detector (Ocean Optics, USA). The whole system was controlled by AutoAnalysis 5.0 software (Sciware, Spain).

Metal analysis

Graphite furnace atomic absorption spectrometry (GF-AAS) was used to optimize the extraction/elution steps on TrisKem Pb resin and validate lead amounts in real samples (after filtration on a 0.45 µm polyethersulfone membrane). The measurements were carried out on a Thermo Scientific ICE3500 (USA) atomic absorption spectrometer equipped with a lead hollow-cathode lamp operated at 10 mA (wavelength of 217 nm). Argon flow was 0.2 L.min -1 except during atomisation step (no flow). The furnace settings were: drying at 110 °C, ramp for 9 s, hold for 35 s; cracking at 800 °C, ramp for 5 s, hold for 20 s; atomising at 1200 °C, ramp for 1 s and 3 s hold; cleaning at 2500 °C, no ramp and 3 s hold.

Interfering cations were determined by inductively coupled plasma-atomic emission spectrometry (ICP-AES) with a Jobin YVON JY2000 Ultratrace spectrometer, equipped with a CMA spray chamber and a Meinhard R50-C1 glass nebuliser. Determinations were performed with the following settings: power 1000W, pump speed 20 mL.min -1 , plasma flow rate 12 L.min -1 , coating gas flow rate 0.2 L.min -1 , nebuliser flow rate 0.9 L.min -1 and nebuliser pressure 2.08 bar.

2.3.Flow procedure

The MPFS system was operated according to the procedure in 6 steps given in Table 1: in step 1, 9 mL of nitric acid solution (0.05 mol.L -1 ) was used for conditioning the Pb resin and washing the system before analysis. Then 50 mL of acidified sample were introduced in the system (step 2). In step 3, the resin was washed again with nitric acid solution (0.05 mol.L -1 ) to eliminate potential interfering species extracted by the resin. 5 mL of ammonium oxalate solution (0.025 mol.L -1 ) and PAR (0.01 mmol.L -1 ) were simultaneously pumped through the resin to eluate extracted lead and at the inlet of the mixing coil, respectively (steps 4-5).

Measurements were based on the peak height. The analytical signal was recorded at 520 nm when the elution started. Absorbance spectra were acquired every 0.5 s, with an integration time of 55 ms and an average of 3 spectra.

Results and discussion

3.1.Extraction and washing steps 3.1.1. Influence of concentration of nitric acid

The developed system was based on the spectrophotometric determination of lead in the presence of PAR reagent. However, it is well known that lead determination by this simple method often suffers from interferences of cadmium, copper and zinc [START_REF] Dagnall | Determination of lead with 4-(2-pyridylazo)resorcinol[END_REF]. Therefore, the concentration of these interfering metals should be reduced before PAR detection to allow selective determination of lead. The results in Fig. 2 showed that more than 90% of lead was extracted with HNO3 concentration above 0.02 mol.L -1 . The best lead extraction (97%) was obtained with HNO3 0.1 mol.L -1 .

Lead determination was thus based on extraction properties of the resin

Concerning the other metal cations, iron was also partially extracted (20-50%, with the highest extraction at 0.1 mol.L -1 ). Aluminium, copper and zinc were not retained by the resin whatever the acid concentration. For HNO3 concentrations of 0.02 and 0.05 mol.L -1 , a very few amount (≤ 5%) of cobalt, chromium, and cadmium was extracted. It can therefore be concluded that the resin has an excellent selectivity for lead over other metals, except for iron for which extraction is dependent on HNO3 concentration. A nitric acid concentration of 0.05 mol.L -1 appeared to be sufficient to extract lead quantitatively and to limit the extraction rate of iron below 20%.

After metal extraction on the resin, a washing step can improve the selectivity by removing potentially interfering metals without eluting the analyte of interest. This washing step has also been studied with different HNO3 concentrations. Based on the previous results, only iron and lead were studied. The washing solutions were collected and analysed by ICP-AES. As can be seen in Fig. 3 the amount of both iron and lead washed out decreased with increasing HNO3 concentration. With HNO3 0.02 mol.L -1 , all the iron retained within the resin was eliminated, but an important lead elution rate of 42% was also observed. A nitric acid concentration of 0.05 mol.L -1 seems to be the best compromise between iron elimination (75%) and limited lead elution rate (about 10%).

To summarize 0.05 mol.L -1 HNO3 may be used for both extraction and washing steps, thereby limiting the number of pumps to be used in the system.

Sample flow-rate

Sample flow rate was studied from 2 to 5 mL.min -1 for a sample volume of 30 mL. Samples were acidified with 0.05 mol.L -1 HNO 3 before extraction. As displayed in supplementary material (Fig. S1), results showed an important increase of lead extraction from 3 mL.min -1 to 5 mL.min -1 . At the highest flow rate (5 mL.min -1 ) the extraction percentage of lead reached 94 %.

To avoid overpressure problem of solenoid pumps or clogging of resin column, the column inner diameter and height has been set respectively at 7.8 mm and 4.5 mm, which are higher than in a traditional packed column. None of the previously mentioned problems were observed during analysis. Although solenoid pumps are not usually suitable for resin solid phase extraction, it seems that high sample flow rates enable fluidizing of the resin bed and thus good lead extraction. Therefore, a 5mL.min -1 flow rate was chosen for the system.

Elution step

Ammonium oxalate is recommended by the resin manufacturer for lead elution. The lead elution has thus been studied at different concentrations of ammonium oxalate. The elution flow rate has been adjusted at 4 mL.min -1 , according to the resin manufacturer instructions.

The results shown in Fig. 4 highlighted that elution efficiency increased with ammonium oxalate concentration up to 0.05 mol.L -1 .

For concentration of ammonium oxalate between 0.025 and 0.1 mol.L -1 , the results obtained were not significantly different. To limit reagent consumption, the concentration of ammonium oxalate was adjusted at 0.025 mol.L -1 for further experiments.

Detection

Lead elution profile has been studied by collecting small fractions (0.5 mL) of ammonium oxalate used for elution with a flow rate of 5 mL.min -1 . As can be seen on Fig. S2 (Supplementary material), elution was complete within one minute (elution rate of 99% of total extracted lead after 5 mL).

Lead was detected by spectrophotometry using a PAR solution pumped simultaneously and mixed with ammonium oxalate in the mixing coil. PAR solution was buffered with borate solution at pH=12 in order to obtain a mixture at pH=9 (optimised for UV-Vis detection) when mixed with ammonium oxalate which pH is around 6.6.

PAR reagent and ammonium oxalate flow rates were studied from 1 to 5 mL.min -1 .

Optimization was carried out with a standard lead solution of 100 µg.L -1 . PAR flow rate was optimized using a fixed ammonium oxalate flow rate of 5 mL.min -1 and inversely. Fig. 5 showed that the best absorbance was obtained with the maximum flow rate of the pumps useable for the two reagents (5 mL.min -1 ). The absorbance increased regularly with increasing PAR flow rate. The absorbance values increased slightly with an ammonium oxalate flow rate between 1 and 4 mL.min -1 and, strongly increased for a flow rate of 5 mL.min -1 . This observation was consistent with results obtained for sample flow rate with better efficiency at high flow rates.

However, eluent flow rate seemed to have less influence than sample flow rate. Flow rates of PAR and ammonium oxalate were thus fixed at 5 mL.min -1 .

Analytical features

Calibration curves have been constructed for various sample volumes (10, 25 and 50 mL). From these data, limits of detection (LOD; 3; n=10) and coefficient of variation (CV; n=6) were determined and summarized in Table S1 (Supplementary material). The analytical features obtained for a 50 mL sample volume seemed adapted for typical lead concentrations in natural waters: LOD was calculated at 2.7 µg.L -1 , linear domain range was between 3 and 120 µg.L -1 .

The LOD was acceptable for environmental sample analysis, but sample volume can potentially be increased if lower LOD needs to be reached for water samples with very small amounts of lead. An additional experiment has been carried out and a volume of 200 mL of sample has been passed through the TrisKem Pb resin column. The results obtained showed that no traces of lead were detected at the column outlet and that the breakthrough volume of the resin was not reached. Coefficient of variation obtained with the optimal conditions (50 mL sample volume), for a lead concentration of 50 µg.L -1 was 5.4 %.

A brief comparison of previously reported flow methods for lead determination with proposed 3D-printed flow system is given in Table 2. Some of these methods require complicated or expensive equipment in particular for detection step [START_REF] Mitani | On-line liquid phase micro-extraction based on drop-in-plug sequential injection lab-at-valve platform for metal determination[END_REF][START_REF] Beltran | Determination of lead by atomic fluorescence spectrometry using an automated extraction/pre-concentration flow system[END_REF][START_REF] Anthemidis | On-line sequential injection dispersive liquid-liquid microextraction system for flame atomic absorption spectrometric determination of copper and lead in water samples[END_REF][START_REF] Ampan | Exploiting sequential injection analysis with bead injection and lab-onvalve for determination of lead using electrothermal atomic absorption spectrometry[END_REF]. Compared to other simpler flow procedures using spectrophotometric detection [START_REF] Di Nezio | A sensitive spectrophotometric method for lead determination by flow injection analysis with on-line preconcentration[END_REF][START_REF] Mesquita | A flow system for the spectrophotometric determination of lead in different types of waters using ion-exchange for pre-concentration and elimination of interferences[END_REF], the proposed method has a lower detection limit better suited for analysis of natural water samples.

Validation

In order to validate the 3D-printed optimized system, five samples of freshwater have been collected at purposely chosen points in a coastal river ("The Arc", south of France) in areas close to anthropogenic activities. The samples were UV-photooxidized at 254 nm during 30 minutes and filtered at 0.45 µm before analysis. These samples have been analyzed in duplicate by ICP-AES and developed 3D-printed system with a sample volume of 50 mL, and the results were summarized in Table 3.

The values obtained by the two methods were consistent. The mean difference between the two methods was 5.8%, (min -11.2%; max 11.7%). The results obtained by the proposed system were compared (t test) with the reference method values (ICP-AES) and no significant differences at the 95% confidence level were found.

Conclusion

A 3D-printed system was developed for determination of lead in natural waters. Lead quantification was based on the selective solid phase extraction of lead on TrisKem Pb resin followed by elution with ammonium oxalate and spectrophotometric detection using 4-(2pyridylazo)-resorcinol as chromophoric reagent. Interferences were eliminated by optimisation of extracting and washing steps on TrisKem Pb resin. Detection limit obtained (2.7 µg.L -1 ) was consistent with environmental samples analysis but sample volume may be increased if lower detection limits are needed. The proposed flow system was compared to a reference method (ICP-AES) and was satisfactorily applied to natural waters samples. The optimized 3D-printed MPFS flow system could be controlled by an open-source microcontroller board to design a low cost portable on-line analyzer [START_REF] González | Low cost analyzer for the determination of phosphorus based on open-source hardware and pulsed flows[END_REF][START_REF] González | Multi-pumping flow system for the determination of boron in eye drops, drinking water and ocean water[END_REF]. (V = 30 mL at 3 mL.min-1), nitric acid at 5 mL.min -1 (V = 9 mL for conditioning step, V = 6

3 mL for washing step), ammonium oxalate at 0.1 mol.L -1 (V = 10 mL at 4 mL.min -1 ). µg.L -1 (V = 30 mL at 5 mL.min -1 ), nitric acid at 0.05 mol.L -1 (V = 9 mL for conditioning step, V = 6 mL for washing step, at 5 mL.min -1 ), ammonium oxalate at 4 mL.min -1 (V = 10 mL)]. (V = 30 mL at 5 mL.min -1 ), nitric acid at 0.05 mol.L -1 (V = 9 mL for conditioning step, V = 6 mL for washing step, at 5 mL.min -1 ), ammonium oxalate at 0.025 mol.L -1 (V = 10 mL), PAR at 0.01 mmol.L -1 (V= 10 mL)].

0.17

  TrisKem Pb resin. This resin was initially developed for the separation of 210 Pb and 210 Po and is constituted by crownethers diluted in isodecanol and impregnated on an inert support. The length of the carbon chain of isodecanol facilitates lead elution. Its retention capacity is 29 mg Pb/g of resin. The manufacturer's instructions recommend acidification of samples with 1 mol.L -1 nitric acid before extraction of lead on TrisKem Pb resin. However, the first experiments carried out with 1 mol.L -1 nitric acid led, after a few injections, to yellowing and cracking of the inlet of 3Dprinted column part. The effect of HNO3 concentration on the extraction of lead and potentially interfering metals was thus tested in the range 0-0.1 mol.L -1 in order to preserve 3D-printed parts.

Fig. 1 .

 1 Fig. 1. Schema of the MPFS system with the three 3D printed units. C: resin column, R1: nitric

Fig. 2 .Fig. 3 .

 23 Fig. 2. Extraction of several metals on TrisKem Pb resin for different concentrations of nitricacid [multi-metal solution at 0.037 mmol.L -1 (V = 30 mL at 3 mL.min -1 ), nitric acid at 5 mL.min - 1 (V = 9 mL for conditioning step, V = 6 mL for washing step), ammonium oxalate at 0.1 mol.L - 1 (V = 10 mL at 4 mL.min -1 )].

Fig. 4 .

 4 Fig. 4. Elution of lead at different concentrations of ammonium oxalate [lead solution at 100

Fig. 5 .

 5 Fig. 5. Absorbance versus PAR and ammonium oxalate flow rates [lead solution at 100 µg.L -1

Table 1

 1 Flow procedure for lead determination in water

				Volume	Flow-rate
	Step	Description	Action			P1	P2	P3	P4
				(mL)	(mL.min -1 )
	1	Conditioning	Pumping nitric acid	9	5	ON OFF OFF OFF
	2	Pb extraction	Pumping sample	50	5	OFF ON OFF OFF
	3	Washing	Pumping nitric acid	6	5	ON OFF OFF OFF
	4	Start data acquisition	-	-	-
	5	Elution & detection	Pumping ammonium	5	4	OFF OFF ON	ON
			oxalate and PAR
	6	Stop data acquisition	-	-	-

Acknowledgment

This work was included in the project "Lab-on-Ship" funded by the French Research Agency (ANR-14-CE04-0004).