Assyr Abdulle
email: assyr.abdulle@epfl.ch

David Cohen
email: david.cohen@kit.edu

Gilles Vilmart
email: gilles.vilmart@bretagne.ens-cachan.fr

AND Konstantinos C Zygalakis
email: k.zygalakis@soton.ac.uk

HIGH WEAK ORDER METHODS FOR STOCHASTIC DIFFERENTIAL EQUATIONS BASED ON MODIFIED EQUATIONS

Keywords: weak convergence, modified equations, backward error analysis, stiff integrator, invariant preserving integrator AMS subject classifications. 65C30, 60H35, 65L20

Inspired by recent advances in the theory of modified differential equations, we propose a new methodology for constructing numerical integrators with high weak order for the time integration of stochastic differential equations. This approach is illustrated with the constructions of new methods of weak order two, in particular, semi-implicit integrators well suited for stiff (meansquare stable) stochastic problems, and implicit integrators that exactly conserve all quadratic first integrals of a stochastic dynamical system. Numerical examples confirm the theoretical results and show the versatility of our methodology.

1. Introduction. The problem of computing the expectation of some functional of a random process appears in many practical situations, for example: in finance [START_REF] Rachev (e.D | Handbook of computational and numerical methods in finance[END_REF], in random mechanics [START_REF] Spanos | Computational stochastic mechanics[END_REF], in nonlinear filtering [START_REF] Crisan | The Oxford handbook of nonlinear filtering[END_REF] or bio-chemical processes [START_REF] Gardiner | Handbook of stochastic processes for physics, chemistry and natural sciences[END_REF], to mention a few examples. Here, we are interested in the situation where the random process is the solution of an Itô stochastic system differential equations (SDEs)

dX = f (X)dt + g(X)dW (t), X(0) = X 0 , (1.1)
where X(t), X 0 are random variables with values in R d , f : R d → R d is the drift term, g : R d → R d×m is the diffusion term with d × m matrix values, and the components W [j] (t), j = 1, . . . , m of W (t) = (W [START_REF] Abdulle | Stabilized methods for stiff stochastic systems[END_REF] (t), . . . , W [m] (t)) T are independent Wiener processes. We assume that the drift and diffusion terms are smooth enough, Lipschitz continuous and satisfy a growth bound, to ensure a unique (mean-square bounded) solution of (1.1) [START_REF] Arnold | Stochastic Differential Equations, Theory and applications[END_REF][START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF]. Analytic solutions of SDEs are rarely known and their practical computation are usually done numerically. A one-step numerical method for the approximation of (1.1) is given by

X n+1 = Ψ(f, g, X n , h, ξ n ), (1.2)
where Ψ(f, g, •, h, ξ n ) : R d → R d , X n ∈ R d for n ≥ 0, h denotes the timestep size, and ξ n denotes a random vector. Of interest in this paper is the approximation of E(φ(X(τ ))), where φ : R d → R is a smooth function, by E(φ(X N )), N = τ /h. For a practical computation, E(φ(X N )) is further approximated by a Monte-Carlo method [START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF]. The efficiency of this later approximation, which is not addressed in the present paper, is very important in practice and is still an active research topic. In particular, the methods developed in this paper could be combined with the recently proposed Multilevel Monte-Carlo method [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF].

The accuracy of the approximation can be measured by the weak order of convergence of the numerical method. We recall that a numerical approximation (1.2), starting from the exact initial condition X 0 of (1.1) is said to have weak order p if for τ > 0, we have

|E(φ(X N )) -E(φ(X(t N )))| ≤ Ch p , (1.3)
for any fixed t N = N h ∈ [0, τ ], for all h sufficiently small, and all functions φ : R d → R ∈ C

2(p+1) P (R d , R), with a constant C independent of h. Here and in what follows, C P (R d , R) denotes the space of times continuously differentiable functions R d → R with all partial derivatives with polynomial growth.

Remark 1.1. A well-known theorem of Milstein [START_REF] Milstein | Weak approximation of solutions of systems of stochastic differential equations[END_REF] allows to infer the weak order from the error after one step. Assuming that f, g are Lipschitz continuous and satisfy

f ∈ C 2(p+1) P (R d , R d ), g ∈ C 2(p+1) P (R d , R d×m
) that the moments of the exact solution of the SDE (1.1) exist and are bounded (up to a sufficiently high order) and that φ ∈ C 2(p+1) P (R d , R), then, the local error bound

|E(φ(X 1 )) -E(φ(X(t 1 )))| ≤ Ch p+1 (1.4)
for all initial values X(0) = X 0 and for all h sufficiently small, implies the global error bound (1.3). Here the constant C is again independent of h.

The simplest method to approximate solutions to (1.1) is the so-called Euler-Maruyama method [START_REF] Maruyama | Continuous markov processes and stochastic equations[END_REF], which has weak order one. In many applications, it is of interest to approximate the moments of the solution of an SDE (or expectations of functionals of it) with a better accuracy. The construction of higher order schemes has been pursued by many authors. Classical approaches for getting high weak order numerical schemes for stochastic differential equations are based on weak Taylor approximation or Runge-Kutta type methods [START_REF] Burrage | General order conditions for stochastic Runge-Kutta methods for both commuting and non-commuting stochastic ordinary differential equation systems[END_REF][START_REF] Debrabant | Classification of stochastic Runge-Kutta methods for the weak approximation of stochastic differential equations[END_REF]. For example, weak second order methods were proposed by Milstein [START_REF] Milstein | A method of second order accuracy integration of stochastic differential equation[END_REF][START_REF] Milstein | Weak approximation of solutions of systems of stochastic differential equations[END_REF], Platen [START_REF] Platen | High-order weak approximation of ito diffusions by markov chains[END_REF], Mackevicius [START_REF] Mackevicius | Second order weak approximation for Stratonovich stochastic differential equations[END_REF], Talay [START_REF] Talay | Efficient numerical schemes for the approximation of expectations of functionals of the solution of a SDE and applications[END_REF] (see also [START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF][START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF]) and Tocino and Vigo-Aguiar [START_REF] Tocino | Weak second order conditions for stochastic Runge-Kutta methods[END_REF]. We mention also the extrapolation methods of Talay and Tubaro [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF] and of [START_REF] Kloeden | Extrapolation methods for the weak approximation of Itô diffusions[END_REF] that combines methods with different stepsizes to achieve higher weak order convergence.

In this paper we propose yet another approach inspired by the construction of high order numerical integrators for deterministic problems proposed in [START_REF] Chartier | Numerical integrators based on modified differential equations[END_REF] and the newly developed theory of modified equations for stochastic differential equations [START_REF] Debussche | Weak backward error analysis for SDEs[END_REF][START_REF] Zygalakis | On the existence and the applications of modified equations for stochastic differential equations[END_REF][START_REF] Shardlow | Modified equations for stochastic differential equations[END_REF]. The basic idea of our new approach can be summarized as follows. Instead of applying the numerical method (1.2) to the SDE (1.1), we apply it to a suitably modified differential equation (a perturbation of (1.1)) so that the resulting numerical scheme yields a higher order approximation of the original SDE. This permits to fulfill automatically the order conditions, which can be very numerous for SDEs (for instance, 59 weak order two conditions have been listed for a class of stochastic Runge-Kutta type methods in [ [START_REF] Rößler | Second order Runge-Kutta methods for Itô stochastic differential equations[END_REF], Thm. 5.1]). We present a criterion (see Theorem 2.1) to construct weak methods of arbitrary order. Classical methods (Milstein or Talay methods) can be derived in a new way with our methodology. New methods will also be derived.

As an example, we propose a weak second order mean-square stable method suitable for the integration of so-called stiff problems. By stiff stochastic problems, we refer to mean-square stable problems with multiple scales for which classical explicit methods face a severe step size restriction. We also show how the methodology can be used to construct high weak order methods for random mechanical problems. In particular, we derive new weak second order methods preserving exactly all quadratic first integrals of the underlying SDE. As an illustration, we study the stochastic rigid body problem.

The paper is organized as follows. In Section 2, we present our new methodology and give a criterion for the construction of high weak order methods. In Section 3, we give explicit constructions of weak second order methods with emphasis on the numerical integration of stiff problems and random mechanical problems. Numerical examples illustrate the behavior of our new methods and corroborate the claimed weak orders of convergence.

2. Integrators based on modified equations. The general idea of constructing high order integrators based on modifying equation for SDEs can be summarized as follows. Consider a numerical method (1.2) for problem (1.1) and assume that its weak order of convergence (1.3) is p ≥ 1. We then consider (1.1) with suitably modified drift and diffusion functions

d X = f h ( X)dt + g h ( X)dW (t), X(0) = X 0 , (2.1)
where

1 f h (x) = f (x) + hf 1 (x) + h 2 f 2 (x) + . . . , (2.2)
g h (x) = g(x) + hg 1 (x) + h 2 g 2 (x) + . . . , (2.3)
and apply the numerical method (1.2) to (2.1), i.e.,

X n+1 = Ψ(f h , g h , X n , h, ξ n ).
The goal is to choose f h , g h in such a way that ( X n ) n≥0 is a better weak approximation to the solution of the original SDE (1.1), i.e.,

|E(φ( X

N )) -E(φ(X(t N )))| ≤ Ch p+r ,
with r > 0.

Remark 2.1. The above procedure should not be confused with a procedure called backward error analysis for SDEs [START_REF] Debussche | Weak backward error analysis for SDEs[END_REF][START_REF] Zygalakis | On the existence and the applications of modified equations for stochastic differential equations[END_REF][START_REF] Shardlow | Modified equations for stochastic differential equations[END_REF] developed to study the long time behavior of numerical methods for SDEs. There, one tries to find a modified equation

d X = a h ( X)dt + b h ( X)dW (t), X(0) = X 0 , (2.4)
such that its exact solution is closer to the numerical solution (1.2), i.e.,

|E(φ(X N )) -E(φ( X(t N )))| ≤ Ch p+q ,
with q > 0. In general, the modified SDEs (2.4) and (2.1) are different (see Remark

below).

A natural way of looking at expectations of functionals of diffusion processes is by using the backward Kolmogorov equation associated to (1.1), which is the (deterministic) partial differential equation

∂u ∂t = Lu, u(x, 0) = φ(x), (2.5)
where φ : R d → R is a smooth function, and the differential operator L, called the generator of the SDE (1.1), is given by

L := f • ∇ x + 1 2 (gg T ) : ∇ 2 x . (2.6)
In (2.6), ∇ x and ∇ 2 x denote respectively the gradient and the Hessian matrix operator 2 with respect to x. In the case m = d = 1, the generator reduces to

L = f ∂ ∂x + 1 2 g 2 ∂ 2 ∂x 2 .
The probabilistic interpretation (see for example [START_REF]Stochastic Differential Equations[END_REF][START_REF] Pavliotis | Multiscale methods: averaging and homogenization[END_REF][START_REF] Risken | The Fokker-Planck equation[END_REF]) of the solution u = u f,g (φ, x, t) to (2.5) is that

u f,g (φ, x, t) = E (φ(X(t))|X(0) = x) ,
where X(t) solves (1.1). For the rest of the paper we assume for simplicity that the initial condition X(0) is deterministic. We emphasize here that the results are still valid with random initial conditions, provided obvious notational changes. Using (2.5) one can easily derive the following formal Taylor expansion [START_REF] Debussche | Weak backward error analysis for SDEs[END_REF][START_REF] Zygalakis | On the existence and the applications of modified equations for stochastic differential equations[END_REF]

u f,g (φ, x, h) -φ(x) = ∞ j=1 h j j! L j φ(x).
Under appropriate smoothness assumptions on f, g and φ one can prove that

u f,g (φ, x, h) -φ(x) = k j=1 h j j! L j φ(x) + O(h k+1 ), (2.7)
for all integer k. By defining

U f,g (φ, x, h) = E (φ (Ψ(f, g, X 0 , h, ξ 0 )|X 0 = x)) , (2.8)
for the numerical integrator (1.2), we see that the local weak error of the numerical integrator applied to (1.1) after one step is given by

E(φ(X 1 )) -E(φ(X(t 1 ))) = U f,g (φ, x, h) -u f,g (φ, x, h). (2.9)
Notice that (2.9) is the reformulation of the left-hand side of the local error bound (1.4) in terms of the solution of the backward Kolmogorov equation (2.5) associated to (1.1). Motivated by an expansion of (2.8) in Taylor series, we assume Assumption 2.1. The numerical solution (2.8) has the following expansion

U f,g (φ, x, h) = φ(x) + hA 0 (f, g)φ(x) + h 2 A 1 (f, g)φ(x) + . . . , (2.10)
where A i (f, g), i = 0, 1, 2, . . . are linear differential operators depending on the drift and diffusion functions of the SDE to which the numerical integrator is applied to. We further assume that these differential operators A i (f, g), i = 0, 1, 2, . . . satisfy for all f, f , g, ĝ and ε → 0,

A i (f + ε f , g + εĝ) = A i (f, g) + ε Âi (f, f , g, ĝ) + O(ε 2 ),
2 Here, we consider the usual scalar product on matrices defined by A : B = trace(A T B).

where Âi (f, f , g, ĝ), i = 0, 1, 2, . . . are again differential operators.

The above smoothness assumption is usually satisfied by numerical integrators. For the expansion in integer powers of the stepsize h, special care has to be taken, as explained in the following remark.

Remark 2.2. The assumption that the expansion (2.10) holds with integer powers of the timestep h is essential to avoid non integer powers of h in the modified equation (2.1). For instance, let us consider the scalar θ-Milstein method

X n+1 = X n +(1-θ)hf (X n )+θhf (X n+1 )+g(X n )∆W n + 1 2 g (X n )g(X n )((∆W n ) 2 -h),
where ∆W n are independent N (0, h) distributed random variables and θ is a fixed parameter. More details on this scheme will be given in Section 3. The assumption (2.10) is not satisfied if one uses the Platen [START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF] approximation

1 2 g (X n )g(X n ) 1 2 √ h g(X n + √ hg(X n )) -g(X n ) ,
for approaching the derivative of the diffusion function in the θ-Milstein method, because (2.10) would contain a term of size O(h 5/2 ). However, if one considers instead the approximation used by Rößler [START_REF] Rößler | Second order Runge-Kutta methods for Itô stochastic differential equations[END_REF] in which the noise part is evaluated as

1 2 g (X n )g(X n ) 1 4 √ h g(X n + √ hg(X n )) -g(X n - √ hg(X n )) ,
then, the assumption (2.10) is satisfied. This can be checked by observing that the substitution √ h ↔ -√ h leaves the definition of the method unchanged. Construction of modified equations. For a numerical method (1.2) with an expansion (2.10) satisfying (see (2.9)) U f,g (φ, x, h) -u f,g (φ, x, h) = O(h p+1 ), i.e., of weak order p in view of Remark 1.1, the task now is to find a modified SDE (2.1) such that

U f h ,g h (φ, x, h) -u f,g (φ, x, h) = O(h p+r+1 ), (2.11)
i.e., a numerical method ( X n ) n≥0 of weak order p + r with r > 0 for the original problem (1.1). A second assumption that we make on the numerical integrator is that it is consistent, i.e. of weak order at least one. This assumption implies A 0 (f, g)φ = Lφ and A 0 (f h , g h )φ = Lφ, where

Lφ := f h • ∇ x φ + 1 2 (g h g T h ) : ∇ 2 x φ, (2.12)
for all function φ. Substituting f h , g h given by (2.2),(2.3), respectively, in (2.12) yields the following expansion for

L L = L + hL 1 + h 2 L 2 + . . . ,
where for j = 0, 1, 2, . . . , L j is given by

L j = f j • ∇ x + 1 2 j k=0 (g k g T j-k ) : ∇ 2 x , (2.13)
where we used the notations f 0 := f and g 0 := g. We will also sometimes write L j = L j (f j , g, g 1 , . . . , g j ) to emphasize the dependence of those operators on the functions f j ,g,g 1 ,. . . , g j .

We may now state in this section the main result of this paper. We show that under suitable assumptions, the weak order p of the numerical integrator (1.2) can be increased to p + r with r ≥ 1 by applying it to a suitably modified SDE (2.1), with modified drift and diffusion of the form

f h,s (x) = f (x) + hf 1 (x) + . . . + h s f s (x),
(2.14)

g h,s (x) = g(x) + hg 1 (x) + . . . + h s g s (x), (2.15)
where s = p + r -1. The integrator with improved weak order r can be written as

X n+1 = Ψ(f h,p+r-1 , g h,p+r-1 , X n , h, ξ n ).
(2.16)

Theorem 2.1. Assume that the numerical scheme (1.2) has order p ≥ 1 and that Assumption 2.1 holds. Let r ≥ 1 and assume that the functions f j and g j for j = 1, . . . , p+r-2 have been constructed such that X n+1 = Ψ(f h,p+r-2 , g h,p+r-2 , X n , h, ξ n ) has weak order p + r -1. Consider the differential operator defined as

L p+r-1 φ := lim h→0 u f,g (φ, x, h) -U f h,p+r-1 ,g h,p+r-1 (φ, x, h) h p+r , (2.17
)

where u f,g (φ, x, h) is expanded in (2.7) and U f,g (φ, x, h) is defined in (2.8).
If there exist functions f p+r-1 : R d → R d and g p+r-1 : R d → R d×m such that the differential operator (2.17) can be written in the form (2.13), then the numerical integrator (2.16) applied to the SDE with the modified drift and diffusion (2.14),(2.15) has weak order of accuracy p + r for the original system of SDEs (1.1) provided

f h,p+r-1 ∈ C 2(p+r+1) P (R d , R d ) g h,p+r-1 ∈ C 2(p+r+1) P (R d , R d×m ). The error bound |E(φ( X N )) -E(φ(X(t N )))| ≤ Ch p+r ,
holds for any fixed t N = N h ∈ [0, τ ] with h sufficiently small and for all functions φ ∈ C

2(p+r+1) P (R d , R). Proof. By induction hypothesis, X n+1 = Ψ(f h,p+r-2 , g h,p+r-2 , X n , h, ξ n
) is an integrator of weak order p + r -1. Thus, it has a weak expansion of the form

U f h,p+r-2 ,g h,p+r-2 (φ, x, h) = φ(x) + hA 0 (f h,p+r-2 , g h,p+r-2 )φ(x) + . . . + h p+r A p+r-1 (f h,p+r-2 , g h,p+r-2 )φ(x) + O(h p+r+1 ) = φ(x) + hL(f, g)φ(x) + . . . + h p+r-1 (p + r -1)! L p+r-1 (f, g)φ(x) + h p+r B p+r (f, g)φ(x) + O(h p+r+1 ),
where B p+r (f, g) is a certain differential operator. Using Assumption 2.1, we have the relation

A i (f h,p+r-1 , g h,p+r-1 ) = A i (f h,p+r-2 , g h,p+r-2 ) + O(h p+r-1
). The weak expansion of the modified integrator (2.16) can then be written as

U f h,p+r-1 ,g h,p+r-1 (φ, x, h) = φ(x) + hL(f, g)φ(x) + . . . + h p+r-1 (p + r -1)! L p+r-1 (f, g)φ(x) + h p+r L p+r (f p+r-1 , g, g 1 , . . . , g p+r-1 ) + B p+r (f, g) φ(x) + O(h p+r+1 ),
where L p+r-1 is defined in (2.13). If f p+r-1 and g p+r-1 are such that

L p+r-1 (f p+r-1 , g, g 1 , . . . , g p+r-1 ) = L p+r (p + r)! -B p+r (f, g),
then (2.11) holds for U f h,p+r-1 ,g h,p+r-1 (φ, x, h). Now observing that the right-hand side of the above equality is equal to the right-hand side of (2.17) together with Remark 1.1 proves the theorem.

Relation with backward error analysis.. We close this section by relating the previous construction of modified integrators with the backward error analysis for SDEs [START_REF] Debussche | Weak backward error analysis for SDEs[END_REF][START_REF] Zygalakis | On the existence and the applications of modified equations for stochastic differential equations[END_REF][START_REF] Shardlow | Modified equations for stochastic differential equations[END_REF] mentioned in Remark 2.1. Applying the numerical integrator (1.2) to the original SDE (1.1), we search for a modified differential equation (2.4) such that

U f,g (φ, x, h) -u fh ,ĝ h (φ, x, h) = O(h p+q+1 ) (2.18)
with q > 0. The aim in such a procedure is to better understand the behavior of the numerical method (1.2) (applied to (1.1)) by studying the modified SDE (2.4). The modified SDE (2.4), with fh , ĝh given by an expansion

fh (x) = f (x) + h f1 (x) + h 2 f2 (x) + . . . , ĝh (x) = g(x) + hĝ 1 (x) + h 2 ĝ2 (x) + . . . ,
has an associated backward Kolmogorov equation (the formula (2.5) with L replaced by L) in (2.6), where

L = L0 + h L1 + h 2 L2 + . . . , (2.19)
where L0 = L and for j = 1, 2, . . . Lj is given by

Lj = fj • ∇ x + 1 2 j k=0 (ĝ k ĝT j-k ) : ∇ 2 x ,
where ĝ0 = g. The Taylor expansion (2.7) becomes

u fh ,ĝ h (φ, x, h) -φ(x) = k j=1 h j j! Lj φ(x) + O(h k+1 ),
which gives in terms of the expansion (2.19) (see [START_REF] Zygalakis | On the existence and the applications of modified equations for stochastic differential equations[END_REF])

u fh ,ĝ h (φ, x, h) -φ(x) = k j=1 h j i1+i2+...+i l +l=j 1 l! ( Li1 • • • Li l φ)(x) + O(h k+1 ).
The task in this approach is to find fh , ĝh such that for U f,g (φ, x, h) given by (2.10) holds

A j (f, g) = i1+i2+...+i l +l=j 1 l! Li1 • • • Li l ,
for j = p + q. The above relation permits to define by induction the differential operators Lj used to construct the modified equation for backward error analysis. We emphasize once more that the aim and the theory for integrators based on modified equations and backward error analysis are different. In the former approach, the modified SDE constitutes only a surrogate to obtain a better numerical approximation of the solution of the original SDE, in the latter approach, the modified SDE is a tool to study a numerical integrator applied to the original SDE. Remark 2.3. In the case p = r = 1, the above procedure yields for backward error analysis and for modified integrators the operators L1 = A 1 -1 2 L 2 and L 1 = 1 2 L 2 -A 1 , respectively. Thus, the perturbations f1 , ĝ1 in the modified equations for backward error analysis and f 1 , g 1 for modified integrators are identical up to the multiplicative factor -1.

3. High weak order methods with application to stiff problems and geometric integration. In this section we show two applications of the methodology developed in Section 2. We first derive a class of weak second order methods based on first order methods. Classical methods (Milstein or Talay methods) will be recovered, but new methods will also be derived. In particular, we derive a new weak second order method which is mean-square stable, suitable for the integration of so-called stiff problems. This method belongs to a general class of weak second order methods derived by Milstein [START_REF] Milstein | Weak approximation of solutions of systems of stochastic differential equations[END_REF], but seems not to have appeared explicitly in the literature. Secondly, we show how our methodology can be applied to structure preserving integrators and derive weak second order methods preserving quadratic invariants. As an example, we consider the stochastic rigid body problem.

3.1.

Weak second order methods with application to stiff stochastic problems. To illustrate our methodology based on modifying equations, we derive here a family of weak second order methods. For that, we pick a weak first order method

X 1 = Ψ(f, g, X 0 , h, ξ 0 ), consider the modified equation dX = [f (X) + hf 1 (X)] dt + [g(X) + hg 1 (X)] dW (t), X(0) = X 0 , (3.1)
and apply Theorem 2.1. Accordingly, we have to find f 1 , g 1 such that

L 1 = L 2 2 -A 1 (f, g), (3.2)
where

L 1 := f 1 ∇ x + 1 2 (gg T 1 + g 1 g T ) : ∇ 2
x , and where the differential operator A 1 depends on the choice of the weak first order method.

3.1.1. One-dimensional case. For the sake of simplicity let us first consider a one dimensional SDE with one dimensional noise. The simplest weak first order method is the Euler-Maruyama method. However, for reasons explained in Remark 3.2 below, this is not a suitable method to start with. A fairly general class of weak first order methods that can be used for our purpose is the θ-Milstein method [START_REF] Higham | A-stability and stochastic mean-square stability[END_REF] (that we denote θ-M method in what follows)

X n+1 = X n +(1-θ)hf (X n )+θhf (X n+1 )+g(X n )∆W n + 1 2 g (X n )g(X n )((∆W n ) 2 -h), (3.3
) where ∆W n are independent N (0, h) distributed random variables and X 0 = x. For θ = 0, (3.3) reduces to the classical explicit Milstein method, while for θ > 0 it yields semi-implicit methods (implicit in the drift function and explicit in the diffusion function). We expand in Taylor series the function φ up to the 4th order, φ(X 1 ) =

4 i=0 1 i! φ (i) (x)F i + ..., where F = (1 -θ)hf (x) + θhf (X 1 ) + g(x)∆W 0 + 1 2 g (x)g(x)((∆W 0 ) 2 -h),
and obtain 3

U f,g (φ, x, h) = E(φ(X 1 )|X 0 = x) = φ(x) + hLφ(x) + h 2 A 1 (f, g)φ(x) + O(h 3 ),
where

A 1 (f, g)φ(x) = θ f (x)f (x) + 1 2 f (x)g 2 (x) φ (x) + 1 2 f 2 (x) + 2θf (x)g 2 (x) + 1 2 (g (x)g(x)) 2 φ (x) + 1 2 g (x)g 3 (x) + g 2 f (x) φ (x) + h 2 8 g 4 (x)φ (4) (x).
Applying the method (3.3) to the modified equation (3.1) we obtain U f h ,g h (φ, x, h) which is a second order approximation of u f,g (φ, x, h) if we can find f 1 , g 1 such that (3.2) holds. A simple computation reveals that

1 2 L 2 φ -A 1 (f, g)φ (x) = 1 2 -θ f (x)f (x) + 1 2 f (x)g 2 (x) φ (x) + 1 2 -θ f (x)g(x) + 1 2 g (x)f (x) + 1 4 g 2 (x)g (x) g(x)φ (x).
We see from the above formula that we can define the appropriate operator L 1 with

f 1 (x) = 1 2 -θ f (x)f (x) + 1 2 1 2 -θ f (x)g 2 (x), g 1 (x) = 1 2 -θ f (x)g(x) + 1 2 g (x)f (x) + 1 4 g 2 (x)g (x).
Now setting f h,1 = f + hf 1 and g h,1 = g + hg 1 , we obtain the following new integrator

X n+1 = X n + (1 -θ)hf h,1 (X n ) + θhf h,1 (X n+1 ) + g h,1 (X n )∆W n + 1 2 g (X n )g(X n )((∆W n ) 2 -h), (3.4)
which has weak order two for the SDE (1.1) in dimension one. Remark 3.1. In principle one should also replace g by g h,1 in the last term of (3.4), but omitting hg 1 for this term does not affect the weak order two of accuracy. Indeed,

g h,1 (X n )g h,1 (X n ) = g (X n )g(X n ) + C(X n )h + O(h 2 ), where C(x) is a smooth function. Using E(C(X n )h((∆W n ) 2 -h)) = 0, we deduce E g h,1 (X n )g h,1 (X n )((∆W n ) 2 -h) = E(g (X n )g(X n )((∆W n ) 2 -h)) + O(h 3 ), 3 Recall that E(∆W 0 ) = E(∆W 3 0 ) = 0 and E(∆W 2 0 ) = h, E(∆W 4 0 ) = 3h 2 .
and thus it does not influence the accuracy of the method because it induces a O(h 3 ) perturbation of E(φ(X n+1 )). Notice that the integrator (3.4) belongs to a sub-class of a general family of weak second order methods introduced by Milstein [START_REF] Milstein | Weak approximation of solutions of systems of stochastic differential equations[END_REF]. For θ = 0 it has also been considered by Talay who proved its order of convergence [START_REF] Talay | Efficient numerical schemes for the approximation of expectations of functionals of the solution of a SDE and applications[END_REF].

For θ = 1/2 the method was considered by Milstein who showed its good stability behavior for scalar SDEs with additive noise. For θ = 1, the method does not seem to have appeared explicitly in the literature. We will show below that it has favorable stability properties for scalar SDEs with multiplicative noise (mean-square stability). Remark 3.2. Notice that L 2 is a differential operator of order four in general. Thus, the difference 1 2 L 2 -A 1 (f, g) is a differential operator of order two of the same form as L only if A 1 (f, g) contains the same third and fourth order derivatives of φ as 1 2 L 2 . As explained further, this is true for the Milstein method. However, this would not be the case for the Euler-Maruyama method where the term 1 2 g (x)g(x) 3 φ (x) involving the third derivative of φ in 1 2 L 2 φ is not cancelled in general (unless g = 0, i.e. for additive noise). Therefore, as observed in [START_REF] Zygalakis | On the existence and the applications of modified equations for stochastic differential equations[END_REF], a modified SDE cannot be constructed for the Euler-Maruyama method.

3.1.2. Multi-dimensional case. The formula derived for the one-dimensional case can easily be extended to the multi-dimensional case. Consider the multidimensional SDE (1.1), where f is a column vector of size d and g is a d × m matrix (below, we denote by • [i] the ith component of a vector in R d and by • [i,j] the coefficients of a d × m matrix). For a fixed parameter θ, consider the θ-M method

X n+1 = X n + (1 -θ)hf (X n ) + θhf (X n+1 ) + g(X n )∆W n + M (X n , W ), (3.5)
where the Milstein term M (X n , W ) is defined for i = 1, . . . , d by

M [i] = Ξ i (X n ) : I = m j1,j2=1 Ξ i[j 1 ,j2] I [j1,j2] .
The coefficients of the m × m matrix Ξ i are defined for i = 1, . . . , d by

Ξ i[j 1,j2] = d k=1 ∂g [i,j2] ∂x k g [k,j1] ,
and the coefficients of the m × m matrix I of multiple integrals are given by

I [j1,j2] = tn+1 tn s tn dW j1 (t) dW j2 (s). (3.6)
Following the procedure for the one-dimensional case, we obtain the modified θ-M method of weak order two

X n+1 = X n + (1 -θ)hf h,1 (X n ) + θhf h,1 (X n+1 ) + g h,1 ∆W n + M (X n , W ), (3.7)
where f h,1 = f + hf 1 , g h,1 = g + hg 1 , and f 1 , g 1 are given (componentwise) by

f 1,[i] = 1 2 -θ (f f ) [i] + 1 2 1 2 -θ gg T : f [i] , g 1,[i,j] = 1 2 -θ (f g) [i,j] + 1 2 g [i,j] f + 1 4 gg T : g [i,j] , (3.8)
for all i = 1, . . . , d and j = 1, . . . , m. We observe that the above method contains derivatives of the drift and diffusion functions. This is a general feature of the methods obtained using modified equations. In some cases, these derivatives are easy and cheap to compute (see for example the stochastic mechanical problem in Section 3.2.1). In general, these derivatives can be approximated. In particular, one can use formulas based on finite differences. Some care is however required for an efficient implementation (i.e., a low number of function evaluations in dependence on the number of Wiener processes [START_REF] Debrabant | Families of efficient second order Runge-Kutta methods for the weak approximation of Itô stochastic differential equations[END_REF]). Remark 3.3. The multiple integral matrix I in (3.6) is difficult to evaluate in general and needs to be approximated. One can use for instance the following weak approximation for the matrix I in the definition (3.6)(see for instance [22, eq. (5.12.9)])

J = 1 2 ∆W n ∆W T n + E n ,
where E n is a random skew-symmetric matrix whose coefficients E n,[j1,j2] are independent two-point distributed random variables,

P(E n,[j1,j2] = ±h) = 1/2, for all 1 ≤ j 1 < j 2 ≤ d and E n,[j1,j2] = -E n,[j2,j1
] for all j 1 , j 2 = 1, . . . , m. Using J instead of I does not alter the weak order two of accuracy of the method (3.7) (it does however decrease the strong order of the method from 1 to 1/2). The independent Gaussian variables ∆W n,[j] can also be replaced by independent three point random variables with

P(∆W n,[j] = ± √ 3h) = 1 6 , P(∆W n,[j] = 0) = 2 3 , (3.9)
without decreasing the weak order two of the method.

3.1.3.

A mean-square stable modified θ-M method. In this section we show that we can construct a second order modified θ-M method with favorable mean-square stability. To study the stability in the mean-square sense of numerical integrators, a widely used test equation introduced in [START_REF] Saito | Stability analysis of numerical schemes for stochastic differential equations[END_REF] for SDEs is the following scalar SDE with multiplicative noise

dX = λXdt + µXdW (t), (3.10)
where the parameters λ, µ ∈ C. We notice that other test equations have been considered recently in [START_REF] Buckwar | Towards a systematic linear stability analysis of numerical methods for systems of stochastic differential equations[END_REF][START_REF] Buckwar | A comparative linear mean-square stability analysis of Maruyama-and Milstein-type methods[END_REF], to better account for the stability behavior of numerical integrators when applied to systems of SDEs or scalar equations with several multiplicative noise terms. The mean-square stability domain of (3.10) is given by

S = {(λ, µ) ∈ C 2 ; λ + 1 2 |µ| 2 < 0}. (3.11)
The set of (λ, µ) that fulfill condition (3.11) can be visualized, for λ, µ ∈ R, as the shaded area with a boundary given by the dotted parabolas in Figure 3.1. The θ-M method applied to the linear test equation (3.10) yields where V n are independent Gaussian variables with a N (0, 1) distribution and p = λh and q = µ √ h. Squaring the result and taking the expectation, we obtain the relation

X n+1 = (1 + p(1 -θ) + qV n + 1 2 q 2 (V 2 n -1)) 1 -θp X n ,
E(|X n+1 | 2 ) = R θ,M (p, q)E(|X n | 2 ), where R θ,M (p, q) = |1 + p(1 -θ)| 2 + |q| 2 + |q| 4 /2 |1 -θp| 2 . (3.12)
We next define the set

S θ,M = {(p, q) ∈ C 2 ; R θ,M (p, q) < 1}.
The method is called mean-square (MS) stable if R θ,M (p, q) ≤ 1, for all (p, q) ∈ S, or alternatively if S ⊂ S θ,M . It is readily seen that there does not exist a value of θ ∈ [0, 1] such that the θ-M method is MS stable. Furthermore, MS stability is recovered for θ = 3/2 [START_REF] Higham | A-stability and stochastic mean-square stability[END_REF]. Remark 3.4. In contrast, the θ-methods (the methods (3.5) with M ≡ 0), whose stability function reads

R θ (p, q) = |1 + p(1 -θ)| 2 + |q| 2 |1 -θp| 2 ,
can be shown to be MS stable if and only if θ ≥ 1/2 as reported in [START_REF] Higham | Mean-square and asymptotic stability of the stochastic theta method[END_REF]. We come now to study the stability properties of the modified θ-M methods (3.4) whose stability functions can be easily deduced from (3.12) and read

R θ,M (p, q) = |1 + p(1 -θ)| 2 + | q| 2 + |q| 4 /2 |1 -θ p| 2 ,
where p = p + ( 1 2 -θ)p 2 , q = q + (1 -θ)pq. A simple calculation shows that this method is MS stable if and only if θ = 1. Thus, for θ = 1 we have constructed the weak second order method (3.4) which is MS stable. This method is thus suitable for the numerical integration of stiff systems of SDEs as illustrated in the numerical example below.

In Figure 3.1 we plot the mean-square stability domain for the standard θmethods, the standard θ-M methods and the modified θ-M methods (the light-dark region which lies inside the dotted parabola is the stability domain S of the exact solution of the test problem). For θ ∈ [0, 1], it can be seen that the θ-M methods are never MS stable, and that only for θ = 1 is the modified θ-M method MS stable. Numerical experiments. We illustrate the numerical behavior of the modified θ-M methods previously constructed. We consider an economy model for asset prices proposed in [START_REF] Hofmann | Option pricing under incompleteness and stochastic volatility[END_REF], see also [START_REF] Higham | MAPLE and MATLAB for stochastic differential equations in finance[END_REF]. It is an Itô system of SDEs in dimension d = 3, with m = 2 non-commutative noises, given by

dX [1] = β 1 X [1] X [2] dW [1] (t), dX [2] = -(X [2] -X [3] )dt + β 2 X [2] dW [2] (t), (3.13
)

dX [3] = α(X [2] -X [3] )dt,
where X [START_REF] Abdulle | Stabilized methods for stiff stochastic systems[END_REF] (t), X [START_REF] Abdulle | Chebyshev methods for stiff stochastic differential equations[END_REF] (t) and X [START_REF] Abdulle | S-ROCK methods for stiff Itô SDEs[END_REF] (t) represent the asset price, the instantaneous volatility and the average volatility, respectively. We take the parameters

β 1 = 1, β 2 = 0.3, the initial value X(0) = (X [1] (0), X [2] (0), X [3] (0)) T = (1, 0.1, 0.1) T
, and consider the time interval [0, 1] as in [START_REF] Higham | MAPLE and MATLAB for stochastic differential equations in finance[END_REF]. The parameter α > 0 corresponds to the volatility parameter, and is related to the strength of the past dependence of the average volatility. We refer to [START_REF] Hofmann | Option pricing under incompleteness and stochastic volatility[END_REF] for details in the context of economy modeling. For large values of α, the largest eigenvalue in modulus of the Jacobian of the drift function is -|α| and the SDE problem (3.13) becomes stiff. We shall consider various values of the volatility parameter α in the numerical experiments.

Since the drift vector field in (3.13) is linear, the modified θ-M methods (3.14) for (3.13) are linearly implicit and using the formulas (3.7) and (3.8), it can be written as (a) Error for E(X 2 [START_REF] Abdulle | Stabilized methods for stiff stochastic systems[END_REF] ). Nonstiff case α = 1. ). Very stiff case α = 100. Comparison of weak convergence rates for the modified θ-M with θ = 1 (solid lines), the θ-method with θ = 1 (dashed-dotted lines), the θ-method with θ = 0 (Euler-Maruyama method, dotted line), and the modified θ-M with θ = 0 (Milstein-Talay method, dashed line).

(Id -θhA)X n+1 = (Id + (1 -θ)hA)X n + g h,1 (X n )∆W n + M (X n , ∆W n ), (3.14)
where A denotes the matrix

A = 1 -h 1 2 -θ 1 + α    0 0 0 0 -1 1 0 α -α    ,
and where

g h,1 (X)∆W =    β 1 X [1] X [2] + h 2 (X [3] -X [2] ))∆W [1] β 2 1 -h( 1 2 -θ) X [2] + h 2 (X [3] -X [2] ) ∆W [2] αβ 2 h 1 2 -θ X [2] ∆W [2]    , M (X, ∆W ) =    1 2 β 2 1 X [1] X 2 [2] (∆W 2 [1] -h) + β 1 β 2 X [1] X [2] I [2,1] 1 2 β 2 2 X [2] (∆W 2 [2] -h) 0    .
We take for the random variables ∆W n = (∆W n, [START_REF] Abdulle | Stabilized methods for stiff stochastic systems[END_REF] , ∆W n, [START_REF] Abdulle | Chebyshev methods for stiff stochastic differential equations[END_REF] ) T independent Gaussian variables with mean zero and variance h. Notice that similar results have been obtained when considering instead discrete random variables satisfying (3.9). The above multiple integral I [START_REF] Abdulle | Chebyshev methods for stiff stochastic differential equations[END_REF][START_REF] Abdulle | Stabilized methods for stiff stochastic systems[END_REF] in the Milstein term M (X, ∆W ) is approximated by

I [2,1] ≈ (∆W [1] ∆W [2] + ξ n h)/2
, where ξ n are independent random variables satisfying P(ξ n = ±1) = 1/2 as detailed in Remark 3.3.

To confirm the weak order two of convergence of the modified θ-M method (3.7), we compute the relative errors in the quantities E(X 2 [START_REF] Abdulle | Stabilized methods for stiff stochastic systems[END_REF] ) in Figures 3.2(a)-(c) and E(X 2 [START_REF] Abdulle | Chebyshev methods for stiff stochastic differential equations[END_REF] ) in Figure 3.2(d) at the final time t = 1 for the stepsizes h = 2 -i , i = 0, . . . , 7. The reference solutions are computed using the small timestep h = 2 -14 .

To check carefully the accuracy of the methods up to small time steps, we need to drastically reduce the Monte-Carlo error. We thus approximate the required moments of the numerical solutions by averages over 500 millions of trajectories computed in fortran, using the random number generator [START_REF] Petersen | Lagged Fibonacci series random number generators for the NEC SX-3[END_REF]. For a fair comparison, notice that we use the same set of random numbers for each numerical integrator. We observe in Figure 3.2 the expected lines of slope two (solid lines) both in the nonstiff case (α = 1) and the stiff cases (α = 25 and α = 100). Notice that for small timesteps (h < 0.25) the zigzag that we observe is due to the Monte-Carlo error, which could be further reduced by increasing the number of samples. For comparison, we also plot the results for the classical semi-implicit θ-method (θ = 1) (weak order one), obtained from (3.14) by removing the Milstein term M (X n , ∆W n ) and setting h = 0 in the definitions of A and g h,1 (X)∆W . We also compare with two classical explicit integrators, the Euler-Maruyama method (weak order one) and the Talay method (weak order two), obtained by taking θ = 0 in (3.14). Notice these two explicit methods are not (unconditionally) mean-square stable. Indeed, since for large α, the largest eigenvalue in the drift function of (3.13) has size α, the mean-square stability constraint for these explicit methods has the form αh ≤ C, where C is a constant of moderate size, independent of h and α. We observe in Figures 3. This numerical experiment shows that the modified θ-M method with θ = 1 has the same (unconditional) mean-square stability as the standard θ-method (θ ≥ 1/2), but with an improved accuracy by several orders of magnitude due to the improved weak order two of convergence.

High weak order integrators preserving quadratic invariants.

In this section, we construct numerical integrators for Stratonovich SDEs of high weak order two which exactly conserve all quadratic first integrals (up to machine precision). We consider the SDE (1.1) in Stratonovich form with a one-dimensional noise

dX = f (X)dt + g(X) • dW (t), X(0) = X 0 , (3.15)
where the notation •dW (t) emphasizes that the Stratonovich stochastic integrals are considered for (3.15). As a basic numerical integrator to apply our methodology of modified equation, we choose the (fully) implicit midpoint rule, as first introduced in [START_REF] Milstein | Numerical methods for stochastic systems preserving symplectic structure[END_REF],

X n+1 = X n + hf X n + X n+1 2 + g X n + X n+1 2 ∆W n , (3.16)
where ∆W n is a scalar random variable. It is shown in [START_REF] Milstein | Numerical methods for stochastic systems preserving symplectic structure[END_REF] that (3.16) has weak and strong orders one in the case of a one-dimensional or commutative multi-dimensional noise. Notice however that for general SDEs with multi-dimensional noise, the strong order is 1/2 and the weak order is 1. Remark 3.5. The method (3.16) is implicit with respect to both the drift and the noise terms. In the case where ∆W n is a standard Gaussian variable, the unboundedness of ∆W n for arbitrarily small h leads to non-uniqueness of solutions to the non-linear system (3.16) and the integrator is not well defined. One way to address this problem, is to replace ∆W n , with a suitable chosen bounded random variable [START_REF] Milstein | Numerical methods for stochastic systems preserving symplectic structure[END_REF] (see also [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF]Sect. 1.3]). Here we shall simply consider discrete random variable, e.g. (3.9), as in Remark 3.3, which are obviously bounded.

First integral conservation. A smooth quantity C(x) is called a first integral of the system (3.15) if it is exactly conserved along time for all realizations of the Wiener process W (t), i.e. C(X(t)) = C(X 0 ) for all time t and all initial condition X(0) = X 0 . Given a smooth function C(x), the identity

4 dC(X) = ∇C(X) • f (X)dt + ∇C(X) • g(X) • dW (t) shows that C(X) is a first integral of (3.15) if and only if ∇C(x) • f (x) = ∇C(x) • g(x) = 0 for all x ∈ R d .
(3.17)

Proposition 3.1. The implicit midpoint rule (3.16) exactly preserves all quadratic first integrals C(x) of (3.15), i.e. C(x n+1 ) = C(x n ) for all h and all realizations of ∆W n .

Proof. Since C(x) is quadratic, we write C(x) = x T Sx, where S is a constant symmetric matrix. Using (3.17), with

X = (X n+1 + X n )/2 and R(X) = hf (X) + g(X)∆W n , a short computation yields C(X n+1 )-C(X n ) = X T n SR(X)+R(X) T SX n+1 = 2X T SR(X) = ∇C(X)•R(X) = 0.
3.2.1. New invariant preserving integrators of high weak order. Using the framework of integrators based on modified equations, we introduce the following new numerical integrator

X n+1 = X n + hf h,1 X n + X n+1 2 + g h,1 X n + X n+1 2 ∆W n , (3.18)
where f h,1 = f + hf 1 and g h,1 = g + hg 1 and show below that for

f 1 = 1 4 1 2 f (g, g) -g f g g 1 = 1 4 1 2 g (g, g) -g g g , (3.19)
the numerical integrator is a weak second order method for the SDE (3.15) which preserves all quadratic first integrals. We notice that if we consider the modified

Stratonovich SDE dX = [f (X) + hf 1 (X)] dt + [g(X) + hg 1 (X)] • dW (t), (3.20
= f (X) + 1 2 g (X)g(X) dt + g(X)dW (t), (3.21)
where compared to the Itô system of SDEs (1.1), the vector field f is replaced by f + 1 2 g g. This permits to deduce an expansion analogue to (2.7) associated to the Itô SDE (3.21). The weak expansion of (3.18) (applied to (3.15), equivalent to the Itô SDE (3.21)) is computed as follows. First we have (for X 0 = x) X 1 = x + F = X 0 + hf (x + F/2) + g(x + F/2)∆W 0 , where

F = hf + h 2 f F + h 8 f (F, F ) + g∆W 0 + 1 2 g F ∆W 0 + 1 8 g (F, F )∆W 0 + 1 48 g (F, F, F )∆W 0 + O(h 5/2 ).
For the computation of A 1 (f, g)φ we consider the expansion

φ(X 1 ) = φ(x + F ) = φ(X 0 ) + k F [k] ∂ k φ + 1 2 kl F [k] F [l] ∂ kl φ + 1 6 klm F [k] F [l] F [m] ∂ klm φ + 1 24 klmi F [k] F [l] F [m] F [i] ∂ klmi φ + . . .
We then compute E φ(X 1 )|X 0 = x) = E φ(x + F ) , identify the differential operator multiplying the term h 2 , and obtain after some tedious but straightforward computations,

1 2 L 2 φ -A 1 (f, g) φ = 1 4 1 2 f (g, g) -g f g + 1 4 g (g, g, g) - 1 4 g g (g, g) -g g g g • ∇ x φ + 1 8 g 1 2 g (g, g) -g g g T + 1 2 g (g, g) -g g g g T : ∇ 2 x φ = f 1 + 1 2 (g 1 g + g g 1 ) • ∇ x φ + 1 2 gg T 1 + g 1 g T : ∇ 2 x φ,
where we define f 1 = 1 4 1 2 f (g, g) -g f g and g 1 = 1 4 1 2 g (g, g) -g g g . The modified Itô SDE of Theorem 2.1 then reads

dX = f h,1 + 1 2 g g + h 2 (g 1 g + g g 1 ) (X)dt + g h,1 (X)dW (t),
where f h,1 = f + hf 1 and g h,1 = g + hg 1 . Using g h,1 g h,1 = g g + h(g 1 g + g g Proof. Consider the original midpoint rule (3.16) applied to (3.15). Using Remark 2.3, we obtain that the modified SDE up to second order for backward error analysis associated to (3.16) is given by (3.20) with h replaced by -h,

d X = f ( X) -hf 1 ( X) dt + g( X) -hg 1 ( X) • dW (t), X(0) = X 0 ,
Example: a stochastic rigid body model. To illustrate that the integrators previously introduced conserve quadratic first integrals and to compare the performance of the various methods proposed (see Table 3.1), we consider a randomly perturbed rigid body problem that is, the motion of a rigid body in R 3 subject to a scalar white noise perturbation. The equations of motion of an asymmetric rigid body with Stratonovich noise in dimension m = 1 are given by 5

dX = XI -1 Xdt + µ Xe 1 • dW (t), dQ = Q I -1 Xdt + µQ e 1 • dW (t), (3.25)
where e 1 = (1, 0, 0) T , µ ≥ 0 is a parameter and I = diag (I 1 , I 2 , I 3 ). A generalization of equation (3.25) for a 3-dimensional noise is presented in [25, Eq. (6.9)-(6.10)], where one can also find a physical justification for these equations. This model is a variant of the model proposed in [START_REF] Liao | Random motion of a rigid body[END_REF] with the additional feature that it preserves the spatial angular momentum QX, as detailed below. In the case where µ = 0, we recover the standard deterministic equations of motion of an asymmetric rigid body. We refer to [17, Sect. VI.5] for a survey of geometric and invariant preserving integrators for the rigid body problem in the context of ODEs. The constants I 1 , I 2 , I 3 > 0 are the moments of inertia which characterize the rigid body. The function X(t) represents the angular momentum in R 3 in the body frame. The matrix Q(t) is a rotation matrix in R 3 which gives the orientation of the rigid body in a fixed frame. Notice that the first line in (3.25) can be rewritten simply as

dX [1] = 1 I 3 - 1 I 2 X [2] X [3] dt, dX [2] = 1 I 1 - 1 I 3 X [3] X [1] dt + µX [3] • dW (t), dX [3] = 1 I 2 - 1 I 1 X [1] X [2] dt -µX [2] • dW (t).
The system of SDEs (3.25) has several first integrals, all of which are quadratic. It has QX as first integral, which represents the spatial momentum in R 3 with respect to the body frame. It also has

Q T Q = Id as first integral because Q is an orthogonal matrix. Since Q is orthogonal, the Casimir C(X) = 1 2 X 2 [1] + X 2 [2] + X 2 [3] is also conserved. Considering the Hamiltonian H(X) = 1 2 X 2 [1] /I 1 + X 2 [2] /I 2 + X 2 [3] /I 3 , we have dH(X) = µ X [2] X [3] 2 1 I 2 - 1 I 3 • dW (t),
which shows that H(X) is also a first integral if and only if I 2 = I 3 (symmetric body) or µ = 0 (the noise is zero). Using formulas (3.19) where the functions f and g correspond to the right-hand side of (3.25), a straightforward computation yields the modified SDE associated to 5 We use the standard hat notation for the correspondence between 3×3 skew-symmetric matrices and size 3 vectors, X =

  0 -X [3] X [2] X [3] 0 -X [1] -X [2] X [1] 0   , for all X =   X [1] X [2] X [3]   . (3.25), dX = X(I -1 + hµ 2 4 J -1 )Xdt + µ 1 + hµ 2 4 Xe 1 • dW (t), (3.26) dQ = Q( I -1 X + hµ 2 4 J -1 X)dt + µ 1 + hµ 2 4 Q e 1 • dW (t),
where we define J = diag (I 1 , I 3 , I 2 ). We obtain from Theorem 3.2 that applying the implicit midpoint rule (3.16) to the Statonovitch SDE (3.26) yields a weak order two approximation of the solution of (3.25) which exactly conserves all quadratic first integrals, i.e. C(X

n+1 ) = C(X n ), Q n+1 X n+1 = Q n X n and Q T n Q n =
I 1 = 1 I 1 1 + hµ 2 4 , 1
I 2 = 1 I 2 + hµ 2 4I 3 , 1
I 3 = 1 I 3 + hµ 2 4I 2 , 1
yields the modified SDE (3.26). Thus, our modification to high weak order reduces to a perturbation of the parameters and has a negligible overcost. Implementation. We now detail the implementation of the standard implicit midpoint rule (3.16) for the stochastic rigid body problem (3.25). The implementation of the modified implicit midpoint rule (3.18) (and similarly for the method (3.22)) is deduced using Remark 3.6 by modifying the moments of inertia I 1 , I 2 , I 3 . We refer to [START_REF] Chartier | Numerical integrators based on modified differential equations[END_REF] for the implementation of the corrector f 2 in (3.22) and (3.24).

It is a standard approach to use quaternions q n to represent the orthogonal matrices Q n (see [START_REF] Hairer | Geometric numerical integration. Structure-preserving algorithms for ordinary differential equations[END_REF] in the context of rigid body integrator implementations). The implicit midpoint rule (3.16) for the angular momentum X(t) can be written as

X n+1 = X n + h Y I -1 Y + µ Y e 1 ∆W n ,
where we denote Y = (X n+1 + X n )/2. This implicit system can be solved by a few fixed point iterations. Next, the configuration update

Q n+1 = Q n + h Q n + Q n+1 2 I -1 Y + µ Q n + Q n+1 2 e 1 ∆W n , is equivalent to Q n+1 = Q n Ω,
where Ω is the orthogonal matrix defined by the Cayley transform Ω

= Id + Z Id -Z -1 with Z = h 2 I -1 Y + ∆Wnµ 2 e 1 .
Thus, the configuration update for the rotation matrix Q n reduces to a multiplication of quaternions 6

q n+1 = q n • ω ω , with ω = 1 + h 2 i Y [1] I 1 + j Y [2] I 2 + k Y [3] I 3 + i ∆W n 2 µ,
where the matrices Q n , Q n+1 are represented by the quaternions q n+1 , q n , respectively.

6 Notice that ω denotes the norm of the quaternion ω, thus ω/ ω is a quaternion of norm 1. Convergence rates. We consider the system (3.25) on the time interval [0, 10]. We take the moments of inertia I 1 = 0.345, I 2 = 0.653, I 3 = 1.0, which correspond to the water molecule (nearly flat body). Initial values are X(0) = (0.8, 0.6, 0) T and Q(0) is the identity matrix. We have carefully implemented the above integrators in fortran, using quaternions for the rotation matrices. In Figure 3.3, we plot the errors for E(X 2 [START_REF] Abdulle | Stabilized methods for stiff stochastic systems[END_REF] ), E(sin(e T 1 Qe 1 )) at final time t = 10 versus the timestep h = 2 -i , i = 1, . . . , 8. The reference solution is computed using the small timestep h = 2 -14 . To check carefully the accuracy of the methods, we compute numerically E(X 2 [START_REF] Abdulle | Stabilized methods for stiff stochastic systems[END_REF] ) and E(sin(e T 1 Qe 1 )) using the averages over 300 millions of trajectories. We consider two values of the noise parameter: µ = 0.1 and µ = 0.01. We observe in all cases lines of slope two for the modified midpoint rule imr2 (3.18) which confirms its weak order two of accuracy. For the standard midpoint rule imr in (3.16) and the modified version imr(4) in (3.24) which both have weak order one, we observe for large stepsizes h, lines of slope four and two respectively in the case where the deterministic error (h 2 or h 4 ) is dominant compared to µ 2 h. For smaller timesteps, we retrieve lines of slope one, the weak order of these two methods. Similarly, for the improved modified midpoint rule imr2(4) in (3.22), we observe lines of slope four for large timesteps and only two for small timestep. This corroborates the theoretical results collected in Table 3.1.

Conclusion.

In this paper, we introduced a new framework for increasing the weak order of accuracy of a given numerical method for SDE by considering the numerical integration of a suitably modified problem. Our methodology, which uses tools developed for backward error analysis for stochastic problems [START_REF] Zygalakis | On the existence and the applications of modified equations for stochastic differential equations[END_REF][START_REF] Debussche | Weak backward error analysis for SDEs[END_REF][START_REF] Shardlow | Modified equations for stochastic differential equations[END_REF], generalizes to SDEs the framework of numerical integrators based on modified equations introduced in [START_REF] Chartier | Numerical integrators based on modified differential equations[END_REF] for deterministic problems. This approach permits to fulfill automatically the numerous order conditions for high weak order schemes. We illustrated our approach at the example of the θ-Milstein-Talay method, and obtained for θ = 1 a mean-square stable integrator of weak order two. The numerical experiments conducted for a stiff problem in economy show an improvement in accuracy of two orders of magnitude over the classical θ-method of weak order one.

In the spirit of backward error analysis for the study of geometric integrators for ODEs, where the modified equations inherit the geometric properties of the integrators, we also derived new high weak order integrators based on the implicit midpoint rule, that automatically conserve all quadratic first integrals. The efficiency of the approach is illustrated at the example of a stochastic rigid body model which possesses several quadratic first integrals. A natural extension of this work would be to search for modified equations to construct new integrators of weak order three or more with good stability or geometric properties.

We note that this new approach also allows to construct higher order Chebyshev methods for stiff SDEs. An attempt to generalize such methods, introduced in [1, 2, 3], to higher order has been proposed in [START_REF] Burrage | Explicit stochastic Runge-Kutta methods with large stability regions[END_REF][START_REF] Komori | Weak second order s-rock methods for Stratonovich stochastic differential equations[END_REF]. This generalization involves the solution of a large number of order conditions and the resulting methods appear to have less favorable stability properties than the method proposed in [START_REF] Abdulle | Stabilized methods for stiff stochastic systems[END_REF][START_REF] Abdulle | Chebyshev methods for stiff stochastic differential equations[END_REF][START_REF] Abdulle | S-ROCK methods for stiff Itô SDEs[END_REF]. In [START_REF] Abdulle | Second weak order explicit stabilized methods for stiff stochastic differential equations[END_REF], we show that using techniques based on modifying equations as proposed in this paper, it is possible to construct high weak order Chebyshev method in an efficient way with better stability properties than the method given in [START_REF] Burrage | Explicit stochastic Runge-Kutta methods with large stability regions[END_REF][START_REF] Komori | Weak second order s-rock methods for Stratonovich stochastic differential equations[END_REF].

Fig. 3 . 1 .

 31 Fig. 3.1. Comparison of mean-square stability domains (dark gray area) of the θ-methods, the θ-M methods, and the modified θ-M methods. The modified θ-M methods with θ = 0, 1/2 have been proposed by Milstein and Talay. Horizontal axis: p, vertical axis: q.

2 Slope 1 (

 21 (θ=1) θ-method (θ=0) Modified θ-M (θ=0) Modified θ-M (θ=1) Slope c) Error for E(X 2 [1] ). Stiff case α = 25.

2 Slope 1 (

 21 (θ=1) θ-method (θ=0) Modified θ-M (θ=0) Modified θ-M (θ=1) Slope d) Error for E(X 2 [2]). Stiff case α = 25.

Fig. 3 . 2 .

 32 Fig.3.2. Finance model(3.13). Comparison of weak convergence rates for the modified θ-M with θ = 1 (solid lines), the θ-method with θ = 1 (dashed-dotted lines), the θ-method with θ = 0 (Euler-Maruyama method, dotted line), and the modified θ-M with θ = 0 (Milstein-Talay method, dashed line).

2

 2

  (c)-(d) that these methods are indeed unstable for h > 2 -4 for the (moderately) stiff case α = 25. For the very stiff case α = 100 (Figure 3.2(b)), these methods show too much instability to fit in the scales of the figure and are thus omitted.

3 . 6 .

 36 Id for all n, and in the case I 2 = I 3 (symmetric body), we have alsoH(X n+1 ) = H(X n ).Remark Notice that the modified SDE (3.26) is of the same form as the original rigid body equations (3.25) with modified data parameters. Indeed, replacing µ by µ = µ(1 + hµ 2 /4), and replacing I = diag (I 1 , I 2 , I 3 ) in the original SDE (3.25) by I = diag ( I 1 , I 2 , I 3 ), where

  Error for E(X 2[START_REF] Abdulle | Stabilized methods for stiff stochastic systems[END_REF] ). Noise size µ = 0.1.

  Error for E(X 2[START_REF] Abdulle | Stabilized methods for stiff stochastic systems[END_REF] ). Noise size µ = 0.01.

  Error for E(sin(e T 1 Qe 1 )). Noise size µ = 0.1.

  Error for E(sin(e T 1 Qe 1 )). Noise size µ = 0.01.

Fig. 3 . 3 .

 33 Fig. 3.3. Rigid body problem (3.25). Comparison of weak convergence rates for imr, see (3.16) (dotted lines), imr2, see (3.18) (solid lines), imr2(4), see (3.22) (dashed lines), and imr(4), see (3.24) (dashed-dotted lines).

  Proof. The Stratonovich SDE (3.15) is equivalent to the Itô SDE dX

	)
	then (3.18) is equivalent to applying the original midpoint rule (3.16) to the modified
	Stratonovich SDE (3.20).
	Theorem 3.2. The integrator (3.18) for a system of Stratonovich SDEs (3.15)
	with m = 1 noise has weak order 2. It exactly conserves all quadratic first integrals
	of (3.15).

Here, h is the timestep size of the numerical method (1.2).

Notice that Stratonovich calculus is used here.

Acknowledgments. The research of A. A. is partially supported under Swiss National Foundation Grant 200021 140692. K. C. Z. was supported by Award No. KUK-C1-013-04 of the King Abdullah University of Science and Technology (KAUST).

and we have from (2.18) with p = 1, q = 0, E(φ(X 1 )) -E(φ( X(h))) = O(h 2 ). Using Proposition 3.1, we have C(X 1 ) = C(X 0 ). On the one hand, replacing φ(x) by φ(C(x)), we obtain E(φ(C(X 0 ))) -E(φ(C( X(h)))) = O(h 2 ). On the other hand, using (3.17), we have

where dφ(C( X)) has size O(h). This yields E(φ(C(X 0 ))) -E(φ(C( X(h)))) = 0 for all test functions φ, and thus C( X(h)) = C(X 0 ). We obtain

We close this section by indicating that a further modification of the integrator (3.18) allows yet an even better accuracy. Consider

where f h,2 = f + hf 1 + h 2 f 2 and g = g + hg 1 (as previously) with f 1 , g 1 as defined in (3.19) and f 2 given by

The above term f 2 corrects the deterministic error of size O(h 2 ), but the weak order of the integrator (3.22) remains 2. Notice that it would be interesting to search for modified fields f 2,h , g 2,h to achieve the weak order 3. For comparison, we also consider the integrator of weak order one

where (3.23). Notice that the integrators (3.22) and (3.24) are equivalent to the modifying implicit midpoint rule of order 4 for ODEs (deterministic case) introduced in [START_REF] Chartier | Numerical integrators based on modified differential equations[END_REF] in the case where the diffusion function g is zero.

We summarize in the following table our theoretical findings. To study the convergence rates in dependence of the noise size, we consider (3.15), where •dW (t) is replaced by •µ dW (t), with µ ∈ R + a given scaling parameter.

weak order order for ODEs weak error method of accuracy (no noise: µ = 0) E(φ(X N )) -E(φ(X(t N )))


Guillaume Ch Èze
email: guillaume.cheze@math.univ-toulouse.fr

DARBOUX THEORY OF INTEGRABILITY IN THE SPARSE CASE

Darboux's theorem and Jouanolou's theorem deal with the existence of first integrals and rational first integrals of a polynomial vector field. These results are given in terms of the degree of the polynomial vector field.

Here we show that we can get the same kind of results if we consider the size of a Newton polytope associated to the vector field. Furthermore, we show that in this context the bound is optimal.

Introduction

In this paper we study the following polynomial differential system in C n :

dX 1 dt = A 1 (X 1 , . . . , X n ), . . . , dX n dt = A n (X 1 , . . . , X n ),
where A i ∈ C[X 1 , . . . , X n ] and deg A i ≤ d. We associate to this polynomial differential system the polynomial derivation

D = n i=1 A i (X 1 , . . . , X n )∂ Xi .
The computation of first integrals of such polynomial differential systems is an old and classical problem. The situation is the following: we want to compute a function F such that the hypersurfaces F (X 1 , . . . , X n ) = c, where c are constants, give orbits of the differential system. Thus we want to find a function F such that D(F ) = 0.

In 1878, G. Darboux [START_REF] Darboux | Mémoire sur les équations différentielles du premier ordre et du premier degré[END_REF] has given a strategy to find first integrals. One of the tools developed by G. Darboux is now called Darboux polynomials. A polynomial f is said to be a Darboux polynomial, if D(f ) = g.f , where g is a polynomial. The polynomial g is called the cofactor. A lot of properties of a polynomial differential system are related to Darboux polynomials of the corresponding derivation D, see e.g. [START_REF] Goriely | Integrability and nonintegrability of dynamical systems[END_REF][START_REF] Dumortier | Qualitative theory of planar differential systems[END_REF]. There exists a lot of different names in the literature for Darboux polynomials, for example we can find: special integrals, eigenpolynomials, algebraic invariant hypersurfaces, or special polynomials. G. Darboux shows in [START_REF] Darboux | De l'emploi des solutions particulières algébriques dans l'intégration des systèmes d'équations différentielles algébriques[END_REF] that if the derivation D has at least n+d-1 n + 1 Darboux polynomials then D has a first integral which can be expressed by means of these polynomials. More precisely the first integral has the following form: i f λi i where f i are Darboux polynomials and λ i are complex numbers. This kind of integral is called today a Darboux first integral.

Date: October 29, 2012.

In 1979, J.-P. Jouanolou shows in his book [START_REF] Jouanolou | Équations de Pfaff algébriques[END_REF], that if a derivation has at least n+d-1 n + n Darboux polynomials then the derivation has a rational first integral. We recall that a rational first integral is a first integral which belongs to C(X 1 , . . . , X n ). Several authors have given simplified proof for this result. M. Singer proves this result in C 2 , see [START_REF] Michael | Liouvillian first integrals of differential equations[END_REF]. This approach is based on a work of Rosenlicht [START_REF] Rosenlicht | On Liouville's theory of elementary functions[END_REF]. J.-A. Weil generalizes this strategy and gives a proof for derivations in C n , see [START_REF] Weil | Constantes et polynômes de Darboux en algèbre différentielle: applications aux systèmes diffrentiels linéaires[END_REF]. J. Llibre and X. Zhang gives a direct proof of Jouanolou's result in [START_REF] Llibre | Rational first integrals in the Darboux theory of integrability in C n[END_REF].

Darboux and Jouanolou theorem are improved in [START_REF] Llibre | Darboux theory of integrability for polynomial vector fields in R n taking into account the multiplicity at infinity[END_REF][START_REF] Llibre | Darboux theory of integrability in C n taking into account the multiplicity[END_REF]. The authors show that we get the same kind of result if we take into account the multiplicity of Darboux polynomials. The multiplicity of a Darboux polynomial is defined and studied in [START_REF] Christopher | Multiplicity of invariant algebraic curves in polynomial vector fields[END_REF].

The Darboux theory of integrability has been successfully used in the study of some physical problems, see e.g. [START_REF] Valls | Rikitake system: analytic and Darbouxian integrals[END_REF][START_REF] Llibre | Darboux integrability and algebraic invariant surfaces for the Rikitake system[END_REF], and in the study of limit cycles and centers, see e.g. [START_REF] Christopher | Invariant algebraic curves and conditions for a centre[END_REF][START_REF] Schlomiuk | Algebraic particular integrals, integrability and the problem of the center[END_REF][START_REF] Llibre | Configurations of limit cycles and planar polynomial vector fields[END_REF]. Unfortunately, to the author knowledge, there do not exist example showing if these bounds are optimal. In this note we are going to study the situation in the sparse case. This means that we are going to consider polynomials A i with some coefficients equals to zero. In this situation, the size of the polynomials A i is not measured by the degree but by the size of its Newton polytope. We recall that the Newton polytope of a Laurent polynomial f (X) = α c α X α , where X = X 1 , . . . , X n and α is a multi-index (α 1 , . . . , α n ) ∈ Z n , is the convex hull in R n of the exponent α of all nonzero terms of f . We denote this polytope by N (f ).

In this note we prove a result improving Darboux and Jouanolou theorem. Our result depends on the size of a Newton polytope associated to the derivation and not on the degree d. Furthermore, in this context we can give example showing that the bound is optimal. This is our result:

Theorem 1. Let D = n i=1 A i (X 1 , . . . , X n )∂ Xi be a derivation. Consider generic values (x 1 , . . . , x n ) in C n and the polytope N D = N n i=1 x i A i X i .
Let B be the number of integer points in We can remark that if we consider dense polynomials A i with degree d, that is to say each coefficient of A i is nonzero, then B = n+d-1 n . Thus Theorem 1 gives the classical bounds in the dense case. Now, we illustrate why these bounds can be better than the classical ones. We give an example with n = 2 in order to give a picture. If A i has the following form:

N D ∩ N n , then (
A i (X 1 , X 2 ) = c e,e X e
1 X e 2 + c e-1,e X e-1 1 X e 2 + c e,e-1 X e 1 X e-1 2 + c 0,0 , then B = 3e + 2, and d = deg(A i ) = 2e. In this situation we have n+d-1 2 = 2e(2e + 1)/2. Thus for such examples Theorem 1 gives a linear bound instead of a quadratic bound.

Figure 1 shows the Newton polygon of A i (X 1 , X 2 ), when e = 3. The triangle corresponds to the Newton polygon of dense polynomials with total degree 6. In this situation, Jouanolou's theorem says that if we have 23 Darboux polynomials then we have a rational first integral. Here, our bound improves this result and means that 13 Darboux polynomials are sufficient to construct a rational first integral.

N A i (X 1 , X 2 ) .
The Newton polygon of

A i (X 1 , X 2 ) X 1 corresponds to a translation of the Newton polygon of A i (X 1 , X 2 ). In Figure 2 we show N A i (X 1 , X 2 ) X 1 . 0 - 6 X 2
A i (X 1 , X 2 ) X 1 .
Figure 3 shows the part of the Newton polygon of

x 1 A 1 (X 1 , X 2 ) X 1 + x 2 A 2 (X 1 , X 2 ) X 2 in N 2 , when x 1 , x 2 are generic. 0 - 6 X 2 X 1 q q q q q p q p q p p p p q q p p q q p p p p p p p p p p d d ¨¡ ¡ ¡ Figure 3. Newton polygon N x 1 A 1 (X 1 , X 2 ) X 1 + x 2 A 2 (X 1 , X 2 ) X 2 ∩ N 2 .
Structure of the paper. In Section 1 we give some results about Newton polytopes and weighted degree. In Section 2 we prove Theorem 1 and we show with an example that the bound is optimal.

1. ToolBox

In this section we introduce some notations and results that will be useful in Section 2. These kinds of tools are already present in the work of Ostrowski, see [START_REF] Ostrowski | On multiplication and factorization of polynomials. I. Lexicographic orderings and extreme aggregates of terms[END_REF]. For more results about sparse polynomials, see e.g. [START_REF] Gelfand | Discriminants, resultants and multidimensional determinants[END_REF][START_REF] Cox | Using algebraic geometry[END_REF].

Definition 2. Let P be a polytope, then H is a supporting hyperplane of P if

(1)

H ∩ P = ∅, (2)
P is fully contained in one of the two halfspaces defined by H.

In our situation, as we consider Newton polytopes, the equation of H has integer coefficients. More precisely, the equation of H is ν.m = a, where . denotes the usual scalar product, ν is a vector with integer coefficients, and a is an integer.

We can represent a Newton polytope with the equations of its supporting hyperplanes:

N (f ) = {m ∈ Z n | ν j .m ≤ a j , for j = 1, . . . , k},
where k is the number of supporting hyperplanes. Now we define a degree related to a vector ν ∈ Z n .

Definition 3. Let ν ∈ Z n , we set deg ν (f ) = max m∈N (f ) ν.m.
Now, we explain why we can call deg ν (f ) a degree.

Proposition 4. Let f and g be two polynomials in C[X 1 , . . . , X n ], ν = (ν 1 , . . . , ν n ) in Z n and x 1 , x 2 two generic elements in C 2 .

(1

) deg ν (f + g) ≤ max deg ν (f ), deg ν (g) , (2) deg ν (x 1 f + x 2 g) = max deg ν (f ), deg ν (g) , (3) deg ν (f.g) = deg ν (f ) + deg ν (g), (4) deg ν (∂ Xi f ) ≤ deg ν (f ) -ν i = deg ν (f /X i ).
Proof.

(1) As N (f + g) is included in the convex hull of N (f ) ∪ N (g), we have

deg ν (f + g) = max m∈N (f +g) ν.m ≤ max m∈Conv(N (f )∪N (g)) ν.m,
where Conv(.) denotes the convex hull. Furthemore, as we consider a convex set, we deduce that this maximum is reached on a point in N (f ) ∪ N (g). Thus

deg ν (f + g) ≤ max m∈N (f )∪N (g) ν.m ≤ max max m∈N (f ) ν.m, max m∈N (g) ν.m ≤ max deg ν (f ), deg ν (g) .
(2) As x 1 , x 2 are generic then N (x 1 f + x 2 g) is equal to the convex hull of N (f ) ∪ N (g). Indeed, with generic x 1 and x 2 we avoid simplifications in the sum x 1 f + x 2 g. Then the proof in this case proceeds as before.

(3) This result comes from the well-known result by Ostrowski, [START_REF] Ostrowski | On multiplication and factorization of polynomials. I. Lexicographic orderings and extreme aggregates of terms[END_REF], which gives: N (f.g) = N (f ) + N (g), where + in this situation is the Minkowski sum. ( 4) Let e i be the i-th vector of the canonical basis of R n , then we have max

m∈N (∂X i f ) ν.m ≤ max m∈N (f ) ν.(m -e i ) = max m∈N (f ) ν.m -ν i = deg ν (f ) -ν i .
We also have max

m∈N (f ) ν.(m -e i ) = max m∈N (f /Xi) ν.m
this completes the proof.

Newton polytopes and degree deg ν are related by the following proposition.

Proposition 5. Let f ∈ C[X ±1 1 , . . . , X ±1
n ] be a Laurent polynomial with corresponding Newton polytope:

N (f ) = {m ∈ Z n | ν j .m ≤ a j , for j = 1, . . . , k},
where ν j .m = a j , whith a j ∈ Z, and ν j ∈ Z n are the equations of the k supporting hyperplanes of

N (f ). Let g ∈ C[X ±1 1 , . . . , X ±1 n ] such that for j = 1, . . . , k, deg νj (g) ≤ deg νj (f ) then N (g) is included in N (f ).
Proof. We just have to remark that max n∈N (g)

ν j .n = deg νj (g) ≤ deg νj (f ) = max m∈N (f ) ν j .m = a j .
Thus each element in N (g) satisfies the equations of N (f ).

Proof of Theorem 1

2.1. Newton polytope and cofactors. We are going to show that if A i are sparse then the cofactors are sparse. This property will be the main tool of the proof of Theorem 1. Proposition 6. Let D = n i=1 A i (X 1 , . . . , X n )∂ Xi be a derivation. Let f be a Darboux polynomial with corresponding cofactor g. Consider generic values (x 1 , . . . , x n ) in C n and let N D be the convex set

N n i=1 x i A i X i then N (g) ⊂ N D ∩ N n . Proof. First, obviously N (g) ∈ N n since g is a polynomial. Second, we have g.f = n i=1 A i ∂ Xi f , thus for all ν ∈ Z n we have deg ν (g.f ) = deg ν ( n i=1 A i ∂ Xi f ).
Thanks to Proposition 4, we deduce these inequalities deg ν (g

) + deg ν (f ) ≤ max i (deg ν (A i )∂ Xi f ) ≤ max i (deg ν (A i ) + deg ν ∂ Xi f ) ≤ max i (deg ν (A i ) + deg ν (f ) -ν i ).
Thus we have

deg ν (g) ≤ max i (deg ν (A i ) -ν i ) deg ν (g) ≤ max i deg ν ( A i X i ) deg ν (g) ≤ deg ν n i=1 x i . A i X i .
Now, we apply Proposition 5 with ν ∈ Z n corresponding to supporting hyperplanes of N n i=1 x i .

A i X i and we get the desired result.

We can now prove easily Theorem 1. Indeed, we use the classical strategy to prove Darboux Theorem in our situation. As for all cofactors g fi associated to a Darboux polynomial f i , we have, by Proposition 6, N (g fi ) ⊂ N D ∩ N n , then all cofactors belong to a C-vector space of dimension B, where B is the number of integer points in N D ∩ N n . Thus if we have B + 1 cofactors, then there exists a relation between them: (⋆) i∈I λ i g fi = 0, where λ i are complex numbers. Now, we recall a fundamental and straightforward result on Darboux polynomials: g f1.f2 = g f1 + g f2 . Thus relation (⋆) gives the Darboux first integral i∈I f λi i . This proves the first part of Theorem 1. Now, in order to prove the second part of our theorem, we can use the strategy proposed in [START_REF] Llibre | Rational first integrals in the Darboux theory of integrability in C n[END_REF]. In [START_REF] Llibre | Rational first integrals in the Darboux theory of integrability in C n[END_REF] the authors show that if we have n relations of the type (⋆) then we can deduce a relation with integer coefficients, i.e. λ i ∈ Z. This gives a first integral of this kind: i∈I f λi i with λ i ∈ Z, thus this first integral belongs to C(X 1 , . . . , X n ) and we have a rational first integral. In the sparse case as the cofactors belong to a C-vector space of dimension B, if we have B + n Darboux polynomials then we have B + n cofactors and then we deduce n relations between the cofactors. With the strategy used in [START_REF] Llibre | Rational first integrals in the Darboux theory of integrability in C n[END_REF] we obtain that the derivation has a rational first integral.

2.2. The bound is optimal. Consider a polynomial p(X 1 ) ∈ C[X 1 ] with degree d. Let α be a root of p and ξ 2 , . . . , ξ n be distinct complex numbers such that p ′ (α), ξ 2 , . . . , ξ n are Z-independent. We denote by D the following derivation:

D = p(X 1 )∂ X1 + ξ 2 X 2 ∂ X2 + • • • + ξ n X n ∂ Xn .
This derivation has no non-trivial rational first integrals, by Corollary 5.3 in [START_REF] Goriely | Integrability and nonintegrability of dynamical systems[END_REF]. Indeed, (α, 0, . . . , 0) is a fixed point of the polynomial vector field, and the corresponding eigenvalues are distinct and Z independent.

By Proposition 6, if g is a cofactor then N (g) ⊂ N D ∩ N n . Here, N D ∩ N n is the set of univariate polynomials in X 1 with degree smaller than d-1. Thus B = d.

Furthermore, (X 1α), (X 1α 2 ), . . . , (X 1α d ) where α i are roots of p, and X 2 , . . . , X n are Darboux polynomials. Thus we have d+n-1 Darboux polynomials.

In conclusion, we cannot improve the bound given in Theorem 1, since there exists a derivation without rational first integrals which has B + n -1 Darboux polynomials.

Figure 1 .

 1 Figure 1. Newton polygon N A i (X 1 , X 2 ) .

Figure 2 .

 2 Figure 2. Newton polygon NA i (X 1 , X 2 ) X 1 .