# Multirevolution Integrators for Differential Equations with Fast Stochastic Oscillations

Adrien Laurent, Gilles Vilmart

► To cite this version:

Adrien Laurent, Gilles Vilmart. Multirevolution Integrators for Differential Equations with Fast Stochastic Oscillations. SIAM Journal on Scientific Computing, 2020, 42 (1), pp.A115-A139. 10.1137/19M1243075. hal-04347677

HAL Id: hal-04347677

https://hal.science/hal-04347677v1

Submitted on 15 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Multirevolution integrators for differential equations with fast stochastic oscillations

Adrien Laurent<fnref n="1" /> and Gilles Vilmart<fnref n="1" />

December 16, 2019

## Abstract

We introduce a new methodology based on the multirevolution idea for constructing integrators for stochastic differential equations in the situation where the fast oscillations themselves are driven by a Stratonovich noise. Applications include in particular highly-oscillatory Kubo oscillators and spatial discretizations of the nonlinear Schrödinger equation with fast white noise dispersion. We construct a method of weak order two with computational cost and accuracy both independent of the stiffness of the oscillations. A geometric modification that conserves exactly quadratic invariants is also presented.

**Keywords:** highly-oscillatory stochastic differential equations, nonlinear Schrödinger equation, white noise dispersion, geometric integration, quadratic first integral.

**AMS subject classification (2010):** 60H35, 35Q55, 34E13.

## 1 Introduction

This article aims at developing invariant-preserving integrators of second weak order that are robust with respect to the stiffness $\varepsilon$ both in accuracy and cost for the following class of highly-oscillatory $d$-dimensional SDEs driven by a one-dimensional Stratonovich noise

$$dX(t) = \frac{1}{\sqrt{\varepsilon}}AX(t) \circ dW(t) + F(X(t))dt, \quad t > 0, \quad X(0) = X_0, \quad (1.1)$$

where $W$ is a standard one-dimensional Wiener process, the function $F: \mathbb{R}^d \to \mathbb{R}^d$ is a smooth non-linear map, the stiff parameter $\varepsilon > 0$ is fixed and assumed small, and $A \in \mathbb{R}^{d \times d}$ is a given matrix satisfying $e^A = \text{Id}$ (equivalently $A$ is diagonalizable and has all its eigenvalues in $2i\pi\mathbb{Z}$). In the deterministic setting, this last property yields that the solution $x(t) = \exp(\varepsilon At)x_0$ of $\frac{dx}{dt} = \varepsilon^{-1}Ax$ is $\varepsilon$-periodic. For stochastic oscillations, it means that the solution $X(t) = \exp(\varepsilon^{-1/2}AW(t))X_0$ of $dX = \varepsilon^{-1/2}AX \circ dW$ satisfies $X(T) = X(0)$ for a random time $T = \inf\{t > 0, |\varepsilon^{-1/2}W(t)| = 1\}$ of mean $\varepsilon$. The class of SDEs (1.1) includes in particular highly-oscillatory Kubo oscillators (see [8])

$$dX = \frac{2\pi}{\sqrt{\varepsilon}} \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} X \circ dW + \begin{pmatrix} 0 & -a \\ a & 0 \end{pmatrix} X dt, \quad a \in \mathbb{R}, \quad (1.2)$$

<footnote n="1">¹Université de Genève, Section de mathématiques, 2-4 rue du Lièvre, CP 64, CH-1211 Genève 4, Switzerland. Adrien.Laurent@unige.ch, Gilles.Vilmart@unige.ch.</footnote>
or equivalently, $dY = 2i\pi\epsilon^{-1/2}Y \circ dW + iaYdt$ in the complex setting where $Y = X_1 + iX_2$.

Applying standard SDE integrators to solve equation (1.1) requires in general a time stepsize $h \le \epsilon$ to be accurate, which makes these methods dramatically expensive when $\epsilon$ is small. The goal of this paper is to create robust numerical methods, i.e., numerical integrators whose cost and accuracy do not deteriorate when $\epsilon$ becomes small. Several classes of methods have already been developed for highly-oscillatory SDEs with a deterministic fast oscillation (see for instance [10, 25]), but not in the case where the stiff oscillatory part is applied to the noise itself. To numerically face this challenge, we introduce in this paper a new methodology to develop robust methods of any high weak order to approximate the solution of equation (1.1). In particular, we propose a method of weak order two, and a geometric modification of this algorithm that preserves quadratic invariants.

Stochastic oscillations as defined in (1.1) typically arise in fiber optics models (see [2, 3, 15]) with a spatial discretization of the highly-oscillatory nonlinear Schrödinger equation (NLS) with white noise dispersion

$$du(t) = \frac{i}{\sqrt{\epsilon}} \Delta u(t) \circ dW(t) + F(u(t))dt, \quad u(t=0) = u_0. \qquad (1.3)$$

As described for instance in [15], in the case $\epsilon = 1$, the NLS equation (1.3) with a cubic nonlinearity $F(u) = |u|^2u$ is a model in dimension $d = 1$ describing the propagation of a signal in optical fibers where $x$ corresponds to the retarded time, while $t$ corresponds to the distance along the fiber. Taking into account the inevitable chromatic dispersion effects of the signal, modeled by a random centered stationary process $m$ with a coefficient $\nu > 0$, yields the following random PDE,

$$\frac{\partial v}{\partial x}(x,t) = \nu im(x) \frac{\partial^2 v}{\partial t^2}(x,t) + \nu^2 F(v(x,t)), \quad v(x=0,t) = u_0(t).$$

The perfect fiber would satisfy $m = 0$, but in practice, engineers build fibers with a small varying dispersion coefficient. To limit the pulse broadening induced by random dispersion, specialists use a wide range of dispersion management techniques (see for instance [15] and references therein). In [20, 11], the authors show that if we denote $u^{\nu}(x,t) = v(x/\nu^2, t)$, then as $\nu$ tends to 0 and under some ergodicity assumptions on $m$, $u^{\nu}$ converges to the solution $u$ of equation (1.3) with $\epsilon = 1$. The non-stiff counterpart of equation (1.3), i.e., for $\epsilon = 1$, has also been studied theoretically in [12] for a particular nonlinearity. The highly-oscillatory behavior ($\epsilon \ll 1$) appears naturally when observing the propagation in long time with a small nonlinearity (via the change of variable $t \leftarrow \epsilon t$) or the propagation of a small initial data in an optical fiber with a polynomial nonlinearity (via the change of variable $u \leftarrow u/\epsilon$). A goal of this article is to develop efficient and cheap numerical methods that can model the propagation of pulses in this context, in order to observe some specific behaviors and, ultimately, to build enhanced fibers. Models of the form (1.3) also appear in the recent work [14] in the context of stochastic three-wave semi-linear systems. We emphasize that there is a growing interest in the recent literature for stochastic models involving a fast Stratonovitch noise in the context of ergodic stochastic dynamics. In [1], it is shown for a class of overdamped Langevin equations that adding an appropriate fast Stratonovitch noise permits to increase the convergence rate to equilibrium, while reducing the asymptotic variance at infinity. This suggests that new efficient samplers for the invariant distribution of Langevin type models in context of large dimensional molecular dynamics models could be developed. We also mention the recent homogenization results on stochastic dynamics
<figure><img src="image_2.png" /><figcaption>Figure 1: Exact solution evaluated at revolution times for the deterministic oscillator (2.1) with $F(y) = iy$ and $\varepsilon = 10^{-1}$.</figcaption></figure>

with fast Stratonovitch noises in [18] where our periodicity assumption is replaced by an ergodicity assumption on the fast component of the dynamics posed on manifolds.

Numerous possibilities exist for numerically integrating equations (1.1) or (1.3). We highlight in particular the exponential integrators [8, 13] for the SDE (1.1), and the exponential integrators [9], the Fourier split-step method [20] or the Crank-Nicholson scheme [4] for the SPDE (1.3). These methods have the advantage that they preserve the $L^2$ invariant of the equation (that is, $\|u(t)\|_{L^2} = \|u_0\|_{L^2}$ for all $t \ge 0$) for a class of polynomial nonlinearities. However, they face a severe timestep restriction $h \le \varepsilon$ when the stiff parameter $\varepsilon$ is small. Even in the case of deterministic oscillations, there are restrictions in general, though some robust algorithms exist (see [10] for instance). The methods presented in this paper solve this issue of stepsize restriction. The idea is to approximate the solution of equation (1.1) at random times called revolution times because they correspond to complete revolutions of the oscillatory part $dX = \varepsilon^{-1/2}AX \circ dW$. This is in the spirit of [17] which also approximates the solution of SDEs at random times.

The article is organized as follows. Section 2 is devoted to the presentation of the new integrators. In section 3, we build an asymptotic expansion of the solution of (1.1) and evaluate it at revolution times to derive the new integrators and a limit model for equation (1.1). Section 4 is devoted to the weak convergence theorems and their proofs. In section 5, we present numerical experiments to confirm our theoretical error estimates, and we apply the new methods to solve numerically the Schrödinger equation (1.3).

## 2 Multirevolution integrators for stochastic oscillators

Initially created in [21, 5] in the context of celestial mechanics and later extended using geometric integration (see for instance [23, 6, 7]), multirevolution methods represent a class of numerical methods used for solving highly-oscillatory differential equations while reducing the cost of computation. In particular, they can approximate the solution of highly-oscillatory ODEs of the following form at stroboscopic times $\varepsilon NT$, where $T = 1$ is the period of $\frac{dx}{dt} = Ax$, and $N$ is an integer,

$$ \frac{dx}{dt} = \frac{1}{\varepsilon}Ax + F(x), \quad x(0) = x_0. \qquad (2.1) $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Revolution times (2.2) of a Brownian path (top) and exact solution evaluated at revolution times for the Kubo oscillator (1.2) with $a = 1$ and $\varepsilon = 10^{-1}$ (bottom).</figcaption></figure>

The solution $x$ of this equation at times $\varepsilon NT$ is a perturbation of identity, that is, $x$ satisfies $x(\varepsilon t) = x_0 + \mathcal{O}(\varepsilon t)$, thus the solution loses its highly-oscillatory feature when evaluated at stroboscopic times, as shown in Figure 1 for the first component of the solution of equation (2.1) with $F(x) = ix$ (respectively $F(y) = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} y$ in the real setting). The idea of multirevolution is to approximate $x(\varepsilon N)$ with $N = \mathcal{O}(\varepsilon^{-1})$ with a cost independent of $\varepsilon$.

For stochastic oscillations, the solution $X(t) = e^{\varepsilon^{-1/2}AW(t)}X_0$ of $dX = \varepsilon^{-1/2}AX \circ dW$ is not periodic, but satisfies $X(\varepsilon T_N) = X_0$ where the $T_N$ are random variables called revolution times and defined by

$$ T_0 = 0, \tag{2.2} \\ T_{N+1} = \inf \left\{ t > T_N, \varepsilon^{-1/2} |W(\varepsilon t) - W(\varepsilon T_N)| \ge 1 \right\}, \quad N = 0, 1, 2, \dots $$

If $X$ is the solution of (1.1), we show in section 3.1 that $X$ evaluated at times $\varepsilon T_N$ is a perturbation of identity (in a strong and weak sense). Figure 2 illustrates the definition of revolution times and shows the perturbation of identity property on the first component of a Kubo oscillator (1.2) with $a = 1$. We highlight that the revolutions times $T_N$ can be simulated without simulating the exact path $W$. Also we emphasize that the proposed algorithms do not require to simulate $W$ thanks to the use of appropriate discrete random variables. This will be detailed in section 3.4.
We show in section 3.3 that the solution $X$ of (1.1) evaluated at times $\varepsilon T_{t\varepsilon^{-1}}$ (when $t/\varepsilon \in \mathbb{N}$ is an integer) converges weakly when $\varepsilon \to 0$ to the solution $y_t$ of the deterministic ODE

$$ \frac{dy_t}{dt} = \langle g^0 \rangle(y_t), \quad y_0 = X_0, $$

where $g_{\theta}^{0}(y) = e^{-A\theta}F(e^{A\theta}y)$ and $\langle g^0 \rangle := \int_{0}^{1} g_{\theta}^{0}d\theta$. This ODE is exactly the same one as the asymptotic model for deterministic oscillators of the form (2.1). This asymptotic model naturally yields a weak order 1 deterministic integrator. We propose the two following new multirevolution methods of second weak order for integrating equation (1.1) at the revolution times $\varepsilon T_{Nm}$ for $m = 0, 1, 2, \dots$ with cost in $H = N\varepsilon = O(1)$ independent of $\varepsilon$. Method B is a geometric modification of Method A to preserve quadratic invariants of the form $Q(y) = \frac{1}{2}y^T S y$ where $S \in \mathbb{R}^{d \times d}$ is a given symmetric matrix. Methods A and B involve a Fourier decomposition of the following functions that are 1-periodic with respect to $\theta$,

$$ g_{\theta}^{0}(y) = e^{-A\theta}F(e^{A\theta}y) = \sum_{k \in \mathbb{Z}} c_{k}^{0}(y)e^{2i\pi k\theta} \quad (2.3) $$

$$ g_{\theta}^{1}(y)(z) = e^{-A\theta}F'(e^{A\theta}y)(e^{A\theta}z) = \sum_{p \in \mathbb{Z}} c_{p}^{1}(y)(z)e^{2i\pi p\theta} $$

with Fourier coefficients $(c_k^0(y))_{k \in \mathbb{Z}}$ and $(c_p^1(y))_{p \in \mathbb{Z}}$, respectively. The series appearing in (2.3) have an infinite number of terms in general. For a practical implementation of the new methods, we truncate these series up to an even number of modes $K_t$, while inducing an exponentially small error (see Remark 4.3). For each timestep, we also introduce the bounded discrete random variables $(\hat{\alpha}_k^N)_{k}$, and deterministic sequences $(\hat{\beta}_{p,k}^N)_{p,k}$ and $(\hat{\beta}_{p,k}^N)_{p,k}$ that satisfy

$$
\begin{aligned}
\mathbb{E}[\hat{\alpha}_k^N] &= \begin{cases} 1 & \text{if } k = 0 \\ 0 & \text{else} \end{cases} &
\mathbb{E}[\hat{\alpha}_p^N \hat{\alpha}_k^N] &= \begin{cases} 1 + \frac{2}{3N} & \text{if } p = k = 0 \\ \frac{1}{\pi^2 p^2 N} & \text{if } p + k = 0, p, k \neq 0 \\ 0 & \text{else} \end{cases} \\
\hat{\beta}_{p,k}^N &= \begin{cases} \frac{1}{2} + \frac{1}{3N} & \text{if } p = k = 0 \\ \frac{1}{2\pi^2 k^2 N} & \text{if } p = 0, k \neq 0 \\ \frac{-1}{2\pi^2 p^2 N} & \text{if } p \neq 0, k = 0 \\ \frac{1}{2\pi^2 p^2 N} & \text{if } p + k = 0, p, k \neq 0 \\ 0 & \text{else} \end{cases} &
\hat{\beta}_{p,k}^N &= \begin{cases} \frac{1}{2\pi^2 k^2 N} & \text{if } p = 0, k \neq 0 \\ \frac{-1}{2\pi^2 p^2 N} & \text{if } p \neq 0, k = 0 \\ 0 & \text{else} \end{cases}
\end{aligned}
$$

The definition and construction of these random variables are further discussed in sections 3.2 and 3.4.

**Method A** (Explicit integrator of weak order two in $H = N\varepsilon$ to approximate the solution of equation (1.1) at times $\varepsilon T_{Nm}$ for $m = 0, 1, 2, \dots$)

$$ Y_0 = X_0 \\ \text{for } m \ge 0 \text{ do} $$

$$ Y_{m+1} = Y_m + H \sum_{k=-K_t/2}^{K_t/2-1} c_k^0(Y_m) \hat{\alpha}_k^N + H^2 \sum_{p,k=-K_t/2}^{K_t/2-1} c_p^1(Y_m)(c_k^0(Y_m)) \hat{\beta}_{p,k}^N \quad (2.4) $$

end for
**Method B (Geometric integrator of weak order two in $H = N\epsilon$ to approximate the solution of equation (1.1) at times $\epsilon T_{Nm}$ for $m = 0, 1, 2, \dots$ while preserving quadratic invariants)**

$Y_0 = X_0$ 
for $m \ge 0$ do

$$ 
\begin{aligned}
Y_{m+1} = Y_m + H \sum_{k=-K_t/2}^{K_t/2-1} c_k^0 \left( \frac{Y_m + Y_{m+1}}{2} \right) \hat{\alpha}_k^N & \quad (2.5) \\
& + H^2 \sum_{p,k=-K_t/2}^{K_t/2-1} c_p^1 \left( \frac{Y_m + Y_{m+1}}{2} \right) \left( c_k^0 \left( \frac{Y_m + Y_{m+1}}{2} \right) \right) \hat{\beta}_{p,k}^N
\end{aligned}
$$

end for

**Remark 2.1.** One could apply a Newton iteration to solve the implicit equation (2.5) in Method B. However, a few fixed point iterations are sufficient (see discussion in [16, Chap. VIII] for non-stiff implicit methods). Indeed, since the Lipschitz constant of the iterated map has size $O(H)$, the convergence rate of the fixed point iterations is independent of the smallness of $\epsilon$.

**Remark 2.2.** We observe that $\hat{\beta}_{p,k}^N$ and $\hat{\beta}_{p,k}^N$ are always zero except when $p = 0, k = 0$ or $p + k = 0$. Thus the computational cost of one step of Methods A and B grows linearly in the number of modes in (2.3).

# 3 Analysis and asymptotic expansion of the exact solution

In this section, we first obtain a local expansion of the solution of (1.1) and then evaluate it at particular random times to deal with the highly-oscillatory patterns of the exact solution. Finally we derive from this expansion an asymptotic limit for equation (1.1) when $\epsilon \to 0$.

## 3.1 Asymptotic expansion of the exact solution

Instead of studying directly equation (1.1), we apply the change of variable $t \leftarrow \epsilon^{-1}t$ to obtain the following equation, whose solution satisfies $Y(t) = X(\epsilon t)$ with $X$ solution of (1.1),

$$ dY(t) = AY(t) \circ d\widetilde{W}(t) + \epsilon F(Y(t))dt, \quad Y(0) = X_0, \quad (3.1) $$

where we denote for simplicity the Brownian motion $\widetilde{W}(t) = \epsilon^{-1/2}W(\epsilon t)$ again by $W$. We introduce the following assumption which guarantees in particular global existence and uniqueness of the solution.

**Assumption 3.1.** The function $F$ is globally Lipschitz continuous and lies in $C_p^3$, i.e., there exists constants $L, C, K > 0$ such that for all $y, y_1, y_2 \in \mathbb{R}^d$

$$ |F(y_1) - F(y_2)| \le L |y_1 - y_2| \quad |F^{(i)}(y)| \le C(1 + |y|^K), \quad i \in \{0, 1, 2, 3\}. \quad (3.2) $$

Also the initial condition $X_0$ has bounded moments, that is, $E[|X_0|^p] < \infty$ for $p \ge 0$.
Therefore we denote $\varphi_{\varepsilon,t}(X_0) = Y(t)$ the solution of equation (3.1) and focus in the rest of the paper on the approximation of $\varphi_{\varepsilon,t}(y)$ at times $t = O(\varepsilon^{-1})$. The variation of constants formula yields

$$ \varphi_{\varepsilon,t}(y) = e^{AW(t)}y + \varepsilon \int_{0}^{t} e^{A(W(t)-W(s))} F(\varphi_{\varepsilon,s}(y))ds. \quad (3.3) $$

We deduce the following regularity properties.

**Lemma 3.2.** Under Assumption 3.1, the following estimates hold for all $y, y_1, y_2 \in \mathbb{R}^d$, where C and K are independent of $\varepsilon$ and $t$,

$$ 1. \quad |\varphi_{\varepsilon,t}(y_1) - \varphi_{\varepsilon,t}(y_2)| \le C |y_1 - y_2| e^{C\varepsilon t}, $$

$$ 2. \quad |\varphi_{\varepsilon,t}(y)| \le C(1 + |y|)e^{C\varepsilon t}, $$

$$ 3. \quad \varphi_{\varepsilon,t}(y) \text{ is } \mathcal{C}^3 \text{ in } y \text{ and } \left|\varphi_{\varepsilon,t}^{(i)}(y)\right| \le C(\varepsilon t)^{i-1}(1 + |y|^K)e^{C\varepsilon t} \text{ for } i \in 1, 2, 3. $$

The proof is postponed to the appendix. It mainly relies on the Gronwall theorem and the boundedness of the one-periodic function $\theta \mapsto e^{\theta A}$. Using a local expansion of the solution of (3.1) in $\varepsilon$, we define the following first and second order approximations of $\varphi_{\varepsilon,t}(y)$,

$$ \psi_{\varepsilon,t}^{1}(y) = e^{AW(t)}y + \varepsilon e^{AW(t)} \int_{0}^{t} e^{-AW(s)} F(e^{AW(s)}y)ds \quad (3.4) $$

$$ \psi_{\varepsilon,t}^{2}(y) = \psi_{\varepsilon,t}^{1}(y) + \varepsilon^{2}e^{AW(t)} \int_{0}^{t} e^{-AW(s)} \cdot F'(e^{AW(s)}y) \left( e^{AW(s)} \int_{0}^{s} e^{-AW(r)} F(e^{AW(r)}y)dr \right) ds. $$

**Proposition 3.3** (local expansion). Under Assumption 3.1, for all $y \in \mathbb{R}^d$, $j \in \{1, 2\}$ and $t \ge 0$, there exist C and K two positive constants independent of $\varepsilon$ and $t$ such that

$$ \left|\varphi_{\varepsilon,t}(y) - \psi_{\varepsilon,t}^{j}(y)\right| \le C(1 + |y|^K)e^{C\varepsilon t}(\varepsilon t)^{j+1}. $$

The functions $\psi_{\varepsilon,t}^j$ satisfy the following straightforward inequalities proved with similar arguments as for Lemma 3.2.

**Lemma 3.4.** With the assumptions and notations of Proposition 3.3, the following estimates hold for all $y \in \mathbb{R}^d$, where C and K are independent of $\varepsilon$ and $t$,

$$ |\psi_{\varepsilon,t}^{1}(y)| \le C(1 + |y|)e^{C\varepsilon t}, \quad (3.5) $$

$$ |\psi_{\varepsilon,t}^{2}(y)| \le C(1 + |y|^{K})e^{C\varepsilon t}, \quad (3.6) $$

$$ \left|\psi_{\varepsilon,t}^{2}(y) - e^{AW(t)}y\right| \le C(1 + |y|^{K})(\varepsilon t)e^{C\varepsilon t}. \quad (3.7) $$

*Proof of Proposition 3.3.* Using Assumption 3.1, we get

$$ \left|\varphi_{\varepsilon,t}(y) - \psi_{\varepsilon,t}^{1}(y)\right| \le L\varepsilon \int_{0}^{t} \left|\varphi_{\varepsilon,s}(y) - e^{AW(s)}y\right| ds. $$
Then Lemma 3.2 yields

$$
\begin{aligned}
|\varphi_{\varepsilon,s}(y) - e^{AW(s)}y| &\le C\varepsilon \int_0^s |F(\varphi_{\varepsilon,r}(y))| dr \le C\varepsilon \int_0^s (1 + |\varphi_{\varepsilon,r}(y)|) dr \\
&\le C\varepsilon \int_0^s (1 + C(1 + |y|)e^{C\varepsilon r}) dr \le C(1 + |y|)e^{C\varepsilon s}(\varepsilon s).
\end{aligned}
$$

We deduce $|\varphi_{\varepsilon,t}(y) - \psi_{\varepsilon,t}^1(y)| \le C(1 + |y|)e^{C\varepsilon t}(\varepsilon t)^2$.

For $j = 2$, we first denote

$$
\tilde{\psi}_{\varepsilon,t}^2(y) = e^{AW(t)}y + \varepsilon e^{AW(t)} \int_0^t e^{-AW(s)} F(\psi_{\varepsilon,s}^1(y)) ds.
$$

With the same arguments we used for $j = 1$ and inequality (3.5), we have

$$
|\varphi_{\varepsilon,t}(y) - \tilde{\psi}_{\varepsilon,t}^2(y)| \le C(1 + |y|^K) e^{C\varepsilon t} (\varepsilon t)^3.
$$

It is sufficient to prove that $|\psi_{\varepsilon,t}^2(y) - \tilde{\psi}_{\varepsilon,t}^2(y)| \le C(1 + |y|^K)e^{C\varepsilon t}(\varepsilon t)^3$. A Taylor expansion of $F(\psi_{\varepsilon,s}^1(y))$ in $\varepsilon$ gives

$$
F(\psi_{\varepsilon,s}^1(y)) = F(e^{AW(s)}y) + \varepsilon F'(e^{AW(s)}y) \left( e^{AW(s)} \int_0^s e^{-AW(r)} F(e^{AW(r)}y) dr \right) + R_{\varepsilon,s}.
$$

The remainder $R_{\varepsilon,s}$ satisfies

$$
|R_{\varepsilon,s}| \le C\varepsilon^2 \sup_{x\in[e^{AW(s)}y, \psi_{\varepsilon,s}^1(y)]} \|F''(x)\| e^{AW(s)} \int_0^s e^{-AW(r)} F(e^{AW(r)}y) dr^2.
$$

Then, using the polynomial growth of $F''$ and inequality (3.5), we get

$$
|R_{\varepsilon,s}| \le C(1 + |e^{AW(s)}y|^K + |\psi_{\varepsilon,s}^1(y)|^K)(\varepsilon s)^2 e^{C\varepsilon s} \le C(1 + |y|^K)(\varepsilon s)^2 e^{C\varepsilon s}.
$$

Hence the result.

I did not find any text in the image.

We shall prove in section 3.3 that the function $\psi_{\varepsilon,t}^2$ in (3.4) evaluated at the revolution times $T_N$ (defined in (2.2)) yields a strong order 2 approximation in $H = \varepsilon N$.

**Remark 3.5.** If we replace the Brownian motion $W$ in (3.4) by a piecewise linear function $W_\tau$ defined by

$$
W_{\tau} = \left(1 - \frac{t}{\tau} + i\right) W_i + \left(\frac{t}{\tau} - i\right) W_{i+1} \quad \text{for } i\tau \le t \le (i+1)\tau, \quad (3.8)
$$

where $W_0 = 0$ and $W_{i+1} = Wi + \sqrt{\tau}\xi_i$ with $(\xi_i)_i$ a family of independent standard Gaussian random variables, then it can be shown that we obtain an integrator of strong order two in $\epsilon t$. However, the cost of a standard method computing an approximation of the integrals of equation (3.4) by replacing $W$ with $W_\tau$ is in $O(t^2/\tau^2)$, which makes this method tremen- dously expensive for $t = O(\epsilon^{-1})$. This is why we develop in section 3.4 weak integrators based on a weak approximation of equation (3.4) with a cost independent of $t$. We shall re- place stochastic integrals by appropriate discrete random variables in order not to simulate any expensive Brownian path $W$.
<figure><img src="image_4.png" /><figcaption>Figure 3: Convergence in law of $\frac{\sqrt{N}}{\mathrm{Var}(T_1)}(\frac{T_N}{N} - \mathbb{E}[T_1])$ to a standard Gaussian random variable.</figcaption></figure>

## 3.2 Properties of the revolution times

In this section, we study some properties linked to the revolution times $T_N$ that will be useful for the analysis.

**Proposition 3.6.** The revolution times $T_N$ defined in (2.2) are positive and finite almost surely. Their differences $(T_{N+1} - T_N)_{N \ge 0}$ are independent identically distributed random variables (same law as $T_1$). The Laplace transform $\mathbb{E}[e^{zT_1}]$ of $T_1$ exists and is analytic for $\mathrm{Re}(z) < \frac{\pi^2}{8}$. In addition, for $x \in [0, \frac{\pi^2}{8}]$, $\mathbb{E}[e^{xT_1}] = \frac{1}{\cos(\sqrt{2x})}$. The variable $T_1$ has bounded moments, and they are given by

$$ \mathbb{E}[T_1^k] = \frac{(-2)^k k!}{(2k)!} \sum_{j=1}^{p} (-1)^j \sum_{\substack{n_1+\cdots+n_j=p \\ n_i \in \mathbb{N}^*}} \binom{2p}{2n_1, \ldots, 2n_j}. \quad (3.9) $$

In particular, $\mathbb{E}[T_1] = 1$, $\mathbb{E}[T_1^2] = \frac{5}{3}$, and $\mathrm{Var}(T_1) = \frac{2}{3}$. Finally, for a fixed $\varepsilon_0 \in ]0, \frac{\pi^2}{16}]$, for all $\varepsilon \in ]0, \varepsilon_0]$ and $p \ge 0$, we have the estimate

$$ \mathbb{E}[e^{\varepsilon T_N} (\varepsilon T_N)^p] \le C e^{C \varepsilon N} (\varepsilon N)^p. \quad (3.10) $$

The law of the first revolution time $T_1$ has an analytic density, but there is no closed formula for it. It can be numerically approximated accurately by inverting the Laplace transform. In Figure 3, we observe the convergence in law of $T_N$ to a Gaussian variable according to the central limit theorem.

**Proof of Proposition 3.6.** The first properties can be deduced from [24, Chap. 2.3], where the Laplace transform formula is obtained with an analytic continuation of the equality $\mathbb{E}[e^{-xT_1}] = \frac{1}{\cosh(\sqrt{2x})}$ for $x > 0$. Comparing the Taylor expansions of $\mathbb{E}[e^{xT_1}]$ and $\frac{1}{\cos(\sqrt{2x})}$ yields (3.9). The estimate (3.10) is proved as follows

$$ \mathbb{E}[e^{\varepsilon T_N} (\varepsilon T_N)^p] \le \mathbb{E}[e^{2\varepsilon T_N}]^{1/2} \mathbb{E}[(\varepsilon T_N)^{2p}]^{1/2} = \mathbb{E}[e^{2\varepsilon T_1}]^{N/2} \varepsilon^p \mathbb{E}[T_N^{2p}]^{1/2} $$
$$ \leq \mathbb{E}[e^{2\epsilon_0 T_1}]^{\epsilon N/2\epsilon_0} (\epsilon N)^p \mathbb{E}[T_1^{2p}]^{1/2} \leq C e^{C\epsilon N} (\epsilon N)^p, $$

where we used first the Cauchy-Schwarz inequality and then twice the Jensen inequality. □

For developing an algorithm for the weak error, it is useful to know the moments of the random variables appearing in the discretization, that are costly to simulate numerically, in order to replace them with cheap discrete approximations with the same first and second moments. This is the goal of the following proposition.

**Proposition 3.7.** The random variables

$$ \begin{aligned} \alpha_k^N &= \frac{1}{N} \int_0^{T_N} e^{2i\pi k W(s)} ds, \\ \beta_{p,k}^N &= \frac{1}{N^2} \int_0^{T_N} e^{2i\pi p W(s)} \int_0^s e^{2i\pi k W(r)} dr ds, \\ \tilde{\beta}_{p,k}^N &= \frac{1}{N^2} \int_0^{T_N} e^{2i\pi p W(s)} \left( \int_0^s e^{2i\pi k W(r)} - \frac{1}{2} \int_0^{T_N} e^{2i\pi k W(r)} dr \right) ds = \beta_{p,k}^N - \frac{\alpha_p^N \alpha_k^N}{2} \end{aligned} $$

satisfy

$$ \begin{aligned} \mathbb{E}[\alpha_k^N] &= \delta_k = \begin{cases} 1 & \text{if } k = 0 \\ 0 & \text{else} \end{cases}, & \mathbb{E}[\alpha_p^N \alpha_k^N] &= \begin{cases} 1 + \frac{2}{3N} & \text{if } p = k = 0 \\ \frac{1}{\pi^2 p^2 N} & \text{if } p + k = 0, p, k \neq 0 \\ 0 & \text{else} \end{cases} \\ \mathbb{E}[\beta_{p,k}^N] &= \begin{cases} \frac{1}{2} + \frac{1}{3N} & \text{if } p = k = 0 \\ \frac{1}{2\pi^2 k^2 N} & \text{if } p = 0, k \neq 0 \\ \frac{1}{2\pi^2 p^2 N} & \text{if } p \neq 0, k = 0 \\ \frac{1}{2\pi^2 p^2 N} & \text{if } p + k = 0, p, k \neq 0 \\ 0 & \text{else} \end{cases}, & \mathbb{E}[\tilde{\beta}_{p,k}^N] &= \begin{cases} \frac{1}{2\pi^2 k^2 N} & \text{if } p = 0, k \neq 0 \\ \frac{-1}{2\pi^2 p^2 N} & \text{if } p \neq 0, k = 0 \\ 0 & \text{else} \end{cases} \end{aligned} $$

*Proof.* Let $k \neq 0$ (the case $k = 0$ is straightforward using Proposition 3.6), then the Itô formula applied to $e^{2i\pi k W(s)}$ gives

$$ \frac{1}{2\pi^2 k^2 N} (e^{2i\pi k W(t)} - 1) = \frac{i}{\pi k N} \int_0^t e^{2i\pi k W(s)} dW(s) - \frac{1}{N} \int_0^t e^{2i\pi k W(s)} ds, $$

which yields at time $t = T_N$,

$$ \alpha_k^N = \frac{i}{\pi k N} \int_0^{T_N} e^{2i\pi k W(s)} dW(s). $$

Then $t \mapsto \int_0^t e^{2i\pi k W(s)} dW(s)$ is a martingale, so by the Doob theorem, as $t \wedge T_N$ is finite, $\mathbb{E}[\int_0^{t \wedge T_N} e^{2i\pi k W(s)} dW(s)] = 0$ for all $t$. The dominated convergence theorem for stochastic integrals allows us to take the limit $t \to \infty$ and yields $\mathbb{E}[\alpha_k^N] = 0$.

For the coefficients $\beta_{p,k}^N$, let $(p, k) \neq (0, 0)$ (the case $p = k = 0$ is obtained straightforwardly using Proposition 3.6), we use the Itô formula on $e^{2i\pi p W(s)}$ and we integrate from $s$ to $T_N$,

$$ 1 - e^{2i\pi p W(s)} = 2i\pi p \int_s^{T_N} e^{2i\pi p W(r)} dW(r) - 2\pi^2 p^2 \int_s^{T_N} e^{2i\pi p W(r)} dr. $$

Then, multiplying by $\frac{1}{N^2} e^{2i\pi k W(s)}$ and integrating from 0 to $T_N$ yields

$$ \frac{\alpha_k^N - \alpha_{p+k}^N}{N} = \frac{2i\pi p}{N^2} \int_0^{T_N} e^{2i\pi k W(s)} \int_s^{T_N} e^{2i\pi p W(r)} dW(r) ds $$
$$ -\frac{2\pi^2 p^2}{N^2} \int_0^{T_N} e^{2i\pi k W(s)} \int_s^{T_N} e^{2i\pi p W(r)} dr ds. $$

Using the stochastic Fubini theorem, we deduce

$$ \int_0^{T_N} e^{2i\pi k W(s)} \int_s^{T_N} e^{2i\pi p W(r)} dW(r) ds = \int_0^{T_N} e^{2i\pi p W(r)} \int_0^r e^{2i\pi k W(s)} ds dW(r), $$

which has zero average by the same arguments as before. Also by the Fubini theorem for stochastic integrals, $\beta_{p,k}^N = \frac{1}{N^2} \int_0^{T_N} e^{2i\pi k W(s)} \int_s^{T_N} e^{2i\pi p W(r)} dr ds$, so that we get if $p \neq 0$,

$$ \mathbb{E}[\beta_{p,k}^N] = \frac{\delta_{p+k} - \delta_k}{2\pi^2 p^2 N}. $$

The case $p = 0$ is obtained by integrating by parts and using the same arguments. Indeed, we find

$$ \beta_{0,k}^N = \frac{T_N}{N} \alpha_k^N - \beta_{k,0}^N = \frac{i}{\pi k N^2} \int_0^{T_N} T_N e^{2i\pi k W(s)} dW(s) - \beta_{k,0}^N, $$

and $\mathbb{E}[\beta_{0,k}^N] = -\mathbb{E}[\beta_{k,0}^N]$. Finally, the moments $\mathbb{E}[\alpha_p^N \alpha_k^N]$ are computed via the equality $\beta_{p,k}^N + \beta_{k,p}^N = \alpha_p^N \alpha_k^N$. Then we obtain $\mathbb{E}[\tilde{\beta}_{p,k}^N]$ from the formula $\tilde{\beta}_{p,k}^N = \beta_{p,k}^N - \frac{\alpha_p^N \alpha_k^N}{2}$. $\square$

**Remark 3.8** (stochastic Fourier series). Let $f$ be a $L^2$ function on $]0, 1[$ extended on $\mathbb{R}$ by 1-periodicity, whose Fourier coefficients are denoted as $(c_k)_{k \in \mathbb{Z}}$. Then we deduce from Proposition 3.7 the following equalities, where the second one is the stochastic version of the Bessel-Parseval theorem,

$$ \begin{aligned} \mathbb{E} \left[ \left[ \int_0^{T_1} f(W(s)) ds \right] \right] &= c_0 = \int_0^1 f(\theta) d\theta, \\ \mathbb{E} \left[ \left[ \int_0^{T_1} |f(W(s))|^2 ds \right] \right] &= \sum_{k \in \mathbb{Z}} |c_k|^2, \\ \mathbb{E} \left[ \left[ \left| \int_0^{T_1} f(W(s)) ds \right|^2 \right] \right] &= \frac{5 |c_0|^2}{3} + \sum_{k \in \mathbb{Z}*} \frac{|c_k|^2}{\pi^2 k^2}. \end{aligned} $$

## 3.3 Asymptotic expansion at revolution times and limit model

With the results of subsection 3.2, it is now possible to evaluate the local expansions (3.4) at revolution times. To approximate numerically the integrals appearing in equation (3.4) without evaluating $F$ and $F'$ too many times, we first replace the 1-periodic functions $g_\theta^0(y)$ and $g_\theta^1(y)(z)$ defined in (2.3) by their associated Fourier series with Fourier coefficients $(c_k^0(y))_{k \in \mathbb{Z}}$ and $(c_p^1(y))_{p \in \mathbb{Z}}$. We define the following approximation of $\varphi_{\varepsilon,t}(y)$,

$$ \begin{aligned} \psi_{\varepsilon,t}(y) = e^{AW(t)} y &+ \varepsilon \sum_{k \in \mathbb{Z}} e^{AW(t)} c_k^0(y) \int_0^t e^{2i\pi k W(s)} ds \quad (3.11) \\ &+ \varepsilon^2 \sum_{p,k \in \mathbb{Z}} e^{AW(t)} c_p^1(y) \left( c_k^0(y) \int_0^t e^{2i\pi p W(s)} e^{2i\pi k W(r)} dr ds \right). \end{aligned} $$

Notice that $c_k^0(y) \in \mathbb{C}^d$ and $c_p^1(y) = (c_k^0)'(y) \in \mathbb{C}^{d \times d}$ but $\psi_{\varepsilon,t}(y) \in \mathbb{R}^d$. We now evaluate this function $\psi_{\varepsilon,t}(y)$ at time $t = T_N$ to get a second order strong approximation.
**Proposition 3.9.** We define the following quantity

$$ \psi_{\varepsilon,N}(y) = y + H \sum_{k \in \mathbb{Z}} c_k^0(y) \alpha_k^N + H^2 \sum_{p,k \in \mathbb{Z}} c_p^1(y) (c_k^0(y)) \beta_{p,k}^N, $$

where $(c_k^0(y))_{k \in \mathbb{Z}}$ and $(c_p^1(y))_{p \in \mathbb{Z}}$ are the Fourier coefficients of the 1-periodic functions $g_\theta^0(y)$ and $g_\theta^1(y)$ defined in (2.3), $\alpha_k^N$, $\beta_{p,k}^N$ are the random variables defined in Proposition 3.7 and $y \in \mathbb{R}^d$ is deterministic. Under Assumption 3.1, for all test functions $\phi \in C_P^3$, there exists $H_0 > 0$ such that for all $H = N\varepsilon \le H_0$, the following estimates hold, where C and K are independent of $\varepsilon$ and N,

$$ \mathbb{E}\left[\left|\varphi_{\varepsilon,T_N}(y) - \psi_{\varepsilon,N}(y)\right|^2\right]^{1/2} \le C(1 + |y|^K)H^3, \quad (3.12) $$

$$ |\mathbb{E}[\phi(\varphi_{\varepsilon,T_N}(y))] - \mathbb{E}[\phi(\psi_{\varepsilon,N}(y))]| \le C(1 + |y|^K)H^3, \quad (3.13) $$

that is, $\psi_{\varepsilon,N}(y)$ is a numerical approximation of $\varphi_{\varepsilon,T_N}(y)$ of strong/weak local order two.

*Proof.* Inequality (3.12) is a straightforward consequence of Proposition 3.6 when evaluating the estimates of Proposition 3.3 at time $T_N$. For the weak local estimate (3.13), using inequality (3.12), the mean value inequality, Lemma 3.2 and equations (3.5) and (3.6), we get

$$ \begin{aligned} & \mathbb{E}\left[\left|\phi(\varphi_{\varepsilon,T_N}(y)) - \phi(\psi_{\varepsilon,T_N}^j(y))\right|\right] \\ & \le \mathbb{E}\left[\sup_{x \in [\varphi_{\varepsilon,T_N}(y), \psi_{\varepsilon,T_N}^j(y)]} |\phi'(x)| |\varphi_{\varepsilon,T_N}(y) - \psi_{\varepsilon,T_N}^j(y)|\right] \\ & \le C(1 + |y|^K)\mathbb{E}\left[e^{C\varepsilon T_N} (\varepsilon T_N)^{j+1} \sup_{x \in [\varphi_{\varepsilon,T_N}(y), \psi_{\varepsilon,T_N}^j(y)]} (1 + |x|^p)\right] \\ & \le C(1 + |y|^K)\mathbb{E}\left[(\varepsilon T_N)^{j+1} e^{C\varepsilon T_N} \left(1 + |\varphi_{\varepsilon,T_N}(y)|^p + |\psi_{\varepsilon,T_N}^j(y)|^p\right)\right] \\ & \le C(1 + |y|^K)\mathbb{E}\left[(\varepsilon T_N)^{j+1} e^{C\varepsilon T_N}\right]. \end{aligned} $$

Finally we obtain inequality (3.13) by taking $H$ small enough so that we can apply Proposition 3.6. $\square$

For a fixed $T = N\varepsilon$, when $\varepsilon \to 0$ (or equivalently $N \to \infty$), the solution of (3.1) evaluated at stroboscopic times $T_N = T_{T\varepsilon^{-1}}$ converges weakly to the solution of a deterministic ODE that involves only the first mode $c_0^0 = \langle g^0 \rangle = \int_0^1 g_\theta^0 d\theta$ of $g^0$. This asymptotic model is the same one as for the deterministic equation (2.1). The proof is postponed to subsection 4.3.

**Proposition 3.10 (asymptotic model).** Under Assumption 3.1, for $T > 0$, the solution $\varphi_{\varepsilon,T_{T\varepsilon^{-1}}}(X_0)$ (for $\varepsilon$ such that $T\varepsilon^{-1}$ is an integer) of equation (1.1) converges weakly when $\varepsilon \to 0$ to the solution at time T of

$$ \frac{dy_t}{dt} = \langle g^0 \rangle(y_t), \quad y_0 = X_0, \quad (3.14) $$

that is, for all test function $\phi \in C_P^3$,

$$ \lim_{\varepsilon \to 0} |\mathbb{E}[\phi(\varphi_{\varepsilon,T_{T\varepsilon^{-1}}}(X_0))] - \mathbb{E}[\phi(y_T)]| = 0. $$
**Remark 3.11.** It can be proven using the results of section 4 that the solution of the asymptotic model (Proposition 3.10) is an order one weak approximation of $X(\varepsilon T_{Nm})$ for $m \ge 0$ and $X$ solution of equation (1.1). We deduce the following simple one-step explicit deterministic integrator that corresponds to the Euler method applied to equation (3.14),

$$y_0 = X_0, \quad y_{m+1} = y_m + Hc_0^0(y_m). \qquad (3.15)$$

Its cost is independent of $\varepsilon$ and $N$, and it has weak order one w.r.t. $H$, that is, for all $m \ge 0$, $\mathbb{E}[\phi(\varphi_{\varepsilon,T_{Nm}}(X_0))] - \mathbb{E}[\phi(y_m)] = \mathcal{O}(H)$.

## 3.4 Construction of the second order integrators

To obtain an integrator of weak order two with a cost independent of $\varepsilon$ and $N$, we truncate the local expansion of Proposition 3.9. We also replace the involved random variables with cheap discrete random variables with the same first and second moments. To simulate the random variable $\alpha_k^N$ with discrete random variables $\hat{\alpha}_k^N$ with the same first and second moments, we introduce a set $(\xi_k)_{k \in \mathbb{N}}$ of independent random variables such that $\mathbb{P}(\xi_k = \pm 1) = \frac{1}{2}$, the covariance matrix $(C_\alpha^N)_{p,k}$ such that

$$(C_{\alpha}^{N})_{2p-1:2p,2k-1:2k} = \begin{pmatrix} \text{Cov}(\text{Re}(\alpha_p^N), \text{Re}(\alpha_k^N)) & \text{Cov}(\text{Re}(\alpha_p^N), \text{Im}(\alpha_k^N)) \\ \text{Cov}(\text{Im}(\alpha_p^N), \text{Re}(\alpha_k^N)) & \text{Cov}(\text{Im}(\alpha_p^N), \text{Im}(\alpha_k^N)) \end{pmatrix},$$

and $\Gamma^N$ its square root. Then, $\hat{\alpha}_k^N$ is defined for $k \ge 0$ as

$$\hat{\alpha}_k^N = \delta_k + \sum_{l \in \mathbb{N}} (\Gamma_{2k-1,l}^N + i\Gamma_{2k,l}^N)\xi_l \quad \text{with} \quad \delta_k = \begin{cases} 1 & \text{if } k=0 \\ 0 & \text{else} \end{cases}$$

and we fix $\hat{\alpha}_k^N = \overline{\hat{\alpha}_k^N}$ for $k < 0$ (so that the solution stays real while still having the good moments). We also define $\hat{\beta}_{p,k}^N = \mathbb{E}[\hat{\beta}_{p,k}^N]$ with the values of Proposition 3.7. Doing so yields Method A.

For Method B, we adapt Method A in the spirit of the middle point scheme for ODEs (see [16, Chap. IV]) so that it preserves any quadratic invariant. We also replace $\hat{\beta}_{p,k}^N$ by $\tilde{\beta}_{p,k}^N = \mathbb{E}[\tilde{\beta}_{p,k}^N]$, using the values of Proposition 3.7.

**Remark 3.12.** The methodology presented in section 3 can be generalized to any order. Thus, under more regularity assumptions on $F$, it is possible to build algorithms similar to Method A of any weak order and that are still robust with respect to the stiffness $\varepsilon$. For order 3, Method A becomes

$$Y_{m+1} = Y_m + H \sum_{k \in \mathbb{Z}} c_k^0(Y_m) \hat{\alpha}_k^N + H^2 \sum_{p,k \in \mathbb{Z}} c_p^1(Y_m)(c_k^0(Y_m)) \hat{\beta}_{p,k}^N \\ + H^3 \sum_{l,p,k \in \mathbb{Z}} c_l^1(c_p^1(c_k^0))(Y_m) \hat{\gamma}_{l,p,k}^{(1),N} + c_l^2(c_p^0, c_k^0)(Y_m) \hat{\gamma}_{l,p,k}^{(2),N}$$

with the new random variables

$$\gamma_{l,p,k}^{(1),N} = \frac{1}{N^2} \int_0^{T_N} e^{2i\pi qW(s)} \int_0^s e^{2i\pi pW(r)} \int_0^r e^{2i\pi kW(q)} dqdrds,$$
$$ \gamma_{l,p,k}^{(2),N} = \frac{1}{2N^2} \int_0^{T_N} e^{2i\pi qW(s)} \int_0^s e^{2i\pi pW(r)} \int_0^s e^{2i\pi kW(q)} dq dr ds, $$

and where the discrete random variables share the same moments up to order 3 for the $\hat{\alpha}_k^N$, order 2 for the $\hat{\beta}_{p,k}^N$, and order 1 for the $\hat{\gamma}_{l,p,k}^{(i),N}$. It is also possible to generalize Method B up to any order in the spirit of the middle-point scheme, but the construction of discrete random variables allowing the preservation of quadratic invariants is not obvious for higher orders (although backward error analysis guarantees the preservation of quadratic invariants for the exact random variables based on W).

# 4 Weak convergence analysis

This section focuses on the proofs of the following two theorems, showing the order two convergence of Methods A and B.

**Theorem 4.1.** Assume that the Fourier coefficients $c_k^0$, $c_p^1$ of $g_\theta^0$ and $g_\theta^1$ in (2.3) are non-zero only for $-K_t/2 \le k, p < K_t/2$. Then, under Assumption 3.1, Method A has weak order two, that is, for all $T > 0$, for all test functions $\phi \in C_p^3$, there exists $H_0 > 0$ such that for all $H \le H_0$, for all $m \ge 0$ such that $mN\epsilon = mH \le T$, there exist two positive constants $K$ and $C$ both independent of $\epsilon, N$ and $K_t$ such that

$$ |\mathbb{E}[\phi(\varphi_{\epsilon, T_{Nm}}(X_0))] - \mathbb{E}[\phi(Y_m)]| \le CH^2(1 + \mathbb{E}[|X_0|^K]). \quad (4.1) $$

**Theorem 4.2.** Assume that the Fourier coefficients $c_k^0$, $c_p^1$ of $g_\theta^0$ and $g_\theta^1$ in (2.3) are non-zero only for $-K_t/2 \le k, p < K_t/2$. Then, under Assumption 3.1, if $c_0^1(c_k^0)$ and $c_p^1(c_0^0)$ are Lipschitz continuous uniformly in $k$ and $p$, Method B is well defined and has weak order two (i.e., it satisfies an estimate of the form (4.1)). In addition, if for a fixed symmetric matrix $S \in \mathbb{R}^{d \times d}$, the quantity $Q(y) = \frac{1}{2}y^T Sy$ is preserved by equation (1.1), then Method B also preserves the invariant $Q(y) = \frac{1}{2}y^T Sy$, that is, the solution $Y_{m+1}$ of equation (2.5) satisfies $Q(Y_{m+1}) = Q(Y_m)$.

These two theorems focus on approximating the exact solution of equation (1.1) at the revolution times $\epsilon T_{Nm}$, $m = 0, 1, \dots$, but one could compute an approximation at different times by composing with other methods at the end of the integration.

**Remark 4.3.** Since the error constant $C$ in (4.1) is independent of the number $K_t$ of Fourier modes, we emphasize that Theorems 4.1 and 4.2 remain valid for infinitely many modes ($K_t \to \infty$). In addition, assuming that $F$ is of class $C_p^{s+1}$ yields a truncation error of the Fourier series in (2.3) of size $\mathcal{O}((1 + |y|^K)K_t^{-s})$ (see, e.g., [19, Sect. III.1.3]), and if $g_\theta^0$ is assumed analytic in $\theta$ (for example if $F$ is a polynomial), this error becomes exponentially small as $\mathcal{O}((1 + |y|)e^{-cK_t})$. For simplicity of the analysis, we thus assumed in Theorems 4.1 and 4.2 that $g_\theta^0$ and $g_\theta^1$ have a finite number $K_t$ of non-zero Fourier modes in (2.3). If this assumption does not hold, the truncation errors $\mathcal{O}((1 + |y|^K)K_t^{-s})$ or $\mathcal{O}((1 + |y|)e^{-cK_t})$ should be added in the right-hand side of the error estimate (4.1). Let us prove it in the analytic case. We first apply the change of variable $\tilde{\varphi}_{\epsilon,t}(y) = e^{-AW(t)}\varphi_{\epsilon,t}(y)$ that has no effect at time $t = T_{Nm}$. We now have to compare the two solutions of the following integral formulations

$$ \tilde{\varphi}_{\epsilon,t}(y) = y + \epsilon \int_0^t g_{W(s)}^0(\tilde{\varphi}_{\epsilon,s}(y))ds, $$
$$ \tilde{\varphi}_{\varepsilon,t}^{(K_t)}(y) = y + \varepsilon \int_0^t \sum_{k=-K_t/2}^{K_t/2-1} c_k^0(\tilde{\varphi}_{\varepsilon,s}^{(K_t)}(y)) e^{2i\pi k W(s)} ds. $$

Using the truncation estimates that we previously discussed and the Lipschitz property of $g_\theta^0$, one gets

$$ \begin{aligned} \left| \tilde{\varphi}_{\varepsilon,t} - \tilde{\varphi}_{\varepsilon,t}^{(K_t)} \right| (y) &\le \varepsilon \int_0^t \left| g_{W(s)}^0(\tilde{\varphi}_{\varepsilon,s}(y)) - g_{W(s)}^0(\tilde{\varphi}_{\varepsilon,s}^{(K_t)}(y)) \right| ds + C \varepsilon t e^{-cK_t} \sup_{[0,1]} |g^0(y)| \\ &\le C \varepsilon \int_0^t \left| \tilde{\varphi}_{\varepsilon,s}(y) - \tilde{\varphi}_{\varepsilon,s}^{(K_t)}(y) \right| ds + C \varepsilon t (1 + |y|) e^{-cK_t}. \end{aligned} $$

The Gronwall lemma and Proposition 3.6 yield for $mN\varepsilon \le T$,

$$ \mathbb{E} \left[ \left| \tilde{\varphi}_{\varepsilon, T_{Nm}}(y) - \tilde{\varphi}_{\varepsilon, T_{Nm}}^{(K_t)}(y) \right|^2 \right]^{1/2} \le C(1 + |y|)e^{-cK_t}. $$

The structure of the convergence proof is similar to the one for standard SDE integrators, see e.g. [22, Chap. 2], but one has to be cautious because our solution is evaluated at stochastic times and the error constants should not depend on $\varepsilon$ or $N$.

## 4.1 Boundedness of the numerical moments

**Proposition 4.4** (bounded moments for the integrator (2.4)). Assume that for $y \in \mathbb{R}^d$, the numerical integrator $\hat{\psi}_{\varepsilon,N}(y)$ is given by

$$ \hat{\psi}_{\varepsilon,N}(y) = y + H \sum_{k \in \mathbb{Z}} c_k^0(y) \hat{\alpha}_k^N + H^2 \sum_{p,k \in \mathbb{Z}} c_p^1(y) (c_k^0(y)) \hat{\beta}_{p,k}^N, \quad (4.2) $$

where $\hat{\alpha}_k^N$, $\hat{\beta}_{p,k}^N$ are random variables defined such that for all $q > 0$, $\mathbb{E} \left[ \left( \sum_k \frac{|\hat{\alpha}_k^N|^2}{k^2} \right)^q \right]$ and $\mathbb{E} \left[ \left( \sum_{p,k} \frac{|\hat{\beta}_{p,k}^N|^2}{k^2} \right)^q \right]$ are bounded uniformly in $N$. Then, under Assumption 3.1 and if $|y|$ has bounded moments, for any $T > 0$, for all $m$, $H$ such that $m\varepsilon N = mH \le T$, for all $q > 0$, we have $\mathbb{E}[|\hat{\psi}_{\varepsilon,N}^m(y)|^{2q}] \le C_q(1 + \mathbb{E}[|y|^{2q}])$, where $C_q$ is independent of $m$, $\varepsilon$ and $N$.

*Proof.* We first prove

$$ |\hat{\psi}_{\varepsilon,N}(y) - y| \le CH(1 + |y|)M_N, \quad (4.3) $$

where $\mathbb{E}[(M_N)^{2q}] \le C_q$ for all $q > 0$. We have

$$ \begin{aligned} |\hat{\psi}_{\varepsilon,N}(y) - y| &= H \left| \sum_{k \in \mathbb{Z}} c_k^0(y) \hat{\alpha}_k^N + H \sum_{p,k \in \mathbb{Z}} c_p^1(y) (c_k^0(y)) \hat{\beta}_{p,k}^N \right| \\ &\le CH \left( M_N^{(0)} \sqrt{\sum_{k \in \mathbb{Z}} k^2 |c_k^0(y)|^2} + M_N^{(1)} \sqrt{\sum_{p \in \mathbb{Z}} |c_p^1(y)|^2} \sqrt{\sum_{k \in \mathbb{Z}} k^2 |c_k^0(y)|^2} \right) \end{aligned} $$

where $M_N^{(0)} = \sqrt{\sum_k \frac{|\hat{\alpha}_k^N|^2}{k^2}}$ and $M_N^{(l)} = \sqrt{\sum_{p,k} \frac{|\hat{\beta}_{p,k}^N|^2}{k^2}}$ have moments bounded uniformly in $N$. Then using the Bessel-Parseval theorem, we get $\sum_k k^2 |c_k^0(y)|^2 = \int_0^1 |\partial_\theta g_\theta(y)^0|^2 d\theta$.
Assumption 3.1 yields $|\partial_{\theta}g_{\theta}^{0}(y)| \le C(1 + |y|)$. Then, the Bessel-Parseval theorem applied on $g_{\theta}^{1}$ gives $\sqrt{\sum_{p} |c_{p}^{1}(y)|^{2}} \le C$, hence the result.

We define $\Delta\psi_{m} = \hat{\psi}_{\varepsilon,N}^{m+1}(y) - \hat{\psi}_{\varepsilon,N}^{m}(y) = (\hat{\psi}_{\varepsilon,N} - \text{Id})(\hat{\psi}_{\varepsilon,N}^{m}(y))$, then

$$ (\hat{\psi}_{\varepsilon,N}^{m+1}(y))^{2q} = (\hat{\psi}_{\varepsilon,N}^{m}(y))^{2q} + \sum_{j=1}^{2q} \binom{2l}{j} (\hat{\psi}_{\varepsilon,N}^{m}(y))^{2q-j} \Delta\psi_{m}^{j}. $$

Equation (4.3) and the bounded moments of $M_{N}$ give

$$ \begin{aligned} \left| \mathbb{E} \left[ (\hat{\psi}_{\varepsilon,N}^{m}(y))^{2q-j} \Delta \psi_{m}^{j} \right] \right| &\le \mathbb{E} \left[ \left| \hat{\psi}_{\varepsilon,N}^{m}(y) \right|^{2q-j} C H^{j} \left( 1 + \left| \hat{\psi}_{\varepsilon,N}^{m}(y) \right|^{j} \right) M_{N}^{j} \right] \\ &\le C_{q} H \left( 1 + \mathbb{E} \left[ \left| \hat{\psi}_{\varepsilon,N}^{m}(y) \right|^{2q} \right] \right). \end{aligned} $$

We deduce

$$ 1 + \mathbb{E} \left[ \left| \hat{\psi}_{\varepsilon,N}^{m+1}(y) \right|^{2q} \right] \le e^{C_q H} \left( 1 + \mathbb{E} \left[ \left| \hat{\psi}_{\varepsilon,N}^{m}(y) \right|^{2q} \right] \right) $$

$$ \text{and by induction } \mathbb{E} \left[ \left| \hat{\psi}_{\varepsilon,N}^{m}(y) \right|^{2q} \right] \le e^{C_q m H} (1 + \mathbb{E}[|y|^{2q}]) \le e^{C_q T} (1 + \mathbb{E}[|y|^{2q}]). \quad \square $$

**Proposition 4.5** (bounded moments for the integrator (2.5)). Assume that for $y \in \mathbb{R}^{d}$, the numerical scheme $\hat{\psi}_{\varepsilon,N}(y)$ satisfies

$$ \begin{aligned} \hat{\psi}_{\varepsilon,N}(y) = y + H \sum_{k \in \mathbb{Z}} c_k^0 \left( \frac{y + \hat{\psi}_{\varepsilon,N}(y)}{2} \right) \hat{\alpha}_k^N & \qquad (4.4) \\ & + H^2 \sum_{p,k \in \mathbb{Z}} c_p^1 \left( \frac{y + \hat{\psi}_{\varepsilon,N}(y)}{2} \right) \left( c_k^0 \left( \frac{y + \hat{\psi}_{\varepsilon,N}(y)}{2} \right) \right) \hat{\beta}_{p,k}^N, \end{aligned} $$

where $\hat{\alpha}_{k}^{N}$, $\hat{\beta}_{p,k}^{N}$ are random variables defined such that for all $q > 0$, $\sum_{k} |\hat{\alpha}_{k}^{N}|$, $\sum_{p,k} |\hat{\beta}_{p,k}^{N}|$, $\mathbb{E}\left[\left(\sum_{k} \frac{|\hat{\alpha}_{k}^{N}|^2}{k^2}\right)^q\right]$ and $\mathbb{E}\left[\left(\sum_{p,k} \frac{|\hat{\beta}_{p,k}^{N}|^2}{k^2}\right)^q\right]$ are bounded uniformly in $N$. Then, under Assumption 3.1 and if $|y|$ has bounded moments, for $H_0$ small enough and any $T > 0$, for all $m$, $H$ such that $m \varepsilon N = mH \le T$ and $H \le H_0$, for all $q > 0$, we have $\mathbb{E}[|\hat{\psi}_{\varepsilon,N}^{m}(y)|^{2q}] \le C_q(1 + \mathbb{E}[|y|^{2q}])$, where $C_q$ is independent of $m$, $\varepsilon$ and $N$.

*Proof.* We prove an equivalent of the estimate (4.3) for $\hat{\psi}_{\varepsilon,N}(y)$. The growth properties of the Fourier coefficients yield

$$ \begin{aligned} \left| \hat{\psi}_{\varepsilon,N}(y) - y \right| &\le CH \left( \sum_{k} |c_k^0|(y) |\hat{\alpha}_k^N| + \sum_{p,k} |c_k^0|(y) |\hat{\beta}_{p,k}^N| \right) \\ &\quad + CH \left( \sum_{k} |\hat{\alpha}_k^N| + \sum_{p,k} |\hat{\beta}_{p,k}^N| \right) \left| \hat{\psi}_{\varepsilon,N}(y) - y \right|, \end{aligned} $$
As $\sum_k |\hat{\alpha}_k^N| + \sum_{p,k} |\hat{\beta}_{p,k}^N|$ is bounded, using the same estimates as in the proof of Proposition 4.4, we get for all $H \le H_0$ small enough,

$$|\hat{\psi}_{\varepsilon,N}(y) - y| \le CH(1 + |y|)M_N, \quad (4.5)$$

where $M_N$ has bounded moments. The remaining of the proof is the same as in the proof of Proposition 4.4. □

## 4.2 Local weak error

**Proposition 4.6 (local error estimate).** Assume that for $y \in \mathbb{R}^d$ deterministic, the numerical scheme can be written as

$$\hat{\psi}_{\varepsilon,N}(y) = y + H \sum_{k \in \mathbb{Z}} c_k^0(y) \hat{\alpha}_k^N + H^2 \sum_{p,k \in \mathbb{Z}} c_p^1(y) (c_k^0(y)) \hat{\beta}_{p,k}^N + R,$$

where $\mathbb{E}[|R|] \le C(1 + |y|^K)H^3$ and $\hat{\alpha}_k^N \in \mathbb{C}$, $\hat{\beta}_{p,k}^N \in \mathbb{R}$ are random variables such that $\hat{\alpha}_k^N = \hat{\alpha}_{-k}^N$ and

$$\mathbb{E}[\hat{\alpha}_k^N] = \mathbb{E}[\alpha_k^N], \quad \mathbb{E}[\hat{\beta}_{p,k}^N] = \mathbb{E}[\beta_{p,k}^N], \quad \mathbb{E}[\hat{\alpha}_{k_1}^N \hat{\alpha}_{k_2}^N] = \mathbb{E}[\alpha_{k_1}^N \alpha_{k_2}^N].$$

Under Assumption 3.1, if $\hat{\psi}_{\varepsilon,N}(y)$ satisfies the assumptions of Proposition 4.4 (or Proposition 4.5), for all test function $\phi \in C_P^3$, there exists $H_0 > 0$ such that for all $H = N\varepsilon \le H_0$, the following estimate holds, where C and K are independent of $\varepsilon$ and N,

$$|\mathbb{E}[\phi(\varphi_{\varepsilon,T_N}(y))] - \mathbb{E}[\phi(\hat{\psi}_{\varepsilon,N}(y))]| \le C(1 + |y|^K)H^3,$$

that is, the numerical scheme has weak local order two.

*Proof.* Using Proposition 3.9 and its notation $\psi_{\varepsilon,N}(y)$, it is enough to prove that

$$|\mathbb{E}[\phi(\psi_{\varepsilon,N}(y))] - \mathbb{E}[\phi(\hat{\psi}_{\varepsilon,N}(y))]| \le C(1 + |y|^K)H^3.$$

A local expansion gives

$$\phi(\psi_{\varepsilon,N}(y)) = \phi(y) + \phi'(y)(\psi_{\varepsilon,N}(y) - y) + \phi''(y)(\psi_{\varepsilon,N}(y) - y, \psi_{\varepsilon,N}(y) - y) + R_1.$$

As $\psi_{\varepsilon,N}(y) = \psi_{\varepsilon,T_N}^2(y)$ (see equation (3.4)), using inequalities (3.6), (3.7) evaluated at $T_N$ and Proposition 3.6 yields

$$\begin{aligned}
\mathbb{E}[|R_1|] & \le \mathbb{E} \left[ \sup_{x \in [y, \psi_{\varepsilon,N}(y)]} |\phi^{(3)}(x)| |\psi_{\varepsilon,N}(y) - y|^3 \right] \\
& \le \mathbb{E} \left[ C(1 + |y|^K + |\psi_{\varepsilon,N}(y)|^K)(1 + |y|^K)(\varepsilon T_N)^3 e^{C \varepsilon T_N} \right] \\
& \le \mathbb{E} \left[ C(1 + |y|^K)(\varepsilon T_N)^3 e^{C \varepsilon T_N} \right] \\
& \le C(1 + |y|^K)H^3.
\end{aligned}$$
We obtain a similar expansion for $\phi(\hat{\psi}_{\varepsilon,N}(y))$:

$$ \phi(\hat{\psi}_{\varepsilon,N}(y)) = \phi(y) + \phi'(y)(\hat{\psi}_{\varepsilon,N}(y) - y) + \phi''(y)(\hat{\psi}_{\varepsilon,N}(y) - y, \hat{\psi}_{\varepsilon,N}(y) - y) + \hat{R}_1, $$

where, using inequality (4.3) (or (4.5)),

$$ \begin{aligned} \mathbb{E}[|\hat{R}_1|] &\le \mathbb{E} \left[ \sup_{x \in [y, \hat{\psi}_{\varepsilon,N}(y)]} |\phi^{(3)}(x)| |\hat{\psi}_{\varepsilon,N}(y) - y|^3 \right] \\ &\le C \mathbb{E} \left[ (1 + |y|^K + |\hat{\psi}_{\varepsilon,N}(y)|^K)(1 + |y|^K) H^3 M_N^3 \right] \\ &\le C(1 + |y|^K) H^3 \mathbb{E}[(1 + M_N^K)] \\ &\le C(1 + |y|^K) H^3. \end{aligned} $$

Making the difference of both equations gives

$$ \begin{aligned} \phi(\psi_{\varepsilon,N}(y)) - \phi(\hat{\psi}_{\varepsilon,N}(y)) = & \phi'(y)(\psi_{\varepsilon,N}(y) - \hat{\psi}_{\varepsilon,N}(y)) - \phi''(y)(\psi_{\varepsilon,N}(y) - y)^2 & (4.6) \\ & + \phi''(y)(\hat{\psi}_{\varepsilon,N}(y) - y)^2 + R, \end{aligned} $$

where $\mathbb{E}[|R|] \le C(1 + |y|^K)H^3$. For the first term of (4.6), we have

$$ \begin{aligned} \mathbb{E}[\phi'(y)(\psi_{\varepsilon,N}(y) - \hat{\psi}_{\varepsilon,N}(y))] = & H \sum_{k \in \mathbb{Z}} \mathbb{E}[\phi'(y)(c_k^0(y)(\alpha_k^N - \hat{\alpha}_k^N))] \\ & + H^2 \sum_{p,k \in \mathbb{Z}} \mathbb{E}[\phi'(y)(c_p^1(y)(c_k^0(y))(\beta_{p,k}^N - \hat{\beta}_{p,k}^N))]. \end{aligned} $$

Then, we get

$$ \mathbb{E}[\phi'(y)(c_k^0(y)(\alpha_k^N - \hat{\alpha}_k^N))] = \mathbb{E}[\alpha_k^N - \hat{\alpha}_k^N]\phi'(y)(c_k^0(y)) = 0. $$

We can do the same thing with the term in $\beta_{p,k}^N$ and obtain

$$ \mathbb{E}[\phi'(y)(\psi_{\varepsilon,N}(y) - \hat{\psi}_{\varepsilon,N}(y))] = 0. $$

Let us now study the second order term $Z = \phi''(y)(\hat{\psi}_{\varepsilon,N}(y) - y)^2 - \phi''(y)(\psi_{\varepsilon,N}(y) - y)^2$ that appears in (4.6). We develop this expression and keep only the order one and two terms to obtain $Z = H^2Y + R$ where $\mathbb{E}[|R|] \le C(1 + |y|^K)H^3$ (by the same arguments as before) and

$$ \begin{aligned} Y &= \sum_{k_1, k_2} [\phi''(y)(c_{k_1}^0(y)\hat{\alpha}_{k_1}^N, c_{k_2}^0(y)\hat{\alpha}_{k_2}^N) - \phi''(y)(c_{k_1}^0(y)\alpha_{k_1}^N, c_{k_2}^0(y)\alpha_{k_2}^N)] \\ &= \sum_{k_1, k_2} (\hat{\alpha}_{k_1}^N \hat{\alpha}_{k_2}^N - \alpha_{k_1}^N \alpha_{k_2}^N)\phi''(y)(c_{k_1}^0(y), c_{k_2}^0(y)) \end{aligned} $$

The condition on the moments of the $\hat{\alpha}_k^N$ yields $\mathbb{E}[Y] = 0$.

Putting all these arguments together in (4.6), we finally get that

$$ |\mathbb{E}[\phi(\psi_{\varepsilon,N}(y))] - \mathbb{E}[\phi(\hat{\psi}_{\varepsilon,N}(y))]| \le C(1 + |y|^K)H^3. $$

We deduce the local order two of the proposed numerical scheme. □
**Remark.** The constant $H_0$ in Proposition 4.6 depends on $F$, but also on the polynomial growth power of $\phi$ and its first three derivatives. This dependence is expected when trying to evaluate the solution of SDEs at random times. To make $H_0$ independent of the test functions, one can consider the following sets of test functions

$$C_{P,K}^3 = \{\phi \in C^3, \exists C > 0, \exists k \le K, \forall y, |\phi^{(i)}(y)| \le C(1 + |y|^k), i \in \{0, 1, 2, 3\}\}.$$

## 4.3 Global error

**Theorem 4.7 (global convergence).** Assume that the numerical scheme $\hat{\psi}_{\varepsilon,N}$ satisfies equation (4.2) (respectively equation (4.4)) where $\hat{\alpha}_k^N \in \mathbb{C}$, $\hat{\beta}_{p,k}^N \in \mathbb{R}$ (respectively $\hat{\beta}_{p,k}^N \in \mathbb{R}$) are random variables such that $\hat{\alpha}_k^N = \hat{\alpha}_{-k}^N$ and

$$\mathbb{E}[\hat{\alpha}_k^N] = \mathbb{E}[\alpha_k^N], \quad \mathbb{E}[\hat{\beta}_{p,k}^N] = \mathbb{E}[\beta_{p,k}^N], \quad \mathbb{E}[\hat{\alpha}_{k_1}^N \hat{\alpha}_{k_2}^N] = \mathbb{E}[\alpha_{k_1}^N \alpha_{k_2}^N].$$

(respectively $\hat{\alpha}_k^N$ satisfies the same conditions and $\hat{\beta}_{p,k}^N$ satisfies $\mathbb{E}[\hat{\beta}_{p,k}^N] = \mathbb{E}[\tilde{\beta}_{p,k}^N]$). Under Assumption 3.1, if for all $q > 0$, $\mathbb{E}\left[\left(\sum_k \frac{|\hat{\alpha}_k^N|^2}{k^2}\right)^q\right]$ and $\mathbb{E}\left[\left(\sum_{p,k} \frac{|\hat{\beta}_{p,k}^N|^2}{k^2}\right)^q\right]$ are bounded uniformly in $N$ (respectively $\sum_k |\hat{\alpha}_k^N|$, $\sum_{p,k} |\hat{\beta}_{p,k}^N|$, $\mathbb{E}\left[\left(\sum_k \frac{|\hat{\alpha}_k^N|^2}{k^2}\right)^q\right]$ and $\mathbb{E}\left[\left(\sum_{p,k} \frac{|\hat{\beta}_{p,k}^N|^2}{k^2}\right)^q\right]$ are bounded uniformly in $N$), for all $T > 0$, for all test function $\phi \in C_P^3$, there exists $H_0 > 0$ such that for all $H \le H_0$, for all $M \ge 0$ such that $MN\varepsilon = MH \le T$, there exist two positive constants $K$ and $C$ both independent of $\varepsilon$ and $N$ such that

$$\left|\mathbb{E}[\phi(\varphi_{\varepsilon,T_{NM}}(X_0))] - \mathbb{E}[\phi(\hat{\psi}_{\varepsilon,N}^M(X_0))]\right| \le C H^2(1 + \mathbb{E}[|X_0|^K]).$$

*Proof.* We denote

$$e_M = \mathbb{E}[\phi(\varphi_{\varepsilon,T_{NM}}(X_0))] - \mathbb{E}[\phi(\hat{\psi}_{\varepsilon,N}^M(X_0))]$$

and rewrite it with a telescopic sum

$$e_M = \sum_{m=1}^{M} \mathbb{E}[\phi(\varphi_{\varepsilon, T_{N(m-1)}}(\hat{\psi}_{\varepsilon,N}^{M-m+1}(X_0)))] - \mathbb{E}[\phi(\varphi_{\varepsilon, T_{Nm}}(\hat{\psi}_{\varepsilon,N}^{M-m}(X_0)))] \\ = \sum_{m=1}^{M} \mathbb{E}[\tilde{\phi}_{m-1}(\hat{\psi}_{\varepsilon,N}(\hat{\psi}_{\varepsilon,N}^{M-m}(X_0)))] - \mathbb{E}[\tilde{\phi}_{m-1}(\varphi_{\varepsilon, T_N}(\hat{\psi}_{\varepsilon,N}^{M-m}(X_0)))]$$

where $\tilde{\phi}_{m-1} = \phi \circ \varphi_{\varepsilon, T_{N(m-1)}}$. Using Lemma 3.2 and $\phi \in C_P^3$, we obtain for $0 \le i \le 3$,

$$|\tilde{\phi}_m^{(i)}(y)| \le C e^{C\varepsilon T_{Nm}} (1 + |y|^K).$$

Thus, knowing the hitting times involved, $\tilde{\phi}_m \in C_P^3$. Using Assumption 3.1, $(c_p^0)' = c_p^1$ and $\beta_{p,k}^N = \tilde{\beta}_{p,k}^N + \frac{\alpha_p \alpha_k}{2}$, we deduce that $\hat{\psi}_{\varepsilon,N}$ satisfies the assumptions of Proposition 4.6. Applying Proposition 4.6 to each term of $e_M$ gives

$$|e_M| \le \sum_{m=1}^{M} C \mathbb{E}[e^{C\varepsilon T_{Nm}}] H^3 \left(1 + \mathbb{E}\left[\left|\hat{\psi}_{\varepsilon,N}^{M-m}(X_0)\right|^K\right]\right).$$
Finally, the moments of $\hat{\psi}_{\varepsilon,N}^m(X_0)$ are all bounded uniformly in $\varepsilon$, $N$ and $m$ according to Proposition 4.4 (respectively 4.5). Thus

$$|e_M| \le \sum_{m=1}^{M} CH^3(1 + \mathbb{E}[|X_0|^K]) \le CH^2(1 + \mathbb{E}[|X_0|^K]).$$

We deduce the global weak order two. □

With the help of Theorem 4.7, we prove Proposition 3.10 and the convergence of Methods A and B.

*Proof of Proposition 3.10.* Rewriting Theorem 4.7 for order one yields for all $H = N\varepsilon$ small enough and all $M \ge 0$,

$$|\mathbb{E}[\phi(\varphi_{\varepsilon,T_{NM}}(X_0))] - \mathbb{E}[\phi(y_{\varepsilon NM})]| \le C(\varepsilon N)^2(1 + \mathbb{E}[|X_0|^K]).$$

Evaluating in $N = 1$, $M = \frac{T}{\varepsilon}$ and taking the limit $\varepsilon \to 0$ yield the result. □

*Proof of Theorem 4.1.* As $\hat{\alpha}_k^N \le C$ and $\sum_{p,k} \frac{|\mathbb{E}[\beta_{p,k}^N]|^2}{k^2}$ converges by Proposition 3.7, Theorem 4.7 applies and concludes the proof. □

*Proof of Theorem 4.2.* The regularity assumptions yield the Lipschitzness of the $c_k^0(y)$ and the involved $c_p^1(y)(c_k^0(y))$ with constants independent of $k$ and $p$. As $\sum_k |\hat{\alpha}_k^N|$ and $\sum_{p,k} |\hat{\beta}_{p,k}^N|$ are bounded, the right-hand side of equation (2.5) is a contraction for all $H \le H_0$ small enough and the constant does not depend on $Y_m$, so $H_0$ depends only on $F$ and $F'$. Thus, the integrator is well-posed for all $H \le H_0$.

The weak order two is obtained using Theorem 4.7. Indeed the use of discrete random variables and Proposition 3.7 gives the convergence of the involved series.

For showing that Method B preserves quadratic invariants, it is sufficient to prove that $Q'(y)(\sum_k c_k^0(y)\hat{\alpha}_k^N) = 0$ and $Q'(y)(\sum_{p,k} c_p^1(y)(c_k^0(y))\hat{\beta}_{p,k}^N) = 0$ (see [16, Chap. IV]). The preservation of $Q$ by equation (1.1) yields $Q'(y)(Ay) = 0$ and $Q'(y)(F(y)) = 0$. We deduce the following two equations, valid for all $y \in \mathbb{R}^d$,

$$y^T S g_{\theta}^0(y) = 0, \qquad (4.7)$$

$$y^T S g_{\theta}^1(y)(g_{\nu}^0(y)) = -(g_{\nu}^0(y))^T S g_{\theta}^0(y), \qquad (4.8)$$

where equation (4.8) is obtained by differentiating equation (4.7) in the direction $g_{\nu}^0$. Using equation (4.7), we have

$$Q'(y) \left( \sum_k c_k^0(y) \hat{\alpha}_k^N \right) = \int_0^1 Q'(y)(g_{\theta}^0(y)) \sum_k e^{-2i\pi k \theta} \hat{\alpha}_k^N d\theta = 0.$$

For the second order term, equation (4.8) and the values of Proposition 3.7 yield
$$
\begin{align*}
Q'(y) \left( \sum_{p,k} c_p^1(y) (c_k^0(y)) \hat{\beta}_{p,k}^N \right) \\
&= \int_0^1 \int_0^1 y^T S g_\theta^1(y) (g_\nu^0(y)) \sum_{p,k} e^{-2i\pi p\theta} e^{-2i\pi k\nu} \hat{\beta}_{p,k}^N d\nu d\theta \\
&= - \int_0^1 \int_0^1 (g_\nu^0(y))^T S g_\theta^0(y) \sum_{p,k} e^{-2i\pi p\theta} e^{-2i\pi k\nu} \hat{\beta}_{p,k}^N d\nu d\theta \\
&= -\frac{1}{2} \int_0^1 \int_0^1 (g_\nu^0(y))^T S g_\theta^0(y) \sum_{p,k} [e^{-2i\pi p\theta} e^{-2i\pi k\nu} + e^{-2i\pi p\nu} e^{-2i\pi k\theta}] \hat{\beta}_{p,k}^N d\nu d\theta \\
&= 0.
\end{align*}
$$

Hence Method B is well-posed, has weak order 2 and preserves the invariant $Q$. $\square$

<figure><img src="image_5.png" /></figure>

# 5 Numerical experiments

In this section, we first illustrate numerically the weak order two of Methods A and B with convergence curves. Then, we apply the new algorithms to solve the nonlinear Schrödinger equation with highly-oscillatory white noise dispersion (1.3).

## 5.1 Weak order of convergence

To confirm the results of Theorem 4.1 and Theorem 4.2, we check numerically if Methods A and B have weak order two of accuracy w.r.t. $H$ uniformly in $\varepsilon$ and $N$. As the Euler-Maruyama method and the algorithms presented in [8, 4, 9] are completely inaccurate if they do not satisfy the severe timestep restriction $h \ll \varepsilon$, we compare the performance of Methods A and B to the performance of the Euler method (3.15). We first apply the algorithms on equation (1.1) with the linearity $F(y) = iy$, $A = 2i\pi$, $X_0 = 1$ and $\varepsilon = 10^{-3}$. Equivalently we can write it in the real setting as

$$dX = \frac{2\pi}{\sqrt{\varepsilon}} \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} X \circ dW + \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} Xdt, \quad X_0 = \begin{pmatrix} 1 \\ 0 \end{pmatrix}.$$

We plot on a logarithmic scale an estimate of the weak error for approximating $X$ at time $T = 10^{-3}T_2s$ where $E[T] = 0.256$. The exact solution $X(T)$ is approximated by the output of Method B for $H = \varepsilon$. The parameters $N$ and $m$ are varying under the condition that $Nm = 2^8$. The test function is $\phi(y) = 2y_1 + 4y_2$ and the average is taken over $10^7$ trajectories. We choose the tolerance $10^{-13}$ for the fixed point. On the right picture of Figure 4, we use a modification of a Kubo oscillator introduced in [8] with the nonlinearity $F(y) = i(1 + \mathrm{Re}(y)^3 + \mathrm{Im}(y)^5)y$. In the real setting, it yields the following two-dimensional SDE

$$dX = \frac{2\pi}{\sqrt{\varepsilon}} \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} X \circ dW + \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} (1 + X_1^3 + X_2^5) Xdt, \quad X_0 = \begin{pmatrix} 1 \\ 0 \end{pmatrix}.$$

We take 8 modes for the Fourier decomposition and the same other parameters as before. The average is taken over $10^6$ trajectories.
<figure><img src="image_6.png" /><figcaption>Figure 4: Weak error versus the stepsize $H = N\epsilon$ for approximating the solution of equation (1.1) at time $\epsilon T_2^s$ for the linear $F(y) = iy$ (left) and the non-linear $F(y) = i(1+\mathrm{Re}(y)^3+\mathrm{Im}(y)^5)y$ (right) with $A = 2i\pi$, $X_0 = 1$, $\epsilon = 10^{-3}$ and the test function $\phi(y) = 2\mathrm{Re}(y) + 4\mathrm{Im}(y)$.</figcaption></figure>

In both cases, we observe the weak order two of Methods A and B. The irregularities of the curve for a small $H$ come from Monte-Carlo errors. We repeated the same experiment on many other examples and we always observe the desired order two as long as $H$ is small enough.

## 5.2 Numerical experiments on NLS equation with white noise dispersion

We now apply the algorithms to solve on the torus $\mathbb{T} = [-\pi, \pi]$ the following SPDE of the form (1.3), with a polynomial linearity and the stiffness parameter $\epsilon = 10^{-2}$,

$$du = \frac{2\pi}{\sqrt{\epsilon}} \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} \Delta u \circ dW + \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} |u|^{2\sigma} u dt, \quad x \in \mathbb{T}, t > 0, \qquad (5.1)$$

where the unknown $u$ is a random process depending on $x \in \mathbb{T}$ and $t \ge 0$. We consider a spectral discretization in space of this equation with $K_x = 2^7$ modes $u(x,t) \approx \sum_{|l| \le K_x} Y_l(t) e^{ilx}$. We obtain an equation of the desired form (1.1) with a truncated nonlinearity and the block-diagonal matrix

$$A = \mathrm{diag}(-2\pi l^2 \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}, |l| \le K_x).$$

Beginning with the initial condition $u_0(x) = \exp(-3x^4 + x^2)$ on $\mathbb{T}$ that decreases fast enough, we apply Methods A and B in the two cases $\sigma = 2$ and $\sigma = 4$ with $K_t = 2^6$ modes, $N = 10$ revolutions, $m = 150$ iterations and a tolerance of $10^{-13}$ for the fixed point iteration. Figure 5 shows the evolution in time of one trajectory given by Method B (with a 300 points evaluation grid in space).

In Figure 6, we observe the discrete $L^2$ and $H^1$ norms behavior of one trajectory given by our two algorithms and the Euler method (3.15) (the simulated $(\alpha_k)_k$ are the same for Methods A and B). The Euler method quickly blows up in both norms. The $L^2$ norm of Method A is not conserved. In contrast, Method B preserves the $L^2$ norm according to Theorem 4.2. When $\sigma = 4$, numerical simulations hint that a blow-up in the $H^1$ norm always happens for all considered methods at a certain time that increases as $\epsilon$ goes to zero. We recall that in the optic fiber model (1.3), $t$ represents the distance along the optic fiber
<figure><img src="image_7.png" /><figcaption>Figure 5: Approximation by Method B of $|u|$ and $|\partial_x u|$ with $u$ solution of a spatial discretization with $K_x = 2^7$ modes of the nonlinear Schrödinger equation with white noise dispersion (5.1) on the torus $T = [-\pi, \pi]$ with the parameters $\epsilon = 10^{-2}$, $\sigma = 2$ (top) and $\sigma = 4$ (bottom).</figcaption></figure>

and a cubic nonlinearity ($\sigma = 2$) is typically considered [15]. For $\sigma = 2$, we do not observe any blow-up in the $H^1$ norm in Figure 6, suggesting the well-posedness of the model for all optic fiber distance. Also, the larger $\sigma$ is, the sooner the blow-up happens. These behaviors agree with the blow-up conjecture for $\epsilon = 1$ and $\sigma \ge 4$ presented in [4], and suggest that the conjecture persists in the highly-oscillatory regime $\epsilon \ll 1$.

## Acknowledgments

The authors would like to thank Georg Gottwald for helpful and stimulating discussions. This work was partially supported by the Swiss National Science Foundation, grants No. 200020_184614, No. 200021_162404 and No. 200020_178752. The computations were performed at the University of Geneva on the Baobab cluster using the Julia programming language.
<figure><img src="image_8.png" /><figcaption>Figure 6: Evolution in long time of the quantities $||U_t||_{L^2} - ||u_0||_{L^2}$ (left) and $||U_t||_{H^1} - ||u_0||_{H^1}$ (right) with $U_t$ the approximation computed with Euler method and Methods A and B for one trajectory of equation (5.1) with $\epsilon = 10^{-2}$, $\sigma = 2$ (top) and $\sigma = 4$ (bottom).</figcaption></figure>

## References

[1] A. Abdulle, G. A. Pavliotis, and G. Vilmart. Accelerated convergence to equilibrium and reduced asymptotic variance for Langevin dynamics using Stratonovich perturbations. *Comptes Rendus Mathematique*, 357(4):349-354, 2019.

[2] G. Agrawal. *Nonlinear Fiber Optics*. Electronics & Electrical. Academic Press, 2007.

[3] G. Agrawal. *Applications of Nonlinear Fiber Optics*. Number vol. 10 in Applications of nonlinear fiber optics. Elsevier, 2008.

[4] R. Belaouar, A. de Bouard, and A. Debussche. Numerical analysis of the nonlinear Schrödinger equation with white noise dispersion. *Stoch. Partial Differ. Equ. Anal. Comput.*, 3(1):103-132, 2015.

[5] M. Calvo, L. O. Jay, J. I. Montijano, and L. Rández. Approximate compositions of a near identity map by multi-revolution Runge-Kutta methods. *Numer. Math.*, 97(4):635-666, 2004.
[6] M. Calvo, J. I. Montijano, and L. Rández. On explicit multi-revolution Runge-Kutta schemes. *Adv. Comput. Math.*, 26(1-3):105-120, 2007.

[7] P. Chartier, J. Makazaga, A. Murua, and G. Vilmart. Multi-revolution composition methods for highly oscillatory differential equations. *Numer. Math.*, 128(1):167-192, 2014.

[8] D. Cohen. On the numerical discretisation of stochastic oscillators. *Math. Comput. Simulation*, 82(8):1478-1495, 2012.

[9] D. Cohen and G. Dujardin. Exponential integrators for nonlinear Schrödinger equations with white noise dispersion. *Stoch. Partial Differ. Equ. Anal. Comput.*, 5(4):592-613, 2017.

[10] D. Cohen and M. Sigg. Convergence analysis of trigonometric methods for stiff second-order stochastic differential equations. *Numer. Math.*, 121(1):1-29, 2012.

[11] A. de Bouard and A. Debussche. The nonlinear Schrödinger equation with white noise dispersion. *J. Funct. Anal.*, 259(5):1300-1321, 2010.

[12] A. Debussche and Y. Tsutsumi. 1D quintic nonlinear Schrödinger equation with white noise dispersion. *J. Math. Pures Appl. (9)*, 96(4):363-376, 2011.

[13] U. Erdoğan and G. J. Lord. A new class of exponential integrators for sdes with multiplicative noise. *IMA Journal of Numerical Analysis*, 2018.

[14] E. Faou. Linearized wave turbulence convergence results for three-wave systems. *arXiv e-prints*, page arXiv:1805.11269, May 2018.

[15] J. Garnier. Stabilization of dispersion-managed solitons in random optical fibers by strong dispersion management. *Optics communications*, 206(4-6):411-438, 2002.

[16] E. Hairer, C. Lubich, and G. Wanner. *Geometric numerical integration*, volume 31 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second edition, 2006. Structure-preserving algorithms for ordinary differential equations.

[17] N. Hofmann, T. Müller-Gronbach, and K. Ritter. Optimal approximation of stochastic differential equations by adaptive step-size control. *Math. Comp.*, 69(231):1017-1034, 2000.

[18] X.-M. Li. Homogenisation on homogeneous spaces. *J. Math. Soc. Japan*, 70(2):519-572, 2018. With an appendix by Dmitriy Rumynin.

[19] C. Lubich. *From quantum to classical molecular dynamics: reduced models and numerical analysis*. Zurich Lectures in Advanced Mathematics. European Mathematical Society (EMS), Zürich, 2008.

[20] R. Marty. On a splitting scheme for the nonlinear Schrödinger equation in a random medium. *Commun. Math. Sci.*, 4(4):679-705, 2006.

[21] B. n. Melendo and M. Palacios. A new approach to the construction of multirevolution methods and their implementation. *Appl. Numer. Math.*, 23(2):259-274, 1997.
[22] G. N. Milstein and M. V. Tretyakov. *Stochastic numerics for mathematical physics*. Scientific Computation. Springer-Verlag, Berlin, 2004.

[23] A. Murua and J. M. Sanz-Serna. Order conditions for numerical integrators obtained by composing simpler integrators. *R. Soc. Lond. Philos. Trans. Ser. A Math. Phys. Eng. Sci.*, 357(1754):1079–1100, 1999.

[24] D. Revuz and M. Yor. *Continuous martingales and Brownian motion*, volume 293 of *Grundlehren der Mathematischen Wissenschaften* [Fundamental Principles of Mathematical Sciences]. Springer-Verlag, Berlin, second edition, 1994.

[25] G. Vilmart. Weak second order multirevolution composition methods for highly oscillatory stochastic differential equations with additive or multiplicative noise. *SIAM J. Sci. Comput.*, 36(4):A1770–A1796, 2014.

# Appendix

*Proof of Lemma 3.2.*

1. First, $\varphi_{\varepsilon,t}(y)$ is the solution of

$$ \varphi_{\varepsilon,t}(y) = e^{AW(t)}y + \varepsilon e^{AW(t)} \int_{0}^{t} e^{-AW(s)} F(\varphi_{\varepsilon,s}(y))ds. $$

Using the boundedness of the continuous periodic function $\theta \to e^{\theta A}$ and Assumption 3.1, we get

$$ |\varphi_{\varepsilon,t}(y_1) - \varphi_{\varepsilon,t}(y_2)| \le |y| + L\varepsilon \int_{0}^{t} |\varphi_{\varepsilon,s}(y_1) - \varphi_{\varepsilon,s}(y_2)| ds. $$

The Gronwall lemma yields the desired bound.

2. Straightforward using previous statement.

3. Differentiating the integral formulation defining $\varphi_{\varepsilon,t}(y)$ gives

$$ \partial_y \varphi_{\varepsilon,t}(y)(h) = e^{AW(t)}h + \varepsilon e^{AW(t)} \int_{0}^{t} e^{-AW(s)} F'(\varphi_{\varepsilon,s}(y))(\partial_y \varphi_{\varepsilon,s}(y)(h))ds. $$

Then Assumption 3.1 yields

$$ |\partial_y \varphi_{\varepsilon,t}(y)(h)| \le |h| + L\varepsilon \int_{0}^{t} |\partial_y \varphi_{\varepsilon,t}(y)(h)| ds. $$

The Gronwall lemma allows us to obtain

$$ |\partial_y \varphi_{\varepsilon,t}(y)| \le e^{L\varepsilon t}. $$

For the second derivative, we get

$$ \partial_y^2 \varphi_{\varepsilon,t}(y)(h,k) = \varepsilon e^{AW(t)} \int_{0}^{t} e^{-AW(s)} [F'(\varphi_{\varepsilon,s}(y))(\partial_y^2 \varphi_{\varepsilon,s}(y)(h,k))] $$
$$ + F''(\varphi_{\varepsilon,s}(y))(\partial_y\varphi_{\varepsilon,s}(y)(h), \partial_y\varphi_{\varepsilon,s}(y)(k))]ds. $$

Then

$$ \begin{aligned} |\partial_y^2 \varphi_{\varepsilon,t}(y)(h,k)| &\le C\varepsilon \int_0^t [(1 + |\varphi_{\varepsilon,s}(y)|^K) |\partial_y \varphi_{\varepsilon,s}(y)(h)| |\partial_y \varphi_{\varepsilon,s}(y)(k)| \\ &\quad + |\partial_y^2 \varphi_{\varepsilon,s}(y)(h,k)|]ds \\ &\le C\varepsilon t(1 + |y|^K) e^{C\varepsilon t} |h| |k| + C\varepsilon \int_0^t |\partial_y^2 \varphi_{\varepsilon,s}(y)(h,k)|. \end{aligned} $$

Then the Gronwall lemma allows to conclude. The proof is similar for the third derivative.

□