¹

1 Introduction↩︎

In this article, we study the temporal approximation of semilinear stochastic evolution equations of the form \[\label{eq:introSEE} \begin{cases} \mathrm{d}U + AU \,\mathrm{d}t&= F(U)\,\mathrm{d}t+ G(U) \,\mathrm{d}W\quad \text{on }(0,T], \\ U_0&=\xi\in L^p(\Omega;X). \end{cases}\tag{1}\] Here, \(-A\) is the generator of a contractive but in general not analytic \(C_0\)-semigroup on a Hilbert space \(X\), the nonlinearity \(F\) and the multiplicative noise \(G\) are globally Lipschitz, \(W\) is a cylindrical Brownian motion, \(\xi\) is the initial data, and \(p\in [2,\infty)\). In general, \(F\) and \(G\) can be time-dependent and random but for simplicity we restrict this introductory discussion to autonomous, deterministic coefficients.

The aim of this paper is to obtain strong convergence rates exceeding \(1/2\) for a class of time discretization schemes, the so-called Milstein schemes, in the hyperbolic setting. We present regularity assumptions on \(F\), \(G\), and \(\xi\) under which the optimal rate \(1\) is attained, such as linear growth of \(F\) and \(G\) on \(D(A)\) for the classical Milstein scheme or on \(D(A^2)\) for rational versions thereof.

Hyperbolic SPDEs such as Schrödinger or Maxwell’s equations have attracted considerable attention in recent years (see [1]–[10] and references therein). Strong convergence rates have been studied in this setting for the exponential Euler method for some non-parabolic equations [1], [4], which was generalised to rational schemes with a unifying approach by one of the authors in [11]. However, it is known that the rate \(1/2\) up to a logarithmic correction factor is optimal among all schemes solely using the Wiener increments as information from the Brownian motion. Indeed, already for SDEs, the error of such schemes has to grow at least like \(\sqrt{\log(T/h)}h^{1/2}\) for step size \(h\to 0\) [12].

The Milstein scheme for SPDEs↩︎

The Milstein scheme overcomes this limitation by incorporating more information from the Brownian motion. Originally developed for SDEs [13] (see also [14]) via an Itô–Taylor expansion of first order, Jentzen and Röckner extended it to SPDEs in their seminal paper [15].

The central idea behind the Milstein scheme is as follows. Assume for simplicity that the mild solution of 1 is even a strong one, i.e.it can be written as an Itô process \(U_{h}=\xi+\int_0^h (-AU_{s}+ F(U_{s}))\,\mathrm{d}s+ \int_0^h G(U_{s})\,\mathrm{d}W_s\), where \(h>0\) is the time step. To construct a first-order approximation, it suffices to approximate the temporal integrand \(-AU_{s}+F(U_{s})\approx -AU_{0}+F(U_{0})\) to order \(0\). In contrast, since the stochastic integral is only of order \(1/2\), the stochastic integrand should be approximated to order \(1/2\) rather than \(0\). From Itô’s formula, we deduce \(G(U_{s}) \approx G(U_{0})+\int_0^s G'(U_{r})\,\mathrm{d}U_{r}\) and thus, using 1 and omitting all terms of order greater than \(1/2\), \[\begin{align} G(U_{s}) &\approx G(U_{0})+\int_0^s G'(U_{r})\p[\big]{-AU_{r}+F(U_{r})}\,\mathrm{d}r+\int_0^s G'(U_{r})G(U_{r})\,\mathrm{d}W_r\\ &\approx G(U_{0})+\int_0^s G'(U_{0})G(U_{0})\,\mathrm{d}W_r. \end{align}\] Inserting these approximations in the solution formula results in the expansion \[\begin{align} U_{h} \approx \xi-hAU_{0}+h F(U_{0}) + \int_0^h G(U_{0}) \,\mathrm{d}W_s+ \int_0^h \int_0^s G'(U_{0})G(U_{0})\,\mathrm{d}W_r\,\mathrm{d}W_s \end{align}\] with an iterated stochastic integral characteristic of Milstein schemes. An analogous reasoning applies to mild (rather than strong) solutions via mild Itô formulas, cf.[16], where convolutions with the semigroup \((S(t))_{t\ge 0}\) arise. Motivated by this observation, we define the rational Milstein scheme \(u=(u_{j})_{j=0,\ldots,M}\) by \(u_{0}\mathrel{\vcenter{:}}=\xi\) and for \(1 \le j \le M\) \[\begin{align} u_{j} &\mathrel{\vcenter{:}}= R_h^j \xi+ h\sum_{i=0}^{j-1} R_h^{j-i}F(u_{i}) + \sum_{i=0}^{j-1} R_h^{j-i}\p[\Big]{\int_{t_{i}}^{t_{i+1}}G(u_{i})\,\mathrm{d}W_s+ \int_{t_{i}}^{t_{i+1}} G'(u_{i})\br[\Big]{\int_{t_{i}}^s G(u_{i}) \,\mathrm{d}W_r}\,\mathrm{d}W_s}. \end{align}\] Here, \(t_{j}=jh\), \(0\le j \le M\), and \(R=(R_h)_{h>0}\) is a time discretization scheme approximating the semigroup \(S\), i.e.\(R_h \approx S(h)\). If \(R=S\), it is called the (exponential) Milstein scheme.

The Milstein scheme has been studied extensively for parabolic SPDEs. Convergence of the root-mean-square error at rate \(1\) was shown in [15] under a commutativity condition on the noise, which has subsequently been lifted in [17]. Further developments include its analysis for SPDEs driven by non-continuous martingale noise [18], its \(L^p\)- and almost sure convergence for advection-diffusion equations [19], and the space-time-discretisation by a Milstein–Galerkin scheme [20]. Rational Milstein schemes have been formulated as mild Itô processes in [16] and a software package for the numerical approximation of the iterated stochastic integrals is available [21]. Finally, derivative-free Milstein-type schemes have been proposed [22], [23], offering an easier implementation while preserving high convergence orders.

However, all of these higher-order convergence rates pertain to the parabolic case, where regularisation phenomena of the underlying analytic semigroup can be used to improve convergence rates. These are not applicable to 1 if the associated semigroup is merely contractive. Establishing higher-order convergence of the Milstein scheme for hyperbolic stochastic evolution equations has already been raised as an open research direction over a decade ago in [15], together with the generalisation from exponential to rational Milstein schemes. Numerical simulations suggest convergence at rate \(1\) of the Milstein scheme for stochastic Schrödinger equations [24]. Recently, convergence rates of up to \(3/2\) have been established for the stochastic wave equation with finite-dimensional noise using a so-called \((\hat{\alpha},\beta)\)-scheme (with \((\hat{\alpha},\beta)=(1,0)\)) [25]. Further recent developments for hyperbolic SPDEs include [3], where rate \(1\) was shown for a splitting scheme for the nonlinear Schrödinger equation with additive noise and in [8] with multiplicative noise.

Often in the literature, the error considered is the pointwise strong error \[\begin{align} \max_{0 \le j\le M}\mathbb{E}\br[\big]{\|U_{t_{j}}-u_{j}\|_X^p}, \end{align}\] where \(U\) is the mild solution to 1 and \(u\) its temporal discretisation. However, smallness of the pointwise strong error does not imply convergence of the path of the approximations or that numerical simulations converge to the mild solution. Hence, we investigate convergence rates of the pathwise uniform strong error \[\begin{align} \mathbb{E}\br[\Big]{\max_{0 \le j\le M} \|U_{t_{j}} - u_{j}\|_X^p} \end{align}\] instead, where now the maximum over \(j\) is inside the expectation. In general, the pathwise uniform error cannot be estimated by means of the pointwise strong error without deteriorating the convergence rate in the non-deterministic case. If the pathwise uniform strong error decays at rate \(\alpha>0\) for all \(p\in [2,\infty)\), a Borel–Cantelli argument yields almost sure convergence at rate \(\alpha-\varepsilon\) for all \(\varepsilon\in (0,\alpha)\). That is, asymptotically, one has \(\max_{0\le j \le M}\|U_{t_{j}} - u_{j}\|_X\lesssim h^{\alpha-\varepsilon}\) almost surely.

Main result↩︎

In the spirit of the Kato setting [26], let \(X\) and \(Y\) be Hilbert spaces such that \(Y\hookrightarrow X\) continuously and let \((S(t))_{t\ge 0}\) be a \(C_0\)-semigroup on \(X\). A time discretization scheme \(R\colon [0,\infty) \to \mathscr{L}(X),~h\mapsto R_h \mathrel{\vcenter{:}}= R(h)\) approximates \(S\) to order \(\alpha\in(0,1]\) on \(Y\) if for all \(T>0\) \[\norm{(S(t_{j})-R_h^j)u}_X \lesssim_{T,\alpha} h^\alpha \norm{u}_Y\] for all \(u \in Y\), \(h>0\), and \(j \in \mathbb{N}\) such that \(t_{j}=jh \in [0,T]\), where \(R_h^j=(R_h)^j\). It is called contractive on \(X\) if \(\norm{R_h}_{\mathscr{L}(X)}\le 1\) for all \(h>0\). Let \(H\) be a Hilbert space and denote by \(\mathscr{L}_2(H,X)\) and \(\mathscr{L}_2^{(2)}(H,X)\) the spaces of linear and bilinear Hilbert–Schmidt operators, respectively. Our main result establishing convergence rates for the pathwise uniform error of the rational Milstein scheme is as follows.

Theorem 1 is obtained as a special case of the main results in Theorem 15 and Theorem 16 for the exponential Milstein scheme. They also allow for \(F\) and \(G\) to be time-dependent and random. If the assumptions are satisfied for \(p=2\), the same bound is obtained for the \(r\)-th moment, \(r\in [1,2]\). For a suitable continuous-time extension \((\bar{u}_t)_{t\in [0,T]}\) of the rational Milstein scheme, an analogous estimate to ?? holds for the supremum of \(U_{t}-\bar{u}_t\) over the full time interval \([0,T]\), cf. Theorem 20. Common choices for \(Y\) are suitable intermediate spaces between \(X\) and \(\mathop{\mathrm{D}}(A^2)\) such as domains of fractional powers of \(A\). In the special case that \(Y = \mathop{\mathrm{D}}(A)\), the optimal rate of convergence \(1\) is achieved for the exponential Milstein scheme, and, if one can choose \(Y=\mathop{\mathrm{D}}(A^2)\), also for most commonly used rational Milstein schemes. Table 1 illustrates this dependence of the convergence rate on the choice of \(Y\) in more detail.

The implicit Euler and the Crank–Nicolson scheme are two common choices of rational schemes to approximate the semigroup, but Theorem 1 is not limited to these. From functional calculus, it follows that any scheme \(R_h=r(-hA)\) is admissible if it is induced by a rational, holomorphic function \(r\colon\mathbb{C}_-\to\mathbb{C}\) satisfying \(|r|\le 1\) on the left open half-plane \(\mathbb{C}_-\) and approximating the exponential function (cf.Proposition 6). In particular, this includes A-acceptable and consistent rational A-stable schemes such as the higher-order implicit Runge–Kutta methods Lobatto IIIA, IIIB, and IIIC, Radau methods, and certain DIRK schemes (see [27] for definitions and properties of these schemes).

The error estimate ?? is optimal in the sense that it matches the convergence rate of the initial-value term on its own, up to a logarithmic factor for rational Milstein schemes. Moreover, in the case of sufficient regularity, we achieve the optimal rate \(1\), which coincides with the order of the Itô–Taylor expansion in the SDE case.

To the best of the authors’ knowledge, the present work provides the first rigorous error analysis of the Milstein scheme for hyperbolic SPDEs, both in terms of an abstract framework and for concrete model equations listed below. Our primary contributions are:

• first optimal convergence rates for Milstein schemes for hyperbolic SPDEs, partially addressing an open problem raised in [15] and [11],

• treatment of rational Milstein schemes based on rational semigroup approximations \(R_h\neq S(h)\) such as the Crank–Nicolson Milstein scheme, answering the open problem posed for parabolic SPDEs in [15] in the case of hyperbolic SPDEs,

• maximal estimates in \(p\)-th moment, \(p\in [2,\infty)\), leading to pathwise uniform convergence rates rather than pointwise root-mean-square estimates,

• error estimates on the full time interval for a suitable continuous-time extension of Milstein-type schemes for hyperbolic SPDEs.

For concrete equations like Schrödinger or Maxwell’s equations, our results improve results from the literature for rational schemes to higher rates \(\alpha\in (1/2,1]\) for Milstein schemes. Numerical simulations confirm these rates.

While our results also apply to additive noise, they are not novel in this case. This can be attributed to the observation that the Milstein scheme then reduces to the exponential Euler method or the corresponding rational scheme, for which convergence at rate \(1\) was shown in [11]. Likewise, we do not achieve improved rates for wave equations compared to [11], but our result covers nonlinearities and noise depending on both position and velocity components of the solution as opposed to just the position. In the latter case, rate \(3/2\) was achieved for finite-dimensional noise in [25] for a so-called \((\hat{\alpha},\beta)\)-scheme (with \((\hat{\alpha},\beta)=(1,0)\)), which seems out of reach in our setting, since the decay of the semigroup difference \(\|[S(t)-S(s)]x\|_X\) is limited by \((t-s)^1\) regardless of the smoothness of \(x\).

As shown in [28], pathwise uniform error bounds can be derived from pointwise ones via Hölder continuity of the \(p\)-th moment and the Kolmogorov–Chentsov theorem. However, for fixed integrability \(p\) of the initial values, this decreases the convergence rate by \(1/p\) [11].

To render the above results applicable to an implementable numerical scheme for SPDEs, a spatial discretization is required. As the main findings of the present work concern the temporal discretization, we refer the interested reader to the literature on space discretization, e.g.via spectral Galerkin methods as in our numerical simulation or [15], [29], finite differences [7], [30], finite elements [6], [10] or a discontinuous Galerkin approach [2], [9].

For nonlinear Nemytskii operators \(F\) and \(G\), the linear growth condition on \(Y\) limits one’s choice of \(Y\) to \(H^1\), prohibiting optimal convergence rates for second order equations with Nemytskii-type nonlinearities. However, optimal rates can be achieved for some non-Nemytskii nonlinearities (cf. Subsection 6.2). Relaxing the framework to allow for polynomial growth on \(Y\) would overcome this restriction and is left for future work. The analysis of the Milstein scheme in the globally Lipschitz setting constitutes a first step towards treating locally Lipschitz nonlinearities, which occur more frequently in model equations. A systematic error analysis of uniform strong errors for rational schemes in the locally Lipschitz setting would then provide a natural starting point for the analysis of Milstein schemes, both of which remain open.

Method of proof↩︎

We extend Kato’s framework [26], which was first used for SPDEs in [11], to systematically treat hyperbolic problems beyond convergence rate \(1/2\). This entails working with a pair of Hilbert spaces \(Y\hookrightarrow X\) such that \(F\) and \(G\) in 1 are, as in Theorem 1, Lipschitz continuous on \(X\) and of linear growth on \(Y\). Lipschitz continuity on higher order Sobolev spaces typically fails for Nemytskii operators and is thus not assumed. Although some parabolic equations are covered by the Kato setting, alternative approaches based on regularisation phenomena yield better results. Consequently, we focus on hyperbolic problems.

Achieving higher-order convergence relies on suitable Taylor expansions, which require Gâteaux differentiability and continuity of the derivative, which in turn imply Fréchet differentiability. Care must be taken in specifying between which spaces this differentiability assumption is imposed. Indeed, assuming Fréchet differentiability of \(F\colon X\to X\) for \(F\) a Nemytskii operator on \(X=L^2\) immediately limits the setting to affine-linear ‘nonlinearities’ \(F\). More precisely, consider the Nemytskii operator \(F\colon L^2(\mathcal{O})\to L^2(\mathcal{O}), F(u)\mathrel{\vcenter{:}}=\phi \circ u\) associated to a Lipschitz continuous \(C^1\)-function \(\phi\colon \mathbb{R}\to\mathbb{R}\), and \(\mathcal{O}\subseteq\mathbb{R}^d\) a bounded domain. If \(F\) is Fréchet differentiable at some \(\tilde{u}\in L^2(\mathcal{O})\), there exist \(a,b\in \mathbb{R}\) such that \(F(u)=au+b\) for all \(u\in L^2(\mathcal{O})\) (see also Proposition 5).

This was recently pointed out as a prevalent problem in the numerical analysis literature [31]. We avoid this pitfall by merely assuming Gâteaux differentiability of \(F\) as a map on \(X\). This is satisfied for any \(\phi\) as above and, together with Lipschitz continuity of \(F\) on \(X\), ensures uniform boundedness of the Gâteaux derivative \(F'\colon X \to \mathscr{L}(X)\). Furthermore, we assume Hölder continuity of the Gâteaux derivative \(F'\colon Y \to \mathscr{L}(Y,X)\), which ensures the existence of a suitable form of Taylor expansion. The latter also implies Fréchet differentiability of \(F\colon Y \to X\). However, since this only requires controlling the norm of the derivative of \(F\) uniformly over a unit ball in \(Y\) rather than \(X\), genuinely nonlinear Nemytskii operators are admissible provided that \(Y\) is sufficiently regular. For instance, in one spatial dimension, \(Y=H^s\) for \(s>\frac{1}{2}\) suffices due to the embedding into \(L^\infty(\mathcal{O})\). One then has the Taylor expansion \[\begin{align} F(u) &= F(v) + F'(v)\br{u-v}+ \int_{0}^{1} \big( F'(v+\zeta(u-v))-F'(v) \big)\br{u-v} \,\mathrm{d}\zeta \end{align}\] for all \(u,v\in Y\). Together with Lipschitz continuity of \(F'\colon Y \to \mathscr{L}(Y,X)\), which implies Fréchet differentiability of \(F\colon Y \to X\), this suffices to prove convergence at rate \(1\). Contrary to [15], we do not assume \(F\) and \(G\) to be twice Fréchet differentiable.

To consider pathwise uniform errors in \(L^p(\Omega)\) for general \(p\in [2,\infty)\) instead of root-mean-square errors (i.e.\(p=2\)) only as in [15], a stochastic Fubini argument is employed in one of the terms arising from the Taylor expansion and maximal inequalities for stochastic convolutions are used. A crucial ingredient that enables the extension from exponential to rational Milstein schemes is a logarithmic square function estimate ([11] and [32]) of the form \[\mathbb{E}\br[\bigg]{\sup_{i\in \{1, \ldots, n\}} \sup_{t\geq 0}\Big\|\int_0^t \Phi_i(s) \,\mathrm{d}W_s\Big\|_X^p} \lesssim \sqrt{\log(n)^p} \|(\Phi_i)_{i=1}^n\|^p\] for a suitable square function norm of \(\Phi\) (see Proposition 3 below).

Furthermore, regularity estimates are needed. It is only possible to show \(1/2\)-Hölder continuity of the \(p\)-th moment of the mild solution \(U\) to 1 in \(X\) but not in \(Y\). A different splitting of the error permits us to circumvent this: Rather than \((\mathbb{E}[\|U_{t}-U_{s}\|_X^p])^{1/p}\) for \(s\le t\), we estimate differences with the semigroup \((\mathbb{E}[\|U_{t}-S(t-s)U_{s}\|_X^p])^{1/p}\) in terms of \((t-s)^{1/2}\). Since some terms then vanish, the remaining terms can be estimated via linear growth and thus also in \(Y\), as illustrated in Lemma 6. On \(X\), this decay can be improved to \((t-s)^1\) by additionally considering the difference with the stochastic integral term \(\int_s^t G(S(t-s)U_{s})\,\mathrm{d}W_r\) (cf.Lemma 7). The proof is finished by an application of a discrete Grönwall inequality.

Overview↩︎

Section 2 recalls some facts from stochastic integration, Fréchet derivatives, and semigroup approximation required subsequently. We introduce our setting and assumptions in full generality in Section 3. Section 4 then recalls the underlying well-posedness results from the literature and contains a proof of pointwise strong stability of the Milstein scheme. Pathwise uniform convergence rates for the Milstein scheme are stated and proven in 5, where the main result in Theorem 15 extends Theorem 1 above. Improvements for the exponential Milstein scheme in Theorem 16 and the linear case in Corollary 2 as well as possible generalisations are discussed. Our results are illustrated for the stochastic Schrödinger, Maxwell’s, and transport equations in Section 6. Numerical simulations for different versions of the stochastic Schrödinger equation in Section 7 validate our theoretical findings.

Acknowledgements↩︎

The authors thank Foivos Evangelopoulos-Ntemiris, Mark Veraar, and Joris van Winden for helpful discussions and comments.

2 Preliminaries↩︎

Notation↩︎

Denote the natural numbers by \(\mathbb{N}\mathrel{\vcenter{:}}=\{1,2,3,\ldots\}\). Throughout the paper, \(H\), \(X\), and \(Y\) are separable Hilbert spaces, \((\Omega, \mathscr{F}, \mathbb{P})\) is a fixed probability space with filtration \((\mathscr{F}_t)_{t \in [0,T]}\) satisfying the usual conditions, and \((W_t)_{t\ge 0}\) denotes an \(H\)-cylindrical Brownian motion. For \(p\in[2,\infty)\), we abbreviate by \(\norm{\cdot}_p\) the canonical norm in \(L^p(\Omega)\) and write \(\mathcal{B}(X)\) for the Borel \(\sigma\)-algebra of \(X\) and \(\mathcal{P}\) for the predictable \(\sigma\)-algebra on \(\Omega\times[0,T]\). We use the subscript \(L_{\mathcal{G}}^p\) to denote \(\mathcal{G}\)-measurable elements or processes in \(L^p\). Denote by \(C^{\alpha}(I;X)\), or simply \(C^{\alpha}(I)\) if \(X=\mathbb{R}\), the space of bounded and \(\alpha\)-Hölder continuous functions \(\phi\colon I\to X\) for \(\alpha \in (0,1]\) and \(I\subseteq\mathbb{R}\). The notation \(f(x) \lesssim_{a,b} g(x)\) is used if there is a constant \(C \ge 0\) depending on \(a,b\) such that for all \(x\) in the respective set \(f(x) \le C g(x)\).

Fix a final time \(T>0\). We consider a uniform time grid \(\{t_{j} = jh:j=0,\ldots,M\}\) on \([0,T]\) with \(M\in \mathbb{N}\) steps and time step size \(h\mathrel{\vcenter{:}}= T/M>0\). For \(t \in [0,T]\), the time grid point before \(t\) is given by \(\lfloor t \rfloor \mathrel{\vcenter{:}}=\max\{t_j:\,t_j\le t\}\). We approximate the exact solution \(U_{}=(U_{t})_{t\in [0,T]}\) of a given evolution equation governed by a \(C_0\)-semigroup \(S=(S(t))_{t\geq 0}\) by a numerical solution \(u_{}=(u_{j})_{j=0,\ldots,M}\) given by the Milstein scheme. The computation of \(u_{}\) relies on a time discretization scheme \(R=(R_h)_{h>0}\) that approximates the semigroup \(S\).

2.1 Stochastic integration in Hilbert spaces↩︎

By \(\mathscr{L}_2(H,X)\), we denote the space of Hilbert–Schmidt operators from \(H\) to \(X\) and by \({\mathscr{L}_2^{(2)}(H,X)}\) the space of bilinear Hilbert–Schmidt operators, which consists of all bilinear bounded operators \(\Phi \colon H \times H \to X\) such that \[\|\Phi\|_{\mathscr{L}_2^{(2)}(H,X)}\mathrel{\vcenter{:}}=\p[\Big]{\sum_{m \in \mathbb{N}} \sum_{n\in\mathbb{N}} \|\Phi(h_m,h_n)\|_X^2}^{1/2}<\infty,\] where \((h_n)_{n\in\mathbb{N}}\) is an orthonormal basis of \(H\). It is easily checked that \({\mathscr{L}_2^{(2)}(H,X)}\cong \mathscr{L}_2(H;{\mathscr{L}_2(H,X)})\) is an isometric isomorphism and thus also \({\mathscr{L}_2^{(2)}(H,X)}\cong \mathscr{L}_2(\overline{H\otimes H},X)\) [33]. Here, \(\overline{H\otimes H}\) denotes the Hilbert space tensor product, which is given by the completion of the algebraic tensor product with respect to the canonical inner product in \(H\otimes H\).

For \(\phi \in {\mathscr{L}_2(H,X)}\) and a sequence \(\gamma = (\gamma_n)_{n\in\mathbb{N}}\) of centred i.i.d.normally distributed random variables we define \[\label{eq:convradonW} \phi \gamma \mathrel{\vcenter{:}}=\sum_{n\in\mathbb{N}} \gamma_n \phi h_n,\tag{2}\] where the convergence is in \(L^p(\Omega;X)\) for every \(p \in [1,\infty)\) and almost surely (see [33]).

We recall some fundamental definitions and statements on stochastic integration in Hilbert spaces from [34]. Given an \(\mathscr{L}_2(H,X)\)-valued integrand, we take stochastic integrals w.r.t.an \(H\)-cylindrical Brownian motion as integrator, which is a bounded linear mapping \(W_H\colon L^2(0,T;H) \to L^2(\Omega)\) such that

1. \(W_H b\) is centered Gaussian for all \(b \in L^2(0,T;H)\),

2. \(\mathbb{E}[W_H b_1 \cdot W_H b_2] = \langle b_1,b_2 \rangle_{L^2(0,T;H)}\) for all \(b_1, b_2 \in L^2(0,T;H)\),

3. \(W_H b\) is \(\mathscr{F}_t\)-measurable for all \(b\in L^2(0,T;H)\) supported in \([0,t]\),

4. \(W_H b\) is independent of \(\mathscr{F}_s\) for all \(b\in L^2(0,T;H)\) supported in \([s,T]\).

A complex \(H\)-cylindrical Brownian motion is defined analogously with a complex conjugate on \(W_H b_2\). For each fixed \(h\in H\) of unit norm, \((W_H(t)h)_{t \in[0,T]}\mathrel{\vcenter{:}}=(W_H(\mathbf{1}_{(0,t)} \otimes h))_{t \in[0,T]}\) is a (standard) Brownian motion. In the special case \(H=\mathbb{R}\), we recover the classical real-valued Brownian motion.

Given a linear, bounded, positive self-adjoint operator \(Q\in \mathscr{L}(H)\) of trace class, coloured noise can be described with the help of a \(Q\)-Wiener process. This is an equivalent description in the sense that \(W_Q\) is a \(Q\)-Wiener process if and only if \(Q^{1/2}W_H(t)\mathrel{\vcenter{:}}=\sum_{n\geq 1} Q^{1/2} h_n W_H(t) h_n = W_Q(t)\) for some \(H\)-cylindrical Brownian motion \(W_H\) for all \(t\in[0,T]\), where we used 2 . One can reduce the study of 1 with a \(Q\)-Wiener process \(W_Q\) to the case of cylindrical Brownian motion by replacing \(G\) by \(G Q^{1/2}\). Subsequently, we will omit the index \(H\) from \(W_H\) for the sake of readability. We now recall a standard property of stochastic integrals from [34].

A central estimate for stochastic integrals is the following maximal inequality for stochastic convolutions with a quasi-contractive semigroup.

If the semigroup considered is the identity, the classical Burkholder–Davis–Gundy inequalities are recovered. In this case, the inequality remains valid in the absence of a supremum, and is referred to as Itô’s isomorphism.

Moreover, the following logarithmic square function estimate is an essential tool in the convergence estimate of the rational Milstein scheme. Its original version [36] with \(\log(M)\)-scaling was improved to \(\sqrt{\log(M)}\)-scaling in [11] and generalised in [32].

2.2 Gâteaux and Fréchet differentiability↩︎

The analysis of the Milstein scheme relies on both Gâteaux and Fréchet derivatives, whose definition we recall. In the following, let \(Z_1\) and \(Z_2\) be two arbitrary Banach spaces.

Both Gâteaux and Fréchet derivatives are unique whenever they exist by linearity and uniqueness of limits in \(Z_2\). We point out that the limit in the definition of Gâteaux differentiability only involves the norm on the codomain \(Z_2\), whereas for the Fréchet derivative, also the norm on the domain \(Z_1\) appears. By setting \(y=\varepsilon z\) in the Fréchet derivative, one immediately sees that every Fréchet differentiable operator is Gâteaux differentiable and the derivatives agree. The converse is false in general, since only directional derivatives are considered for Gâteaux differentiability and no uniform norm control is required. However, under a continuity assumption on the Gâteaux derivatives, Fréchet differentiability is assured.

Nevertheless, Fréchet differentiability is not required to obtain uniform bounds on the Gâteaux derivative. Indeed, Lipschitz continuity is sufficient.

The case of Nemytskii operators \(B\) on \(Z_1=Z_2=L^2(\mathcal{O};\mathbb{R})\) for a bounded domain \(\mathcal{O}\subseteq\mathbb{R}^d\), \(d\in \mathbb{N}\), illustrates that Fréchet differentiability on all of \(L^2\) is a restrictive assumption.

In case \(\phi(x,s)\equiv \phi(s)\), the conditions on \(\phi\) in the proposition above reduce to \(\phi\in C^1(\mathbb{R})\) being Lipschitz continuous. Replacing \(Z_1\) by a subspace with a stronger norm, Fréchet differentiability no longer implies that \(B\) is affine linear.

2.3 Approximation of semigroups↩︎

A key step in solving a stochastic evolution equation numerically is to approximate the associated semigroup describing the behaviour of the linear part.

Common choices of schemes approximating the semigroup \(S\) generated by \(-A\) with their respective approximation orders include (cf.Corollary 4.4 in [38] with \(q=1\) for IE and \(q=2\) for CN):

• exponential Euler (EXE): \(R_h = S(h)\), any order \(\alpha >0\) on \(X\);

• implicit Euler (IE): \(R_h = (1+hA)^{-1}\), order \(\alpha \in (0,1]\) on \(\mathop{\mathrm{D}}(A^{2\alpha})\);

• Crank–Nicolson (CN): \(R_h = (2-hA)(2+hA)^{-1}\), order \(\alpha \in (0,2]\) on \(\mathop{\mathrm{D}}(A^{3\alpha/2})\) provided that \(S\) is contractive.

An essential assumption in our hyperbolic setting consists of contractivity of the semigroup and of the scheme approximating it. For EXE, contractivity of the scheme clearly is equivalent to contractivity of the semigroup. Both IE and CN are instances of rational schemes, which are schemes \(R\) that are induced by a rational function \(r\colon \mathbb{C}\to\mathbb{C}\cup \{\infty\}\) via \(R_h=r(h\cdot(-A))\) for all sufficiently small \(h>0\). This definition uses the bounded \(H^\infty\)-calculus of \(A\) as the negative generator of a \(C_0\)-contraction semigroup to define \(R_h\). Under additional assumptions on \(r\) on the open left half-plane \(\mathbb{C}_-\), contractivity of the scheme \(R\) is recovered as a consequence of [33].

This proposition ensures the contractivity of IE, CN, but also higher-order implicit Runge–Kutta methods including Radau methods, Lobatto IIIA, IIIB, and IIIC as well as some DIRK schemes. More generally, all A-stable schemes [39] and all A-acceptable schemes [40] satisfy the condition above on \(r\), where the latter require it to hold on \(\overline{\mathbb{C}_-}\) and an additional consistency condition. As is the case for these examples, a common choice for the spaces \(Y\) on which the semigroup \(S\) is approximated to some rate depending on \(\alpha>0\) consists of domains \(\mathop{\mathrm{D}}(A^\alpha)\) of fractional powers of \(A\). It is useful to know that the embedding \(\mathop{\mathrm{D}}(A^{\alpha}) \hookrightarrow \mathop{\mathrm{D}}_A(\alpha, \infty) \mathrel{\vcenter{:}}=(X,\mathop{\mathrm{D}}(A))_{\alpha,\infty}\) into the real interpolation space with parameter \(\infty\) holds. For details on interpolation spaces, we refer the interested reader to [41], [42].

2.4 A discrete Grönwall inequality↩︎

The following variant of a discrete Grönwall inequality [11] based on [43] proves to be useful.

3 Setting and Assumptions↩︎

We consider the stochastic evolution equation \[\label{eq:SEE} \begin{cases} \mathrm{d}U + AU \,\mathrm{d}t&= F(t,U)\,\mathrm{d}t+ G(t,U) \,\mathrm{d}W\quad \text{on }(0,T], \\ U_{0}&=\xi\in L_{\mathscr{F}_0}^p(\Omega;X). \end{cases}\tag{3}\] on a Hilbert space \(X\) for some \(p \in [2,\infty)\) and \(T>0\). Here, \(-A\) generates a contractive \(C_0\)-semigroup, \(F\) is a nonlinearity, \(G\) multiplicative noise, and \((W_t)_{t\ge 0}\) is an \(H\)-cylindrical Brownian motion for some separable Hilbert space \(H\). Assumptions 1, 2[assCond:FG95rate95Lipschitz], and 2[assCond:FG95rate95Yinvariance] ensure that 3 is well-posed (see Theorem 11 for the precise statement).

The embedding of \(Y\) into a real interpolation space allows us to obtain decay rates for semigroup differences via interpolation [11].

The conditions [assCond:FG95rate95Lipschitz]–[assCond:FG95rate95lineargrowthY] of the following assumption are a slight modification of [11] for schemes without Milstein terms.

We consider contractive approximations of the semigroup associated with the stochastic evolution equation and recall Definition 3.

Admissible choices for \(R\) include the exponential Euler (EXE) method \(R=S\), the implicit Euler (IE) method \(R_h=(1+hA)^{-1}\), the Crank–Nicolson (CN) method \(R_h=(2-hA)(2+hA)^{-1}\), and other A-stable rational schemes, as discussed in Subsection 2.3. In order to approximate the mild solution \(U\) of 3 in time, the rational and exponential Milstein schemes are considered.

The series in 6 shall be understood as the \(L^p\)-limit of rectangular partial sums. Indeed, they converge in \(L^p(\Omega;X)\) for every \(p\in [1,\infty)\) by Itô’s isometry for \(p=2\), monotonicity for \(p\in[1,2)\), and for \(p>2\) by Nelson’s hypercontractivity estimate for the second Wiener chaos [44].

Recalling Notation 7 and Lemma 1, we can reformulate the noise-related terms of the rational Milstein scheme as \[\begin{align} \label{eq:rewriteGpG} &\phantom{\le }\sum_{i=0}^{j-1} R_h^{j-i}\p[\Big]{\int_{t_{i}}^{t_{i+1}}G(t_{i},u_{i})\,\mathrm{d}W_s+ \int_{t_{i}}^{t_{i+1}} G'(t_{i},u_{i})\br[\Big]{\int_{t_{i}}^s G(t_{i},u_{i}) \,\mathrm{d}W_r}\,\mathrm{d}W_s}\nonumber\\ &= \sum_{i=0}^{j-1} R_h^{j-i}\p[\Big]{\int_{t_{i}}^{t_{i+1}}G(t_{i},u_{i})\,\mathrm{d}W_s+ \int_{t_{i}}^{t_{i+1}} \int_{t_{i}}^s (G'G)(t_{i},u_{i}) \,\mathrm{d}W_r\,\mathrm{d}W_s}. \end{align}\tag{7}\]

We focus on the case \(\alpha>\frac{1}{2}\). Otherwise, the exponential Euler scheme or rational schemes, which are both using only the Wiener increments of the noise, are easier to implement and allow for the same pathwise uniform convergence rate \(\alpha \in (0,\frac{1}{2}]\) under less restrictive assumptions and with a shorter proof, see [11].

4 Well-posedness and Stability↩︎

This section contains the well-posedness and stability results required to prove the convergence result in the next section. Here, well-posedness shall be understood in the sense of existence and uniqueness of mild solutions to 3 satisfying an a priori estimate. Stability refers to pointwise stability of the rational Milstein scheme, i.e.moment bounds uniform in the number of time steps.

As a consequence of the definition of the mild solution to 3 , for all \(s,t\in[0,T]\) with \(s\leq t\), \[\label{eq:mild-solution-consequence} U_{t} = S(t-s)U_{s} + \int_{s}^{t}S(t-r)F(r,U_{r})\,\mathrm{d}r+ \int_{s}^{t}S(t-r)G(r,U_{r}) \,\mathrm{d}W_r.\tag{8}\] The following standard well-posedness theorem for globally Lipschitz nonlinearities on \(X\) (see [34] and [11]) is included for completeness. Additionally assuming merely linear growth on \(Y\) rather than Lipschitz continuity, [11] yields well-posedness on \(Y\).

As an inspection of the proofs shows, the statement of the theorem remains true with the same constants if one is added on the left-hand side. Hence, for \(q\in [2,\infty)\) and \(Z \in \{X,Y\}\), the constants \[\begin{align} \label{eq:defCWPqZ} C_{q,Z}^{\mathrm{WP}} &\mathrel{\vcenter{:}}= C_{q,Z}^{\mathrm{bdd}}(1 + \|\xi\|_{L^q(\Omega;Z)} + \|f\|_{L^q(\Omega;L^1(0,T;Z))} + B_q \norm{g}_{L^q(\Omega;L^2(0,T;\mathscr{L}_2(H,Z)))}) \end{align}\tag{9}\] are finite under the assumptions of Theorem 11. Moreover, the estimate \[\label{eq:WPboundCWPqZ} 1+\sup_{t \in [0,T]} \|U_{t}\|_{L^q(\Omega;Z)} \le C_{q,Z}^{\mathrm{WP}} <\infty\tag{10}\] holds. Note that the additional \(1\) on the left-hand side here and in the stability estimate below is included for notational convenience. The following shorthand notation is used throughout the paper.

5 Convergence Rates↩︎

Our main goal is to establish pathwise uniform convergence rates for the rational and the exponential Milstein scheme as a temporal approximation of stochastic evolution equations of the form \[\label{eq:SEEconv} \begin{cases} \mathrm{d}U + AU \,\mathrm{d}t&= F(t,U)\,\mathrm{d}t+ G(t,U) \,\mathrm{d}W\quad \text{on }(0,T], \\ U_{0}&=\xi\in L_{\mathscr{F}_0}^p(\Omega;X). \end{cases}\tag{11}\] on \([0,T]\), where \(-A\) generates a contractive \(C_0\)-semigroup \((S(t))_{t\ge 0}\) on a Hilbert space \(X\) and \(U\) is the mild solution. The conditions on the nonlinearity \(F\) and the multiplicative noise \(G\) as well as \(A\), \(\xi\), \(X\), and \(W\) are as in Section 3. As we have seen in Section 4, these ensure the well-posedness of 11 both in \(X\) and a more regular continuously embedded subspace \(Y\) as well as pointwise stability of the rational Milstein scheme from Definition 4 in \(Y\). We recall Notation 12.

This section is split into two subsections: The central error estimates are presented, proved, and discussed in Subsection 5.1. Subsequently, the error estimates are generalised to a suitable extension of the rational Milstein scheme to the full time interval \([0,T]\) in Subsection 5.2. Before discussing our main result in Theorem 15, a pathwise uniform convergence rate for the rational Milstein scheme, we first state two useful lemmas on the path regularity of the mild solution to 11 . Two of the error terms are estimated in the auxiliary Lemmas 8 and 9 as part of the proof of Theorem 15. The error estimate is extended to \(p\in [1,2)\) in Corollary 1 and improved for exponential Milstein schemes in our second main result, Theorem 16. Possible generalisations regarding less regular nonlinearities, weaker differentiability assumptions, or \(2\)-smooth Banach spaces are discussed in Remarks 17–19 and simplifications in the linear case in Corollary 2. The error estimate remains valid for a suitable extension of the scheme to \([0,T]\), as shown in Theorem 20.

5.1 Main error estimates at the grid points↩︎

In order to achieve decay of order higher than \(\frac{1}{2}\), an additional stochastic integral term is taken into account in the difference. This yields decay of order \(\min\{\alpha+\frac{1}{2},1\}\) on the space \(X\).

We are now in the position to state and prove the main result of this paper.

The statement remains true if one merely assumes \(f\in L_\mathcal{P}^\infty(0,T;L^{2\alpha p}(\Omega;Y)) \cap C([0,T];L^p(\Omega;Y))\) rather than \(f\in C([0,T];L^{2\alpha p}(\Omega;Y))\) as above, and likewise for \(g\) with values in \({\mathscr{L}_2(H,Y)}\). Moreover, it suffices if Assumption 2[assCond:FG95rate95HolderContTime] holds for \(q=p\) or if the semigroup \((S(t))_{t\ge 0}\) is quasi-contractive (cf.Definition 1), where the latter can be seen by a simple scaling argument.

Since \(\|\phi\|_{L^r(\Omega)} \le \|\phi\|_{L^p(\Omega)}\) for all \(1\le r\le p\) and \(\phi \in L^p(\Omega)\), the error estimate extends to lower moments \(r\in [1,2)\). Note that the regularity assumptions such as integrability of the initial values are required to hold for some \(p\ge 2\), with \(p=2\) corresponding to the weakest assumptions.

To prove the main error estimate, we define the following auxiliary time-discrete stochastic processes, which allow us to split the error suitably and analyse two of the error terms in the subsequent lemmas.

In the following lemma, the error arising from the difference of the nonlinearities evaluated at the solution and at the approximation at the grid points is analysed. It is estimated in terms of sums of the full error at previous grid points, which will finally be dealt with by a discrete Grönwall argument in the proof of Theorem 15.

Next, we estimate the error \(E_3\) caused by the rational approximation of the semigroup. To this end, we can no longer make use of the maximal inequality for stochastic convolutions of Theorem 2, since the semigroup is approximated by a rational scheme, leading to a discrete stochastic convolution. Instead, we apply the logarithmic square function estimate from Proposition 3, which introduces a logarithmic correction factor in the estimate, but in turn allows us to proceed via the approximation order \(\alpha\) of the scheme \(R\) and linear growth of \(G-g\) on \(Y\).

We can now pass to the proof of the main theorem.

The logarithmic correction factor can be omitted for the exponential Milstein scheme, as it is not necessary to use the logarithmic square function estimate for terms not containing the semigroup.

We comment on the case of lower rates \(\alpha \in (0,1/2]\) and make two observations regarding possible generalisations of Theorems 15 and 16.

In the linear case, the necessary assumptions simplify considerably, since for example Lipschitz continuity is immediate and the Gâteaux derivatives \(F'(u)[v]=F(v)\) and \(G'(u)\) are independent of \(u\), and thus trivially Hölder continuous. In the following we also allow for affine linear coefficients.

5.2 Error estimates on the full time interval↩︎

The error estimates of Theorems 15 and 16 hold for the difference between the solution \(U_{t}\) at the grid points \(t=t_j\) and the scheme \(u_{j}\). The convergence rate obtained there remains valid for arbitrary times \(t\in [0,T]\) if one chooses a suitable extension of the Milstein scheme to general times \(t\). For the extension introduced in Definition 8, such an error estimate is contained in Theorem 20, whose proof relies on the maximal estimates for (iterated) stochastic integral increments of Lemma 10.

Lemma 10 also holds for \(\mathcal{X}\) and \(Z\) being \(2\)-smooth Banach spaces if one replaces Hilbert–Schmidt by \(\gamma\)-radonifying operators and \(K\) by some space-dependent constant (see also Remark 19).

Note that \(\bar{u}\) is a stochastic continuous-time extension of the rational Milstein scheme \(u\) but not an interpolation constructed only from its grid values \(u_{j}\). While the nonlinearities \(F,G,G'G\) are only evaluated at the grid points \(t_{j}\) and at \(u_{j}\) for all \(j\) such that \(t_{j}\le t\) to compute \(\bar{u}_t\), the scheme \(\mathcal{R}_h(t)\), the Wiener increments, and the iterated stochastic integral increments are used inside the intervals. By contrast, interpolations of the grid values are easier to compute. However, the error of a piecewise constant extension of the rational Milstein scheme only decays at the lower rate \(h^{1/2}\sqrt{\log(T/h)}\) due to the path regularity of \(U\), which can be seen by combining the proof of [11] with Theorem 15. For the scalar SDE \(\mathrm{d}U_t=\mathrm{d}W_t\), the same holds for a piecewise affine interpolation, since its error is the Brownian bridge (see [48]), whose supremum is of order \(h^{1/2}\sqrt{\log(M)}\) over \(M\) intervals of length \(h\). Moreover, for the SDE the iterated stochastic integral over \([t_j,t_{j+1}]\) does not provide additional information, since it equals \(\frac{1}{2}((W_{t_{j+1}}-W_{t_{j}})^2-h)\).

6 Applications to Hyperbolic Equations↩︎

We consider a linear and a nonlinear stochastic Schrödinger equation, the stochastic Maxwell’s equations, and the stochastic transport equation with a Nemytskii-type nonlinearity in Subsections 6.1, 6.2, 6.3, and 6.4, respectively. Throughout this section, we abbreviate the exponential Euler, Crank–Nicolson, and implicit Euler methods (cf.Subsection 2.3 for their definition) by EXE, CN, and IE. The corresponding Milstein schemes are referred to as EX-Milstein (or exponential Milstein), CN-Milstein, and IE-Milstein.

In the rest of the paper, we consider SPDEs with coloured noise driven by a \(Q\)-Wiener process \(W_Q\) with trace class operator \(Q\) (cf.Subsection 2.1 for the definition) rather than an \(H\)-cylindrical Brownian motion \(W\). This allows for an easier comparison with results from the literature [1], [4], [11], where \(Q\)-Wiener processes were used. Recall from Subsection 2.1 that every \(Q\)-Wiener process can be written in the form \(W_Q(t)=Q^{1/2}W_H(t)\) for some \(H\)-cylindrical Brownian motion \(W_H\) for all \(t\in[0,T]\). Hence, the results from the previous sections are applicable to stochastic evolution equations with noise term \(\mathbb{G}(t,u)\,\mathrm{d}W_Q\) by considering \(G(t,u) \mathrel{\vcenter{:}}=\mathbb{G}(t,u) Q^{1/2}\).

6.1 Linear stochastic Schrödinger equation↩︎

Denote by \(\mathbb{T}^d\) the \(d\)-dimensional torus for \(d\in \mathbb{N}\) and let \(L^q\mathrel{\vcenter{:}}= L^q(\mathbb{T}^d;\mathbb{C})\) for \(q\in[2,\infty]\) as well as \(H^\sigma \mathrel{\vcenter{:}}= H^\sigma(\mathbb{T}^d;\mathbb{C})\) for \(\sigma \ge 0\) in this and the subsequent subsection. Consider the linear stochastic Schrödinger equation \[\label{eq:schroedinger-linear} \begin{cases} \mathrm{d}U &= -{\mathrm{i}}(\Delta U + V\cdot U) \,\mathrm{d}t- {\mathrm{i}}U \,\mathrm{d}W_Q \quad \text{on }(0,T], \\ U_0 &= \xi \end{cases}\tag{38}\] on the Hilbert space \(X\mathrel{\vcenter{:}}= H^\sigma\), \(\sigma\ge 0\), with a potential \(V\colon\mathbb{T}^d\to\mathbb{C}\) and a \(Q\)-Wiener process \(W_Q\) with trace class operator \(Q\). In the abstract setting, we thus consider \(H\mathrel{\vcenter{:}}= L^2\) and the Schrödinger operator \(A\mathrel{\vcenter{:}}={\mathrm{i}}\Delta\) with domain \(\mathop{\mathrm{D}}(A)=\{f\in H^\sigma: {\mathrm{i}}\Delta f \in H^\sigma\}=H^{\sigma+2}\) to account for \(X=H^\sigma\), the linear operator \(F_0=-{\mathrm{i}}M_V\) given by multiplication by the potential \(V\) and the noise \(G_0 = [u \mapsto -{\mathrm{i}}M_u Q^{1/2}]\). In the following cases, we verify that \(F_0\) and \(G_0\) satisfy the conditions of Corollary 2 for the exponential Milstein scheme with \(Y\mathrel{\vcenter{:}}= H^{\sigma+2\alpha}\) for some \(\alpha \in (0,1]\). These cases are obtained by combining all possible cases in which \(X=H^\sigma\) and \(Y=H^{\sigma+2\alpha}\) satisfy [11] and rearranging by dimension.

We could also include \(d=1\), \(\sigma \in [0,\frac{1}{2})\), \(\alpha\in (0,\frac{1}{4}-\frac{\sigma}{2})\), and \(\beta>\frac{1}{2}\) as a fifth case, or as a sixth case \(d\ge 2\), \(\sigma \in [0,1)\), \(\alpha\in (0,\frac{1-\sigma}{2}]\), and \(\beta>\frac{d}{2}\) with analogous arguments to the ones presented below. However, both cases imply \(\alpha\le 1/2\), which is the parameter regime in which Milstein schemes are not advantageous (see Remark 17).

For rational Milstein schemes, the additional assumption that \(R\) approximates \(S\) to rate \(\alpha\in (0,1]\) on \(Y\) enters, leading to higher regularity assumptions on \(Y\). Namely, we can choose \(Y=H^{\sigma+\ell \alpha}\) with \(\ell=4\) for IE-Milstein and \(\ell=3\) for CN-Milstein. To see this, recall the definition of IE and CN from Subsection 2.3. Note that IE approximates \(S\) to order \(\alpha\) on \(\mathop{\mathrm{D}}(A^{2\alpha})\), where the fractional domain is given by \(\mathop{\mathrm{D}}(({\mathrm{i}}\Delta)^{2\alpha}) = H^{\sigma+4\alpha}\) for the Schrödinger equation. Likewise, CN approximates \(S\) to order \(\alpha\) on \(\mathop{\mathrm{D}}(A^{3\alpha/2})=H^{\sigma+3\alpha}\). This observation allows us to generalise Assumption 5 to rational Milstein schemes.

Here, the cases that only allow for suboptimal rates \(\alpha\le \frac{1}{2}\) are included. This can be useful when comparing the rates of convergence in regularity regimes where the exponential Milstein method attains the optimal rate \(1\) but for e.g.IE-Milstein, the regularity of \(V\) and \(Q\) limit us to rate \(\frac{1}{2}\).

Assumption 5 for EX-Milstein corresponds to \(\ell=2\) in the above. Even then, it can be beneficial to consider \(\ell>2\) in order to treat higher dimensions. If the stronger assumptions for \(\ell>2\) are satisfied, the \(L^2\)-error (\(\sigma=0\)) could be considered in dimension four, which was not feasible in the setting of Assumption 5[case:SoLin45D].

We point out that, up to the logarithmic factor, rational Milstein schemes attain the optimal rate of convergence \(1\) for a smoother potential and noise compared to the exponential Milstein scheme. This improves the rates of up to \(1/2\) obtained for the exponential Euler method in [1] as well as [11] and for IE and CN (without Milstein term) in [11].

6.2 Stochastic Schrödinger equation with a nonlocal nonlinearity↩︎

As a second example, we again consider the stochastic Schrödinger equation, now with nonlinear \(F\) and \(G\), which are nonlocal due to convolution. More precisely, consider \[\label{eq:nl-schroedinger} \begin{cases} \mathrm{d}U &= -{\mathrm{i}}\p[\big]{\Delta U + \eta\ast [\phi(U)] } \,\mathrm{d}t- {\mathrm{i}}\kappa \ast [\psi(U)] \,\mathrm{d}W_Q\quad \text{on }(0,T], \\ U_0 &= \xi \end{cases}\tag{40}\] on \(H^\sigma=H^\sigma(\mathbb{T}^d;\mathbb{C})\) for some \(\sigma\ge 0,~d\in \mathbb{N}\), nonlinearities \(\phi,\psi\colon\mathbb{C}\to\mathbb{C}\), convolution kernels \(\eta,\kappa\colon\mathbb{T}^d\to\mathbb{C}\), and a \(Q\)-Wiener process \(W_Q\).

In the abstract setting, this amounts to \(X\mathrel{\vcenter{:}}= H^\sigma\), \(Y\mathrel{\vcenter{:}}= H^{\sigma+\ell\alpha}\), where \(\alpha\in (0,1]\) and \(\ell\in [2,\infty)\), and \(H\mathrel{\vcenter{:}}= L^2\) with non-Nemytskii-type nonlinearities \[\begin{align} F\colon H^\sigma &\to H^\sigma,\quad u \mapsto -{\mathrm{i}}\eta \ast [\phi(u)]=\br[\Big]{ x \mapsto - {\mathrm{i}}\int_{\mathbb{T}^d} \eta(y)\phi(u(x-y)) \,\mathrm{d}y } \end{align}\] and \(G\colon H^\sigma \to \mathscr{L}_2(L^2,H^\sigma)\), \(G(u)\mathrel{\vcenter{:}}=-{\mathrm{i}}M_{\kappa \ast [\psi(u)]}Q^{1/2}\). Under the following assumption, we obtain pathwise uniform convergence at rate \(\alpha\).

6.3 Stochastic Maxwell’s equations↩︎

The stochastic Maxwell’s equations \[\begin{align} \label{eq:Maxwell} \Bigg\{\begin{aligned} \mathrm{d}U &= [AU+F(U)] \,\mathrm{d}t+ G(U)\,\mathrm{d}W_Q\quad \text{on }(0,T],\\ U_{0}&=\xi=(\mathbf{E}_0^\top,\mathbf{H}_0^\top)^\top \end{aligned} \end{align}\tag{41}\] with boundary conditions of a perfect conductor as in [4] describe how the electric and magnetic field \(\mathbf{E}\) and \(\mathbf{H}\) behave. This example is treated in [4] for the exponential Euler method and in [11] additionally for rational schemes, all without Milstein terms. For the sake of readability, we repeat some arguments in our analysis of Milstein schemes for 41 .

For a bounded, simply connected domain \(\mathcal{O}\subseteq\mathbb{R}^3\) with smooth boundary, let \(X\mathrel{\vcenter{:}}= L^2(\mathcal{O})^6=L^2(\mathcal{O})^3\times L^2(\mathcal{O})^3\) be equipped with the weighted scalar product \[\label{eq:defScalarProdMaxwell} \left\langle \begin{pmatrix} \mathbf{E}_1\\\mathbf{H}_1 \end{pmatrix},\begin{pmatrix} \mathbf{E}_2\\\mathbf{H}_2 \end{pmatrix}\right\rangle_X \mathrel{\vcenter{:}}=\int_\mathcal{O}\big(\mu(x) \langle \mathbf{H}_1,\mathbf{H}_2\rangle + \varepsilon(x) \langle \mathbf{E}_1,\mathbf{E}_2\rangle\big)\,\mathrm{d}x,\tag{42}\] where \(\langle\cdot,\cdot\rangle\) denotes the standard scalar product in \(L^2(\mathcal{O})^3\) and the permittivity and permeability \(\varepsilon, \mu \in L^\infty(\mathcal{O})\) are assumed to be uniformly positive, i.e.\(\varepsilon,\mu \ge \kappa >0\) for some constant \(\kappa\in(0,\infty)\). The Maxwell operator \(A\colon \mathop{\mathrm{D}}(A) \to X\) is given by \[A\begin{pmatrix} \mathbf{E}\\ \mathbf{H}\end{pmatrix} \mathrel{\vcenter{:}}=\begin{pmatrix} 0 & \varepsilon^{-1}\nabla\times\\-\mu^{-1}\nabla \times&0 \end{pmatrix}\begin{pmatrix} \mathbf{E}\\ \mathbf{H}\end{pmatrix} = \begin{pmatrix} \varepsilon^{-1}\nabla\times \mathbf{H}\\-\mu^{-1}\nabla \times \mathbf{E}\end{pmatrix}\] on \(\mathop{\mathrm{D}}(A) \mathrel{\vcenter{:}}= H_0(\mathop{\mathrm{curl}},\mathcal{O}) \times H(\mathop{\mathrm{curl}},\mathcal{O})\) with \(H(\mathop{\mathrm{curl}},\mathcal{O}) \mathrel{\vcenter{:}}=\{\mathbf{H}\in (L^2(\mathcal{O}))^3\,:\,\nabla \times \mathbf{H}\in L^2(\mathcal{O})^3\}\) and \(H_0(\mathop{\mathrm{curl}},\mathcal{O})\) the subspace of \(\mathbf{H}\) with vanishing tangential trace \(\mathbf{n} \times \mathbf{H}\vert_{\partial\mathcal{O}}\), where \(\mathbf{n}\) is the unit outward normal vector. Moreover, \(W_Q\) is a \(Q\)-Wiener process for symmetric, non-negative \(Q\) with finite trace such that \(Q^{1/2}\in \mathscr{L}_2(H,X)\) and further regularity specified below. We equip \(H\mathrel{\vcenter{:}}= L^2(\mathcal{O})^6\) with the canonical norm. In 41 , \(F\colon [0,T] \times X \to X\) is the linear drift term given by \[\label{eq:defFMaxwell} F(t,V)(x)=\begin{pmatrix} \sigma_1(x,t) \mathbf{E}_V(x)\\\sigma_2(x,t) \mathbf{H}_V(x)\end{pmatrix},\quad V=\begin{pmatrix}\mathbf{E}_V\\\mathbf{H}_V\end{pmatrix}\in L^2(\mathcal{O})^6,~~ x\in \mathcal{O}\tag{43}\] for \(\sigma_1,\sigma_2\colon \mathcal{O}\times [0,T] \to \mathbb{R}\), whose smoothness is specified later. In the notation of Corollary 2, the error estimate for the linear case, this corresponds to \(f\equiv 0\) and \[F_0\colon \Omega\times [0,T]\to \mathscr{L}(L^2(\mathcal{O})^6),~ F_0(\omega,t)(V)\mathrel{\vcenter{:}}=(\sigma_1(\cdot,t)\mathbf{E}_V^\top,\sigma_2(\cdot,t)\mathbf{H}_V^\top)^\top.\] As noise \(G(V)\) for \(V=(\mathbf{E}_V^\top,\mathbf{H}_V^\top)^\top\in L^2(\mathcal{O})^6\), we consider the Nemytskii map associated to \(\mathop{\mathrm{diag}}((-\varepsilon^{-1}\mathbf{E}_V^\top,-\mu^{-1}\mathbf{H}_V^\top))Q^{1/2}\). We can recast this in the linear setting by choosing \(g\equiv 0\) and \(G_0\colon\Omega \times [0,T]\to \mathscr{L}(X,\mathscr{L}_2(H,X))\) defined via \[\begin{align} \label{eq:defG0Maxwell} (G_0(\omega,t)(V))(h) &\mathrel{\vcenter{:}}=\begin{pmatrix} -\varepsilon^{-1}\mathop{\mathrm{diag}}(\mathbf{E}_V)&0\\ 0&-\mu^{-1}\mathop{\mathrm{diag}}(\mathbf{H}_V) \end{pmatrix} (Q^{1/2}h)\in L^2(\mathcal{O})^6 \end{align}\tag{44}\] for all \(h\in L^2(\mathcal{O})^6\). To verify the assumptions of Corollary 2 for Maxwell’s equations, regularity of \(\sigma_1\) and \(\sigma_2\) is required in addition to the assumption on \(\mu,\varepsilon\), which we repeat.

Indeed, uniform boundedness of \(F_0\) in \(\mathscr{L}(X)\)-norm by \(C_{F_0,X}\mathrel{\vcenter{:}}=\max\{\|\sigma_1\|_\infty,\|\sigma_2\|_\infty\}\) follows via \[\begin{align} \sup_{t\in[0,T]} \|F_0(t)\|_{\mathscr{L}(X)}^2 &= \sup_{t\in[0,T]} \sup_{\|V\|_X=1} \int_\mathcal{O}\left( \mu(x) \|\sigma_2(\cdot,t)\mathbf{H}_V\|_{L^2(\mathcal{O})^3}^2+\varepsilon(x) \|\sigma_1(\cdot,t)\mathbf{E}_V\|_{L^2(\mathcal{O})^3}^2\right)\,\mathrm{d}x\\ &\le \sup_{\|V\|_X=1} \sup_{t\in[0,T]} \max\{\|\sigma_1(\cdot,t)\|_\infty,\|\sigma_2(\cdot,t)\|_\infty\}^2 \|V\|_X^2 = C_{F_0,X}^2. \end{align}\] As the space \(Y\), we choose \(Y\mathrel{\vcenter{:}}=\mathop{\mathrm{D}}(A)\). Since \[\|F_0(t)\|_{\mathscr{L}(Y)}^2 = \sup_{\|V\|_{\mathop{\mathrm{D}}(A)}=1} \|F_0(t)V\|_Y^2 = \sup_{\|V\|_{\mathop{\mathrm{D}}(A)}=1} \p[\big]{\|AF_0(t)V\|_X^2+\|F_0(t)V\|_X^2}\] and the arguments above allow us to estimate the second term by \(C_{F_0,X}^2\|V\|_X^2\), it suffices to estimate the first term. By an explicit calculation of the curl operator, \[\begin{align} \|AF_0(t)V\|_X^2 &\le \kappa^{-2} \int_\mathcal{O}\mu \|\nabla \times (\sigma_1(\cdot,t) \mathbf{E}_V)\|_{L^2(\mathcal{O})^3}^2 + \varepsilon \|\nabla \times(\sigma_2(\cdot,t)\mathbf{H}_V)\|_{L^2(\mathcal{O})^3}^2\,\mathrm{d}x\\ &\le 3\kappa^{-2}\p[\Big]{C_{F_0,X}^2 \|AV\|_X^2+2\p[\Big]{\max_{j=1,2}\max_{i=1,2,3} \|\partial_{x_i}\sigma_j\|_\infty^2}\|V\|_X^2}. \end{align}\] Thus, we conclude uniform boundedness of \(\|F_0(t)\|_{\mathscr{L}(Y)}\) by \[C_{F_0,Y}\mathrel{\vcenter{:}}=\max\Big\{\sqrt{3}\kappa^{-1}C_{F_0,X}, \sqrt{6}\kappa^{-1}\max_{j=1,2}\max_{i=1,2,3} \|\partial_{x_i}\sigma_j\|_\infty+C_{F_0,X}\Big\}.\] Moreover, \(t \mapsto F_0(t)\) is Lipschitz continuous, since by Hölder’s inequality in \(L^2(\mathcal{O})\) and uniform Lipschitz continuity of \(\sigma_1\), \(\sigma_2\), we have \[\begin{align} \|F_0(t)-F_0(s)\|_{q_\alpha,\mathscr{L}(X)}^2 &= \sup_{\|V\|_X=1} \int_\mathcal{O}\mu\|(\sigma_2(\cdot,t)-\sigma_2(\cdot,s))\mathbf{H}_V\|^2+\varepsilon\|(\sigma_1(\cdot,t)-\sigma_1(\cdot,s))\mathbf{E}_V\|^2\,\mathrm{d}x\\ &\le \sup_{\|V\|_X=1} \max_{j=1,2}\|\sigma_j(\cdot,t)-\sigma_j(\cdot,s)\|_{L^\infty(\mathcal{O})}^2 \|V\|_X^2 \le \max_{j=1,2}\, C_{\sigma_j}^2(t-s)^2, \end{align}\] where the norms in the first line are taken in \({L^2(\mathcal{O})^3}\). Uniform boundedness of \(G_0\) w.r.t.the norm in \(\mathscr{L}(X,{\mathscr{L}_2(H,X)})\) by \(C_{G_0,X}\mathrel{\vcenter{:}}=\kappa^{-1}C_{H^\beta \hookrightarrow L^\infty}\|Q^{1/2}\|_{\mathscr{L}_2(L^2(\mathcal{O})^6,H^\beta(\mathcal{O})^6)}\) follows from the definitions, uniform positivity of the coefficients \(\mu,\varepsilon\), Hölder’s inequality, and the embedding \(H^\beta(\mathcal{O}) \hookrightarrow L^\infty(\mathcal{O})\) for any \(\beta>\frac{3}{2}\). W.r.t.the norm in \(\mathscr{L}(Y,{\mathscr{L}_2(H,Y)})\), this is a consequence of \(Q^{1/2} \in \mathscr{L}_2(L^2(\mathcal{O})^6,H^{1+\beta}(\mathcal{O})^6)\) again for \(\beta>\frac{3}{2}\) and the linear growth estimate [4] \[\begin{align} \|G(V)\|_{\mathscr{L}_2(H,\mathop{\mathrm{D}}(A))} &\le C \|Q^{1/2}\|_{\mathscr{L}_2(L^2(\mathcal{O})^6,H^{1+\beta}(\mathcal{O})^6)}(1+\|V\|_{\mathop{\mathrm{D}}(A)}) \end{align}\] for some \(C>0\) as detailed in [11]. Since \(G_0\) is independent of \(t\), temporal Hölder continuity is trivially satisfied. Reasoning as in [11] based on [50] and [4], we see that \(X\), \(Y\), and \((S(t))_{t\ge 0}\) fulfil Assumption 1. We can thus improve the convergence rates \(1/2\) obtained for Maxwell’s equations in [4] for EXE and \(1/2\) up to logarithmic correction in [11] for IE and CN to up to \(1\) for Milstein schemes.

This theorem illustrates that rational Milstein schemes result in worse rates of convergence compared to the exponential Milstein scheme given the same regularity of the coefficients and noise. The IE-Milstein scheme only converges at rate \(1/2\) up to a logarithmic correction factor, just like the IE method without Milstein term [11], which should be preferred in numerical simulations (see Remark 17). Increasing the regularity such that the assumptions of Corollary 2 are satisfied for \(Y=\mathop{\mathrm{D}}(A^2)\) and \(Y=\mathop{\mathrm{D}}(A^{3/2})\), we expect the IE- and CN-Milstein, respectively, to achieve the optimal rate of convergence \(1\) up to a logarithmic correction factor. However, to the authors’ best knowledge, no explicit characterisations of \(\mathop{\mathrm{D}}(A^\alpha)\) or \(\mathop{\mathrm{D}}_A(\alpha,\infty)\) are available for non-integer \(\alpha\) for the Maxwell operator.

6.4 Nonlinear stochastic transport equation↩︎

Finally, we analyse a first order equation for which our framework yields optimal convergence rates of the exponential Milstein scheme, even for nonlinear Nemytskii operators. Consider \[\label{eq:transport} \begin{cases} \mathrm{d}U &= \br[\big]{\nabla U + \phi(U)} \,\mathrm{d}t+ \psi(U)\,\mathrm{d}W_Q \quad \text{on }(0,T], \\ U_0 &= \xi \end{cases}\tag{45}\] in dimension \(d=1\) with nonlinearities \(\phi\colon\mathbb{R}\to\mathbb{R}\) and \(\psi\colon\mathbb{R}\to\mathbb{R}\), and a \(Q\)-Wiener process \(W_Q\). In this subsection, we abbreviate \(L^q\mathrel{\vcenter{:}}= L^q(\mathbb{R};\mathbb{R})\) for \(q\in[2,\infty)\) and \(H^\alpha \mathrel{\vcenter{:}}= H^\alpha(\mathbb{R};\mathbb{R})\) for \(\alpha\in [0,1]\). Since the best known Lipschitz estimate for \(\sigma \in (0,1)\) even with smooth Lipschitz \(\phi\in C_b^2\) is \[\|\phi(u)-\phi(v)\|_{H^\sigma} \lesssim \|u-v\|_{H^\sigma}+ \p[\big]{1+\|u\|_{H^\sigma}+\|v\|_{H^\sigma}}\|u-v\|_{L^\infty},\] cf. [51], which is nonlinear in \(u\) and \(v\), we restrict our considerations to the case \(\sigma=0\). We set \(X=H=L^2\), \(Y=H^\alpha\) for some \(\alpha\in (0,1]\), and consider the linear operator \(Au = -u'\) for \(u\in \mathop{\mathrm{D}}(A) \mathrel{\vcenter{:}}= H^1\). Note that \(\mathop{\mathrm{D}}(A^\alpha)=H^\alpha\) for \(\alpha\in[0,1]\) and \(-A\) generates the left shift semigroup [52], which is contractive on both \(L^2\) and \(H^\alpha\). Furthermore, \(F\) and \(G\) are given as Nemytskii operators \(F\colon L^2 \to L^2\), \(u \mapsto \phi \circ u\) and \(G\colon L^2 \to \mathscr{L}_2(L^2,L^2)\), \(G(u)=[h \mapsto M_{\psi \circ u} Q^{1/2}h]\).

To show convergence, we need to assume sufficient regularity of \(\phi\) and \(\psi\) and the covariance operator. We make this precise in the next theorem. Before, we recall useful estimates of composition and multiplication operators in the one-dimensional case that have been used in the previous subsections.

7 Numerical Simulations of the Stochastic Schrödinger Equation↩︎

To complement the theoretical findings presented above, we include simulations of three variants of the stochastic Schrödinger equation. First, the linear case is considered with a potential and multiplicative trace-class noise in Subsection 7.1. Second, a nonlocal nonlinearity as considered in Subsection 6.2 is simulated in Subsection 7.2. Lastly, we include simulations with a nonlinearity not covered by our results as an outlook in Subsection 7.3.

All three are considered on \(L^2 \mathrel{\vcenter{:}}= L^2(\mathbb{T};\mathbb{C})\) with the orthonormal Fourier basis functions \(e_\ell(x)=(2\pi)^{-1/2}{\mathrm{e}}^{{\mathrm{i}}\ell x}\) for \(x\in \mathbb{T}\), \(\ell\in\mathbb{Z}\) and up to the final time \(T=\frac{1}{2}\). We abbreviate \(H^s \mathrel{\vcenter{:}}= H^s(\mathbb{T};\mathbb{C})\) for \(s\in[0,\infty)\). For the spatial and the noise discretization we used a spectral Galerkin approach utilising the first \(2^{10}\) Fourier coefficients, i.e.we only considered the components with \(\ell\in\{-2^9+1,\dots,2^9\}\). As initial value, we chose \(\xi=\sum_{\ell\in\mathbb{Z}} \mu_\ell e_\ell\) with \(\mu_\ell=(1+\lvert \ell\rvert^{2.51})^{-1}\), \(\ell\in\mathbb{Z}\), being the Fourier coefficients of \(\xi\). Using the Fourier characterisation \(\|\xi\|_{H^s}^2 \simeq \sum_{\ell\in\mathbb{Z}} (1+|\ell|^2)^s|\mu_\ell|^2\) of the \(H^s\)-norm due to \(\norm{e_\ell}_{H^s}^2 \simeq (1+\abs{\ell}^{2})^s\), we see that \(\xi\in H^s\) for all \(s\in[0,\frac{5.02-1}{2})=[0,2.01)\), in particular, \(\xi\in H^2\).

Similarly, we define the covariance operator \(Q = \sum_{\ell\in\mathbb{Z}} \lambda_\ell \langle\cdot,e_\ell\rangle e_\ell\) with \(\lambda_\ell = (1+\abs{\ell}^{5.1})^{-1}\) for \(\ell\in\mathbb{Z}\). The exponent \(5.1\) determines the regularity of the covariance operator and thus of the driving noise. We have \(Q^{\frac{1}{2}}\in \mathscr{L}_2(L^2,H^s)\) for all \(s\in[0,\frac{5.1-1}{2})=[0,2.05)\), in particular, \(Q^{\frac{1}{2}}\in \mathscr{L}_2(L^2,H^2)\). Indeed, reasoning as for \(\xi\), we can bound \[\begin{align} \norm[\big]{Q^{\frac{1}{2}}}_{\mathscr{L}_2(L^2,H^s)}^2 &= \sum_{\ell\in\mathbb{Z}} \norm[\big]{Q^{\frac{1}{2}}e_\ell}_{H^s}^2 \leq \sum_{\ell\in\mathbb{Z}} \lambda_\ell \norm{e_\ell}_{H^s}^2 \lesssim \sum_{\ell\in\mathbb{Z}} \frac{1}{1+\abs{\ell}^{5.1-2s}}, \end{align}\] which is finite in the given range of \(s\).

All simulations were performed with linear multiplicative noise given by \(G(u)=-{\mathrm{i}}M_uQ^{\frac{1}{2}}\). This implies that \((G'G)(u)\) is symmetric as a bilinear Hilbert–Schmidt operator for all \(u\in L^2\) (also referred to as the commutativity condition). Thus, we can use the identity (cf.[15]) \[\begin{align} (G'G)(u_{i}) \Delta_2 W_{i+1}^K &= \frac{1}{2} G'(u_{i})\br[\big]{G(u_{i})\Delta W_{i+1}^K}\Delta W_{i+1}^K - \frac{1}{2}h \sum_{\ell=-K/2+1}^{K/2} G'(u_{i})\br[\big]{G(u_{i})e_\ell}e_\ell \end{align}\] to simulate the iterated integral terms required to implement Milstein schemes. Here, the superscript \(K\) denotes the truncation of the Wiener process to the first \(K\) Fourier coefficients. In general, i.e.if the operator does not fulfil the commutativity condition, it is currently not known how to simulate the iterated integrals exactly. Therefore, an approximation of the iterated integrals becomes necessary. This means that in practice different schemes may have substantially different computational costs per time step and therefore it is instructive to consider the convergence order not with respect to the step size but with respect to the computational cost of an algorithm. This is called the effective order of convergence and was introduced in [54]. In the case of a Milstein scheme in the non-commutative setting, estimates for the computational cost to simulate the iterated integrals with a desired precision are required. Here, we refer to [21] for a recent survey of suitable approximation algorithms in the finite-dimensional case as well as [55] for the best-known estimates in the infinite-dimensional case. For a detailed analysis of the effective order of convergence for Milstein-type schemes for SPDEs, we refer the interested reader to [22] and [56].

We denote by EXE, CN, and IE the exponential Euler, the Crank–Nicolson, and the implicit Euler schemes, respectively, and by EX-Milstein, CN-Milstein, and IE-Milstein (or simply EXM, CNM, IEM) the corresponding exponential or rational Milstein schemes. Approximations computed using EXM, CNM, and IEM are compared with the EXE and the rational schemes CN and IE simulated in [11] for five different step sizes \(h\in\{2^{-5},2^{-6},\dots,2^{-9}\}\). To estimate the pathwise uniform error for \(p=2\), we computed the root-mean-square error over \(100\) samples. As no closed form of the analytical solution is available for the equations considered, the errors are computed w.r.t.a reference solution, which was generated using the exponential Euler scheme with the significantly smaller step size \(h=2^{-16}\).

In Subsections 7.1 and 7.2, we make use of a centred bump function \(\chi\colon\mathbb{T}\to\mathbb{C}\) on the torus defined by \[\begin{align} \label{eq:defBump} f(x)&=\begin{cases} C\exp\p[\big]{\frac{1}{x^2-c^2}}, & \lvert x\rvert < c, \\ 0 & \text{else}, \end{cases} & \chi(x)&=\begin{cases} f(x), & x\in[0,\pi), \\ f(x-2\pi), & x\in[\pi,2\pi) \end{cases} \end{align}\tag{47}\] for \(c\in(0,\pi)\). Here, \(C=\exp(1/c^2)\) is a normalisation constant such that \(\chi(0)=f(0)=1\), and we chose \(c=\frac{\pi}{2}\) for the simulations.

Finally, we mention that the simulations were performed using Julia v1.11.7 [57] on a standard laptop. The code is publicly available under [58].

7.1 A linear example↩︎

We consider the linear Schrödinger equation 38 investigated in Section 6.1 on \(X=L^2\) with the potential \(V=\chi\), i.e.\(d=1\) and \(\sigma=0\). Since \(\xi\in H^2\), \(V\in H^2\), and \(Q^{\frac{1}{2}}\in \mathscr{L}_2(L^2,H^2)\), we can choose \(Y=H^2\) and Assumption 5[case:SoLin1D] is satisfied with \(\alpha=1\) for EXM. Hence, Theorem 21 shows that we should expect a convergence rate of one. For the rational Milstein schemes CNM and IEM, the expected convergence rates are \(2/3\) and \(1/2\), respectively, by Theorem 22. For all schemes without Milstein terms, we expect a rate of \(1/2\) due to [11].

Figure 1 (a) illustrates the numerical errors obtained for different step sizes. The resulting experimental rates of convergence are contrasted with the expected convergence rates discussed above in Table 2. The simulations confirm the higher rates of convergence of some Milstein schemes compared to the exponential Euler scheme and the rational schemes, whose rates are essentially limited by \(1/2\) also in the simulations. In contrast, higher numerical rates of convergence close to \(1\) and \(2/3\) are obtained for EXM and CNM, respectively, while the numerical rate for IEM is indeed close to \(1/2\), thus validating the theory.

7.2 A nonlinear example↩︎

Next, we present a simulation in the setting of Section 6.2 with a nonlocal drift \(F(u)=-{\mathrm{i}}\eta\ast\phi(u)\) that is not of Nemytskii-type but given by the convolution of a Nemytskii-type nonlinearity with a smooth kernel \(\eta\). As convolution kernel, we choose the centred bump function \(\eta=\chi\) from 47 and \(\phi(z)=z(1+\abs{z}^2)^{-1}\). They fulfil Assumption 7 with \(\alpha=\frac{2}{\ell_0}\), \(\ell_0=2,3,4\) for EXM, CNM, and IEM, respectively, and \(\beta=2\) such that Theorem 23 is applicable for the drift part. For simplicity, we again consider linear noise \(G(u)=-{\mathrm{i}}M_uQ^{\frac{1}{2}}\) with \(Q^{1/2}\) as in Subsection 7.1. To estimate the part of the error related to \(G\), we apply Theorems 21 and 22, also with \(\alpha=\frac{2}{\ell_0}\) and \(\beta=2\). An inspection of the proofs of Theorems  23, 21, and 22 shows that we can combine the different error estimates for the different \(F\)- and \(G\)-terms.

The numerical errors as well as the expected and numerical rates of convergence are presented in Figure 1 (b) and Table 2, respectively. Also for this nonlinear example, the advantage of the EXM and CNM can be clearly discerned. The numerical rates conform closely to the expected rates.

7.3 A further nonlinear example↩︎

As an outlook, we consider the nonlinear Nemytskii operator \(F(u)=-{\mathrm{i}}\phi\circ u\) associated with \(\phi(z) =z(1+\abs{z}^2)^{-1}\), again with linear multiplicative noise \(G(u) = -{\mathrm{i}}M_uQ^{\frac{1}{2}}\). On \(Y=H^2\), \(F\) is of quadratic growth, which is not covered by our setting. We could apply Theorem 23 on \(Y=H^1\), as \(F\) is of linear growth on this space, but this limits the expected convergence rates to \(1/2\) for Milstein schemes.

The numerical simulations, whose results can be found in Figure 1 (c) and Table 3, indicate that, nonetheless, the Milstein schemes EXM and CNM converge at rates higher than \(1/2\). It is an interesting open question whether error estimates for Milstein schemes can also be obtained under weaker conditions on \(F\) and \(G\) such as polynomial rather than linear growth assumptions.

References↩︎

---

1. The second author gratefully acknowledges support by the Alexander von Humboldt foundation through a Feodor Lynen Research Fellowship. The authors are also grateful for support from the VICI subsidy VI.C.212.027 of the Netherlands Organisation for Scientific Research (NWO)↩︎