Drift-preserving numerical integrators for stochastic Hamiltonian systems

The paper deals with numerical discretizations of separable nonlinear Hamiltonian systems with additive noise. For such problems, the expected value of the total energy, along the exact solution, drifts linearly with time. We present and analyze a time integrator having the same property for all times. Furthermore, strong and weak convergence of the numerical scheme along with efficient multilevel Monte Carlo estimators are studied. Finally, extensive numerical experiments illustrate the performance of the proposed numerical scheme.


Introduction
Hamiltonian systems are universally used as mathematical models to describe the dynamical evolution of physical systems in science and engineering.If the Hamiltonian is not explicitly time dependent, then its value is constant along exact trajectories of the problem.This constant equals the total energy of the system.
This research has naturally come to the realm of stochastic differential equations (SDEs).Indeed, many publications on energy-preserving numerical integration of stochastic (canonical) Hamiltonian systems have lately appeared.Without being too exhaustive, we mention [8,12,19,28] as well as the works [27,42,9] on invariant preserving schemes.Observe, that such extensions describe Hamiltonian motions perturbed by a multiplicative white noise in the sense of Stratonovich.In some sense, this respects the geometric structure of the phase space.Hence one also has conservation of the Hamiltonian (along any realizations).The derivation of energy-preserving schemes in the Stratonovich setting follows without too much effort from already existing deterministic energy-preserving schemes.This is due to the fact that most of the geometrical features underlying deterministic Hamiltonian mechanics are preserved in the Stratonovich setting.This is not the case in the Itô framework considered here.
An alternative to the above Stratonovich setting is to add a random term to the deterministic Hamiltonian in an additive way.One can then show that the expected value of the Hamiltonian (along the exact trajectories) now drifts linearly with time, leading to a so called trace formula, see for instance the work [40] in the case of a linear stochastic oscillator.To the best of our knowledge, drift-preserving numerical schemes for such a problem have only been theoretically studied in the case of quadratic Hamiltonians [6,5,10,17,15,26,40,39] and in the case of linear stochastic partial differential equations [1,2,11,14,38].For instance, recent contributions in structure-preserving numerics for stochastic problems have addressed conservation properties of stochastic Runge-Kutta methods applied to stochastic Hamiltonian problems of Itô type [5,6].In particular, proper stochastic perturbations of symplectic Runge-Kutta methods have been investigated as drift-preserving schemes.However, these methods seem particularly effective only for linear problems, while, in the nonlinear case, the drift is not accurately preserved in time.Moreover, in the case of additive noise, the more the amplitude of the stochastic term increases, the more the accuracy of the drift preservation deteriorates (see, for instance, Table 1 in [5]).Hence, this evidence reveals a gap in the existing literature that requires ad hoc numerical methods, effective in preserving the drift in the expected Hamiltonian also for the nonlinear case.Closing this gap is the main purpose of this article.
In the present publication, we develop and analyze drift-preserving schemes for stochastic separable Hamiltonian systems, not necessarily quadratic, driven by an additive Brownian motion.In particular, we propose a novel numerical scheme that exactly satisfies the trace formula for the expected value of the Hamiltonian for all times (Section 2).In addition, under general assumptions, we show meansquare and weak convergence of the newly proposed numerical scheme in Sections 3 and 4 to give a complete picture of its properties.We show that both errors converge with order 1 and that the weak order is in general not twice the strong order in this specific case.On top of that, in Section 4 we introduce Monte Carlo and multilevel Monte Carlo estimators for the given numerical scheme and derive their convergence properties along with their computational costs.The main properties of the proposed time integrators are then numerically illustrated in Section 5.

Presentation of the drift-preserving scheme
Let T ą 0, d be a fixed positive integer, pΩ, F , Pq be a probability space with a normal filtration pF t q tPr0,T s , and let W : r0, T s ˆΩ Ñ R d be a standard pF t q tPr0,T s -Brownian motion with continuous sample paths on pΩ, F , Pq.
For a positive integer m and a smooth potential V : R m Ñ R, let us consider the separable Hamiltonian Consider now the following corresponding stochastic Hamiltonian system with additive noise dqptq " ∇ p Hppptq, qptqq dt " pptq dt dpptq " ´∇q Hppptq, qptqq dt `Σ dW ptq " ´V 1 pqptqq dt `Σ dW ptq.
Here, the constant matrix Σ P R mˆd has entries denoted by σ ij .In addition, we assume that the initial values pp 0 , q 0 q of the above SDE have finite energy in expectation, i. e., ErHpp 0 , q 0 qs ă `8.This setting covers, for instance, the following examples: the Hamiltonian considered in [5] (where the matrix Σ is diagonal), the linear stochastic oscillator from [40], various Hamiltonians studied in [34,Chap. 4].
Remark 1.In general, most results from the present publication could be extended to the case of a separable Hamiltonian system with additive martingale Lévy noise.Extension to the case of a multiplicative (Itô) noise would lead to a trace formula for the energy that depends on q t (when the matrix Σ depends only on q t ).This would correspond to the extra term appearing when converting between Stratonovich and Itô stochastic integrals.Using Itô's lemma, one gets the trace formula for the energy of the above problem (see for instance [5]): (2) E rHppptq, qptqqs " E rHpp 0 , q 0 qs `1 2 Tr `ΣJ Σ ˘t, for all times t ą 0. We now want to design a numerical scheme having the same property.The idea we aim to present is inspired by the deterministic energy-preserving schemes from the literature (see [13,22,37] and references therein), able to provide the exact conservation of the energy of a given mechanical system.
Energy-preserving methods represent a follow-up of the classical results on geometric numerical integration relying on the employ of symplectic Runge-Kutta methods, projection methods and nearly preserving integrators (see the comprehensive monograph [23] and references therein).In the general setting of deterministic B-series methods, a general proof for the existence of energy-preserving methods was given in [18] and practical examples of methods were first developed in [37].
We propose the following time integrator for problem (1), which is shown in Theorem 3 to be a drift-preserving integrator for all times, Ψ n`1 " p n `Σ∆W n ´h 2 where h ą 0 is the stepsize of the numerical scheme and ∆W n are Wiener increments.Observe that the deterministic integral above needs to be computed exactly (as in the deterministic setting).This is not a problem for polynomial potentials V for instance.We also observe that, as it happens for deterministic energy-preserving methods, the numerical scheme (3) only requires the evaluation to the vector field of problem (1) in selected points, without requiring additional projection steps for the exact conservation of the trace formula (2), that would inflate the computational cost of the procedure.We would like to remark that the proposed drift-preserving scheme is not the only possibility to get a time integrator that exactly satisfy a trace formula for the energy for all times.Another possibility could be to use a splitting strategy.This idea is currently under investigation in [16], see also Section 5 below.
Remark 2. Since the proposed numerical scheme is implicit, we present two possible ways of showing its solvability.
To show the solvability of the drift-preserving scheme (3), we use the fixed point theorem.Assuming that V 1 pxq is globally Lipschitz continuous, we define the function The solvability of the numerical integrator (3) is thus equivalent to showing that the function F is a contractive mapping.Since F pψ 1 q ´F pψ 2 q " ´h 2 Therefore, there exists an h ˚" b 4 L such that for all h ă h ˚, F is a contractive mapping.If V P C 2 , we could also use the implicit function theorem instead.Indeed, let us define Gpψ, h, ∆W n q " ψ ´pn ´Σ∆W n `h 2 Then the solvability of ψ for sufficiently small h is equivalent in showing that ˇˇ∇ψGpψ, h, ∆W n q ˇˇh "0 ˇˇ‰ 0.
In fact By the implicit function theorem, we can thus conclude that for sufficiently small h, ψ is solvable.
We next show that the numerical scheme (3) satisfies, for all times, the same trace formula for the energy as the exact solution to the SDE (1), i. e., that it is a drift-preserving integrator.
Theorem 3. Assume that V P C 1 .Consider the numerical discretization of the stochastic Hamiltonian system with additive noise (1) by the drift-preserving numerical scheme (3).Then the expected energy satisfies the following trace formula (4) E rHpp n , q n qs " E rHpp 0 , q 0 qs `1 2 Tr `ΣJ Σ ˘tn for all discrete times t n " nh, where n P N. Observe that h is an arbitrary step size which is sufficiently small for the implicit numerical scheme to be solvable.
Proof.Before we start, let us for convenience introduce the notation Int " Using the definitions of the Hamiltonian, of the numerical scheme and a Taylor expansion for the potential V yields By the definition of the numerical scheme and properties of the Wiener increments, one obtains further Due to cancellations and thanks to properties of the Wiener increments, the expression above simplifies to Tr `ΣJ Σ ˘h.
A recursion now completes the proof.
The above result can also be seen as a longtime weak convergence result and provides a longtime stability result (in a certain sense) for the drift-preserving numerical scheme (3).Already in the case of a quadratic Hamiltonian, such trace formulas are not valid for classical numerical schemes for SDEs, as seen in [40] for instance and in the numerical experiments provided in Section 5.

Mean-square convergence analysis
In this section, we assume that V 1 is a globally Lipschitz continuous function and consider a compact time interval r0, T s.Then from the standard analysis of SDEs, see for instance [30, Theorems 4.5.3 and 4.5.4],we known that for all t, s P r0, T s one has where the constant C may depend on the coefficients of the SDE (1) and its initial values.
We now state the result on the mean-square convergence of the drift-preserving integrator (3).
Theorem 4. Consider the stochastic Hamiltonian problem with additive noise (1) on a fixed time interval r0, T s and the drift-preserving integrator (3).Assume that the potential V 1 is globally Lipschitz continuous.Then, there exists h ˚ą 0 such that for all 0 ă h ď h ˚, the numerical scheme converges with order 1 in mean-square, i. e., where the constant C may depend on the initial data, the end time T , and the coefficients of (1), but is independent of h and n for n " 1, . . ., N .Here, N is an integer such that N h " T .
Proof.We base the proof of this result on Milstein's fundamental convergence theorem, see [34, Theorem 1.1] and first consider one step in the approximation of the numerical scheme (3), starting from the point pp, qq at time t " 0 to pp, qq at time t " h: Let us introduce a similar notation for the exact solution of (1) starting from the point pp, qq at time t " 0 to ppphq, qphqq at time t " h: Finally, we define the local errors e q loc " qphq ´q and e p loc " pphq ´p.In order to use Milstein's fundamental convergence theorem in our situation, one has to show that It turns out that, due to the particular form of (1), we actually achieve better rates for the second term in (8).Using ( 5) and ( 6), we obtain e q loc " hpψphq ´ψq " Hence by the property of Itô integral, we have where |¨| F denotes the Frobenius norm of a matrix.For the term e q,2 loc , using Hölder's inequality and the linear growth condition on V 1 , we obtain Here in the last step, we used the classical bounds from the beginning of this section, The last inequality comes from the fact that, from (5), we have This implies that there exists an h ˚ą 0 such that for all 0 ă h ď h ˚, inequality (9) is verified.
For the estimate of e p loc , we get from ( 5) and ( 6) that Then using that V 1 is globally Lipschitz continuous, we obtain In the last step, we use the fact that Therefore we obtain the bounds |E re p loc s| ď An application of Milstein's fundamental convergence theorem completes the proof.

Weak convergence analysis
The weak error analysis for the Hamiltonian functional in the given situation is a direct consequence of the preceding sections.For the considered SDE, the main quantity of interest is given by a specific test function, i. e., we are interested in computing E rHppptq, qptqqs , for t ą 0. Equation (2) and Theorem 3 imply for all discrete times t n " nh that |E rHpppt n q, qpt n qqs ´E rHpp n , q n qs| " 0.
For fixed h ą 0 let us further extend the discrete processes pp n , n " 0, . . ., N q and pq n , n " 0, . . ., N q to an adapted process for all t ą 0 by setting p h ptq " p n , q h ptq " q n , for t P rt n , t n`1 q.Then p h and q h are piecewise continuous adapted processes which satisfy Tr `ΣJ Σ ˘h.
In conclusion using the trace formulas ( 2) and ( 4) we have just shown the following corollary.
Corollary 5. Assume that V P C 1 .Consider the numerical discretization of the stochastic Hamiltonian system with additive noise (1) by the drift-preserving numerical scheme (3), then the weak error in the Hamiltonian is bounded by Tr `ΣJ Σ ˘h, for all t ą 0.More specifically, for all t " t n , with n " 0, . . ., N , the error is zero and the approximation scheme is preserving the quantity of interest.
Let us next consider the weak error for other test functions than H and introduce for convenience the space of square integrable random variables L 2 pΩq " v : Ω Ñ R, v strongly measurable, }v} L 2 pΩq ă `8( , where }v} L 2 pΩq " Er|v| 2 s 1{2 .In the following result on weak convergence, we obtain the same rate as for strong convergence in Section 3 and will see in the examples in Section 5 that this seems optimal.This is caused by our combination of smoothness assumptions and the additive noise.Theorem 6. Assume that V 1 is globally Lipschitz continuous.Then the drift preserving scheme (3) converges weakly with order 1 to the solution of (1) on any finite time interval r0, T s.More specifically, for all differentiable test functions f : R 2d Ñ R with f 1 of polynomial growth there exist C ą 0 and h ˚ą 0 such that for all 0 ă h ď h ˚and corresponding n " 1, . . ., N with hN " T |E rf pppt n q, qpt n qqs ´E rf pp n , q n qs| ď Ch.
Proof.If f is globally Lipschitz continuous, the claim follows directly from Theorem 4.
Else let us use the abbreviations Xpt n q " pppt n q, qpt n qq and X n " pp n , q n q.Then the mean value theorem and the Cauchy-Schwarz inequality imply that |E rf pXpt n qqs ´E rf pX n qs| " ˇˇˇE "ż 1 0 f 1 pX n `spXpt n q ´Xn qq ds ¨pXpt n q ´Xn q ˇˇˇď ż 1 0 }f 1 psXpt n q `p1 ´sqX n q} L 2 pΩq }Xpt n q ´Xn } L 2 pΩq ds.
While the second term in the product is bounded by }Xpt n q ´Xn } L 2 pΩq ď Ch by Theorem 4, we use the polynomial growth assumption on f 1 for the first term to obtain The solution to (1) has bounded 2m-th moment by Theorem 4.5.4 in [30].The corresponding boundedness of the numerical solution follows by Lemma 2.2.2 in [34] for h sufficiently small depending on the polynomial growth constant on f 1 .Choosing h ˚to be the minimum of this restriction and the one from Theorem 4, we conclude the proof.
Remark 7. As we will see in Section 5, there exist combinations of equations and test functions for which the weak order of convergence of the drift-preserving scheme is in fact 2.More specifically, we observe this behaviour in simulations for the identity as test function.ErV 1 pq n `shΨ n`1 qs ds " 2ErΨ n`1 s ´Erp n s.
We can thus conclude that the weak order of convergence in the first moment of the drift-preserving scheme is 2.
In general, we will not be able to compute E rf pp n , q n qs analytically but have to approximate the expectation.This can for example be done by a standard Monte Carlo method.Following closely [32], we denote by the Monte Carlo estimator of a random variable Y based on M independent, identically distributed random variables Ŷ i .It is well-known that E M rY s converges P-almost surely to ErY s for M Ñ `8, by the strong law of large numbers.Furthermore, it converges in mean square if Y is square integrable and satisfies }ErY s ´EM rY s} By Corollary 5 for the Hamiltonian and the splitting of the error into the weak error in Theorem 6 and the Monte Carlo error we therefore obtain the following lemma.
Lemma 8. Assume that V P C 1 .Consider the numerical discretization of the stochastic Hamiltonian system with additive noise (1) by the drift-preserving numerical scheme (3), then the single level Monte Carlo error is bounded by }E rHpppt n q, qpt n qqs ´EM rHpp n , q n qs} L 2 pΩq " M ´1{2 VarrHpp n , q n qs 1{2 ď M ´1{2 }Hpp n , q n q} L 2 pΩq , for all t n " nh.For any test function f under the assumptions of Theorem 6, this result generalizes on any finite time interval to }E rf pppt n q, qpt n qqs ´EM rf pp n , q n qs} L 2 pΩq " Ch `M ´1{2 Varrf pppt n q, qpt n qqs 1{2 ď Ch `M ´1{2 }f pppt n q, qpt n qq} L 2 pΩq with computational work for balanced error contributions For completeness of presentation of the numerical analysis of the drift-preserving scheme (3), although not necessarily relevant for Hamiltonian SDEs, let us look at the savings in computational costs of the multilevel Monte Carlo estimator [25,20].Therefore we assume that pY ℓ , ℓ P N 0 q is a sequence of approximations of Y .For given L P N 0 , it holds that and due to the linearity of the expectation that ErY ℓ ´Yℓ´1 s.
A possible way to approximate ErY L s is to approximate ErY ℓ ´Yℓ´1 s with the corresponding Monte Carlo estimator E M ℓ rY ℓ ´Yℓ´1 s with a number of independent samples M ℓ depending on the level ℓ.
We set and call E L rY L s the multilevel Monte Carlo estimator of ErY L s.
We consider for the remainder of this section for convenience that the approximation scheme is based on a sequence of equidistant nested time discretizations τ " pτ ℓ , ℓ P N 0 q given by τ ℓ " pt ℓ n " T ¨n ¨2´ℓ , n " 0, . . ., 2 ℓ q with sequence of step sizes ph ℓ " T ¨2´ℓ , ℓ P N 0 q.
Let us denote by ppp ℓ , q ℓ q, ℓ P N 0 q the sequence of approximation schemes with respect to the sequence of time discretizations τ .We observe that Theorem 6 implies in the setting of [32, Theorem 1] that a ℓ " 2 ´ℓ, η " 2 ´1 by Theorem 4, and κ " 1. Plugging in these values we obtain the following corollary.
Corollary 9. Consider the numerical discretizations (3) of the stochastic Hamiltonian system (1) satisfying the assumptions in Theorem 4 and Theorem 6.Then for any differentiable test function f with derivative of polynomial growth, the multilevel Monte Carlo estimator on level L ą 0 satisfies for any ǫ ą 0 at the final time T }Erf pppT q, qpT qqs ´EL rf pp L 2 L , q L 2 L qs} L 2 pΩq ď pC 1 `C3 `C2 ζp1 `ǫqq h L with sample sizes given by M ℓ " r2 2pL´ℓ{2q ℓ 2p1`ǫq s, ℓ " 1, . . ., L, ǫ ą 0, and M 0 " r2 2L s, where r¨s denotes the rounding to the next larger integer and ζ the Riemann zeta function.The resulting computational work is bounded by W L " Oph ´2 L L 3`2ǫ q.
In conclusion we have seen that our drift-preserving scheme converges in general weakly with order 1, i. e., the same order as in mean square in Section 3. Approximating quantities of interest with a standard Monte Carlo error with accuracy h requires computational work h ´3.This can be reduced to essentially h ´2 when a multilevel Monte Carlo estimator is applied instead.

Numerical experiments
This section presents various numerical experiments in order to illustrate the main properties of the drift-preserving scheme (3), denoted by DP below.In some numerical experiments, we will compare this numerical scheme with classical ones for SDEs such as the Euler-Maruyama scheme (denoted by EM) and the backward Euler-Maruyama scheme (denoted by BEM).

5.1.
The linear stochastic oscillator.We first consider the SDE (1) with the following Hamiltonian and with Σ " 1 and W scalars.We take the initial values pp 0 , q 0 q " p0, 1q.
In this situation the integral in the drift-preserving scheme (3) can be computed exactly, resulting in the following time integrator Using the stepsize h " 5{2 4 , resp.h " 100{2 8 , and the time interval r0, 5s, resp.r0, 150s, we compute the expected values of the energy Hpp, qq along the numerical solutions.For this problem, we also use the stochastic trigonometric method from [10], denoted by STM, which is know to preserve the trace formula for the energy for this problem.For comparison, the following splitting strategies are also used ‚ composition of the (deterministic) symplectic Euler scheme for the Hamiltonian part with an analytical integration of the Brownian motion (this scheme is denoted by SYMP); ‚ a splitting based on the decomposition dqptq " pptq dt, dpptq " Σ dW ptq and dqptq " 0 dt, dpptq " ´V 1 pqptqq dt (this time integrator is denoted by SPLIT).
The expected values are approximated by computing averages over M " 10 6 samples.The results are presented in Figure 1, where one can clearly observe the excellent behaviour of the drift-preserving scheme as stated in Theorem 3. Observe that it can be shown that the expected value of the Hamiltonian along the Euler-Maruyama scheme drifts exponentially with time.Furthermore, the growth rate of this quantity along the backward Euler-Maruyama scheme is slower than the growth rate of the exact solution to the considered SDE, see [40] for details.These growth rates are qualitatively different from the linear growth rate in the expected value of the Hamiltonian of the exact solution (2), of the STM from [10], and of the drift-preserving scheme (4).Although not having the correct growth rates, the splitting schemes behave much better than the classical Euler-Maruyama schemes.Further splitting strategies are under investigation in [16].We next illustrate numerically the strong convergence rate of the drift-preserving scheme stated in Theorem 4. Using the same parameters as above, we discretize the SDE on the interval r0, 1s using step sizes ranging from 2 ´5 to 2 ´10 .The loglog plots of the errors are presented in Figure 2, where mean-square convergence of order 1 for the proposed integrator is observed.The reference solution is computed with the stochastic trigonometric method using h ref " 2 ´12 .The expected values are approximated by computing averages over M " 10 5 samples.
To conclude this subsection, we numerically illustrate the weak rates of convergence of the above time integrators.In order to avoid Monte Carlo approximations, we focus on weak errors in the first and second moments only, where all the expectations can be computed exactly.We use the same parameters as above except for Σ " 0.1, T " 1, and step sizes ranging from 2 ´4 to 2 ´16 .The results are presented in Figure 3, where one can observe weak order 2 in the first moment and weak order 1 in the second moment for the drift-preserving scheme.This is in accordance with the results from the preceding section.
For this problem, the proposed time integrator reads Ψ n`1 " p n `Σ∆W n ´h 2 ˆcospq n q ´cospq n `hΨ n`1 q hΨ n`1 ˙, q n`1 " q n `hΨ n`1 , p n`1 " p n `Σ∆W n ´h ˆcospq n q ´cospq n `hΨ n`1 q hΨ n`1

˙.
Using the stepsize h " 5{2 8 , resp.h " 10{2 10 , and the time interval r0, 5s, resp.r0, 10s, we compute the expected values of the energy Hpp, qq along the numerical solutions.Newton's iterations are used to solve the nonlinear systems in the drift-preserving scheme (3) as well as in the BEM scheme.The expected values are approximated by computing averages over M " 10 6 samples, resp.M " 10 5 samples.The results are presented in Figure 4, where one can clearly observe the perfect behaviour of the drift-preserving scheme as stated in Theorem 3 as well as the excellent behaviour of the splitting strategies.We also illustrate numerically the strong convergence rate of the drift-preserving scheme stated in Theorem 4. For this, we discretize the stochastic mathematical pendulum with Σ " 0.1 on the interval r0, 0.5s using step sizes ranging from 2 ´5 to 2 ´10 .The loglog plots of the errors are presented in Figure 5, where mean-square convergence of order 1 for the proposed integrator is observed.The reference solution is computed with the drift-preserving scheme using h ref " 2 ´12 .The expected values are approximated by computing averages over M " 10 5 samples.
We conclude this subsection by reporting numerical experiments illustrating the weak rates of convergence of the above time integrators.In order to reduce the Monte Carlo error and thus produce nice plots, we had to multiply the nonlinearity with a small coefficient of 0.2 and considered the interval r0, 1s.The other parameters are the same as above.The step sizes range from 2 ´1 to 2 ´6.The expected values are approximated by computing averages over M " 10 8 samples.The results are presented in Figure 6, where one can observe weak order 2, resp. 1, of convergence for the driftpreserving scheme in the first, resp.second, moment in the variable q.This confirms the theoretical results from the previous section.5.3.Double well potential.We consider the Hamiltonian with double well potential from [5].The SDE (1) is thus given by the Hamiltonian Hpp, qq " 1 2 p 2 `1 4 q 4 ´1 2 q 2 Figure 6.Weak rates of convergence for the drift-preserving scheme (DP), the backward Euler-Maruyama scheme (BEM), and the Euler-Maruyama scheme (EM) when applied to the stochastic mathematical pendulum.Errors in the first moment of q (left) and in the second moment of q (right).
and with Σ " 0.5 and W scalars.We take the initial values pp 0 , q 0 q " p ? 2, ?2q.When applied to the Hamiltonian with double well potential, the time integrator (3) takes the form Ψ n`1 " p n `Σ∆W n ´h 2 Using 2 16 stepsizes of the drift-preserving scheme (3) on the time interval r0, 50s (using fixed-point iterations for solving the implicit systems), and approximating the expected value with M " 10 5 samples, we obtain the result displayed in Figure 7.This again numerically confirms the long time behaviour of the drift-preserving scheme stated in Theorem 3 and shows its superiority compared to the numerical schemes from [5], see the numerical results in [5, Table 1].5.4.Hénon-Heiles problem with two additive noises.Finally, we consider the Hénon-Heiles problem with two additive noises from [5,24].This SDE is given by the Hamiltonian Hpp, qq " 1 2 and with Σ " diagpσ 1 , σ 2 q and W " pW 1 , W 2 q J in (1).We take the following parameters σ 1 " σ 2 " 0.2, α " 1{16, and initial values p 0 " p1, 1q and q 0 " p ? 3, 1q.For this system of SDEs, the drift-preserving scheme reads
Using 2 11 stepsizes of the drift-preserving scheme (3) on the time interval r0, 50s (using fixed-point iterations for solving the implicit systems), and approximating the expected values with M " 10 5 samples, we obtain the result displayed in Figure 8.This figure again clearly illustrates the excellent long time behaviour of the proposed numerical scheme as stated in the above theorem.

2 Figure 5 .
Figure 5. Strong rates of convergence for the drift-preserving scheme (DP), the backward Euler-Maruyama scheme (BEM), and the Euler-Maruyama scheme (EM), when applied to the stochastic mathematical pendulum.

Figure 7 .
Figure 7. Numerical trace formula for the stochastic Hamiltonian with double well potential on r0, 50s.
The above numerical scheme is nothing else than the expectation of the drift-preserving scheme (3)ErV 1 pq n `shΨ n`1 qs ds Erq n`1 s " Erq n s `hErΨ n`1 s n s `Erp n`1 sq Erp n`1 s " Erp n s ´h ż 1 0ErV 1 pp1 ´sqq n `sq n`1 qs ds.