Multirevolution integrators for differential equations with fast stochastic oscillations

TL;DR

This paper presents a novel multirevolution method for efficiently solving stochastic differential equations with fast oscillations driven by noise, achieving high accuracy and invariance preservation.

Contribution

It introduces a new multirevolution integrator for stochastic oscillatory equations, with weak order two and invariant-preserving modifications, independent of oscillation stiffness.

Findings

Achieves weak order two accuracy

Computational cost independent of oscillation stiffness

Exact preservation of quadratic invariants

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\"odinger 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.

Equations

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

$$dX=\frac{2\pi }{\sqrt{\epsilon }}\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)X\circ dW+\left(\begin{array}{cc}0& -a\\ a& 0\end{array}\right)Xdt,a\in \mathbb{R},$$

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

$$\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)),v(x=0,t)=u_0(t).$$

$$\frac{dx}{dt}=\frac{1}{\epsilon }Ax+F(x),x(0)=x_0.$$

$$T_0$$

$$T_{N+1}$$

$$\frac{d{y}_t}{dt}=⟨g^0⟩(y_t),y_0=X_0,$$

$$g_{\theta }^0(y)$$

$$g_{\theta }^1(y)(z)$$

$$\begin{array}[]{rl}\mathbb{E}[\widehat{\alpha}_{k}^{N}]&=\left\{\begin{array}[]{l}1\text{ if }k=0\\ 0\text{ else}\end{array}\right.\\ \widehat{\beta}_{p,k}^{N}&=\left\{\begin{array}[]{l}\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{array}\right.\end{array}\hskip-14.22636pt\begin{array}[]{rl}\mathbb{E}[\widehat{\alpha}_{p}^{N}\widehat{\alpha}_{k}^{N}]&=\left\{\begin{array}[]{l}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{array}\right.\\ \widehat{\widetilde{\beta}}_{p,k}^{N}&=\left\{\begin{array}[]{l}\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{array}\right.\\ \end{array}$$

$$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$$

$$Y_{m+1}$$

$$+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{\tilde{\beta }}}_{p,k}^N$$

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

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

$${\phi }_{\epsilon ,t}(y)=e^{AW(t)}y+\epsilon {\int }_0^t{e}^{A(W(t)-W(s))}F({\phi }_{\epsilon ,s}(y))ds.$$

$${\psi }_{\epsilon ,t}^1(y)$$

$${\psi }_{\epsilon ,t}^2(y)$$

$$\cdot F^{\prime }(e^{AW(s)}y)\left(e^{AW(s)}{\int }_0^s{e}^{-AW(r)}F(e^{AW(r)}y)dr\right)ds.$$

$$\left|{\phi }_{\epsilon ,t}(y)-{\psi }_{\epsilon ,t}^j(y)\right|\le C(1+{\left|y\right|}^K)e^{C\epsilon t}(\epsilon t{)}^{j+1}.$$

$$\left|{\psi }_{\epsilon ,t}^1(y)\right|$$

$$\left|{\psi }_{\epsilon ,t}^2(y)\right|$$

$$\left|{\psi }_{\epsilon ,t}^2(y)-e^{AW(t)}y\right|$$

$$\left|{\phi }_{\epsilon ,t}(y)-{\psi }_{\epsilon ,t}^1(y)\right|\le L\epsilon {\int }_0^t\left|{\phi }_{\epsilon ,s}(y)-e^{AW(s)}y\right|ds.$$

$$\left|{\phi }_{\epsilon ,s}(y)-e^{AW(s)}y\right|$$

$$\le C\epsilon {\int }_0^s(1+C(1+\left|y\right|)e^{C\epsilon r})dr\le C(1+\left|y\right|)e^{C\epsilon s}(\epsilon s).$$

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

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

$$F({\psi }_{\epsilon ,s}^1(y))=F(e^{AW(s)}y)+\epsilon F^{\prime }(e^{AW(s)}y)\left(e^{AW(s)}{\int }_0^s{e}^{-AW(r)}F(e^{AW(r)}y)dr\right)+R_{\epsilon ,s}.$$

$$\left|R_{\epsilon ,s}\right|\le C{\epsilon }^2\underset{x\in [e^{AW(s)}y,{\psi }_{\epsilon ,s}^1(y)]}{\sup }\Vert F^{\prime \prime }(x)\Vert {\left|e^{AW(s)}{\int }_0^s{e}^{-AW(r)}F(e^{AW(r)}y)dr\right|}^2.$$

$$\left|R_{\epsilon ,s}\right|\le C(1+{\left|e^{AW(s)}y\right|}^K+{\left|{\psi }_{\epsilon ,s}^1(y)\right|}^K)(\epsilon s{)}^2{e}^{C\epsilon s}\le C(1+{\left|y\right|}^K)(\epsilon s{)}^2{e}^{C\epsilon s}.$$

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

$$\mathbb{E}[T_1^k]=\frac{(-2{)}^{k}k!}{(2k)!}\sum _{j=1}^p(-1{)}^j\sum _{\underset{n_i\in {\mathbb{N}}^{\ast }}{n_1+\cdots +n_j=p}}(\genfrac{}{}{0ex}{}{2p}{2n_1,\dots ,2n_j}).$$

Full text

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

[TABLE]

where $W$ is a standard one-dimensional Wiener process, the function $F:\mathbb{R}^{d}\rightarrow\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}=\operatorname{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,\left|\varepsilon^{-1/2}W(t)\right|=1\}$ of mean $\varepsilon$. The class of SDEs (1.1) includes in particular highly-oscillatory Kubo oscillators (see [8])

[TABLE]

or equivalently, $dY=2i\pi\varepsilon^{-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\leq\varepsilon$ to be accurate, which makes these methods dramatically expensive when $\varepsilon$ 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 $\varepsilon$ 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 discretizations of the highly-oscillatory nonlinear Schrödinger equation (NLS) with white noise dispersion

[TABLE]

As described for instance in [15], in the case $\varepsilon=1$, the NLS equation (1.3) with a cubic nonlinearity $F(u)=|u|^{2}u$ 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,

[TABLE]

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 [math] and under some ergodicity assumptions on $m$, $u^{\nu}$ converges to the solution $u$ of equation (1.3) with $\varepsilon=1$. The non-stiff counterpart of equation (1.3), i.e. for $\varepsilon=1$, has also been studied theorically in [12] for a particular nonlinearity. The highly-oscillatory behaviour ($\varepsilon\ll 1$) appears naturally when observing the propagation in long time with a small nonlinearity (via the change of variable $t\leftarrow\varepsilon 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/\varepsilon$). 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 litterature 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 Stratnovitch 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 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 $\left\|u(t)\right\|_{L^{2}}=\left\|u_{0}\right\|_{L^{2}}$ for all $t\geq 0$) for a class of polynomial nonlinearities. However they face a severe timestep restriction $h\leq\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,

[TABLE]

The solution $x$ of this equation at times $\varepsilon NT$ is a perturbation of identity, that is $x$ satifies $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 variable called revolution times and defined by

[TABLE]

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\rightarrow 0$ to the solution $y_{t}$ of the deterministic ODE

[TABLE]

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,\ldots$ with cost in $H=N\varepsilon=\mathcal{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}Sy$ 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$,

[TABLE]

with respective Fourier coefficients $(c_{k}^{0}(y))_{k\in\mathbb{Z}}$ and $(c_{p}^{1}(y))_{p\in\mathbb{Z}}$. 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 $(\widehat{\alpha}_{k}^{N})_{k}$, and deterministic sequences $(\widehat{\beta}_{p,k}^{N})_{p,k}$ and $(\widehat{\widetilde{\beta}}_{p,k}^{N})_{p,k}$ that satisfy

[TABLE]

The definition and construction of these random variables is further discussed in Section 3.2 and Section 3.4.

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 $\mathcal{O}(H)$, the convergence rate of the fixed point iterations is independent of the smallness of $\varepsilon$. *

Remark 2.2**.**

We observe that $\widehat{\beta}_{p,k}^{N}$ and $\widehat{\widetilde{\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 $\varepsilon\rightarrow 0$.

3.1 Asymptotic expansion of the exact solution

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

[TABLE]

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

Assumption 3.1**.**

The function $F$ is globally Lipschitz continuous and lies in $\mathcal{C}^{3}_{P}$, i.e. there exists constants $L$, $C$, $K>0$ such that for all $y$, $y_{1}$, $y_{2}\in\mathbb{R}^{d}$

[TABLE]

Also the initial condition $X_{0}$ has bounded moments, that is $\mathbb{E}[\left|X_{0}\right|^{p}]<\infty$ for $p\geq 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=\mathcal{O}(\varepsilon^{-1})$. The variation of constants formula yields

[TABLE]

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$,

$\left|\varphi_{\varepsilon,t}(y_{1})-\varphi_{\varepsilon,t}(y_{2})\right|\leq C\left|y_{1}-y_{2}\right|e^{C\varepsilon t}$, 2. 2.

$\left|\varphi_{\varepsilon,t}(y)\right|\leq C(1+\left|y\right|)e^{C\varepsilon t}$, 3. 3.

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

The proof is postponed to the appendices. 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)$,

[TABLE]

Proposition 3.3** (Local expansion).**

Under Assumption 3.1, for all $y\in\mathbb{R}^{d}$, $j\in\{1,2\}$ and $t\geq 0$, there exists $C$ and $K$ two positive constants independent of $\varepsilon$ and $t$ such that

[TABLE]

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$,

[TABLE]

Proof of Proposition 3.3.

Using Assumption 3.1, we get

[TABLE]

Then Lemma 3.2 yields

[TABLE]

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

For $j=2$, we first denote

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Then, using the polynomial growth of $F^{\prime\prime}$ and inequality (3.5), we get

[TABLE]

Hence the result. $\square$

We shall prove in Section 3.3 that the function $\psi^{2}_{\varepsilon,t}$ 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

[TABLE]

where $W_{0}=0$ and $W_{i+1}=W{i}+\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 $\varepsilon 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 $\mathcal{O}(t^{2}/\tau^{2})$, which makes this method tremendously expensive for $t=\mathcal{O}(\varepsilon^{-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 replace stochastic integrals by appropriate discrete random variables in order not to simulate any expensive Brownian path $W$.

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\geq 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 $\operatorname{Re}(z)<\frac{\pi^{2}}{8}$. In addition, for $x\in\left[0,\frac{\pi^{2}}{8}\right[$, $\mathbb{E}[e^{xT_{1}}]=\frac{1}{\cos(\sqrt{2x})}$. The variable $T_{1}$ has bounded moments and they are given by

[TABLE]

In particular, $\mathbb{E}[T_{1}]=1$, $\mathbb{E}[T_{1}^{2}]=\frac{5}{3}$ and $\operatorname{Var}(T_{1})=\frac{2}{3}$. Finally, for a fixed $\varepsilon_{0}\in\left]0,\frac{\pi^{2}}{16}\right[$, for all $\varepsilon\in]0,\varepsilon_{0}]$ and $p\geq 0$, we have the estimate

[TABLE]

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 expansion of $\mathbb{E}[e^{xT_{1}}]$ and $\frac{1}{\cos(\sqrt{2x})}$ yields (3.9). The estimate (3.10) is proved as follows

[TABLE]

where we used first the Cauchy-Schwarz inequality and then twice the Jensen inequality. $\square$

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 following random variables

[TABLE]

satisfy

[TABLE]

Proof.

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

[TABLE]

which yields at time $t=T_{N}$,

[TABLE]

Then $t\mapsto\int_{0}^{t}e^{2i\pi kW(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 kW(s)}dW(s)]=0$ for all $t$. The dominated convergence theorem for stochastic integrals allows to take the limit $t\rightarrow\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 pW(s)}$ and we integrate from $s$ to $T_{N}$,

[TABLE]

Then, multiplying by $\frac{1}{N^{2}}e^{2i\pi kW(s)}$ and integrating from [math] to $T_{N}$ yields

[TABLE]

Using the stochastic Fubini theorem, we deduce

[TABLE]

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 kW(s)}\int_{s}^{T_{N}}e^{2i\pi pW(r)}drds$, so that we get if $p\neq 0$,

[TABLE]

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

[TABLE]

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}[\widetilde{\beta}_{p,k}^{N}]$ from the formula $\widetilde{\beta}_{p,k}^{N}=\beta_{p,k}^{N}-\frac{\alpha_{p}^{N}\alpha_{k}^{N}}{2}$. $\square$

Remark** (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,

[TABLE]

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^{\prime}$ 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)$,

[TABLE]

Notice that $c_{k}^{0}(y)\in\mathbb{C}^{d}$ and $c_{p}^{1}(y)=(c_{k}^{0})^{\prime}(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.8**.**

We define the following quantity

[TABLE]

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 function $\phi\in\mathcal{C}^{3}_{P}$, there exists $H_{0}>0$ such that for all $H=N\varepsilon\leq H_{0}$, the following estimates hold, where $C$ and $K$ are independent of $\varepsilon$ and $N$,

[TABLE]

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

[TABLE]

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\rightarrow 0$ (or equivalently $N\rightarrow\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^{0}_{\theta}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.9** (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\rightarrow 0$ to the solution at time $T$ of

[TABLE]

that is, for all test function $\phi\in\mathcal{C}^{3}_{P}$,

[TABLE]

Remark 3.10**.**

It can be proven using the results of Section 4 that the solution of the asymptotic model (Proposition 3.9) is an order one weak approximation of $X(\varepsilon T_{Nm})$ for $m\geq 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),

[TABLE]

Its cost is independent of $\varepsilon$ and $N$, and it has weak order one w.r.t. $H$, that is for all $m\geq 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.8. 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 $\widehat{\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

[TABLE]

and $\Gamma^{N}$ its square root. Then, $\widehat{\alpha}_{k}^{N}$ is defined for $k\geq 0$ as

[TABLE]

and we fix $\widehat{\alpha}_{k}^{N}=\overline{\widehat{\alpha}_{-k}^{N}}$ for $k<0$ (so that the solution stays real while still having the good moments). We also define $\widehat{\beta}_{p,k}^{N}=\mathbb{E}[\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 $\widehat{\beta}_{p,k}^{N}$ by $\widehat{\widetilde{\beta}}_{p,k}^{N}=\mathbb{E}[\widetilde{\beta}_{p,k}^{N}]$, using the values of Proposition 3.7.

Remark 3.11**.**

The methodology presented in Section 3 can be generalised 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

[TABLE]

with the new random variables

[TABLE]

and where the discrete random variables share the same moments up to order 3 for the $\widehat{\alpha}_{k}^{N}$, order 2 for the $\widehat{\beta}_{p,k}^{N}$, and order 1 for the $\widehat{\gamma}_{l,p,k}^{(i),N}$. It is also possible to generalise 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\leq 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 function $\phi\in\mathcal{C}^{3}_{P}$, there exists $H_{0}>0$ such that for all $H\leq H_{0}$, for all $m\geq 0$ such that $mN\varepsilon=mH\leq T$, there exists two positive constants $K$ and $C$ both independent of $\varepsilon$, $N$ and $K_{t}$ such that

[TABLE]

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\leq 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 $\varepsilon 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 Theorem 4.1 and Theorem 4.2 remain valid for infinitely many modes ($K_{t}\rightarrow\infty$). In addition, assuming that $F$ is of class $C^{s+1}_{P}$ yields a truncation error of the Fourier series in (2.3) of size $\mathcal{O}((1+\left|y\right|^{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+\left|y\right|)e^{-cK_{t}})$. For simplicity of the analysis, we thus assumed in Theorem 4.1 and Theorem 4.2 that $g_{\theta}^{0}$ and $g_{\theta}^{1}$ have a finite number $K_{t}$ on non-zero Fourier modes in (2.3). If this assumption does not hold, the truncation errors $\mathcal{O}((1+\left|y\right|^{K})K_{t}^{-s})$ or $\mathcal{O}((1+\left|y\right|)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 $\widetilde{\varphi}_{\varepsilon,t}(y)=e^{-AW(t)}\varphi_{\varepsilon,t}(y)$ that has no effect at time $t=T_{Nm}$. We now have to compare the two solutions of the following integral formulations

[TABLE]

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

[TABLE]

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

[TABLE]

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 $\widehat{\psi}_{\varepsilon,N}(y)$ is given by

[TABLE]

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

Proof.

We first prove

[TABLE]

where $\mathbb{E}[(M_{N})^{2q}]\leq C_{q}$ for all $q>0$. We have

[TABLE]

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

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

[TABLE]

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

[TABLE]

We deduce

[TABLE]

and by induction $\mathbb{E}\left[\left|\widehat{\psi}_{\varepsilon,N}^{m}(y)\right|^{2q}\right]\leq e^{C_{q}mH}(1+\mathbb{E}[\left|y\right|^{2q}])\leq e^{C_{q}T}(1+\mathbb{E}[\left|y\right|^{2q}])$. $\square$

Proposition 4.5** (Bounded moments for the integrator (2.5)).**

Assume that for $y\in\mathbb{R}^{d}$, the numerical scheme $\widehat{\psi}_{\varepsilon,N}(y)$ satisfies

[TABLE]

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

Proof.

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

[TABLE]

As $\sum_{k}\left|\widehat{\alpha}_{k}^{N}\right|+\sum_{p,k}\left|\widehat{\widetilde{\beta}}_{p,k}^{N}\right|$ is bounded, using the same estimates as in the proof of Proposition 4.4, we get for all $H\leq H_{0}$ small enough,

[TABLE]

where $M_{N}$ has bounded moments. The remaining of the proof is the same as in the proof of Proposition 4.4. $\square$

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

Proof.

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

[TABLE]

A local expansion gives

[TABLE]

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 yield

[TABLE]

We obtain a similar expansion for $\phi(\widehat{\psi}_{\varepsilon,N}(y))$:

[TABLE]

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

[TABLE]

Making the difference of both equations gives

[TABLE]

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

[TABLE]

Then, we get

[TABLE]

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

[TABLE]

Let us now study the second order term $Z=\phi^{\prime\prime}(y)(\widehat{\psi}_{\varepsilon,N}(y)-y)^{2}-\phi^{\prime\prime}(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^{2}Y+R$ where $\mathbb{E}[\left|R\right|]\leq C(1+\left|y\right|^{K})H^{3}$ (by the same arguments as before) and

[TABLE]

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

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

[TABLE]

We deduce the local order two of the proposed numerical scheme. $\square$

Remark**.**

The constant $H_{0}$ in Proposition 4.6 depends on $F$, but also of 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

[TABLE]

4.3 Global error

Theorem 4.7** (Global convergence).**

Assume that the numerical scheme $\widehat{\psi}_{\varepsilon,N}$ satisfies equation (4.2) (respectively equation (4.4)) where $\widehat{\alpha}_{k}^{N}\in\mathbb{C}$, $\widehat{\beta}_{p,k}^{N}\in\mathbb{R}$ (respectively $\widehat{\widetilde{\beta}}_{p,k}^{N}\in\mathbb{R}$) are random variables such that $\widehat{\alpha}_{k}^{N}=\overline{\widehat{\alpha}_{-k}^{N}}$ and

[TABLE]

(respectively $\widehat{\alpha}_{k}^{N}$ satisfies the same conditions and $\widehat{\widetilde{\beta}}_{p,k}^{N}$ satisfies $\mathbb{E}[\widehat{\widetilde{\beta}}_{p,k}^{N}]=\mathbb{E}[\widetilde{\beta}_{p,k}^{N}]$). Under Assumption 3.1, if for all $q>0$, $\mathbb{E}\left[\left(\sum_{k}\frac{\left|\widehat{\alpha}_{k}^{N}\right|^{2}}{k^{2}}\right)^{q}\right]$ and $\mathbb{E}\left[\left(\sum_{p,k}\frac{\left|\widehat{\beta}_{p,k}^{N}\right|^{2}}{k^{2}}\right)^{q}\right]$ are bounded uniformly in $N$ (respectively $\sum_{k}\left|\widehat{\alpha}_{k}^{N}\right|$, $\sum_{p,k}\left|\widehat{\widetilde{\beta}}_{p,k}^{N}\right|$, $\mathbb{E}\left[\left(\sum_{k}\frac{\left|\widehat{\alpha}_{k}^{N}\right|^{2}}{k^{2}}\right)^{q}\right]$ and \mathbb{E}\bigg{[}\bigg{(}\sum_{p,k}\frac{\left|\widehat{\widetilde{\beta}}_{p,k}^{N}\right|^{2}}{k^{2}}\bigg{)}^{q}\bigg{]} are bounded uniformly in $N$), for all $T>0$, for all test function $\phi\in\mathcal{C}^{3}_{P}$, there exists $H_{0}>0$ such that for all $H\leq H_{0}$, for all $M\geq 0$ such that $MN\varepsilon=MH\leq T$, there exists two positive constants $K$ and $C$ both independent of $\varepsilon$ and $N$ such that

[TABLE]

Proof.

We denote

[TABLE]

and rewrite it with a telescopic sum

[TABLE]

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

[TABLE]

Thus, knowing the hitting times involved, $\widetilde{\phi}_{m}\in\mathcal{C}^{3}_{P}$. Using Assumption 3.1, $(c_{p}^{0})^{\prime}=c_{p}^{1}$ and $\beta_{p,k}^{N}=\widetilde{\beta}_{p,k}^{N}+\frac{\alpha_{p}\alpha_{k}}{2}$, we deduce that $\widehat{\psi}_{\varepsilon,N}$ satisfies the assumptions of Proposition 4.6. Applying Proposition 4.6 to each term of $e_{M}$ gives

[TABLE]

Finally, the moments of $\widehat{\psi}_{\varepsilon,N}^{m}(X_{0})$ are all bounded uniformly in $\varepsilon$, $N$ and $m$ according to Proposition 4.4 (respectively 4.5). Thus

[TABLE]

We deduce the global weak order two. $\square$

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

Proof of Proposition 3.9.

Rewriting Theorem 4.7 for order one yields for all $H=N\varepsilon$ small enough and all $M\geq 0$,

[TABLE]

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

Proof of Theorem 4.1.

As $\widehat{\alpha}_{k}^{N}\leq C$ and $\sum_{p,k}\frac{\left|\mathbb{E}[\beta_{p,k}^{N}]\right|^{2}}{k^{2}}$ converges by Proposition 3.7, Theorem 4.7 applies and concludes the proof. $\square$

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}\left|\widehat{\alpha}_{k}^{N}\right|$ and $\sum_{p,k}\left|\widehat{\widetilde{\beta}}_{p,k}^{N}\right|$ are bounded, the right hand-side of equation (2.5) is a contraction for all $H\leq H_{0}$ small enough and the constant does not depend on $Y_{m}$, so $H_{0}$ depends only of $F$ and $F^{\prime}$. Thus, the integrator is well-posed for all $H\leq H_{0}$.

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

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

[TABLE]

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

[TABLE]

For the second order term, equation (4.8) and the values of Proposition 3.7 yield

[TABLE]

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

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 innacurate 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

[TABLE]

We plot on a logarithmic scale an estimate of the weak error for approximating $X$ at time $T=10^{-3}T_{2^{8}}$ where $\mathbb{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+\operatorname{Re}(y)^{3}+\operatorname{Im}(y)^{5})y$. In the real setting, it yields the following two-dimensional SDE

[TABLE]

We take 8 modes for the Fourier decomposition and the same other parameters as before. The average is taken over $10^{6}$ trajectories.

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 $\varepsilon=10^{-2}$,

[TABLE]

where the unknown $u$ is a random process depending on $x\in\mathbb{T}$ and $t\geq 0$. We consider a spectral discretization in space of this equation with $K_{x}=2^{7}$ modes $u(x,t)\approx\sum_{\left|l\right|\leq 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

[TABLE]

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 behaviour 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 $\varepsilon$ goes to zero. We recall that in the optic fiber model (1.3), $t$ represents the distance along the optic fiber 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 $\varepsilon=1$ and $\sigma\geq 4$ presented in [4], and suggest that the conjecture persists in the highly-oscillatory regime $\varepsilon\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.

Appendix

Proof of Lemma 3.2.

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

[TABLE]

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

[TABLE]

The Gronwall lemma yields the desired bound. 2. 2.

Straightforward using previous statement. 3. 3.

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

[TABLE]

Then Assumption 3.1 yields

[TABLE]

The Gronwall lemma allows to obtain

[TABLE]

For the second derivative, we get

[TABLE]

Then

[TABLE]

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

$\square$