Gilles Vilmart
email: gilles.vilmart@ens-rennes.fr

Weak second order multi-revolution composition methods for highly oscillatory stochastic differential equations with additive or multiplicative noise

Keywords: highly-oscillatory stochastic differential equation, composition method, quadratic first integral conservation, multiplicative noise, time-dependent stochastic Schrödinger equation. MSC numbers: 60H10, 35Q55

HAL is

Introduction

Consider a nonlinear system of (Itô 1 ) stochastic differential equations (SDEs) of the form

dX(t) = (ε -1 AX(t) + f (X(t)))dt + m r=1 g r (X(t))dW r (t), 0 ≤ t ≤ T, (1)
where ε > 0 is fixed, A ∈ R d×d is a constant matrix, f, g r : R d → R d are smooth and Lipschitz vector fields, X(0) = X 0 is the initial condition, and W r , r = 1, . . . , m are independent standard Wiener processes. On the highly oscillatory part of the problem, we assume that exp(A) = I.

(

) 2
This means that the flow of the stiff oscillatory part dX dt = ε -1 AX, given by x → exp(ε -1 tA)x, is periodic with respect to t, with period ε. Thus, multiple oscillatory frequencies in ε -1 A are allowed if they are integer multiples of 2π. Such class of problems includes second order SDEs of the form

Ẍ(t) = -ε -2 K 2 X(t) + a(X(t)) + m r=1 b r (X(t)) Ẇr (t) (3)
where K ∈ R d×d is a constant positive symmetric matrix with eigenvalues in 2πN and a, b r : R d → R d are smooth vector fields. An interesting situation is the case where a = -∇U derives from a potential U : R d → R and where additive noise is considered, i.e. the functions g r are constant and there exists a constant matrix B ∈ R d×m such that (g 1 (q), . . . , g m (q)) = B.

Consider the Hamiltonian, which represents the energy of the system (3),

H(p, q) = 1 2 p T p + ε -2 q T K 2 q + U (q) .
A standard application of the Itô formula to (3) yields that the average of Hamiltonian grows linearly with time due to the additive noise perturbation. Precisely, setting Q(t) = X(t) and P (t) = Ẋ(t), we have

E(H(P (t), Q(t))) = H(P (0), Q(0)) + t 2 trace(BB T ) (5)
and this energy is exactly conserved along time only in the deterministic case (B = 0). The linear growth [START_REF] Castella | Splitting methods with complex times for parabolic equations[END_REF] is not recovered in general by standard explicit integrators, e.g. the Euler-Maruyama method where a super linear growth can be observed [START_REF] Cohen | On the numerical discretisation of stochastic oscillators[END_REF][START_REF] Cohen | A trigonometric method for the linear stochastic wave equation[END_REF].

The class of problems (1) also includes spectral spatial discretizations of the nonlinear Schrödinger equation with a stochastic space-time noise. Consider for instance i ∂u(x, t) ∂t = -∆u(x, t) + εV (x)|u(x, t)| 2q u(x, t) + √ εg(u(x, t)) Ẇ (x, t), t ≥ 0, x ∈ (0, 1), [START_REF] Castella | Stroboscopic averaging for the nonlinear Schrödinger equation[END_REF] where we consider periodic boundary conditions in dimension one of space for simplicity. Here, W (x, t) denotes a real-valued white noise which is white in time and correlated in space. 2 We refer for more details to [START_REF] Debussche | Numerical simulation of focusing stochastic nonlinear Schrödinger equations[END_REF] where numerical simulations of the stochastic Schrödinger equation are presented to investigate the influence of space-time noise over the stochastic blow up time of the solutions, and to [START_REF] De Bouard | Weak and strong order of convergence of a semidiscrete scheme for the stochastic nonlinear Schrödinger equation[END_REF] where the strong and weak convergence rates of the Euler-Maryama method applied to a class of stochastic nonlinear Schrödinger equations are shown to be the same as for SDEs in finite dimension (i.e. weak order 1) in contrast to the parabolic case. In the case of multiplicative noise in the Stratonovitch sense where g(u) = σu (with σ a constant), we have that the L 2 norm arrive at the following system of SDEs with Stratonovitch noise for the Fourier modes

ξ k (t) ∈ C, -ℓ < k ≤ ℓ, dξ k = -ik 2 ξ k + ε f k (. . . , ξ -1 , ξ 0 , ξ 1 , . . .) dt + iσ √ ε ℓ<j≤ℓ ξ k-j • dW j , (7)
where f k corresponds to the nonlinear term V |u| 2q u, W j (-ℓ < j ≤ ℓ) are independent standard Wiener processes, and we set ξ j±2ℓ = ξ j for all j in the above sum (the convolution product of the solution frequencies and the Wiener processes). The Stratonovitch system of SDEs [START_REF] Chartier | Multi-revolution composition methods for highly oscillatory differential equations[END_REF] can be recast into the format (1) by rescaling time (that is, by rewriting the system in terms of the new time variable t = 2π ε t) and using the Stratonovitch to Itô conversion formula.

For the numerical approximation of highly oscillatory problems of the form [START_REF] Abdulle | Stabilized multilevel Monte Carlo method for stiff stochastic differential equations[END_REF], standard SDE integrators require in general a time stepsize of the same magnitude as ε to achieve stability and accuracy, which can be prohibitively expansive for small values of ε, already in the deterministic setting (g r = 0, r = 1, . . . , m in (1)), and thus appropriate integrators are needed. For the numerical study of highly oscillatory stochastic integrators and not assuming [START_REF] Abdulle | The heterogeneous multiscale method[END_REF], we mention the recent work [START_REF] Tao | Structure preserving stochastic impulse methods for stiff Langevin systems with a uniform global error of order 1 or 1/2 on position[END_REF], where splitting methods for the Langevin equations are analyzed, and the works [START_REF] Cohen | On the numerical discretisation of stochastic oscillators[END_REF][START_REF] Cohen | Convergence analysis of trigonometric methods for stiff secondorder stochastic differential equations[END_REF] where stochastic trigonometric methods are proposed and analyzed, and further extended to the linear stochastic wave equation [START_REF] Cohen | A trigonometric method for the linear stochastic wave equation[END_REF]. The trigonometric methods are a modification of the classical leap-frog method for [START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF] where filter functions are introduced to avoid numerical resonances (see [START_REF] Hairer | Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF]Chap. XIII] in the context of deterministic problems). A common aspect of the aforementioned works is the study of the strong convergence (with orders 1/2 or 1), i.e. the error E(|X(t k )-X k |) (where t k = kh and h is the stepsize) for approximating individual trajectories X(t) themselves in the case of additive noise (i.e. (4) holds). In contrast, we consider here the general case of a general nonlinear noncommutative noise, and we focus on the weak error of convergence, i.e. the error |E(φ(X(t k ))) -E(φ(X k ))| where φ is a smooth test function.

In [START_REF] Chartier | Multi-revolution composition methods for highly oscillatory differential equations[END_REF], the class of deterministic multi-revolution composition methods was recently introduced for systems of the form (1) with g r = 0, r = 1, . . . , m. Observe that the exact flow of this problem after one period ε is a smooth perturbation of the identity. The multirevolution approach, as first proposed in [START_REF] Melendo | A new approach to the construction of multirevolution methods and their implementation[END_REF][START_REF] Petzold | Numerical solution of highly oscillatory ordinary differential equations[END_REF], permits the approximation of the N th iterate of a smooth deterministic near identity map ϕ ε (y) = y +O(ε) from R d into itself, but calculating only a few compositions (hence the name multi-revolution). Setting H := N ε, a example of such composition method is the order two approximation of the N th iterate of the map ϕ ε ,

ϕ (H+ε)/2 • ϕ -1 -(H-ε)/2 (y) = ϕ N ε (y) + O(H 3 ).
It is shown in [START_REF] Chartier | Multi-revolution composition methods for highly oscillatory differential equations[END_REF] that the constant symbolized by O is independent of N ≥ 1 and ε ≤ ε 0 if ϕ ε (y) is a C 3 function of (y, ε). Next, considering for ϕ ε the flow of an highly oscillatory system after one highly oscillatory period permits to design large time step geometric integrators. The aim of this paper is to extend and analyze this approach to the stochastic context with general nonlinear multi-dimensional drift and diffusion functions.

It is known (see [START_REF] Sanders | Averaging methods in nonlinear dynamical systems[END_REF] and the recent work [START_REF] Chartier | Higher-order averaging, formal series and numerical integration I: B-series[END_REF]) in a deterministic context (g r = 0) that the solution X(t) of ( 1) is asymptotically close to an effective nonstiff problem of the form

dX(t) dt = a 1 + εa 2 + +ε 2 a 3 + . . . X(t) , X(0) = X 0 (8)
where

a 1 (x) = 1 0 e sA f (e -sA x)ds a 2 (x) = - 1 2 1 0 s 2 0 e s 2 A f (e (s 1 -s 2 )A f (e -s 1 A x)) -e s 1 A f (e (s 2 -s 1 )A f (e -s 2 A x)) ds 1 ds 2 , . . . (9)
at times which are integer multiples of the oscillatory period. Precisely, truncating the above series (8) (which diverges for general nonlinear problems) after the ε p-1 term yields X(t) = X(t) + O(ε p ) for all t = kε ≤ T with k ∈ N, and this remainder can be made exponentially small for analytic data. We observe that the vector field of the effective problem [START_REF] Chartier | Convergence analysis of multi-revolution composition time-splitting pseudo-spectral methods for linear Schrödinger equations[END_REF] involving multiple integrals can be difficult to simulate in general. Although an analogous effective problem can be constructed in some cases (see Remark 3.5 for additive noise), we highlight that the proposed methods do not involve the pre-calculation of such an effective problem but applies directly to the original problem (1), with micro and macro stepsizes in the spirit of the so-called Heterogeneous Multiscale Method [START_REF] Abdulle | The heterogeneous multiscale method[END_REF][START_REF] Engquist | The heterogeneous multiscale methods[END_REF][START_REF] Engquist | Heterogeneous multiscale methods: a review[END_REF]. This paper is organized as follows. In Section 2, we introduce the stochastic multirevolution composition methods and state our main weak convergence results. In Section 3, we perform the error analysis of the proposed methods for general nonlinear systems of SDEs. In Section 4, we present several numerical examples and comparisons with other oscillatory integrators that corroborate the theoretical orders of convergence and illustrate the good qualitative behavior of the proposed methods for various problems over long times, including the stochastic nonlinear Schrödinger equation with multiplicative noise.

Multi-revolution composition methods for oscillatory SDEs

We first introduce semi-discrete multi-revolution methods involving a macro time stepsize H which can be possibly much larger than the oscillatory period ε. We next introduce the fully-discrete methods by coupling the multi-revolution approach with a micro integrator involving in addition a micro stepsize h. We present weak error estimates with respect to H, h with error constants independent of ε. The proofs are postponed to Section 3.

Semi-discrete multi-revolution composition methods

We introduce the following integrator which permits to integrate (1) with oscillatory period ε on a time interval of size O(1) at a computational cost independent of ε by considering appropriate auxiliary non-stiff SDE problems. Algorithm 2.1 (Semi-discrete S-MRCM). For the approximation of the flow after time H = N ε of (1) where N ∈ N, we consider the scheme X k → X k+1 defined by

dK k,1 = (-AK k,1 + α N Hf (K k,1 ))dt + α N H m r=1 g r (K k,1 )dW k,r,1 , K k,1 (0) = X k ( 10
)
dK k,2 = (AK k,2 + β N Hf (K k,2 ))dt + β N H m r=1 g r (K k,2 )dW k,r,2 , K k,2 (0) = K k,1 (1) (11)
X k+1 = K k,2 (1)
where W k,r,1 , W k,r,2 , r = 1, . . . , m are independent Wiener processes, and we define

α N := 1 2 - 1 2N , β N := 1 2 + 1 2N . ( 12
)
For notational brevity, we shall sometimes drop the dependence on k of K k,j and W k,r,j . The advantage of the above Algorithm 2.1 is that the stiff problem (1) on a time interval H = N ε including N fast oscillations (due to the term ε -1 AX in (1)) is approached by the resolution of two non-stiff problems posed on the time interval (0, 1) involving each only one oscillation, and thus more convenient to solve. In other words, in this multirevolution approach, N fast oscillations of (1) are approached by only two oscillations of analogous problems with appropriate coefficients. This is advantageous compared to standard integrators for large values of N , i.e. for small values of ε.

Remark 2.2. In the context of deterministic oscillatory problems, order conditions for MRCMs up to arbitrary high order are derived in [START_REF] Chartier | Multi-revolution composition methods for highly oscillatory differential equations[END_REF] using the algebraic framework of labelled rooted trees, and integrators up to order 4 are exhibited. However, it is known [START_REF] Sheng | Solving linear partial differential equations by exponential splitting[END_REF][START_REF] Sussman | Chaotic evolution of the solar system[END_REF] that a composition method with real coefficients of order strictly larger than 2 necessarily involves negative coefficients, see [START_REF] Blanes | On the necessity of negative coefficients for operator splitting schemes of order higher than two[END_REF] for an elegant geometric proof. Since stochastic problems such as (1) are not reversible in time, such high order methods cannot straightforwardly be applied. A possibility to circumvent this order barrier for non-reversible problems would be to use complex coefficients, as proposed in [START_REF] Castella | Splitting methods with complex times for parabolic equations[END_REF][START_REF] Hansen | High order splitting methods for analytic semigroups exist[END_REF] in the context of deterministic diffusion problems. The use of complex coefficients in this context is however beyond the scope of the present paper.

The next task is to prove accuracy estimates between the exact solution X(t) of (1) and the approximation X k from Algorithm 2.1 where t = kH and H = N ε with error constants independent of H ≤ H 0 and ε ≤ ε 0 . We make the following smoothness assumption of the data.

(H) The functions f ,g r , r = 1, . . . , m are C 6 -functions with all partial derivatives bounded. This implies that f, g r are Lipschitz continuous (but not necessarily bounded), and thus the solution of (1) exists and is unique. In order to state our weak error estimates, we denote C p P (R d , R) the set of functions of class C p where all the partial derivatives have a polynomial growth, i.e. for each partial derivative φ up to order p, there exist C > 0 and

r ∈ N such that |φ(x)| ≤ C(1 + |x| r ). ( 13
)
We shall prove in Section 3 the following second order weak error estimate for the semi-discrete methods.

Theorem 2.3. Let T > 0. Assume (2) and (H). Consider X k the numerical solution of Algorithm 2.1 and X(t) the exact solution of (1). Then, for all φ ∈ C 6 P (R d , R), and all H = N ε with N ∈ N, and k ∈ N with kH ≤ T ,

|E(φ(X k )) -E(φ(X(kH)))| ≤ CH 2 ,
where the constant C is independent of ε, H, k, N .

We highlight that the weak accuracy estimate of Theorem 2.3 holds uniformly with respect to ε where both ε ≪ 1 and ε ≃ 1 are allowed. Observe in addition that Algorithm 2.1 is exact for N = 1 (i.e. the error is zero) because K k,1 (1) = X k in [START_REF] Cohen | On the numerical discretisation of stochastic oscillators[END_REF] and [START_REF] Cohen | A trigonometric method for the linear stochastic wave equation[END_REF] reduces to a time transformation t = t/ε of (1).

Remark 2.4. Notice that setting α N = 0, β N = 1 in [START_REF] Cohen | Convergence analysis of trigonometric methods for stiff secondorder stochastic differential equations[END_REF] in the definition of Algorithm 2.1 would yield the following multi-revolution method X k → X k+1 of weak order 1,

dK k = (AK k + Hf (K k ))dt + √ H m r=1 g r (K k )dW k,r , K k (0) = X k , X k+1 = K k (1).
Indeed, the error estimate of Theorem 2.3 with H 2 replaced by H remains valid for the above scheme for all φ ∈ C 4

P (R d , R) (observe in this case K k,2 (0) = K k,1 (1) = X k ).
In this paper, we shall however focus on the more accurate weak second order methods. The error analysis of the first order method can be performed analogously.

Fully-discrete multi-revolution composition methods

Algorithm 2.1 is called semi-discrete because each step requires the resolution of two systems of SDEs, whose solution is not known in general and needs to be approximated numerically. It was noted in the deterministic context [START_REF] Chartier | Multi-revolution composition methods for highly oscillatory differential equations[END_REF] that in principle, any nonstiff integrator could be used. However, a natural choice for such approximation is to use a splitting method where the oscillatory and non oscillatory parts of the problem are solved separately in an efficient way (sometimes exactly), as proposed and studied recently in [START_REF] Tao | Structure preserving stochastic impulse methods for stiff Langevin systems with a uniform global error of order 1 or 1/2 on position[END_REF] in the stochastic context of the Langevin equation.

We now formulate the fully-discrete multi-revolution stochastic methods for the class of problems [START_REF] Abdulle | Stabilized multilevel Monte Carlo method for stiff stochastic differential equations[END_REF]. It involves a micro stepsize h and a macro stepsize H where H ≫ ε is allowed. We highlight again that both H and h are of moderate size, independently of the smallness of the stiff parameter ε. The approach involves a weak integrator Y k+1 = Φ h (Y k ) with stepsize h for the nonstiff SDE problem

dY (t) = f (Y (t))dt + m r=1 g r (Y (t))dW r (t) (14)
where compared to (1), the stiff term ε -1 A has been removed. The integrator Φ h can be the exact solution if computationally available, or an efficient weak approximation.

Algorithm 2.5 (Fully-discrete S-MRCM). Consider a macro stepsize H = N ε and a micro stepsize h = 1/n with N, n ∈ N * . For the approximation of the flow after time H of (1), we consider the scheme X k+1 → X k defined by the composition

X k+1 = (e hA/2 • Φ β N Hh • e hA/2 ) n • (e -hA/2 • Φ α N Hh • e -hA/2 ) n (X k )
where the micro integrator Φ h is a weak integrator for [START_REF] Debrabant | Diagonally drift-implicit runge-kutta methods of weak order one and two for itô sdes and stability analysis[END_REF] with stepsize h, and α N , β N are defined in [START_REF] Cohen | Convergence analysis of trigonometric methods for stiff secondorder stochastic differential equations[END_REF]. The exponents n indicate that the maps are composed n times with independent random variables.

Notice that Algorithm 2.5 requires for each time step 2n applications of the nonstiff integrator Φ h , applied with independent random variables, and 2n evaluations of exponentials (using e hA/2 • e hA/2 = e hA ). The following theorem with proof postponed to Section 3 states that it has weak second order of accuracy with respect to the micro and macro stepsizes H, h, uniformly with respect to ε. To this aim, we assume that the integrator Φ h satisfies for all x ∈ R d , and all h ≤ h 0 ,

|E(Y k+1 -Y k |Y k = x)| ≤ C(1 + |x|)h, |Y k+1 -Y k | ≤ M k (1 + |Y k |) √ h, (15)
where C, M n are independent of h, x and M k is a random variable with finite moments of all orders. We further assume the following local weak order two estimate, for all φ ∈

C 6 P (R d , R), |E(φ(Φ h (Y 0 ))) -E(φ(Y (h)))| ≤ C(x)h 3 , (16)
where C(x) is independent h and has a polynomial growth [START_REF] De Bouard | Weak and strong order of convergence of a semidiscrete scheme for the stochastic nonlinear Schrödinger equation[END_REF].

Theorem 2.6. Let T > 0. Assume (2) and (H). Consider X k the numerical solution of Algorithm 2.5 and X(t) the exact solution of (1). Assume that the integrator Φ h for (14) satisfies (15) and [START_REF] Engquist | The heterogeneous multiscale methods[END_REF]. Then, for all φ ∈ C 6 P (R d , R), and all h = 1/n and H = N ε small enough with n, N, k ∈ N with kH ≤ T ,

|E(φ(X k )) -E(φ(X(kH)))| ≤ C(H 2 + h 2 ), ( 17
)
where

C is independent of ε, H, n, k, N, h.
The proofs of the above weak convergence estimates (Theorem 2.3 and Theorem 2.6) are provided in Section 3.

Remark 2.7. We highlight that (15) and ( 16) are natural assumptions for a weak order two integrator. Indeed, a classical theorem of Milstein [START_REF] Milstein | Weak approximation of solutions of systems of stochastic differential equations[END_REF] (see [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF]Chap. 2.2]) permits to deduce from the local error (16) the global error |E(φ(Y k )) -E(φ(Y (hk)))| ≤ Ch 2 for all hk ≤ T . Crucial is the assumption [START_REF] Debussche | Numerical simulation of focusing stochastic nonlinear Schrödinger equations[END_REF], which is easily satisfied by any reasonable integrator for [START_REF] Debrabant | Diagonally drift-implicit runge-kutta methods of weak order one and two for itô sdes and stability analysis[END_REF], and that automatically yield that Y k has bounded moments of any order for all kh ≤ T (see [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF]Lemma 2.2,p. 102]). Examples of such integrators with additional favorable geometric properties are presented in Section 4.

We end this section with the following remark, which shown that the multi-revolution approach has links with weak methods for systems of SDEs with small noise of the form

dZ(s) = a(Z(s))dt + εb(Z(s))dt + |ε| m r=1 c r (Z(s))dW r (s), Z(0) = Z 0 ( 18
)
as proposed and analyzed in [START_REF] Milstein | Numerical methods in the weak sense for stochastic differential equations with small noise[END_REF] (see also [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF]Chap. 3]). Such schemes applied to [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF] have global weak order O(h p + εh q ) where p < q on bounded time intervals, using appropriate smoothness assumptions on the vector fields a, b, c r .

Remark 2.8. Setting a(Z) = AZ, b(Z) = f (Z), σ r = g r , Z 0 = X 0 , the system (18) is equivalent to (1)
via the time transformation s = t/ε. In this case, the multi-revolution approach permits to approximate Z(s) on longer time intervals with s = O(ε -1 ) with a computational cost and an accuracy both independent of ε. Indeed, following the lines of the proof of Theorem 2.6, if Ψ h,ε is an integrator for [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF] with weak order O(h p + εh q ), then, provided the numerical moments remain uniformly bounded, the multi-revolution composition scheme

X k+1 = (Ψ h,β N H ) n • (Ψ -h,-α N H ) n (X k ),
with α N , β N defined in [START_REF] Cohen | Convergence analysis of trigonometric methods for stiff secondorder stochastic differential equations[END_REF] and h = 1/n, can be shown to approximate the solution X(t) of (1) at time t = kH with H = N ε with global weak error O(H 2 + H -1 h p + h q ) for all kH ≤ T (equivalently to approximate (18) at time s = kN ≤ ε -1 T ), where the constant in O is independent of H, h, ε.

Weak convergence analysis

A crucial ingredient for the analysis is the variation of constant formula: the solution X(t) of (1) satisfies

X(t) = e tε -1 A X 0 + t 0 e (t-s)ε -1 A f (X(s))ds + t 0 e (t-s)ε -1 A m r=1 g r (X(s))dW r (s). ( 19
)
In other words, considering the change of variables

Y (t) = e -ε -1 tA X(t), we have that Y (t)
is the solution of the non-autonomous SDE problem

dY (t) = e -ε -1 tA f (e ε -1 tA Y (t))dt + m r=1 e -ε -1 tA g r (e ε -1 tA Y (t))dW r (t). (20)
A feature of the semi-discrete Algorithm 2.1 is that it is exact in the case of additive noise (4) and in the absence of the nonlinearity (f = 0), as stated in the following proposition. This is a consequence of the homogeneity and scaling in time properties of Wiener processes. Notice in contrast that the stochastic trigonometric methods [START_REF] Cohen | On the numerical discretisation of stochastic oscillators[END_REF][START_REF] Cohen | Convergence analysis of trigonometric methods for stiff secondorder stochastic differential equations[END_REF] are not exact in this case.

Proposition 3.1. Consider the numerical solution of (1) with f = 0 and additive noise (4). Then, Algorithm 2.1 is exact in the sense that X(t), t = kN ε and X k have the same law of probability for all ε and all N, k, ∈ N.

Proof. The solutions of ( 10) and ( 11) can be expressed using the variation of constant formula [START_REF] Grébert | On the energy exchange between resonant modes in nonlinear Schrödinger equations[END_REF], which yields, using e A = e -A = I,

X 1 = X 0 + √ H 1 0 e -sA Bd W (s)
where we notice that W (s

) := √ α N W 1 (s)+ √ β N W 2 (s) is a standard m-dimensional Wiener process because α N + β N = 1 (recall that W 1
and W 2 are assumed independent). In comparison, using [START_REF] Grébert | On the energy exchange between resonant modes in nonlinear Schrödinger equations[END_REF] and the change of variable ŝ = ε -1 s, and the periodicity of t → e tA with period 1, the exact solution of (1) with f = 0 and (4) satisfies

X(H) = X 0 + N 0 e -sA BdW (εs) = X 0 + 1 0 e -sA B N -1 k=0 dW (εs + k).
Using standard homogeneity and scaling in time properties of the Wiener process, we have that

H -1 N -1 k=0 (W (εs + k) -W (k)
) is also a standard Wiener process. Using the independence of W k,1 , W k,2 , k = 0, 1, 2, . . . in Algorithm 2.1, we obtain that X(H) has the same law of probability as X 1 , and we conclude the proof that X(t k ) and X k have the same law of probability by induction on k.

Using (2), we observe that for integer multiples of the oscillatory period, the solutions of ( 1) and ( 20) coincide: X(N ε) = Y (N ε) for all N ∈ N and all ε. The advantage of considering the form (20) compared to [START_REF] Abdulle | Stabilized multilevel Monte Carlo method for stiff stochastic differential equations[END_REF] is that the drift and diffusion functions (x, t) → e -ε -1 tA f (e ε -1 tA Y (t)) and (x, t) → e -ε -1 tA g r (e ε -1 tA Y (t)) are C 6 -functions with all partial derivatives with respect to the spatial variable x bounded uniformly with respect to ε. Using these regularity, we may recall the following two classical results taken from [START_REF] Talay | Discrétisation d'une équation différentielle stochastique et calcul approché d'espérances de fonctionnelles de la solution[END_REF]. Lemma 3.2. [33, Lemma 1]. Assume (H) and consider the solution Y (t) of (20). Then, there exist a constant C > 0 and r ∈ N such that for all p ∈ N, t

∈ [0, T ], X 0 ∈ R d , E(|X(t)| p ) ≤ C(1 + |X 0 | r ).
In addition, there exists a version of the process X(t) which is almost surely 6 times continuously differentiable with respect to the initial condition X(0) = X 0 , and the derivatives with respect to the initial condition,

(∂ k X(t))/(∂X k 0 ), k = 1, . . . , 6, have bounded moments uniformly with respect to ε ∈ R, t ∈ [0, T ], X 0 ∈ R d . Theorem 3.3. [33, Theorem 2.2]
. Assume (H) and consider the solution Y (t) of (20). For T > 0 and all φ ∈ C 6

P (R d , R), the function v(x, t) := E(φ(Y (t))|Y (T ) = x), 0 ≤ t ≤ T , is solution of the Backward Kolmogorov partial differential equation ∂v ∂t + L ε (-t)v = 0, v(x, T ) = φ(x),
where the generator of (20) is defined by

3 L ε (-t)φ(x) = e -ε -1 tA f (e ε -1 tA x) • ∇φ(x) + 1 2 m r=1 φ ′′ (x) e -ε -1 tA g r (e ε -1 tA x), e -ε -1 tA g r (e ε -1 tA x) (21)
and the partial derivatives

∂ i 0 t ∂ x i 1 • • • ∂ x i k v(x, t) with 2i 0 + k ≤ 6, 1 ≤ i j ≤ d, i 0 , k ≥ 0, are continuous with polynomial growth (13) where r, C are independent of ε, x ∈ R d , and t ∈ [0, T ]. In particular, u(x, t) := E(φ(Y (t))|Y (0) = x) (22)
is equal to v(x, T -t), and assuming T /ε ∈ N and (2), we obtain that u(x, t) in [START_REF] Hansen | High order splitting methods for analytic semigroups exist[END_REF] is the solution of

∂u ∂t = L ε (t)u, u(x, 0) = φ(x), (23)
for all t ∈ [0, T ].
We deduce from Theorem 3.3 a weak Taylor expansion for u(x, H) = E(φ(X(H))) where X(t) is the solution of (1). Proposition 3.4. Assume the hypotheses of Theorem 3.3. Then, u(x, t) defined in [START_REF] Hansen | High order splitting methods for analytic semigroups exist[END_REF] satisfies for all

H = N ε, with N ∈ N, u(x, H) = φ(x) + H 1 0 L 1 (s 1 )φ(x)ds 1 + H 2 -Hε 2 1 0 1 0 L 1 (s 1 )L 1 (s 2 )φ(x)ds 2 ds 1 + Hε 1 0 s 1 0 L 1 (s 1 )L 1 (s 2 )φ(x)ds 2 ds 1 + O(H 3 )
where the constant in O(H 3 ) is independent of N, ε with a polynomial growth (13) with respect to x, and L 1 is defined in [START_REF] Hairer | Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF] with ε = 1.

Proof. Iterating "á la Picard" the integral relation u(x, t) = φ(x)

+ t 0 L ε (s)u(x, s)ds, we obtain u(x, t) = φ(x) + t 0 L ε (s 1 )φ(x)ds 1 + t 0 s 1 0 L ε (s 1 )L ε (s 2 )φ(x)ds 2 ds 1 + t 0 s 1 0 s 2 0 L ε (s 1 )L ε (s 2 )L ε (s 3 )u(x, s 3 )ds 3 ds 2 ds 1 which yields using Theorem 3.3, u(x, H) = φ(x) + H 0 L ε (s 1 )φ(x)ds 1 + H 0 s 1 0 L ε (s 1 )L ε (s 2 )φ(x)ds 2 ds 1 + O(H 3 )
where the constant in O(H 3 ) satisfies the polynomial growth [START_REF] De Bouard | Weak and strong order of convergence of a semidiscrete scheme for the stochastic nonlinear Schrödinger equation[END_REF] uniformly with respect to N, ε. Using the change of variable ŝi = ε -1 s i , we deduce

u(x, H) = φ(x) + ε N 0 L 1 (s 1 )φ(x)ds 1 + ε 2 N 0 s 1 0 L 1 (s 1 )L 1 (s 2 )φ(x)ds 2 ds 1 + O(H 3 ).
Finally, using the periodicity assumption (2) yields the identities

N 0 L 1 (s 1 )φ(x)ds 1 = N 1 0 L 1 (s 1 )φ(x)ds 1 , N 0 s 1 0 L 1 (s 1 )L 1 (s 2 )φ(x)ds 2 ds 1 = N (N -1) 2 1 0 1 0 L 1 (s 1 )L 1 (s 2 )φ(x)ds 2 ds 1 + N 1 0 s 1 0 L 1 (s 1 )L 1 (s 2 )φ(x)ds 2 ds 1 ,
which permits to conclude the proof.

A consequence of Proposition 3.4 is the following remark, which gives some insight on a stochastic effective problem in the case of additive noise, analogously to (8).

(R d , R), |E(φ(X(t k ))) -E(φ(X(t k )))| ≤ Cε 2 ,
where C is independent of N, ε, where X(t) solves the effective SDE

dX = a 1 + ε(a 2 + σ 2 ∆b 2 ) (X)dt + σ(I + εb ′ 2 (X))dW (t), X(0) = X 0 ( 24
)
where a 1 , a 2 are defined in (9) and

b 2 (x) = 1 4 1 0 s 1 e s 1 A f (e -s 1 A x) - s 1 0 e s 2 A f (e -s 2 A x)d s 2 ds 1 .
The proof of ( 24) is deduced showing first the local estimate |E(φ(X(ε))) -E(φ(X(ε)))| ≤ Cε 3 (using Proposition 3.4 with N = 1 for X(ε) and analogously a weak Taylor expansion for X(ε)). Interestingly, the effective SDE (24) has multiplicative noise in general, although the original SDE (1) with B = I in (4) has additive noise.

We may now derive a local weak error estimate for the semi-discrete S-MRCM.

Lemma 3.6. Assume (2) and (H). Consider X 1 the numerical solution of Algorithm 2.1 after one step and X(t) the exact solution of (1). Then, for all initial condition

X 0 = x, all φ ∈ C 6 P (R d , R), and all H = N ε with N ∈ N, |E(φ(X 1 )) -E(φ(X(H)))| ≤ C(x)H 3 ,
where the constant C(x) is independent of ε, H and has a polynomial growth [START_REF] De Bouard | Weak and strong order of convergence of a semidiscrete scheme for the stochastic nonlinear Schrödinger equation[END_REF].

Proof. We observe that the SDE (10) has the form (1) by setting ε = -1 and replacing f by α N Hf and g r by √ α N Hg r . Applying Proposition 3.4 with N = 1, we deduce

E(φ(K 1 (1))|K 1 (0) = x) = φ(x) + α N H 1 0 L -1 (s 1 )φ(x)ds 1 + (α N H) 2 1 0 s 1 0 L -1 (s 1 )L -1 (s 2 )φ(x)ds 2 ds 1 + O(H 3 ). ( 25
)
Analogously for the solution K 2 (t) of ( 11), we have

E(φ(K 2 (1))|K 2 (0) = x) = φ(x) + β N H 1 0 L 1 (s 1 )φ(x)ds 1 + (β N H) 2 1 0 s 1 0 L 1 (s 1 )L 1 (s 2 )φ(x)ds 2 ds 1 + O(H 3 ), ( 26
)
where the constants in the above remainders O(H 3 ) are independent of N, ε and have a polynomial growth [START_REF] De Bouard | Weak and strong order of convergence of a semidiscrete scheme for the stochastic nonlinear Schrödinger equation[END_REF] with respect to x. We next have by the Fubini theorem and using properties of conditional expectancies (recall that W 1 , W 2 in (10), [START_REF] Cohen | A trigonometric method for the linear stochastic wave equation[END_REF] are independent Wiener processes),

E(φ(X 1 )|X 0 = x) = E W 1 E W 2 φ(K 2 (1))|K 2 (0) = K 1 (1) |K 1 (0) = x (27)
where the notations E W 1 , E W 2 refers to the expectation with respect to the Wiener process W 1 , W 2 , respectively. Applying [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF] with

φ replaced by x = K 2 (0) → E W 2 (φ(K 2 (1))),
which is, by Theorem 3.3, almost surely of class C 6 with derivatives of polynomial growth, we deduce for all X 0 = x the estimate 4

E(φ(X 1 )) = I + α N H 1 0 L -1 (s 1 )ds 1 + (α N H) 2 1 0 s 1 0 L -1 (s 1 )L -1 (s 2 )ds 2 ds 1 • I + β N H 1 0 L 1 (s 1 )ds 1 + (β N H) 2 1 0 s 1 0 L 1 (s 1 )L 1 (s 2 )ds 2 ds 1 φ(x) + O(H 3 )
4 Observe that the ordering of composition is the opposite compared to the ordering [START_REF] Cohen | On the numerical discretisation of stochastic oscillators[END_REF], [START_REF] Cohen | A trigonometric method for the linear stochastic wave equation[END_REF] in Algorithm 2.5. This effect is known as the "Vertauschungssatz" of Gröbner in the deterministic literature, see e.

= I + α N H 1 0 L -1 (s 1 )ds 1 + β N H 1 0 L 1 (s 1 )ds 1 + α N β N H 2 1 0 L -1 (s 1 )ds 1 1 0 L 1 (s 1 )ds 1 + (α N H) 2 1 0 s 1 0 L -1 (s 1 )L -1 (s 2 )ds 2 ds 1 + (β N H) 2 1 0 s 1 0 L 1 (s 1 )L 1 (s 2 )ds 2 ds 1 φ(x) + O(H 3 ).
A consequence of ( 2) is

1 0 L -1 (s 1 )ds 1 = 1 0 L 1 (s 1 )ds 1 and 1 0 s 1 0 L -1 (s 1 )L -1 (s 2 )ds 2 ds 1 = 1 0 1 0 L 1 (s 1 )L 1 (s 2 )ds 2 ds 1 - 1 0 s 1 0 L 1 (s 1 )L 1 (s 2 )ds 2 ds 1 .
We deduce

E(φ(X 1 )) = I + (α N + β N )H 1 0 L 1 (s 1 )ds 1 + (α N β N + α 2 N )H 2 1 0 1 0 L 1 (s 1 )L 1 (s 2 )ds 2 ds 1 + (β 2 N -α 2 N )H 2 1 0 s 1 0 L 1 (s 1 )L 1 (s 2 )ds 2 ds 1 φ(X 0 ) + O(H 3 ).
Using [START_REF] Cohen | Convergence analysis of trigonometric methods for stiff secondorder stochastic differential equations[END_REF], we observe that

α N + β N = 1, (α N β N + α 2 N )H 2 = (H 2 -Hε)/2 and (β 2 N - α 2
N )H 2 = Hε. Comparing the above Taylor expansion with the one in Proposition 3.4 permits to conclude the proof.

Based on the local error estimate of Lemma 3.6, we may now give the proof of Theorem 2.3 for the global order two of convergence in the weak sense. To this aim, the following lemma, which is shown in the proof of [25, Lemma 2.2, p. 102] is a crucial ingredient. Lemma 3.7. [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF] Consider a discrete process {Y k } (k ∈ N) satisfying [START_REF] Debussche | Numerical simulation of focusing stochastic nonlinear Schrödinger equations[END_REF]. Assume further that Y 0 has finite moments of all orders. Then, for all integer p, there exists C p such that

E(|Y k | p ) ≤ e Cpkh (E(|Y 0 | p ) + 1)
for all k ∈ N and h ≤ h 0 .

Proof of Theorem 2.3. We use a well known result of Milstein [START_REF] Milstein | Weak approximation of solutions of systems of stochastic differential equations[END_REF] (see [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF]Chap. 2.2]) which permits to deduce automatically the global error estimate of Theorem 2.3 from the local error estimate of Lemma 3.6. To this aim, we only have to check is that the moments E(|X k | 2r ) of the numerical solution are bounded for all k, H with 0 ≤ kH ≤ T uniformly with respect to k and H sufficiently small. This is a consequence of Lemma 3.7 applied with h = H to the discrete process {Y k } defined by Y 2k = X k and Y 2k+1 = K k,1 (1) where K k,1 is given in [START_REF] Cohen | On the numerical discretisation of stochastic oscillators[END_REF]. Considering the SDEs [START_REF] Cohen | On the numerical discretisation of stochastic oscillators[END_REF] and [START_REF] Cohen | A trigonometric method for the linear stochastic wave equation[END_REF], the estimates [START_REF] Debussche | Numerical simulation of focusing stochastic nonlinear Schrödinger equations[END_REF] with h = H for the discrete flows K k,1 (0) → K k,1 (1) and K k,2 (0) → K k,2 (1) are a straightforward consequence of the Lischitzness of f, g r , r = 1, . . . , m. This permits to conclude the proof.

In order to prove Theorem 2.6 for the the global weak accuracy of the fully-discrete method, we need to check that the numerical solution has uniformly bounded moments of any order. This is the purpose of the following lemma which shows that the numerical moments are uniformly bounded in spite of the oscillatory terms e ±hA/2 involved in the scheme.

Lemma 3.8. Assume the hypotheses of Theorem 2.6 and consider the numerical solution X k of Algorithm 2.5. Then, for all p ∈ N, there exists a constant C p such that for all

H = N ε, k ∈ N with kH ≤ T , E(|X k | p ) ≤ C p .
Proof. We observe that the numerical solution X k+1 is calculated inductively from X k using a composition of 2n stochastic mappings of the form

e δhA/2 • Φ γhH • e δhA/2
where δ = -1, γ = α N or δ = 1, γ = β N . Let {Y κ }, κ = 0, 1, 2, . . . denote the discrete process arising from this composition (note that X k = Y 2nk ). We next define the discrete process Y κ from Y κ by

Y 2nk+j = e jhA Y 2nk+j , Y 2nk+n+j = e -jhA Y 2nk+n+j , , j = 1, . . . , n. For 1 ≤ j ≤ n, we obtain Y 2nk+j+1 = e h(j+1/2)A • Φ α N hH • e -h(j+1/2)A ( Y 2nk+j ), which yields Y 2nk+j+1 -Y 2nk+j = e h(j+1/2)A (Φ α N hH (e -h(j+1/2)A Y 2nk+j ) -e -h(j+1/2)A Y 2nk+j ).
Analogously, for n + 1 ≤ j ≤ 2n, the above identity holds with A replaced by -A and α N replaced by β N . Using the estimates (15) on Φ h , the inequality |e λA ξ| ≤ c|ξ| for all λ ∈ R, ξ ∈ R d , and the bounds α N ≤ 1, β N ≤ 1 we deduce for all κ,

|E( Y κ+1 -Y κ | Y κ = x)| ≤ C(1 + |x|)Hh, | Y κ+1 -Y κ | ≤ M κ (1 + | Y κ |) √ Hh,
where C is independent of N, n, H, h, ε, κ, k. Applying Lemma 3.7 then yields that Y κ has uniformly bounded moments of any order for κHh ≤ 2T . We conclude the proof using h = 1/n and X k = Y 2nk for all k.

We shall also need the following result, whose proof follows standard arguments. For the sake of completeness, a proof is provided in Appendix. Lemma 3.9. Assume the hypotheses of Theorem 2.6. Let δ = -1, γ = α N (resp. δ = 1, γ = β N ). Then, the integrator

Y k+1 = e δhA/2 • Φ γHh • e δhA/2 (Y k ),
applied the SDE (10) (respectively [START_REF] Cohen | A trigonometric method for the linear stochastic wave equation[END_REF]) satisfies for all φ ∈ C 6 P (R d , R) and all H, h = 1/n small enough,

|E(φ(Y n )) -E(φ(K j (1)))| ≤ C(x)h 2 H,
where j = 1 (resp. j = 2) and C(x) is independent of h, H, N, n, k, ε and satisfies [START_REF] De Bouard | Weak and strong order of convergence of a semidiscrete scheme for the stochastic nonlinear Schrödinger equation[END_REF].

Proof of Theorem 2.6. Consider the numerical solution denoted Xk of the semi-discrete Algorithm 2.1. Applying Lemma 3.9 for [START_REF] Cohen | On the numerical discretisation of stochastic oscillators[END_REF] and [START_REF] Cohen | A trigonometric method for the linear stochastic wave equation[END_REF], we deduce

|E(φ( X1 )) -E(φ(X 1 ))| ≤ C(x)h 2 H.
Using Lemma 3.6, we deduce the local error estimate

|E(φ(X 1 )) -E(φ(X(H)))| ≤ C(x)(H 3 + h 2 H). (28)
In addition, by Lemma 3.8, the numerical solution X k as uniformly bounded moments of any order, and the Milstein theorem [START_REF] Milstein | Weak approximation of solutions of systems of stochastic differential equations[END_REF] (see [25, Chap. 2.2]) then yields the global error estimate of Theorem 2.6.

Numerical experiments

In this final section, we consider four problems to illustrate numerically the weak convergence rates of the proposed methods and the versatility of the approach. We investigate not only the accuracy of the methods compared to other oscillatory integrators but also their long time qualitative behavior. We also compare our approach with other known highly oscillatory integrators.

An illustrative example: Kubo oscillators

We first focus on the so-called Kubo oscillators. A nonlinear version with multiplicative Stratonovitch noise, as considered recently in [START_REF] Cohen | On the numerical discretisation of stochastic oscillators[END_REF], is

dQ(t) = (-ε -1 P (t) + P f (P (t), Q(t)))dt + σP (t) • dW (t) dP (t) = (ε -1 Q(t) -Qf (P (t), Q(t)))dt -σQ(t) • dW (t) (29)
where P (t), Q(t) ∈ R and the same one-dimensional Wiener process is considered in (29) for both components. Notice that (29) can be put in the form (1) using the time transformation t = t/2π and the Stratonovitch to Itô conversion formula. For the nonstiff part of the system (29), we choose as micro integrator a Strang splitting method

Φ h = Φ V h/2 • Φ W h • Φ V h/2 , (30)
where, on the one hand, the noise part dQ = σP • dW, dP = -σQ • dW is integrated with

Φ V h : q p → Q P = cos(σ √ hξ) sin(σ √ hξ) -sin(σ √ hξ) cos(σ √ hξ) q p
where ξ are independent random variables with P(ξ = ± √ 3) = 1/6, P(ξ = 0) = 2/3. Notice that replacing above √ hξ by ∆W k = W ((n + 1)k) -W (nk) would make Φ V h exact. The above choice of discrete random variable makes Φ V h a weak second order integrator (using E(ξ 2 ) = 1, E(ξ 4 ) = 3). On the other hand, the nonlinear part dQ = σP f (P, Q)dt, dP = -σQf (P, Q)dt is approximated by the implicit midpoint rule defined as Φ W h :

q p → Q P = q + hP M f (P M , Q M ) p -hQ M f (P M , Q M )
where P M = (p + P )/2, Q M = (q + Q)/2. In the implementation of Φ W h , we use fixed point iterations until convergence up to round-off errors. Notice that the splitting scheme [START_REF] Sheng | Solving linear partial differential equations by exponential splitting[END_REF] has weak second order of accuracy (see Proposition 6.1 in Appendix with ε = 1). Conservation of quadratic first integrals An interesting feature of ( 29) is that the quantity C(p, q) = p 2 + q 2 is exactly conserved along time for all trajectories. Precisely, observing that d C(P (t), Q(t)) = 0 (as shown for instance in [10, Prop. 3.2]), we have almost surely

P (t) 2 + Q(t) 2 = P (0) 2 + Q(0) 2 , for all t > 0. ( 31
)
Since both integrators Φ V h and Φ W h exactly conserve the quadratic first integral, we have that Φ h in (30) exactly conserves P (t) 2 + Q(t) 2 , and thus also the corresponding fullydiscrete S-MRCM of Algorithm 2.5:

P 2 k + Q 2 k = P 2 0 + Q 2 0
for all k ∈ N and all numerical trajectories of the method. Note in contrast that the strong method proposed in [START_REF] Cohen | On the numerical discretisation of stochastic oscillators[END_REF] does not conserve exactly this quadratic first integral in general.

Weak convergence rates Since the scheme exactly conserves the quadratic first integral P 2 +Q 2 , we have that the weak error bound [START_REF] Engquist | Heterogeneous multiscale methods: a review[END_REF] holds provided that f : R 2 → R in ( 29) is of class C 6 in a neighbourhood of {(p, q) ; p 2 +q 2 = P 0 +Q 0 } where we consider a deterministic initial condition (here Q 0 = 1, P 0 = 0). We consider the nonlinearity f (p, q) = p 3 + q 5 (similarly to [START_REF] Cohen | On the numerical discretisation of stochastic oscillators[END_REF]). For ε = 2 -6 ≃ 1.6 • 10 -2 (left pictures) and ε = 2 -8 ≃ 3.9 • 10 -3 (right pictures), we plot for many different macro steps H the error for the second moment of the first component E(|Q(T )| 2 ) at the final time T = 2π as a function of the number of macro steps T /H. The expectation is approximated using the average over M = 10 7 trajectories to make the Monte-Carlo error sufficiently small compared to the weak accuracy of the methods. In the top pictures, we take n = 1024 micro steps in each macro step, so that the micro discretization can be considered as nearly exact (see semi-discrete Algorithm 2.1). In the bottom pictures, the four lines correspond respectively to n = 8, 16, 32, 64 micro steps per macro step (from top to bottom lines). In the top pictures, we observe the expected lines of slope 2, as proved in the semi-discrete error analysis of Theorem 2.3. In the bottom pictures, we observe the expected lines of slope 2 only for a sufficiently fine micro stepsize, as shown in Theorem 2.6. As a reference solution, we consider here the standard Strang splitting e ε -1 hA/2 • Φ h • e ε -1 hA/2 with Φ h defined in [START_REF] Sheng | Solving linear partial differential equations by exponential splitting[END_REF] with small stepsize h = 2 -15 T ≃ 1.9 • 10 -4 .

A test problem with non-commutative noise

In some situations (e.g. the stochastic nonlinear Schrödinger model [START_REF] Chartier | Multi-revolution composition methods for highly oscillatory differential equations[END_REF] studied in Section 4.4), the exact solution of the non oscillatory system ( 14) is available and easy to compute. Otherwise a weak second order approximation is needed in general. Notice that already for nonstiff stochastic problems and without structure assumptions (i.e. for a non-commutative noise), methods with high strong order are generally costly to simulate because they involve multiple stochastic integrals, such as h t 0 s 2 0 dW q (s 1 )dW r (s 2 ), which are costly to approximate strongly for q = r. This is however not true for weak methods where such multiple stochastic integrals can be approximated efficiently in a weak sense using appropriate discrete random variables. The aim of this section is to illustrate that this is again not a difficulty in our highly oscillatory context. We consider the following nonlinear oscillatory problem in dimension d = 2, which is a modification of a scalar test SDE from [START_REF] Debrabant | Diagonally drift-implicit runge-kutta methods of weak order one and two for itô sdes and stability analysis[END_REF], with a non-commutative Itô noise with m = 10 independent driving Wiener processes,

dQ(t) = -ε -1 P (t)dt, Q(0) = 1, ( 32
)
dP (t) = ε -1 Q(t)dt + 10 j=1 a -1 j P (t) 2 + Q(t) 2 + b -1 j (1 -Q(t))dW j (t), P (0) = 0.
where P (t), Q(t) ∈ R. The values of the constants a j , j = 1, . . . , 10 are respectively 5, 5, 10, 15, 30, 15, 10, 5, 10, 15, and the values of b j , j = 1, . . . , 10 are respectively 4, 3, 5, 2, 1, 2, 4, 5, 10, 10. Considering the averaged oscillatory energy E = E P 2 + Q 2 , an application of the Itô formula yields dE(t) dt = 10 j=1 a -2 j (E(t) + b -1 j (1 -cos(t/ε))), which can be solved analytically as

E P (t) 2 + Q(t) 2 = e at + b(a + a 3 ε 2 ) -1 e at + a 2 ε 2 cos(t/ε) -aε sin(t/ε) -a 2 ε 2 -1 ,
where a = 10 j=1 a -2 j = 37/225, b = 10 j=1 a -2 j b -1 j = 257/6000. We shall use this formula to check the accuracy of the fully-discrete S-MRCM (Algorithm 2.5).

For the integrator X 1 = Φ h (X 0 ) needed to integrate the non-stiff part of problem [START_REF] Talay | Efficient numerical schemes for the approximation of expectations of functionals of the solution of a SDE and applications[END_REF], which has the form ( 14) with f = 0, we compare two different schemes which fulfil the assumptions of Theorem 2.6. On the one hand we simply use the Euler-Maruyama method of weak order 1,

X 1 = X 0 + √ h m r=1 g r (X 0 )ξ r ,
where ξ r are independent discrete random variables satisfying P(ξ r = ± √ 3) = 1/6, P(ξ r = 0) = 2/3. On the other hand, we consider the following Runge-Kutta type scheme of weak order two, which is a derivative free version of the Milstein-Talay method [START_REF] Talay | Efficient numerical schemes for the approximation of expectations of functionals of the solution of a SDE and applications[END_REF] derived in Solid lines: method (33) with weak order 2 for the nonstiff part. Dashed lines: Euler-Maruyama method with weak order 1 for the nonstiff part. We use in both cases the multi-revolution parameter N = 128 and N = 256 (from top to bottom).

[3, Lemma 3.1] using an idea in [START_REF] Rößler | Second order Runge-Kutta methods for Itô stochastic differential equations[END_REF] to make the number of evaluations of each diffusion function g r , r = 1, . . . , m independent of the dimensions d, m,

X 1 = X 0 + 1 2 m r=1 g r X 0 + m q=1 g q (X 0 )J q,r -g r X 0 - m q=1 g q (X 0 )J q,r + √ h 2 m r=1 g r (X 0 + h 2 m q=1 g q (X 0 )χ q ) + g(X 0 - h 2 m q=1 g q (X 0 )χ q ) ξ r , (33)
where the quantity J q,r is a weak approximation of the multiple stochastic integral h 0 s 2 0 dW q (s 1 )dW r (s 2 ) given by [25, p. 96, eq. (1.25)]

J q,r =    h(ξ r ξ r -1)/2, if q = r, h(ξ q ξ r -χ q )/2, if r < q, h(ξ q ξ r + χ r )/2, if r > q,
and χ l , ξ l , l = 1 . . . m are independent discrete random variables satisfying P(χ l = ±1) = 1/2, P(ξ l = ± √ 3) = 1/6, P(ξ l = 0) = 2/3. In Figure 2, we consider Algorithm 2.5 applied to problem [START_REF] Talay | Efficient numerical schemes for the approximation of expectations of functionals of the solution of a SDE and applications[END_REF] with final time T = 2π, and oscillatory frequency ε = 2 -8 . We plot the error in the averaged energy E(P (T ) 2 + Q(T ) 2 ) ≃ 3.2816345 as a function of the total number of micro steps (left picture) or as a function of the number of evaluations of the diffusion functions g r , r = 1, . . . , m (right picture), taking into account that method (33) requires 5 diffusion function evaluations per step in contrast to the Euler-Maruyama method with a single evaluation of the diffusion functions per step. Since these evaluations dominate the cost of the method for a large dimension m, this is a natural measure of the computational cost. The averages are obtained using 10 8 independent trajectories. We use successively the multi-revolution parameters N = 128 and N = 256, and obtain in Figure 2 the expected lines of slope 2 and 1, except for large numbers of micro steps where the macro error (and also the Monte-Carlo error) is not negligible compared to the micro error. This confirms the error analysis in Theorem 2.6.

We observe in the right picture of Figure 2 that the S-MRCM using a weak second order nonstiff method is more efficient than the S-MRCM using the Euler-Maruyama method in the case where high accuracy is needed (finer than one percent in this example).

A modification of the Fermi-Pasta-Ulam problem with additive noise

We next consider a problem inspired from [START_REF] Hairer | Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF], which is a single-frequency modification of the Fermi-Pasta-Ulam problem often used to test methods for highly-oscillatory problems. It corresponds to a second order Hamiltonian system of the form (3) with d = 6 degrees of freedom and additive noise (4),

dQ = ∇ P H(P, Q)dt, dP = -∇ Q H(P, Q) + σBdW (t). ( 34
)
Its Hamiltonian function is defined by

E(p, q) = 1 2 6 i=1 p 2 i + 1 2ε 2 6 i=4 q 2 i + V (q), (35)
with the quartic interaction potential

V (q) = 1 4 ((q 1 -q 4 ) 4 + (q 2 -q 5 -q 1 -q 4 ) 4 + (q 3 -q 6 -q 2 -q 5 ) 4 + (q 3 + q 6 ) 4 ),
where ε > 0 is a small parameter. Again, following [START_REF] Chartier | Multi-revolution composition methods for highly oscillatory differential equations[END_REF] in the deterministic setting, this problem can be put in the form (1) up to a time transformation t = t/2π, where A ∈ R 12×12 is given by

A =     0 0 0 0 0 0 0 2πηI 0 0 0 0 0 -2πη -1 I 0 0     ,
where η = ε is fixed, and 0, I denote respectively zero and identity matrices of size 3 × 3. This model describes the motion of a chain composed with three soft nonlinear springs and three stiff linear springs (see [START_REF] Hairer | Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF] in the deterministic context). The components i = 1, 2, 3 are associated to a slow motion, while the components i = 4, 5, 6 oscillate rapidly at the frequency ε -1 . We consider the initial condition Q(0) = (1, 0, 0, ε, 0, 0) T , P (0) = (1, 0, 0, 1, 0, 0) T , and similarly to [START_REF] Cohen | Convergence analysis of trigonometric methods for stiff secondorder stochastic differential equations[END_REF], we consider an additive noise perturbation (4) in dimension m = 1 with

B = σ(0, 1, 1/2, 5, 0, 1/100) T . (36)
Remark 4.1. Our stochastic error analysis is performed for simplicity for a constant matrix A independent of ε. However, analogously to the deterministic case in [START_REF] Chartier | Multi-revolution composition methods for highly oscillatory differential equations[END_REF], it can be extended to the present situation where A depends on η = ε. Indeed, since a, b r in (3) are independent of the highly oscillatory velocities, we have that the unbounded factor η -1 involved in A never appears in e tA f (e -tA y) and e tA g r (e -tA y) where f, g r are the vector fields in the formulation (20) of the problem. Weak convergence rates We emphasize that our convergence analysis applies only to Lipschitz vector fields, which is not the case of problem ( 34)-(35). However, numerical experiments still exhibit the high weak order two of convergence of the method predicted in the Lipschitz case. We consider the fully-discrete Algorithm 2.5 where the micro integrator is defined as a Strang Splitting [START_REF] Sheng | Solving linear partial differential equations by exponential splitting[END_REF] where Φ V h integrates exactly the deterministic problem with Hamiltonian 1 2 3 i=1 p 2 i and Φ W h integrates the position components (see [START_REF] Chartier | Multi-revolution composition methods for highly oscillatory differential equations[END_REF] for details on the implementation). In Figure 3, we plot the convergence curves of the S-MRCM for average Hamiltonian energy E(E) at final time T = 2π for ε = 2 -6 (left picture) and ε = 2 -8 (right picture), respectively. We plot the errors at time T = 2π as a function of the number of macro steps T /H where the macro stepsize is H = 2πN ε (taking into account the time transformation t = t/(2π)). The 3 curves in each plot correspond respectively to n = 8, 16, 32 micro steps per macro steps (from top to bottom). We observe convergence curves with slope 2, which corroborate the fully-discrete error analysis of Theorem 2.6. The expectations are obtained as the averages over 10 8 trajectories of the methods, and the reference solution is computed using a standard Strang splitting with stepsize h = 2 -16 T ≃ 10 -4 .

Long time behavior Although our error analysis applies only to bounded time intervals with respect to ε, we investigate in Figure 4 the long time behavior of the S-MRCM on long time intervals of size T = 2πε -1 with ε = 1/200. We consider the evolution along time of the three stiff spring energies

I j = 1 2 p 2 3+j + 1 2ε 2 q 2 3+j , j = 1, 2, 3
the total oscillatory energy I = I 1 + I 2 + I 3 and the Hamiltonian energy in (35). In the deterministic case, it is classical (see e.g. [START_REF] Hairer | Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF]) that the Hamiltonian of the exact solution is exactly conserved, while I(t) = I(0) + O(ε) is nearly conserved along time and called an adiabatic invariant. In addition energy exchanges at the time scale ε -1 can be observed between the stiff spring energies I 1 , I 2 , I 3 . This can be observed in Figure 4 (left pictures) for a reference solution (Strang splitting method with small stepsize h = 10 -3 ) and the MRCM of deterministic order 2 as first proposed in [START_REF] Chartier | Multi-revolution composition methods for highly oscillatory differential equations[END_REF]. In the additive noise case (where we take σ = 0.1 in (36)), there is still an average energy exchange between the stiff springs, as observed in [START_REF] Cohen | Convergence analysis of trigonometric methods for stiff secondorder stochastic differential equations[END_REF], and a linear drift (5) of the average Hamiltonian energy with slope 1 2 σB 2 is expected for the exact solution. We observe an excellent behavior of the proposed S-MRCM, where the expected energy linear drift with the expected slope (5) can be observed (compare the S-MRCM and reference solutions in right pictures). We took the macro parameter N = 10 and n = 8 micro steps per macro time steps. On the time interval (0, 2πε -1 ), this corresponds to nN -1 ε -2 = 3.2 • 10 4 forces evaluations of the potential force ∇V for each numerical trajectory. In comparison, for the reference solution with h ≃ 10 -3 (Strang splitting), we make 2πε -1 h -1 ≃ 1.3 • 10 6 force evaluations to achieve a satisfactory solution. In this experiment, all expectations are obtained using the averages over 10 3 trajectories. Notice that for a standard explicit integrator, such as the Euler-Maruyama method (not represented here), a super linear growth of the numerical Hamiltonian would be observed (see [START_REF] Cohen | Convergence analysis of trigonometric methods for stiff secondorder stochastic differential equations[END_REF]), yielding a completely wrong qualitative behavior.

The stochastic nonlinear Schrödinger equation

Although our analysis applies only to finite dimensional systems of SDEs, we finally consider the SPDE model [START_REF] Castella | Stroboscopic averaging for the nonlinear Schrödinger equation[END_REF] of the nonlinear Schrödinger equation with a multiplicative space-time noise,

i ∂u(x, t) ∂t = -∆u(x, t) + εV (x)|u(x, t)| 2 u(x, t) + σ √ εu(x, t) • Ẇ (x, t), t > 0, (37)
u(x, 0) = cos x + sin x, where V (x) = 2 cos(2x), σ = 10 -2 is a constant, and W (x, t) is a space-time noise which is white in time and meant in the Stratonovitch sense and spatially correlated in space. Problem (37) is posed for simplicity in dimension one of space on the interval (0, 1) with periodic boundary conditions. In the deterministic setting (σ = 0) this problem is analyzed theoretically in [START_REF] Grébert | On the energy exchange between resonant modes in nonlinear Schrödinger equations[END_REF] and considered in [START_REF] Castella | Stroboscopic averaging for the nonlinear Schrödinger equation[END_REF][START_REF] Chartier | Multi-revolution composition methods for highly oscillatory differential equations[END_REF] to illustrate highly oscillatory integrators, and it has a unique global solution in all Sobolev spaces H s for arbitrary s ≥ 0. We consider the S-MRCM (Algorithm 2.5) applied to the spectral formulation [START_REF] Chartier | Multi-revolution composition methods for highly oscillatory differential equations[END_REF], which yields to a system of SDEs with a multiplicative noise of Stratonovitch type in dimension 2ℓ. We observe that y(t) = exp -itα|y 0 | 2 -iβW (t) y 0 is the the exact solution of the complex scalar Stratonovitch SDE idy = α|y| 2 y + βu • dW, y(0) = y 0 for all real parameters α, β. It is thus natural to define the micro integrator Φ h of the S-MRCM algorithm by

Φ h = F ℓ • ψ h • F ℓ -1
where F ℓ denotes the Fast Fourier Transform and ψ h is given by

ψ h : (u j ) -ℓ<j≤ℓ ∈ C 2ℓ → exp -ihV (x j )|u j | 2 -iσ √ h √ ∆x χ j u j -ℓ<j≤ℓ
where x j = j∆x, ∆x = 1/(2ℓ) and χ j are independent random variables with P(χ j = ± √ 3) = 1/6, P(χ j = 0) = 2/3). Features of the considered scheme are that it is explicit and the complexity with respect to the space dimension parameter ℓ reduces to O(ℓ log ℓ) thanks to the Fast Fourier Transform algorithm. In addition, the mass ℓ<j≤ℓ |ξ j (t)| 2 , which is a quadratic first integral of [START_REF] Chartier | Multi-revolution composition methods for highly oscillatory differential equations[END_REF], is exactly conserved by the numerical scheme (up to round-off errors) because the method is a composition of mass preserving flows.

Convergence rates For ε = 10 -4 and ℓ = 128 spatial Fourier modes we plot in Figure 5 the error in E(|ξ 1 (T )| 2 ) and E(|ξ 3 (T )| 2 ) (second moment of the first and third modes) at time T = 2πε -1 as a function of the number of macro steps, here 1, 2, 4, 8, 16, 20, 40, respectively. The number of micro steps per macro step is n = 100, and we take the averages over 10 5 trajectories. Again we observe the expected lines of slope 2. As reference solution, we take the numerical solution with N = 10, which corresponds to 10 3 macro steps on the time interval (0, T ). Long time behavior Again for ε = 10 -4 , we plot in Figure 6 the solution of S-MRCMs with parameter N = 100 for the first Fourier modes |ξ j (t)| for |j| = 1, 3, 5, . . . , 9 (notice that ξ j (t) = 0 for even indices j if σ = 0) along time in the interval (0, 10πε -1 ). In the top pictures, we consider the deterministic case (σ = 0) for which the Hamiltonian is exactly conserved by the exact solution (see the top right picture where the Hamiltonian error of the S-MRCM remains small, of size 10 -8 ). It was shown in [START_REF] Grébert | On the energy exchange between resonant modes in nonlinear Schrödinger equations[END_REF] that a nonlinear beating effect occurs for the first mode of the exact solution

|ξ 1 (t)| 2 = 1 + sin(2εt) 2 + O(ε 1/8 ), |ξ -1 (t)| 2 = 1 -sin(2εt) 2 + O(ε 1/8 ), t ≤ ε -9/8 .
This behavior of the exact solution is well reproduced by the deterministic MRCM as reported in [START_REF] Chartier | Multi-revolution composition methods for highly oscillatory differential equations[END_REF] (see top left picture). In the stochastic case of a multiplicative space-time noise (σ = 0.1), an analogous beating effect with period π/ε can still be observed in the first mode (see a sample trajectory in bottom pictures). Notice that in this stochastic context, the Hamiltonian energy is no longer a conserved quantity for the exact solution (see a sample trajectory in bottom right picture).

Comparison with other oscillatory integrators

In this section, we compare the efficiency of the proposed multi-revolution composition methods with other existing methods for the Kubo oscillator (Section 4.1) and the FPU type problem (Section 4.3). Precisely, we consider the following five methods:

• s-mrcm2: the multi-revolution methods of weak order two (Algorithm 2.5).

• s-mrcm1: the multi-revolution modification of weak order one discussed in Remark 2.4.

• splt2: the standard Strang splitting method between the highly oscillatory part and the nonstiff part, of weak order two,

X k+1 = e ε -1 hA/2 • Φ W h/2 • Φ V h • Φ W h/2 • e ε -1 hA/2 (X k ) (38)
where Φ V h integrates the deterministic nonlinear part and Φ W h integrates the noise, as described in Sections 4.1 and 4.3 for problems [START_REF] Sanders | Averaging methods in nonlinear dynamical systems[END_REF], [START_REF] Tao | Structure preserving stochastic impulse methods for stiff Langevin systems with a uniform global error of order 1 or 1/2 on position[END_REF], respectively.

• trig: the implicit trigonometric method proposed and studied numerically in [10, eq. ( 19)] for the Kubo oscillator (29),

Q k+1 P k+1 = e ε -1 hA Q k P k + h 2 P k f (P k , Q k ) -Q k f (P k , Q k ) + σ 2 P k -Q k ∆W k + h 2 P k+1 f (P k+1 , Q k+1 ) -Q k+1 f (P k+1 , Q k+1 ) + σ 2 P k+1 -Q k+1 ∆W k (39)
where we have the identity e ε -1 hA = cos(ε -1 h) -sin(ε -1 h) sin(ε -1 h) cos(ε -1 h) and ∆W k = W (t k+1 )-W (t k ) are independent centred Gaussian random vectors with variance hI.

• trig1: the explicit trigonometric method proposed in [START_REF] Cohen | On the numerical discretisation of stochastic oscillators[END_REF][START_REF] Cohen | A trigonometric method for the linear stochastic wave equation[END_REF] for second order SDEs with additive noise (3)-( 4),

Q k+1 P k+1 = cos(ε -1 hK) h sinc(ε -1 hK) -ε -1 K sin(ε -1 hK) cos(ε -1 hK) Q k P k (40) + h 2 hΨa(ΦQ k ) Ψ 0 a(ΦQ k ) + Ψ 1 a(ΦQ k+1 ) + h sinc(ε -1 hK)B∆W k cos(ε -1 hK)B∆W k
where we denote the position and velocity vectors X = Q and Ẋ = P . Here, Φ, Ψ, Ψ 0 , Ψ 1 are filter matrices introduced to avoid numerical resonances arising already in the deterministic setting [21, Chap. XIII] and defined in [START_REF] Cohen | On the numerical discretisation of stochastic oscillators[END_REF][START_REF] Cohen | A trigonometric method for the linear stochastic wave equation[END_REF] following [START_REF] Grimm | Error analysis of exponential integrators for oscillatory second-order differential equations[END_REF],

Φ = sinc(ξ), Ψ = sinc 3 (ξ), Ψ 0 = cos(ξ) sinc 2 (ξ), Ψ 1 = sinc 2 (ξ),
with ξ = ε -1 hK where we note that sinc(ξ) = ξ -1 sin(ξ) is defined for possibly singular matrices ξ. We recall that the scheme (40) does not require the eigenvalues of the highly oscillatory matrix K to be integer multiples in contrast to the multirevolution approach where (2) is assumed. In the analysis [START_REF] Cohen | A trigonometric method for the linear stochastic wave equation[END_REF] of the stiff scheme (40), it is proved a strong convergence of order one for the position,

E( Q k -Q(t k ) 2 ) 1/2 ≤
Ch for all kh ≤ T (which implies weak order one) and h ≥ cε, with constants C, c that are independent of ε, h. A bound E( P k -P (t k ) 2 ) 1/2 ≤ C is proved on the velocities. [START_REF] Sanders | Averaging methods in nonlinear dynamical systems[END_REF]. Right picture: error in E(E(T )) (average energy) for the FPU problem [START_REF] Tao | Structure preserving stochastic impulse methods for stiff Langevin systems with a uniform global error of order 1 or 1/2 on position[END_REF]. Fully-discrete S-MRCM with n = 8 micro steps per macro step.

In Figure 7, we consider the Kubo oscillator [START_REF] Sanders | Averaging methods in nonlinear dynamical systems[END_REF] and the FPU problem (34) as described in Sections 4.1 and 4.3 with highly oscillatory frequency ε = 2 -8 . We plot the weak errors (in E(Q(T ) 2 ) and E(E(T )), respectively) versus the total number of steps for multi-revolution methods s-mrcm1 and s-mrcm2, the standard Strang splitting method splt2, and trigonometric methods from [START_REF] Cohen | On the numerical discretisation of stochastic oscillators[END_REF][START_REF] Cohen | A trigonometric method for the linear stochastic wave equation[END_REF] (method trig in (39) for the Kubo oscillator and method trig1 in (40) for the FPU problem). For the multi-revolution methods, the considered number of steps corresponds to the total number of micro-steps, namely 2nT /H for s-mrcm2 and nT /H for s-mrcm1 where T /H is the number of macro steps (we choose here n = 8 micro steps). Since the most expansive part of all integrators is the integration of the nonlinearity of the problems, this is a fair measure of the computational cost of the methods. The expectancies are computed as the averages over 10 7 independent trajectories. 5For both problems, we observe that the multi-revolution composition methods s-mrcm1 and s-mrcm2 exhibit lines of slopes 1 and 2 which corroborates the uniformly accurate weak convergence rates of Theorem 2.6 and Remark 2.4. Thus, any prescribed accuracy at a fixed time T can be achieved with a computation cost which is independent of the smallness of ε. In contrast, the standard splitting methods for both problems and the trigonometric method trig in (39) for the Kubo oscillator require small stepsizes h ≤ Cε for the methods to converge. This means that these standard schemes have a computation cost which grows as O(ε -1 ) as ε → 0 which makes these nonstiff integrators inefficient for highly oscillatory problems. The stiff trigonometric method trig1 in (40) for the FPU problem has a line with slope 1 which corroborates the strong convergence analysis in [START_REF] Cohen | A trigonometric method for the linear stochastic wave equation[END_REF] and is efficient for large time steps. Notice however that the scheme converges only in the regime h ≥ cε. This is not surprising because this stepsize restriction is needed in the analysis in [START_REF] Cohen | A trigonometric method for the linear stochastic wave equation[END_REF]. We did not apply the stiff trigonometric method trig1 to the Kubo oscillator because this scheme is designed for second order SDEs with additive noise and does not generalizes straightforwardly to the multiplicative noise case.

Conclusion

We have presented and analyzed a class of large time step integrators for highly oscillatory SDEs based on the idea of multi-revolution composition methods originally introduced for deterministic problems [START_REF] Chartier | Multi-revolution composition methods for highly oscillatory differential equations[END_REF]. The versatile micro-macro approach can in principle be coupled with any micro integrator in the spirit of the Heterogeneous Multiscale Method [START_REF] Abdulle | The heterogeneous multiscale method[END_REF][START_REF] Engquist | The heterogeneous multiscale methods[END_REF][START_REF] Engquist | Heterogeneous multiscale methods: a review[END_REF]. This permits to construct methods with favourable geometric properties such as quadratic first integral preserving methods, as illustrated in the numerical experiments. Since the multi-revolution methods allow a (macro) stepsize H ≫ ε and have uniform accuracy with respect to the oscillatory period ε, the proposed approach could be coupled with the multilevel Monte-Carlo method for SDEs [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF], which permits to significantly speed up the standard Monte-Carlo method by reducing the variance of the scheme. Indeed, applied to stiff problems, the multilevel approach requires a scheme allowing both coarse and fine time steps, as shown recently in [START_REF] Abdulle | Stabilized multilevel Monte Carlo method for stiff stochastic differential equations[END_REF] in the context of stiff stochastic diffusion problems. Finally, we mention that the analysis of the multi-revolution approach is studied in [START_REF] Chartier | Convergence analysis of multi-revolution composition time-splitting pseudo-spectral methods for linear Schrödinger equations[END_REF] for the linear Schrödinger equation in the deterministic case. The extension to the stochastic case will be investigated in future works.

Proof. We denote by e t(L 1 +εL 2 ) φ(x) and e tL j φ(x), j = 1, 2 the solutions at time t of the backward Kolmogorov equation variant [START_REF] Melendo | A new approach to the construction of multirevolution methods and their implementation[END_REF] with L ε replaced by L 1 + εL 2 and L j , j = 1, 2, respectively. We have by [START_REF] Talay | Discrétisation d'une équation différentielle stochastique et calcul approché d'espérances de fonctionnelles de la solution[END_REF]Thm. 2.2] (see Theorem 3.3 with A = 0), E(φ(X(h))|X(0) = x) = e h(L 1 +εL 2 ) φ(x)

and analogously E(φ(Φ j h (x))) = e hL j φ(x), for j = 1, 2, where the generator L j is defined by L j φ(x) := f j (x) • ∇φ(x) + 1 2 m r=1 φ ′′ g r j (x), g r j (x) and e h(L 1 +εL 2 ) φ(x) and e hL j φ(x) are C 3 functions with respect to h with derivatives with polynomial growth [START_REF] De Bouard | Weak and strong order of convergence of a semidiscrete scheme for the stochastic nonlinear Schrödinger equation[END_REF]. Analogously to [START_REF] Petzold | Numerical solution of highly oscillatory ordinary differential equations[END_REF] in the proof of Lemma 3.6, we deduce

E(φ(X 1 )|X 0 = x) = e h 2 L 1 • e hεL 2 • e h 2 L 1 φ(x).
We obtain for all initial condition X 0 = x, e(h) := E(φ(X 1 )) -E(φ(X(h))) = e where C(x) is independent of h, ε and satisfies the polynomial growth [START_REF] De Bouard | Weak and strong order of convergence of a semidiscrete scheme for the stochastic nonlinear Schrödinger equation[END_REF] with respect to x. We conclude the proof of the global weak error estimate using the Milstein theorem in [START_REF] Milstein | Weak approximation of solutions of systems of stochastic differential equations[END_REF] (see Remark 2.7). In the case where Φ 2 h a weak second order approximation instead of the exact flow, the above local error estimate remains valid for all h small enough, and using Remark 2.7 permits to conclude the proof of Proposition 6.1.

Remark 3 . 5 .

 35 Consider the additive case (4) with dimensions d = m and B = I. It can be checked that for all t k = kε ≤ T with k ∈ N, and all φ ∈ C 6 P

  g. [21, Sect. III.5.1].

2 Figure 1 :

 21 Figure1: Multi-revolution methods for the Kubo oscillator[START_REF] Sanders | Averaging methods in nonlinear dynamical systems[END_REF] with nonlinearity f (p, q) = p 3 + q 5 . Error in E(Q(T ) 2 ) versus number of macro steps (final time T = 2π). Top pictures: semi-discrete S-MRCM. Bottom pictures: fully-discrete S-MRCM with n = 8, 16, 32, 64 micro steps per macro step (respectively from top to bottom).

f e r e n c e s l o p e 2 r e f e r e n c e s l o p e 1 Figure 2 :

 12 Figure2: Multi-revolution methods for the test problem[START_REF] Talay | Efficient numerical schemes for the approximation of expectations of functionals of the solution of a SDE and applications[END_REF] with non-commutative noise with ε = 2 -8 . Error in the averaged energy E(P (T ) 2 + Q(T ) 2 ) at T = 2π versus number of micro steps (left picture) or the number of diffusion function evaluations (right picture). Solid lines: method (33) with weak order 2 for the nonstiff part. Dashed lines: Euler-Maruyama method with weak order 1 for the nonstiff part. We use in both cases the multi-revolution parameter N = 128 and N = 256 (from top to bottom).

2 Figure 3 :

 23 Figure 3: Multi-revolution methods for the FPU like problem (34) with Hamiltonian (35). Error in the Hamiltonian E(E(P (T ), Q(T ))) (from top to bottom pictures, respectively) versus the number of macro steps (final time T = 2π). The lines corresponds to the fullydiscrete S-MRCM with n = 8, 16, 32 micro steps per macro step, respectively.

Figure 4 :

 4 Figure 4: Energy exchanges in the FPU-like problem (34) with ε = 1/200 on the time interval (0, 2πε -1 ). Comparison of the S-MRCM (N = 10,n = 8) (bottom pictures) and a reference solution (Strang splitting with h = 10 -3 , top pictures). Stochastic case (σ = 0.1): averages over 10 3 numerical trajectories.

2 Figure 5 :

 25 Figure 5: Stochastic nonlinear Schrödinger problem (37) with ε = 10 -4 . Error of S-MRCM at time T = 2πε -1 for E(|ξ 1 | 2 ) (first mode) and E(|ξ 3 | 2 ) (third mode) as a function of the number of macro steps. Number of micro steps per macro step: n = 100.

Figure 6 :

 6 Figure 6: Nonlinear Schrödinger problem (37) with ε = 10 -4 on the time interval (0, 2πε -1 ). Top pictures: deterministic case (σ = 0). Bottom pictures: one sample trajectory in the stochastic case (σ = 1). Left pictures: Plot of the actions |ξ j (t)| along time t, for j = 1, 3, 5, . . . , 9 (solid lines) and for j = -1, -3, -5, . . . , -9 (dotted lines) with colors from red (|j| = 1) to black (|j| = 9). Right pictures: Hamiltonian energy along time (in blue).

Figure 7 :

 7 Figure 7: Precision/work diagrams. Weak errors versus the total number of steps for multi-revolution methods or orders 1 and 2 (solid lines), standard Strang splitting methods (38) (dotted lines), and trigonometric methods (dashed lines). Highly oscillatory frequency ε = 2 -8 . Left picture: error in E(Q(T ) 2 ) (second moment of the first component) for the Kubo oscillator[START_REF] Sanders | Averaging methods in nonlinear dynamical systems[END_REF]. Right picture: error in E(E(T )) (average energy) for the FPU problem[START_REF] Tao | Structure preserving stochastic impulse methods for stiff Langevin systems with a uniform global error of order 1 or 1/2 on position[END_REF]. Fully-discrete S-MRCM with n = 8 micro steps per macro step.

h 2 L 2 L 1 φ and e εhL 2 • e h 2 L 0 e 0 e

 22200 1 • e εhL 2 • e h 2 L 1 φ(x) -e hL 1 +εhL 2 φ(x). Next, a straightforward Taylor series calculation yields d p e(h) dh p h=0 = 0 for p = 0, 1, 2 (see e.g. [21, Chap. III.4]). Using [33, Thm. 2.2] (see Theorem 3.3 with A = 0) applied repeatedly with φ replaced by e h 1 φ, the quantity e(h) is of class C 3 with respect to h and the third derivative has the polynomial growth[START_REF] De Bouard | Weak and strong order of convergence of a semidiscrete scheme for the stochastic nonlinear Schrödinger equation[END_REF] with respect to x. In addition, by the variation of constant formula,e(h) = ε h sL 1 +εsL 2 L 2 e (h-s)L 1 φds -ε h hL 1 /2 e εsL 2 L 2 e hL 1 /2 φds,which yields that the third derivative of e(h) has size O(ε) for all h. By the Taylor formula, we deduce for all initial condition X 0 = x the local error bound |e(h)| ≤ C(x)εh3

Here, the correlation in space arises from the fixed number ℓ of Fourier modes considered in the noise.

We denote ∇φ(x) the gradient with respect to x of φ and φ ′′ (x)(•, •) the second derivative of φ, which is a symmetric bilinear form.

In the results for the Kubo oscillator (left picture in Figure7), notice that the errors appear bounded from below for the finer number of steps 1024, 2048. This is due to the Monte-Carlo error, here with size ≃ 5 • 10 -4 for 10 7 trajectories.

Acknowledgements The author is grateful to Arnaud Debussche for helpful discussions on the topic of SPDEs, and to Assyr Abdulle for hospitality and parallel computing resources at the ANMC chair, École Polytechnique Fédérale de Lausanne. This work was partially supported by the Fonds National Suisse, project No. 200020 144313/1.

Appendix

The proof of Lemma 3.9 used in the proof of the main Theorem 2.6 is an immediate consequence of the following more general proposition applied with f 1 = δA, g r 1 = 0, f 2 = f , g r 2 = g r , and ε = γH. Although this result on the Strang splitting method is well known in deterministic contexts (see for instance [START_REF] Hairer | Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF]), it seems not available in the literature in a stochastic context, and it can be useful by itself also for ε = 1. Proposition 6.1. For a fixed ε, consider the system of SDEs

where

. . , m are smooth vector fields satisfying (H). Consider the stochastic flows Φ 1 h and Φ 2 h with time h of the auxiliary SDE problems dX = f 1 dt + m r=1 g r 1 dW r and dX = f 2 dt + m r=1 g r 2 dW r , respectively. Then, the integrator

where Φ 1 h/2 , Φ 2 εh , and Φ 1 h/2 are applied with independent Wiener processes, satisfies for all φ ∈ C 6 P (R d , R) and all hk ≤ T ,

where C is independent of h, ε and X(t) is the solution of (41). In addition, considering for Φ 2 h a weak second order approximation satisfying (15) and (16), the above error estimates remains valid.