Optimal explicit stabilized integrator of weak order one for stiff and ergodic stochastic differential equations
Assyr Abdulle Ibrahim Almuslimani Gilles Vilmart

	Statistical mechanics of specular reflections from fluctuating membranes and interfaces
    David R. Nelson
    December 30, 2023
=======================================================================================

[1]Mathematics Section, École Polytechnique Fédérale de Lausanne,
		Station 8, CH-1015 Lausanne, Switzerland, Assyr.Abdulle@epfl.ch
	[2]Université de Genève, Section de mathématiques, 2-4 rue du Lièvre, CP 64, CH-1211 Genève 4, Switzerland, Ibrahim.Almuslimani@unige.ch, Gilles.Vilmart@unige.ch
	
A new explicit stabilized scheme of weak order one for stiff  and ergodic stochastic differential equations (SDEs)  is introduced. In the absence of noise, the new method coincides with the classical deterministic stabilized scheme (or Chebyshev method) for diffusion dominated advection-diffusion problems and it inherits its optimal stability domain size that grows quadratically with the number of internal stages of the method.
For mean-square stable stiff stochastic problems, the scheme has an  optimal extended mean-square stability domain that grows at the same quadratic rate as the deterministic stability domain size  in contrast to known existing methods for stiff SDEs [A. Abdulle and T. Li. Commun. Math. Sci., 6(4), 2008, A. Abdulle, G. Vilmart, and K. C. Zygalakis, SIAM J. Sci. Comput., 35(4), 2013].
Combined with postprocessing techniques, the new methods
achieve a convergence rate of order two for sampling the invariant measure of a class of ergodic SDEs, achieving a stabilized version of the non-Markovian scheme introduced in [B. Leimkuhler, C. Matthews, and M. V. Tretyakov, Proc. R. Soc. A, 470, 2014].


		
		
		
		Keywords: 
		explicit stochastic methods, stabilized methods, postprocessor, invariant measure, ergodicity, orthogonal Runge-Kutta Chebyshev, SK-ROCK, PSK-ROCK.
		AMS subject classification (2010): 
		 65C30, 60H35, 65L20, 37M25
	
	
	
	
	
	
	
	

§ INTRODUCTION

	
	We consider Itô systems of stochastic differential equations of the form
	
    dX(t) = f(X(t)) dt + ∑_r=1^mg^r(X(t)) dW_r(t),     X(0)=X_0

	where X(t) is a stochastic process with values in ^d, f:^d→^d is the drift term,
	g^r:^d→^d, r=1,…, m are the diffusion terms, and W_r(t), r=1,…,m, are independent 
	one-dimensional Weiner processes fulfilling the usual assumptions.
	We assume that the drift and diffusion functions are smooth enough and Lipschitz continuous to ensure the existence and uniqueness of a solution of (<ref>)
	on a given time interval (0,T).
	We consider autonomous problems to simplify the presentation, but we emphasise that the scheme can also be extended to non-autonomous SDEs. A one step numerical integrator 
for the approximation of (<ref>) at time t=nh is a discrete dynamical system of the form

    X_n+1 = Ψ(X_n,h,ξ_n)

where h denotes the stepsize and ξ_n are independent random vectors. 
Analogously to the deterministic case, standard explicit numerical schemes for stiff stochastic problems, such as the simplest Euler-Maruyama method defined as
	
    X_n+1 = X_n + hf(X_n) + ∑_r=1^mg^r(X_n) Δ W_n,r,     X(0)=X_0,

	where Δ W_n,r=W_r(t_n+1)-W_r(t_n) are the Brownian increments,
face a severe timestep restriction <cit.>, and one can use an implicit or semi-implicit scheme with favourable
	stability properties. In particular, it is shown in <cit.> that the implicit θ-method of weak order one is mean-square A-stable if and only if θ≥ 1/2, while
	weak order two mean-square A-stable are constructed in <cit.>.
	An alternative approach is to consider explicit stabilized schemes with extended stability domains, as proposed in 
	<cit.>.
	In <cit.> the deterministic Chebyshev method is extended to the context of mean-square stiff stochastic differential equations with Itô noise, while the Stratonovitch noise case
	is treated in <cit.>.
	In place of a standard small damping, the main idea in <cit.> is to use a large damping parameter η optimized for each number s of stages to stabilize
	the noise term. This yields a family of Runge-Kutta type schemes with extended stability domain with size L_s≃ 0.33 s^2.
	This stability domain size was improved to L_s≃ 0.42 s^2 in <cit.> where a family of  weak second order stabilized schemes (and strong order one under suitable assumptions) is constructed based on the deterministic orthogonal Runge-Kutta-Chebyshev method of order 2 (ROCK2) <cit.>.

For ergodic SDEs, 	i.e., when (<ref>) has a unique invariant measure 
μ satisfying 
for each test function ϕ and for any deterministic initial condition X_0=x,

    lim_T →∞1/T∫_0^Tϕ(X(s))ds =∫_ϕ(y)dμ(y),     ,

one is interested in approximating numerically the long-time dynamics and to design numerical scheme with a unique invariant measure such that 

    | lim_N →∞1/N+1∑_n=0^Nϕ(X_n)-∫_ϕ(y)dμ(y)|≤ C h^r,
 
where C is independent of h small enough and X_0. In such a situation, we say that the numerical scheme has order r with respect to the invariant measure.
For instance, the Euler-Maruyama method has order 1 with respect to the invariant measure. In <cit.> the following non-Markovian scheme with the same cost as the Euler-Maruyama method was proposed for Brownian dynamics, i.e where the vector field is a gradient f(x)=-∇ V(x) and the noise is additive (g(x)=σ), 
	
    X_n+1 = X_n + hf(X_n) + σΔ W_n,j+Δ W_n+1,j/2,     X(0)=X_0,

	and it was shown in <cit.>
that (<ref>) has order 2 with respect to the invariant measure for Brownian dynamics.
However, the admissible stepsizes for such an explicit method to be stable may face a severe restriction and alternatively to switching to drift-implicit methods, one may ask if a stabilized version of such an attractive non-Markovian scheme exists.

In this paper we introduce a new family of explicit stabilized schemes with optimal mean-square stability domain of size L_s=Cs^2, where C≥ 2-4/3η and η≥ 0 is a small  parameter. We emphasize that in the deterministic case, L_s=2s^2 is the largest, i.e. optimal, stability domain along the negative real axis for an explicit s-stage Runge-Kutta method <cit.>. 
We note that the Chebyshev method (<ref>) (with η=0) realizes such an optimal stability domain. The new schemes have strong order 1/2 and weak order 1.
The main ingredient for the design of the new schemes is to consider second kind Chebyshev polynomials, in addition to the usual first kind Chebyshev polynomials involved in the deterministic Chebychev method and stochastic extensions <cit.>.
For stiff  stochastic problems, the stability domain sizes are close to the optimal value 2s^2
and in the deterministic setting the method coincide with the optimal first order  explicit stabilized method. Thus these methods are more efficient than previously introduced stochastic stabilized methods <cit.>.
For ergodic dynamical systems, in the context of the ergodic Brownian dynamics, the new family of explicit stabilized schemes allows for a postprocessing <cit.> (see also <cit.> in the context of Runge-Kutta methods) to achieve
order two of accuracy for sampling the invariant measure. In this context, our new methods can be seen as a stabilized version of the non-Markovian scheme (<ref>) introduced in <cit.>.

	This paper is organized as follows.
	In Section <ref>, we introduce the new family of schemes with optimal stability domain and we recall the main tools for the study of stiff integrators in the mean-square sense.
	We then analyze its mean-square stability properties (Section <ref>), and convergence properties (Section <ref>).
	In Section <ref>, using a postprocessor we present a modification  with negligible overcost that yields order two of accuracy for the invariant measure of a class of ergodic overdamped Langevin equation.
	Finally, Section <ref> is dedicated to the numerical experiments that confirm our theoretical analysis and illustrate the efficiency of the new schemes.
	
	
	
	

§ NEW SECOND KIND CHEBYSHEV METHODS

	
In this section we introduce our new stabilized stochastic method. We first briefly recall the concept of stabilized methods.
In the context of ordinary differential equations (ODEs),

    dX(t)/dt = f(X(t)),      X(0)=X_0,

and the Euler method X_1=X_0+hf(X(0)), a stabilization procedure  based on recurrence formula has been introduced
in <cit.>. Its construction relies on Chebyshev polynomials (hence the alternative name “Chebyshev methods"),
T_s(cos x) =cos (s x)  and it is based on the explicit s-stage Runge-Kutta method

    K_0   =    X_0,   
    K_1 = X_0+h μ_1 f(K_0), 
    
    K_i   =   μ_ihf(K_i-1)+ν_iK_i-1+κ_iK_i-2,   j=2,…,s,
    
    X_1   =    K_s,

where 

    ω_0=1+η/s^2,  ω_1=T_s(ω_0)/T'_s(ω_0),  μ_1=ω_1/ω_0,
 
and for all i=2,…,s,
	
    μ_i=2ω_1T_i-1(ω_0)/T_i(ω_0),  ν_i=2ω_0T_i-1(ω_0)/T_i(ω_0),  κ_i=-T_i-2(ω_0)/T_i(ω_0)=1-ν_i.

One can easily check that the (family of) methods (<ref>) has the same first order accuracy as the Euler method (recovered for s=1). In addition, the scheme
(<ref>) has a low memory requirement (only two stages  should be stored when applying the recurrence formula) and it has a good internal stability with respect to round-off errors <cit.>.
The attractive feature of such a scheme comes from its stability behavior. Indeed, the method (<ref>) applied to the 
linear test problem dX(t)/dt=λ X(t) yields, using the recurrence relation 
		
    T_j(p) = 2p T_j-1(p) - T_j-2(p),

		where T_0(p)= 1,T_1(p)= p, 
with p=λ h,

    X_1=R_s,η(p) X_0=T_s(ω_0+ω_1 p)/T_s(ω_0) X_0,

where the dependence of the stability function R_s,η on the parameters s and η is emphasized with a corresponding subscript. 
The real negative interval (-C_s(η)· s^2 ,0) is included in the 
stability domain of the method

    S:={p∈ℂ; |R_s,η(p)|≤ 1}.

 The constant C_s(η) ≃ 2-4/3 η depends on the so-called damping parameter η and for η=0,
 it reaches the maximal value C_s(0)=2. Hence, given the stepsize h, for systems with a Jacobian having large real negative eigenvalues (such as diffusion problems) with spectral radius λ_max at X_n, the parameter s for the next step X_n+1 can be chosen adaptively as[The notation [x] stands for the integer rounding of real numbers.]

    s=[√(hλ_max+1.5/2-4/3 η)+0.5],

see <cit.> in the context of deterministically stabilized schemes of order two with adaptative stepsizes.
The method (<ref>)
 is much more efficient as its stability domain increases quadratically with the number s of function evaluations while a composition of s  explicit Euler steps
 (same cost) has a stability domain that only increases linearly with s.
 
	In Figure <ref>(a) we plot the complex stability domain {p∈ ;  |R_s,η(p)| ≤ 1}
	for s=7 stages and different values η=0, η=0.05 and η=3.98, respectively.
	We also plot in Figure <ref>(b) the corresponding stability function R_s,η(p) as a function of p real, to illustrate that 
	the stability domain along the negative real axis corresponds to the values for which |R_s,η(p)|≤1.
	We observe that in the absence of damping (η=0), the stability domain includes the large real interval [-2· s^2,0] of width 2· 7^2=98. However
	for all p that are a local extrema of the stability function, where 
	|R_s,η(p)|=1, the stability domain is very thin and does not include a neighbourhood close to the negative real axis.
	To make the scheme robust with respect to small perturbations of the eigenvalues, it is therefore needed to add some damping 
	and a typical value is η=0.05, see for instance the reviews 
	<cit.>.
	The advantage is that the stability domain now includes a neighbourhood of the negative real axis portion. The price of this improvement is a slight reduction
	of the stability domain size C_η s^2, where C_η≃ 2-4/3η.
	 Chebyshev methods have been first generalized for Itô SDEs in <cit.> (see <cit.> for Stratonovitch SDEs)  
	with the following stochastic orthogonal Runge-Kutta-Chebyshev method (S-ROCK):[A variant with analogous stability properties is proposed in <cit.>
	with g^r(K_s) replaced by g^r(K_s-1) in (<ref>).]
	
    K_0   =   X_0 
    
    	K_1   =   X_0+μ_1hf(X_0) 
    
    	K_i   =   μ_ihf(K_i-1)+ν_iK_i-1+κ_iK_i-2,   i=2,…,s, 
    
    	X_1   =   K_s + ∑_r=1^mg^r(K_s)Δ W_r,

where the coefficients μ_i,ν_i,κ_i are defined in (<ref>),(<ref>). In contrast to the deterministic method (<ref>), where η is chosen small and fixed (typically η=0.05), in stochastic case for the classical S-ROCK method <cit.>, the damping η=η_s is not small and chosen as an increasing function of s that plays a crucial role in stabilizing the noise and in obtaining an increasing portion of the true stability domain (<ref>) as s increases. 


	
	In the context of stiff SDEs, a relevant stability concept is that of mean-square stability.
	A test problem widely used in the literature is <cit.> ,
	
    dX(t)= λ X(t)dt + μ X(t) dW(t),     X(0)=1,

	in dimensions d=m=1 with fixed complex parameters λ,μ.
	Note that other stability test problem in multiple dimensions are also be considered in <cit.> and references therein.
	The exact solution of (<ref>) is called mean-square stable if lim_t→∞𝔼(|X(t)|^2)=0
	and this holds if and only if (λ,μ) ∈^MS,
	where
	
    ^MS = {(λ,μ)∈ℂ^2 ; (λ)+1/2|μ|^2<0}.

	Indeed, the exact solution of (<ref>) is given by X(t)=exp((λ+1/2μ^2) t+μ W(t)), and an application of the the Itô formula yields
	(|X(t)|^2)=exp(( (λ)+1/2μ^2) t)
	which tends to zero at infinity if and only if (λ)+1/2μ^2<0.
	We say that a numerical scheme {X_n} for the test problem (<ref>) is mean-square stable if and only if lim_n→∞(|X_n|^2) =0.
	For a one-step integrator applied to the test SDE (<ref>), we obtain in general a induction of the form
	
    X_n+1= R(p,q,ξ_n)X_n,

	where p=λ h, q =μ√(h), and ξ_n is a random variable
	(e.g.
	a Gaussian
	ξ_n∼𝒩(0,1) 
	or a discrete random variable).
	Using (|X_n+1|^2)= (|R(p,q,ξ_n)|^2) (|X_n|^2), we obtain the mean-square stability condition
	<cit.>
	
    lim_n →∞(|X_n|^2)=0
    	 (p,q)∈_num,
 
	where we define 𝒮_num={(p,q)∈ℂ^2 ;  |R(p,q,ξ)|^2<1 }.
	The function R(p,q,ξ_n) is called the stability function of the one-step integrator.
	For instance, the stability function of the Euler-Maruyama method (<ref>) reads R(p,q,ξ)=1+p+qξ and we have (|R(p,q,ξ)|^2)=(1+p)^2 + q^2.
	
	We say that a numerical integrator is mean-square A-stable if ^MS⊂_num.
	This means that the numerical scheme applied to (<ref>) is mean-square stable for all 
	all h>0 and all (λ,μ) ∈^MS for which the exact solution of (<ref>) is mean-square stable. An explicit Runge-Kutta type scheme cannot however be mean-square stable because its stability domain _num is necessary bounded along the p-axis.
	Following <cit.>, we consider the following portion of the true mean-square stability domain
	
    𝒮_a  = {(p,q) ∈ (-a,0) × ; p+1/2 |q|^2 < 0},

	and define for a given method 
	
    L=sup{a>0 ; 𝒮_a⊂𝒮_num}.

	We search for explicit schemes for which the length L of the stability domain is large.
	For example, for the classical S-ROCK method <cit.>, the value η=3.98 is the optimal damping maximising L for s=7 stages and we can see in Figure <ref> that this damping reduces significantly the stability domain compared to the optimal deterministic domain.

	
	
	
	
	
	
	
	
	
	
	
	

	

		
	The new S-ROCK method, denoted SK-ROCK   (for stochastic second kind orthogonal Runge-Kutta-Chebyshev method) introduced in this paper is defined as
	
    K_0   =   X_0 
    
    	K_1   =   X_0+μ_1hf(X_0+ν_1 Q) +κ_1Q 
    
    	K_i   =   μ_ihf(K_i-1)+ν_iK_i-1+κ_iK_i-2,   i=2,…,s.
    
    	X_1   =   K_s,

	where
	Q=∑_r=1^mg^r(X_0)Δ W_r,
	and μ_1=ω_1/ω_0,ν_1=s ω_1/2,κ_1=s ω_1/ω_0 and μ_i,ν_i,κ_i, i=2,…,s are given by (<ref>), with a fixed small damping parameter η.
	In the absence of noise (g^r=0,r=1,…, m, deterministic case), this method coincides with the standard deterministic order 1 Chebychev method, see the review <cit.>.
	We observe that the new class of methods (<ref>) is closely related to the
	standard S-ROCK method (<ref>).
	Comparing the two schemes (<ref>) and (<ref>),
	the two differences are on the one hand that the noise term is computed at the first internal stage K_1 for (<ref>), whereas it is computed at the final stage in (<ref>), and on the other hand, for the new method (<ref>) the damping parameter η involved in (<ref>)
	is fixed and small independently of s (typically η=0.05), whereas for the standard method (<ref>), the damping η is an increasing function of s, optimized numerically for each number of stages s.

	
	If we apply the above scheme (<ref>) to  the linear test equation (<ref>), 
we obtain 	
    X_n+1 = R(p,q,ξ_n) X_n,

where
 
    (|R(p,q,ξ)|^2) = A(p)^2 + B(p)^2 q^2,

		and
		
    A(p)=T_s(ω_0+ω_1 p)/T_s(ω_0)    
    		B(p)=U_s-1(ω_0+ω_1 p)/U_s-1(ω_0)(1+ ω_1/2 p)

		correspond to the drift and diffusion contributions, respectively. 
The above stability function  (see Lemma <ref> in Section <ref>) is obtained by using the recurrence relation for the first kind Chebyshev polynomials (<ref>) and the similar recurrence relation for the second kind Chebyshev polynomials
	
    U_j(p) = 2p U_j-1(p) - U_j-2(p),

		where U_0(p)= 1,U_1(p)=2p.
		Notice that the relation T_s'(p)=s U_s-1(p) between first and second kind Chebyshev polynomials will be repeatedly used in our analysis.



	
	In Figure <ref>(b)(d), we plot the mean-square stability domain of the SK-ROCK method for s=7 and s=20 stages, respectively and the same small damping η=0.05 as for the deterministic Chebyshev method. We observe that the stability domain has length L_s≃ (2-4/3η)s^2.
	For comparison, we also include in Figure <ref>(a)(c) the mean-square stability domain of the standard S-ROCK method with smaller stability domain size L_s≃ 0.33· s^2.
	
	
In Figure <ref>, we plot the stability function (|R(p,q,ξ)|^2) in (<ref>) as a function
		of p for various scaling of the noise for s=7 stages and damping η=0.05. 
		We see that it is bounded by 1 for p∈(-2(1-2/3η)s^2,0) which is proved asymptotically in Theorem <ref>.
		The case q=0 corresponds to the deterministic case, and we see in Figure <ref>(a), the polynomial (|R(p,0,ξ)|^2)=A(p)^2.
		Noticing that (|R(p,q,ξ)|^2) is an increasing function of q, the case q^2=-2p represented in Figure <ref>(c) corresponds to the upper border of the stability domain _L defined in (<ref>) (note that this is the stability function value along the dashed boundary in Figure <ref>), while the scaling q^2=-p in Figure <ref>(c) is an intermediate regime. In Figures <ref>(b)(c), we also include the drift function A(p)^2 (red dotted lines) and diffusion function B(p)^2 q^2 (blue dashed lines), and it can be observed that their oscillations alternate, which means that any local maxima of one function is close to a zero of the other function.
		This is not surprising because A(p) and B(p) are related to the first kind and second kind Chebyshev polynomials, respectively, corresponding to the cosine and sine functions.
		This also explains how a large mean-square stability domain can be achieved by the new SK-ROCK method (<ref>) with a small damping parameter η, in contrast to the standard S-ROCK method (<ref>) from <cit.> that uses a large and s-dependent damping parameter η with smaller stability domain size L_s≃ 0.33· s^2(see Figures <ref>(a)(c)).






	
	

	
	
		
	

§ MEAN-SQUARE STABILITY ANALYSIS

	
	
	In this section, we prove asymptotically that the new SK-ROCK methods have an extended mean-square stability domain with size Cs^2	growing quadratically as a function of the number of internal stages s, where the constant C ≥ 2-4/3η is the same as the optimal constant of the standard Chebyshev method in the deterministic case, using a fixed and small damping parameter η.
	
	  
		Let s≥ 1 and η≥ 0.
		Applied to the linear test equation dX=λ X dt + μ X dW, the scheme (<ref>) yields
		
    X_n+1 = R(λ h,μ√(h),ξ_n) X_n

		where p=λ h, q=μ√(h),  ξ_n ∼(0,1) is a Gaussian variable and the stability function given by
		
    R(p,q,ξ) = T_s(ω_0+ω_1 p)/T_s(ω_0) 
    		+ U_s-1(ω_0+ω_1 p)/U_s-1(ω_0)(1+ ω_1/2 p) q ξ.

	
	
	Indeed, we take advantage that T_j and U_j have the same recurrence relations (<ref>),(<ref>),	and only the initialization changes with T_1(x)= x and U_1(x)= 2x,
		we deduce Q=X_0 μ√(h)ξ, and we obtain by induction on i≥ 1,
		
    K_i=T_i(ω_0+ω_1p)/T_i(ω_0)X_0+U_i-1(ω_0+ω_1p)/T_i(ω_0) (1+ω_1p/2)s ω_1Q

		and we use T_s'(x)=xU_s-1(x) and sω_1/T_s(ω_0)=1/U_s-1(ω_0), which yields the result for X_1=K_s.
	
	
	
	For a positive damping η, we prove the following main result of this section, showing that a quadratic growth L≥ (2-4/3 η)s^2 of the mean-square stability domain defined in (<ref>) is achieved for all 	η small enough and all stage number s large enough.
	
	 
		There exists η_0>0 and s_0 such that for all η∈ [0,η_0] and all s≥ s_0,
		
		for all p∈ [-2ω_1^-1,0] and p + 1/2|q|^2 ≤ 0, we have (|R(p,q,ξ)|^2) ≤ 1.
	
	
	
		 
		We deduce from Theorem <ref>, that the mean-square stability domain size 
		(<ref>) of SK-ROCK
		grows as (2-4/3 η)s^2 which is arbitrarily close to the optimal stability domain size 2s^2 for η→ 0.
		Indeed, for s→∞ and all η≤η_0, we have 2ω_1^-1s^-2→ 2tanh(√(2η))/√(2η) = 2-4/3 η + (η^2) and for all s,η, we have 2ω_1^-1≥ (2-4/3 η) s^2.
		In addition, in the special case of a zero damping (η=0), 
		
		
		the stability function (<ref>) reduces to
		
    R(p,q,ξ) = T_s(1+p/s^2) + s^-1U_s-1(1+p/s^2)(1+p/2s^2) q ξ,

		and it holds 
		
    (|R(p,q, ξ)|^2) ≤ 1,

		for all s≥1, for all p∈[-2s^2,0] and all q ∈ such that p+|q|^2/2 ≤ 0.
		Indeed, for p∈[-2s^2,0], we denote cosθ = 1+p/s^2∈[-1,1] and using 
    T_s(cos(θ))=cos(sθ),    sin(θ) U_s-1(cos(θ))=sin(sθ),
 we obtain
		
    (|R(p,q, ξ)|^2) ≤(|R(p,√(-2p))  = cos(sθ)^2+ sin(sθ)^2 1+cosθ/2≤ 1,

		where we used -2p=2s^2(1-cosθ), 1+p/2s^2 = 1+cosθ/2 and sin^2 θ = (1+cosθ)(1-cosθ).
	
	
	
	Before we prove Theorem <ref>, we have the following lemma, 
	see <cit.> for analogous results.
	
		We have the following convergences as s→∞ to analytic functions[
			Note that for z<0, we can use √(2z)=i√(-2z) and obtain T_s(1+z/s^2)→cos(√(-2z)) for s→∞ and similarly α(z)=(√(-2z)).] 
		uniformly for z in any bounded set of the complex plan,
		
    T_s(1+z/s^2)    →cosh√(2z),
    
    		s^-1U_s-1(1+z/s^2)   →α(z):= sinh√(2z)/√(2z),
    ω_1s^2   →Ω(η)^-1,   Ω(η):=tanh√(2η)/√(2η).

	

	
		We prove the uniform convergence of the first limit only, since it will be useful in the proof of the next theorem. The others can be proved in a similar way.
		
		First, let us write the two functions of η in Taylor series,
		
    s^-1U_s-1(ω_0)   =   
    			s^-1∑_n=0^s-1U_s-1^(n)(1)/n!(η/s^2)^n
    			
    			=∑_n=0^s-1( 1/n!∏_k=0^n(1-k^2/s^2)∏_k=0^n1/2k+1)η^n,
    α(η)   =   ∑_n=0^∞2^nη^n/(2n+1)!
    			=∑_n=0^∞(1/n!∏_k=0^n1/2k+1)η^n,
	
		where we used the formula
sU_s-1^(n-1)(1)=T_s^(n)(1)=∏_k=0^n-1s^2-k^2/2k+1.
Subtracting the above two identities, we deduce
		
		
    sup_η∈[-η_0,η_0]|s^-1U_s-1(ω_0)-α(η)|
    			
    			
    			   ≤   ∑_n=0^s-1η_0^n/n!(1-∏_k=0^n(1-k^2/s^2))∏_k=0^n1/2k+1
    			+∑_n=s^∞η_0^n/n!
       ≤   ∑_n=0^s-1η_0^n/n!(1-∏_k=0^n(1-k^2/s^2))1/2s-1
    			+∑_n=s^∞η_0^n/n!

		Noticing that η_0^n/n!(1-∏_k=0^n(1-k^2/s^2))1/2s-1 converges to zero as s→∞ and is bounded
		by η_0^n/n! for all integers s,n, which is the general term of the convergent series of exp(η_0)=∑_n=0^∞η_0^n/n!, the Lebesgue dominated convergence theorem implies that (<ref>)
		converges to zero as s→∞, which concludes the proof.
		
		
			
		
		
	
	
	
	
	
	
		For all η small enough and all s large enough, we have the following estimate:
		
    s^2 ω_1/T_s(w_0)^21-(1-ω_1)^2/1-(ω_0-ω_1)^2≤ 1

	
	
		Using the Lemma <ref> we have for s→∞, uniformly for all η∈[0,η_0],
		
    s^2ω_1/T_s(ω_0)^2→2√(2η)/sinh(2√(2η)) and 1-(1-ω_1)^2/1-(ω_0-ω_1)^2→1/1-Ω(η)η.

		Now if we expand both functions in Taylor series we get:
		
    2√(2η)/sinh(2√(2η))=1-4/3η+(η^2),    1/1-Ω(η)η=1+η+(η^2),

		and this implies that for all s large enough and all η≤η_0,
		
    s^2 ω_1/T_s(w_0)^21-(1-ω_1)^2/1-(ω_0-ω_1)^2≤ (1-4/3η_0+(η_0^2))(1+η_0+(η_0^2))=1-1/3η_0+(η_0^2),

		which is less than 1 for η_0 small enough.
	
	
	
	Numerical evidence suggests that the result of Theorem <ref> holds for all s≥ 1 and all η≥ 0.
		Indeed, it can be checked numerically that (<ref>) holds for all η∈(0,1) and all s≥1.
	
	
	
		Setting x=w_0+w_1p, a calculation yields
		
    (|R(p,q,ξ)|^2)    ≤   (|R(p,√(-2p),ξ)|^2)
       =   T_s(x)^2/T_s(w_0)^2 + U_s-1(x)^2/U_s-1(w_0)^2 (1+w_1/2p)^2 (-2p)

		The proof is conducted in two steps, where we treat separately the cases p∈ [-2ω_1^-1,-1] and p∈[-1,0]. 
		For the first case p∈ [-2ω_1^-1,-1], which corresponds to x∈ [-1+η/s^2,ω_0-ω_1], we have
		
    (|R(p,q,ξ)|^2)     =   T_s(x)^2/T_s(w_0)^2 + U_s-1(x)^2/U_s-1(w_0)^2(1-w_0-x/2)^2 2w_0-x/w_1
       =   T_s(x)^2/T_s(w_0)^2 + U_s-1(x)^2(1-x^2) Q_s(x)

		where we denote
		
    Q_s(x)    =   
    			
    			s^2 ω_1/T_s(w_0)^2(1+x-η/s^2/2) 1-(x-η/s^2)^2/1-x^2

		First, we note that 
		1+x-η/s^2/2∈ [0,1-ω_1/2].
		Next, using η/s^2≤ 2, we deduce d/dx( 1-(x-η/s^2)^2/1-x^2) = 2η/s^21+x^2-η/s^2 x/(1-x^2)^2≥
		2η/s^2(1-x)^2/(1-x^2)^2≥ 0.
		Thus, 1-(x-η/s^2)^2/1-x^2 is an increasing function of x, smaller than its value at x=ω_0-ω_1,
		
    1-(x-η/s^2)^2/1-x^2≤1-(1-ω_1)^2/1-(ω_0-ω_1)^2

		Using Lemma <ref> we obtain |Q_s(x)| ≤ 1. This yields (|R(p,q,ξ)|^2) ≤ T_s(x)^2 + U_s-1(x)^2(1-x^2) =1.
		
		For the second case p∈[-1,0] which corresponds to x∈ [ω_0-ω_1,ω_0], we deduce from T_s(x)^2 + U_s-1(x)^2(1-x^2) =1 that
		
    (|R(p,q,ξ)|^2)
    			   ≤   1/T_s(w_0)^2 + U_s-1(x)^2/U_s-1(w_0)^2( (1+w_1/2p)^2 (-2p)-(1-x^2)U_s-1(ω_0)^2/T_s(w_0)^2)

		Using Lemma <ref>, we get 
		
    (|R(p,q,ξ)|^2)
    			   ≤   1/T_s(w_0)^2 + U_s-1(x)^2/U_s-1(w_0)^2( (1+w_1/2p)^2 (-2p)-(1-x^2)U_s-1(ω_0)^2/T_s(w_0)^2)
       →   l(η,p):=1/cosh^2√(2η)+α(η+p/Ω(η))^2/α(η)^2(-2p(Ω(η)-1)+2Ω(η)^2η),

		for s→∞, where the above convergence is uniform for p∈[0,1],η≤η_0.
		Using the fact that Ω(η)=1-2/3η+O(η^2), we deduce 
		 
    ∂ l/∂η|_η=0=-2+α(p)^2(-4/3p+2).

		By Taylor series in the neighbourhood of zero we have 
		α(p)^2=1+2/3p+8/45p^2+O(p^3),
		and for p∈[-1,0], α(p)^2≤ 1+2/3p+8/45p^2, thus
		for all p∈[-1,0],
		
    ∂ l/∂η|_η=0≤-2+(1+2/3p+8/45p^2)(-4/3p+2)
    			=-8/135p^2(4p+9)≤0.

		Therefore, there exists η_0 small enough such that for all p∈[-1,0],η≤η_0, l(η,p)≤ l(0,p)=1.
		This concludes the proof of Theorem <ref>.
	
	
	
	
	

§ CONVERGENCE ANALYSIS

	
	
	
	We show in this section that the proposed scheme (<ref>) has strong order 1/2 and weak order 1
	for general systems of SDEs of the form (<ref>) with Lipschitz and smooth vector fields, analogously to the simplest Euler-Maruyama method.
	
	We denote by C_P^4(^d,^d) the set of functions from ^d to ^d that are 4 times continuously differentiable
	with all derivatives with at most polynomial growth.
The following theorem shows that the proposed SK-ROCK has strong order 1/2 and weak order 1 for general SDEs.
	 
		Consider the system of SDEs (<ref>) on a time interval of length T>0, with f,g ∈ C_P^4(^d,^d), Lipschitz continuous.
		Then the scheme (<ref>) has strong order 1/2 and weak order 1,
		
    |( X(t_n) - X_n )     ≤ Ch^1/2,     t_n=nh≤ T, 
    
    		| (ϕ(X(t_n))) - (ϕ(X_n)) |     ≤ Ch,     t_n=nh≤ T,

		for all ϕ∈ C_P^4(^d,), where C is independent of n,h.
	
For the proof the Theorem <ref>, the following lemma will be useful.
		It relies on the linear stability analysis of Lemma <ref>.
	 
		The scheme (<ref>) has the following Taylor expansion after one timestep,
		
    X_1 = X_0 + hf(X_0) + ∑_r=1^mg^r(X_0)Δ W_r + h (T_s”(ω_0)ω_1^2/T_s(ω_0) + ω_1/2) f'(X_0) ∑_r=1^mg^r(X_0)Δ W_r + h^2 R_h(X_0),

		where all the moments of R_h(X_0) are bounded uniformly with respect to h assumed small enough, with a polynomial growth with respect to X_0.
	
	
		Using the definition (<ref>) of the scheme and the recurrence relations (<ref>),(<ref>), we obtain by induction on i=1,…,s,
		
    K_i    = X_0 + h T_i'(ω_0)ω_1/T_i(ω_0) f(X_0) + sT_i'(ω_0)ω_1/iT_i(ω_0)∑_r=1^mg^r(X_0)Δ W_r 
       + h (sT_i”(ω_0)ω_1^2/iT_i(ω_0) + sT_i'(ω_0)ω_1^2/2iT_i(ω_0)) f'(X_0) ∑_r=1^mg^r(X_0)Δ W_r + h^2 R_i,h(X_0),

		and R_i,h(X_0) has the properties claimed on R_h(X_0). Using ω_1 = T_s(ω_0)/T_s'(ω_0), this yields the result for X_1=K_s.
	
	
	
	
	
		A well-known theorem of Milstein <cit.> (see <cit.>)
		allows to infer the global orders of convergence from the error after one step.
		We first show that for all r∈ℕ the moments (|X_n|^2r) are bounded for all n,h with 0≤ nh ≤ T uniformly
		with respect to all h sufficiently small.
		Then, it is sufficient to show the local error estimate
		
    |(ϕ(X(t_1))) -(ϕ(X_1))| ≤ C h^2,

		for all initial value X(0)=X_0 and where C has at most polynomial growth with respect to X_0, to deduce the weak convergence estimate (<ref>).
		For the strong convergence (<ref>), using the classical result from <cit.>, it is sufficient to show in addition the local error estimate
		
    (X(t_1) - X_1) ≤ C h

		for all initial value X(0)=X_0 and where C has at most polynomial growth with respect to X_0. 
		These later two local estimates are an immediate consequence of Lemma <ref>.
		
		To conclude the proof of the global error estimates, it remains to check that for all r∈ℕ the moments (|X_n|^2r) are bounded  uniformly with respect to all h small enough for all 0≤ nh ≤ T. 
		We use here the approach of <cit.> which states that
		it is sufficient to show
		
    |(X_n+1-X_n| X_n)| ≤ C (1+| X_n|) h,      | X_n+1- X_n| ≤ M_n (1+| X_n|) √(h),

		where C is independent of h and M_n is a random variable with moments of all orders bounded uniformly
		with respect to all h small enough.
		These estimates are a straightforward consequence of the definition 
		(<ref>) of the scheme and the linear growth of f,g (a consequence of their  Lipschitzness). 
		This concludes the proof of Theorem <ref>.
	
	
		In the case of additive noise, i.e. g^r,r=1,…,m are constant functions, one can show that the 
		order of strong convergence (<ref>) become 1, analogously to the case of the Euler-Maruyama method.
		For a general multiplicative noise, a scheme of strong order one can also be constructed with (|R(p,q,ξ)|^2) ≤ 1
		for all p∈[-2ω_1^-1,0] and all q with p+|q|^2/2≤ 0, as it can be check numerically.
		The idea is to modify the first stages of the scheme such that the stability function (<ref>)   becomes
		
    R(p,q,ξ) = T_s(ω_0+ω_1 p)/T_s(ω_0) 
    		+ U_s-1(ω_0+ω_1 p)^2/U_s-1(ω_0)^2(1+ w_1/2 p-ω_1^4/2 p^2)(q ξ+q^2ξ^2-1/2).

		We refer to <cit.> for details.
	

	

§ LONG TERM ACCURACY FOR BROWNIAN DYNAMICS

	
	In this section we discuss the long-time accuracy of the SK-ROCK for Brownian dynamics (also called overdamped Langevin dynamics). We will
	see that using postprocessing techniques we can derive an SK-ROCK method that captures the invariant measure of Brownian dynamics with second order accuracy. In doing so, we do not need our stabilized method to be of weak order 2 on bounded time intervals and we obtain a method that is cheaper than the stochastic orthogonal Runge-Kutta-Chebyshev method of weak order 2 (S-ROCK2) proposed in <cit.>, as S-ROCK2 uses many more function evaluations per time-step and a smaller stability domain.
	

 §.§ An exact SK-ROCK method for the Orstein-Uhlenbeck process

	
	We consider the 1-dimensional Orstein-Uhlenbeck problem with 1-dimensional noise with constants δ,σ>0,

    dX(t) = -δ X(t)dt  + σ dW(t),

	that is ergodic and has a Gaussian invariant measure with mean zero and variance given by lim_t→∞(X(t)^2)=σ^2/(2δ).
	Applying the SK-ROCK method to the above system we obtain 
	
    X_n+1=A(p)X_n+B(p)σ√(h)ξ_n

	where p=-δ h, ξ_n ∼(0,1) is a Gaussian variable and similarly as for  (<ref>) we have
	
    A(p)=T_s(ω_0+ω_1 p)/T_s(ω_0),   B(p)= U_s-1(ω_0+ω_1 p)/U_s-1(ω_0)(1+ ω_1/2 p).

		A simple calculation (using that |A(p)|<1) gives
	
    lim_n→∞(X_n^2)=σ^2/2δ R(p),    
    	R(p)= 2p B(p)^2/A(p)^2-1.

	From the above equation, we see that the SK-ROCK method has order r for the invariant measure of (<ref>) if and only if R(p)=1+ O(p^r) and a short calculation using (<ref>) reveals that
	R(p)=1+ O(p), it has order one for the invariant measure (this is of course not surprising because the 
	SK-ROCK  has weak order one).
	We next apply the techniques of postprocessed integrators popular in the deterministic literature <cit.> and proposed in the stochastic context in <cit.>.
	The idea is to consider a postprocessed dynamics X_n=G_n(X_n) (of negligible cost) such that the process X_n approximates the invariant measure of the dynamical system with higher order. For the process (<ref>), we consider the postprocessor
	
    X_n =X_n+cσ√(h)ξ_n,

	which yields lim_n→∞(X_n^2)=σ^2/2δ (R(p)-2c^2p).
In the case of the SK-ROCK method with η=0 (zero damping), we have A(p)=T_s(1+p/s^2),
B(p)=U_s-1(1+p/s^2)(1+p/(2s^2))/s. Setting c=1/(2s) and using the identity
(1-x^2)U_s-1^2(x)=1-T_s^2(x)
with x=1+p/s^2 reveals that R(p)-2c^2p=1 and we obtain

    lim_n→∞(X_n^2)=σ^2/2δ.

Hence the postprocessed SK-ROCK method (that will be denoted PSK-ROCK) captures exactly the invariant measure of the 1-dimensional Orstein-Uhlenbeck problem (<ref>). Such a behavior is known for the drift-implicit θ method with θ=1/2 (see <cit.> in the context of the stochastic heat equation) and has recently also been shown for the
non-Markovian Euler scheme <cit.>. In <cit.> an interpretation of the scheme <cit.> as an Euler-Maruyama method with postprocessing (<ref>) with c=1/2
has been proposed and we observe that this is exactly the same postprocessor as for the PSK-ROCK method (with s=1,η=0). As the SK-ROCK method with zero damping is mean-square stable (see Remark <ref> for η=0),
it can be seen as a stabilized version of the scheme <cit.>. However, the PSK-ROCK method with s>1 and zero damping is not robust to use as its stability domain along the drift axis does not allow for any imaginary perturbation at the points where  |T_s(1+p/s^2)|=1 and it is not ergodic (see Remark <ref> below). 



  
Stability analysis for Orstein-Uhlenbeck
	
	Let M ∈^d× d denote a symmetric matrix with eigenvalues -λ_d ≤…≤ -λ_1 < 0, and consider the d-dimensional Orstein-Uhlenbeck problem
	
    dX(t) = MX(t)dt + σ dW(t)

	where W(t) denotes a d-dimensional standard Wiener process.
	The following theorem shows that the damping parameter η>0 plays an essential role to warranty the convergence to the 
	numerical invariant measure ρ^h_∞(x)dx at an exponentially fast rate.
	 
		Let η>0.
		Consider the scheme (<ref>) with postprocessor (<ref>) applied to (<ref>)
		with stepsize h and stage parameter s such that 2ω_1^-1≥ hλ_d.
		Then, for all h≤η/λ_1,ϕ∈ C_P^1(^d,),
		
    |(ϕ(X_n)) - ∫_^dϕ(x)ρ^h_∞(x)dx | ≤ Cexp(-λ_1(1+η)^-1 t_n)

		where
		C is independent of h,n,s,λ_1,…,λ_d.
	
	
		It is sufficient to show the estimate
		
    |A(-λ_j h)| ≤exp(-λ_1(1+η)^-1 h)

		for all h≤ h_0,
		where we denote A(z)=T_s(ω_0 +ω_1 z)/T_s(ω_0).
		Indeed, considering two initial conditions X_0^1,X_0^2 for (<ref>) and the corresponding numerical solutions X_n^1,X_n^2  (obtained for the same realizations of {ξ_n}) with postprocessors X_n^1,X_n^2, we obtain X_n^1 - X_n^2 = A(hM) (X_n-1^1 - X_n-1^2) and using the matrix 2-norm A(hM) =max_j |A(-λ_j h)| and (<ref>), we deduce by induction on n,
		
    X_n^1 - X_n^2 = X_n^1 - X_n^2≤exp(-λ_1(1+η)^-1 t_n) X_0^1 - X_0^2,

		and taking X_0^2 distributed according to the numerical invariant measure yields the result. 
		
		For the proof of (<ref>), let z=-λ_j h. 
		Consider first the case z∈ (-ηω_1^-1s^-2,0).
		Using the convexity of A(z) on [-ηω_1^-1s^-2,0] (note that T_s'(x) is increasing on [1,∞)), we can bound A(z) by the affine function
		passing by the points (x_1,A(x_1)),(x_2,A(x_2)) with 
		x_1=-ηω_1^-1s^-2, x_2=0,
		
    A(z) ≤ 1 + z(1-1/T_s(ω_0))η^-1ω_1s^2

		
		Using 
		ω_1s^2≥ 1 and T_s(ω_0)≥ 1+η, we obtain
		
    (1-1/T_s(ω_0))η^-1ω_1s^2≥ 
    		(1-(1+η)^-1)η^-1 = (1+η)^-1.

		This yields for all z∈ [-ηω_1^-1s^-2,0],
		
    A(z) ≤ 1 + z(1+η)^-1≤exp(z(1+η)^-1)

		where we used the convexity of exp(z(1+η)^-1) bounded from below by its tangent at z=0. We obtain
		
    A(-λ_j h) ≤ e^-λ_j h (1+η)^-1≤ e^-λ_1 h (1+η)^-1.

		We now consider the case z∈
		[-L_s,-ηω_1^-1s^-2]. 
		We have |ω_0 +ω_1 z|≤ 1, thus |T_s(ω_0 +ω_1 z)|≤ 1 and
		
    |A(z)| ≤ 1/T_s(ω_0) ≤exp(-λ_1(1+η)^-1 h)
 
		for all h≤ (1+η)log (T_s(ω_0))/λ_1, and thus also for 
		h≤η/λ_1 where
		we use T_s(ω_0) ≥ 1+η and (1+η)log (1+η) ≥η.
		This concludes the proof.
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
	
	 
	  
		
		Note that η>0 is a crucial assumption in Theorem <ref>.
		Indeed, the estimate of Theorem <ref> is false for η=0 already in dimension d=1 for all s>1: for a stepsize h such that 
		1-hλ_1/s^2=cos(π/s) we obtain A(-λ_1h)=-1 and B(-λ_1h)=0 in (<ref>)  (corresponding to the local extrema p=-λ_1 h of A(p) closest to zero)
		and X_n=(-1)^n X_0 for all n, and the scheme is not ergodic.
		In addition, notice that Theorem <ref> allows to use an h-dependent value 
		of η such as η= hλ_1
		where λ_1≥λ_1 is an upper bound for λ_1.
		In this case, the exponential convergence of Theorem <ref> holds for all stepsize
		h≤ 1.
	
	
	
We end this section by noting that being exact for the invariant measure of Brownian dynamics (<ref>) is only true for the PSK-ROCK method (or the method in <cit.>) 
in the linear case, i.e. for a quadratic potential V.

Second order accuracy for the invariant measure 

has been shown in <cit.> (see also <cit.>) for the method in (<ref>) (equivalent to PSK-ROCK with s=1 stage) for general
nonlinear Brownian dynamics (<ref>). This will also be shown for the nonlinear PSK-ROCK method in the next section.

	
	

 §.§ PSK-ROCK: a second order postprocessed SK-ROCK method for nonlinear Brownian dynamics

	We consider the overdamped Langevin equation,
	
    dX(t) = -∇ V(X(t))dt + σ dW(t),

	where the stochastic process X(t) takes values in ^d and W(t) is a d-dimensional Wiener process.
	We assume that the potential V:^d → has class C^∞ and satisfies the at least quadratic growth assumption
	
    x^T∇ V(x) ≥ C_1x^T x -C_2

	for two constants C_1,C_2>0 independent of x∈^d.
	The above assumptions warranty that the system (<ref>) is ergodic with exponential convergence
	to a unique invariant measure with Gibbs density ρ_∞ = Z exp(-2σ^-2V(x)),
	
    |(ϕ(X(t)) - ∫_^dϕ(x) ρ_∞(x) dx| ≤ C e^-λ t,

	for test function ϕ and all initial condition X_0, where C,λ are independent of t.
	
	We propose to modify the internal stage K_1=X_0+μ_1hf(X_0+ν_1Q)+κ_1Q of the method (<ref>) as follows:
	
    K_1=X_0+μ_1hf(X_0+ν_1Q)+κ_1Q+α h(f(X_0+ν_1Q)-2 f(X_0)+f(X_0-ν_1Q)),

	

	where α is a parameter depending on s and η given in 
Theorem <ref> below.

	
		Notice that for α=0, we recover the original definition from (<ref>). 
		We note that the parameter α does not modify 
		the stability function of Lemma <ref>, and yields a perturbation of order (h^2) in the definition of X_1.
		Thus, the results of Theorem <ref> and Theorem <ref>
		remain valid for any value of α for the scheme (<ref>)
		with modified internal stage (<ref>).
	
	

		
		
	 
	Consider the Brownian dynamics (<ref>), where we assume that V:^d → has class C^∞, with ∇ V globally Lipschitz and satisfying (<ref>).
		Consider the scheme (<ref>) applied to (<ref>) with modified internal stage K_1 defined in (<ref>) with α defined in (<ref>),
	and the postprocessor defined as
		
    X_n = X_n + cσ√(h)ξ,

		where 
				
    c^2 = -1/4 + ω_1/2 +ω_1T_s”(ω_0)/T_s'(ω_0)-ω_1^2T_s”(ω_0)/4T_s(ω_0),     α =2/sω_0ω_1 (c^2+ ω_1^2T_s”(ω_0)/2 T_s(ω_0)-r_s),

		and  r_s is defined by induction as r_0=0, r_1=s^2ω_1^3/4ω_0:=Δ_1
    and
		
    r_i=ν_ir_i-1+κ_ir_i-2+Δ_i,    Δ_i=μ_isT_i-1'(ω_0)ω_1/(i-1)T_i-1(ω_0),     i=2,… s.

		
		
		
		
		
		Then, X_n yields order two for the invariant measure, i.e. (<ref>) holds with r=2, and in addition
		
    |(ϕ(X_n) - ∫_^dϕ(x) ρ_∞(x) dx| ≤ C_1 e^-λ t_n + C_2h^2,

		for all t_n=nh, ϕ∈ C_P^∞(^d,), where C_1,C_2 are independent of h assumed small enough, and C_2 is independent of the initial condition X_0.
	
		The proof of Theorem <ref> relies on the following postprocessing analysis  from <cit.>.
	Consider a scheme (<ref>) with bounded moments and assumed ergodic when applied to (<ref>), where the potential V satisfies the above ergodicity assumptions. 
	Assume that the scheme has a weak Taylor expansion after one time step of the form
	
    (ϕ(X_1)|X_0=x) = ϕ(x) + hϕ(x) + h^2𝒜_1 ϕ(x) + (h^3),

		and consider a postprocessor of the form X_n = G_n(X_n) where 
			
    (ϕ(X_1)|X_1=x) = ϕ(x) + h𝒜_1 ϕ(x) + (h^3),

	where the constants in  in (<ref>),(<ref>) have at most a polynomial growth with respect to x. Here ℒϕ=ϕ'f+σ^2/2 Δϕ denotes generator of the SDE and 𝒜_1,𝒜_1 are linear differential operators with smooth coefficients. 
	Note that 𝒜_1ℒ^2/2 in general (otherwise the scheme has weak order 2).
	If the  condition (A_1 + [,𝒜_1])^* ρ_∞ =0 holds,
	equivalently,
		
    ⟨𝒜_1ϕ +[,𝒜_1] ϕ⟩ = 0

		for all test function ϕ, 
		where we define ⟨ϕ⟩ = ∫_ℝ^dϕρ_∞ dx,
		then it is shown in <cit.> that 
		X_n has order two for the invariant measure, i.e. the convergence estimates (<ref>) with r=2 and (<ref>) hold.
	Before we can apply the above result, the following lemma 
	allow to compute the weak Taylor expansion of the modified scheme.
	 
	Consider the scheme (<ref>) with modified stage (<ref>)
	and assume the hypotheses of Theorem <ref>.
	Then (<ref>) holds where
	the linear differential operator 𝒜_1 is given by
		
    𝒜_1ϕ   =   1/2ϕ”(f,f)  + σ^2/2∑_i=1^dϕ”'(e_i,e_i,f) + σ^4/8∑_i,j=1^d ϕ^(4)(e_ie_i,e_j,e_j)
    		+ c_2 ϕ'f'f 
       +   c_3 σ^2/2ϕ' ∑_i=1^df”(e_i,e_i) + c_4 σ^2 ∑_i=1^d ϕ”(f'e_i,e_i),

		where f=-∇ V(x) and
		
    c_2 = ω_1^2T_s”(ω_0)/2 T_s(ω_0),    
    		c_3 = r_s+ω_0/sω_1α,    
    		c_4 = T_s”(ω_0)ω_1/T_s'(ω_0) + ω_1/2.

	
	
	Adapting the proof of Lemma <ref>, the internal stage K_i defined in (<ref>) (and (<ref>) for i=1)
		satisfies (<ref>) where h^2R_i,h(X_0) can be replaced by
		
    ω_1^2T_i”(ω_0)/2 T_i(ω_0) h^2f'(X_0)f(X_0) 
    		+ r̃_i σ^2/2 f”(X_0)(ξ_n,ξ_n) + h^5/2R_i +  h^3 R_i,h(X_0),

		where (R_i)=0 and all the moments of R_i,R_i,h(X_0) are bounded with polynomial growth with respect to X_0. Here, r̃_i is defined by induction as r̃_0=0, 
		r̃_1=Δ_1+α, and 
		
    r̃_i=ν_ir̃_i-1+κ_ir̃_i-2+Δ_i,   i=2,…,s.

		
			
			
			
			
			
			We have (R_i)=0 because R_i is a linear combination of
			f'(X_0)f'(X_0)ξ_n, f”(X_0)(f(X_0),ξ_n), and f”'(X_0)(ξ_n,ξ_n,ξ_n) with zero mean values (recall that odd moments of ξ_n vanish).
			Next, observing that the difference d_i=r̃_i-r_i satisfies d_0=0,d_1=α,
			and 
			d_i=ν_id_i-1+κ_id_i-2,i=2,…,s,
			we deduce 
    r̃_i=r_i+d_i,     d_i=U_i-1(ω_0)/T_i(ω_0)ω_0α  ∀  i=0,..,s.

			In particular, taking i=s in (<ref>),(<ref>),
			and expanding (<ref>), we deduce that (<ref>) holds
			with c_2,c_3,c_4 defined in (<ref>) where we note that c_3=r_s=r_s+d_s.
	
	
	
	Following the proof of <cit.> 
	(see also <cit.>)
	where we apply repeatedly integration by parts for the integral in (<ref>), 
	using Lemma <ref> for the expression of 
	𝒜_1, we deduce that the quantity in (<ref>)
	satisfies
		
    ⟨𝒜_1ϕ +[,𝒜_1]ϕ⟩ = ∑_i=1^d⟨  (c_3-c_2-c^2) σ^2/2ϕ' f”(e_i,e_i) + (c_4-1/4 - c_2/2 -c^2) σ^2 ϕ”(f'e_i,e_i) ⟩,

		where we use 
		[,𝒜_1]=
	-c^2σ^2(1/2 ϕ' ∑_i=1^df”(e_i,e_i) +  ∑_i=1^d ϕ”(f'e_i,e_i) ) for 𝒜_1ϕ =c^2σ^2/2 Δϕ.
		We see that the above quantity (<ref>) vanished if c_3-c_2-c^2=c_4-1/4 - c_2/2 -c^2=0, equivalently,
		
    c_3-c_2  = c^2 = c_4-1/4 - c_2/2.

		For the values of α,c defined in (<ref>), we obtain that (<ref>) indeed holds and we deduce that the order two condition (<ref>) for the invariant measure is satisfied. This concludes the proof.
	
	
	
		
		
		
		
		
		
		
		
		
		
		
	
	
	
		


§ NUMERICAL EXPERIMENTS

	








	
	In this Section, we illustrate numerically our theoretical analysis and we show the performance of the proposed SK-ROCK method and its postprocessed modification PSK-ROCK.
	
	

 §.§ A nonlinear nonstiff problem

	
	We first consider the following non-stiff nonlinear SDE,
	
    dX=(1/4X+1/2√(X^2+1))dt+√(1/2(X^2+1))dW,   X(0)=0.

	whose exact solution is	X(t)=sinh(t/2+W(t)/√(2)).
	In Figure <ref>, we consider the SK-ROCK method (<ref>) and
	plot the strong error (|X(T)-X_N| and the weak error |(arcsinh(X(T))-(arcsinh(X_N))|
	
	at the final time T=Nh=1 using 10^4 samples and number of stages s=1,5,10,100.
	We obtain convergence slopes 1 and 1/2, respectively, which confirms Theorem <ref> stating the strong order 1/2 and weak order 1 of the proposed scheme. 
	Note that s=1 stage is sufficient for the stability of the scheme in the non-stiff case. The results for s=5,10,100 yield nearly identical curves which illustrates that the error constants of the method are nearly independent of the stage number of the scheme.
	
	
	
	

 §.§ Nonlinear nonglobally Lipschitz stiff problems

	
	
	Consider the following nonlinear SDE in dimensions d=2 with a one-dimensional noise (d=2,m=1). This is a modification of a one-dimensional population dynamics model from <cit.> considered in <cit.> for testing stiff integrator performances,
	
    dX   =(ν(Y-1)-λ_1 X(1-X))dt-μ_1X(1-X)dW,   X(0)=0.95,
    
    	dY   =-λ_2Y(1-Y)dt-μ_2Y(1-Y)dW,   Y(0)=0.95.

	Observe that linearizing (<ref>) close to the equilibrium (X,Y)=(1,1), we recover for ν=0 the scalar test problem (<ref>). In Figure <ref>
	we consider the SK-ROCK method applied to (<ref>) with parameters that are identical to those used in <cit.>. We take
	the initial condition X(0)=Y(0)=0.95 close to this steady state and use the parameters ν=2,μ_2=0.5,λ_2=-1.
	In a nonstiff regime (-λ_1=μ_1=1 in Figure <ref>(a)), we observe a convergence slope 1 for the second moment (X(T)^2)
which illustrates the weak order one of the scheme, although our analysis in Theorem <ref> applies only for globally Lipschitz vector fields. The stage number s=1 is sufficient for stability, but we also include for comparison the results for s=10 (note that the results for s=50,100 not displayed here are nearly identical to the case s=10). The convergence curves are obtained as an average over 10^6 samples. 

	In a stiff regime (-λ_1=μ_1^2 =100 in Figure <ref>(b)), we observe for the standard small damping η=0.05 a stable but not very accurate convergence, due to the severe nonlinear stiffness.
However, considering a slightly larger damping η=4, in the spirit of the S-ROCK method, yields a stable integration for all considered timesteps and all trajectories and we observe a line with slope one for the SK-ROCK method. Here, given the timesteps h, the numbers of stages s are adjusted as proposed in (<ref>) where λ_max=|λ_1|=100.

	
	
	
	
	 
		For severely stiff problems, 
		alternatively to switching to drift-implicit schemes <cit.>,
		one can consider in SK-ROCK a slightly larger damping η and the corresponding stage parameter s below, similar to (<ref>) and chosen such that the mean-square stability domain length (<ref>) satisfies L > hλ_max,
		
    s=[√(hλ_max+1.5/2Ω(η))+0.5],
 
		where Ω(η) is given in Lemma <ref>.
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
	
	
	
	
	

 §.§ Linear case: Orstein-Uhlenbeck process

	We now illustrate numerically in details the role of the postprocessor introduced in Theorem <ref>
	for the linear Orstein-Uhlenbeck process in dimension d=m=1,
	
    dX=-λ Xdt+σ dW,  
    	X(0)=2

	where we choose λ=1 and σ=√(2).
	
	
	
	
In Figure <ref>, we consider the SK-ROCK and PSK-ROCK methods with s=1,5,10,100 stages, respectively.
	For a short time T=0.5 (Fig. <ref>(a)(b)), we observe weak convergence slopes one for both SK-ROCK and PSK-ROCK (second moment (X(T)^2))  as predicted by Theorem <ref>, and the postprocessor has nearly no effect of the errors.
	For a long time T=10 where the solution of this ergodic SDE is close to equilibrium, we observe that the weak order one of SK-ROCK (Fig. <ref>(c)) is improved to order two using the postprocessor in PSK-ROCK (Fig. <ref>(d)), which confirms the statement of Theorem <ref> that the postprocessed scheme has order two of accuracy for the invariant measure.
	
	For comparison, in Figure <ref>, we also include the results of PSK-ROCK without damping (η=0) using M=10^8 samples. We recall that for the scalar linear  Orstein-Uhlenbeck process, the PSK-ROCK method with zero damping is exact for the invariant measure (see Section <ref>).
	We observe only Monte-Carlo errors with size ≃ M^-1/2=10^-4, which confirms that the PSK-ROCK method has no bias at equilibrium for the invariant measure in the absence of damping, as shown in (<ref>).
	We emphasise however that this exactness results holds only for linear problems, and a positive damping parameter η should be used for nonlinear SDEs for stabilization, as shown in Sections <ref> and <ref>.	
	



 §.§ Nonglobally Lipschitz Brownian dynamics



To illustrate the advantage of the PSK-ROCK method applied to nonglobally Lipschitz ergodic Brownian dynamics, we next consider the following double well potential V(x)=(1-x^2)^2/4 and the corresponding one-dimensional Brownian dynamics problem
	
    dX=(-X^3+X)dt+√(2)dW,     X(0)=0,

	In Figure <ref>, we compare the performances of S-ROCK, 
	S-ROCK2 considered in <cit.> (a method with weak order 2 for general SDEs), and the new SK-ROCK and PSK-ROCK methods at short time T=0.5 (Figures <ref>(a)(b)) and long time T=10 (Figures <ref>(c)(d)). 
	As we focus on invariant measure convergence and not on strong convergence, we consider here discrete random increments
	with (ξ_n=±√(3))=1/6,(ξ_n=0)=2/3, which has the correct moments so that 
	Theorem <ref> remains valid. Our numerical tests indicate that it makes
	PSK-ROCK with modified stage (<ref>) more stable. For a fair comparison, we use
	the same discrete random increments for all schemes.
	We plot the second moment error versus the time stepsize h and versus the average cost which is the total number of function evaluations during the time integration divided by the total number number of samples. Indeed, the number of function evaluations depends on the trajectories because the stage parameter s is adaptive at each time step. 
	For short time, we can see that the S-ROCK and the SK-ROCK method have order 1 (Figure <ref>(a)) and exhibit similar performance with nearly identical error versus cost curves in Figure <ref>(b), while PSK-ROCK 
is less advantageous for short time. This illustrates that the postprocessing has no advantage for short times. The S-ROCK2 method is the most accurate for small time steps, and it has order 2 as shown in Figures <ref>(a)(c), but at the same time it has a larger average cost as observed in Figures <ref>(b)(d) due to its smaller stability domain with size ≃ 0.42 · s^2.

For long time, the SK-ROCK and S-ROCK both exhibit order 1 of accuracy (Figure <ref>(c)), with an advantage in terms of error versus cost for the
SK-ROCK method that is about 10 times more accurate for large time steps. In contrast, the postprocessed scheme PSK-ROCK exhibits order 2 of convergence (Figure <ref>(c)) which corroborates Theorem <ref>. Since the postprocessing overcost is negligible (two additional vector field evaluations per timestep due to the modified stage K_1 in (<ref>)), this makes PSK-ROCK the most efficient in terms of error versus cost, as shown in Figure <ref>(d). The S-ROCK2 method has order 2 here but with poor accuracy compared to the PSK-ROCK method with approximately the same cost. Note that typically the SK-ROCK method used s=1,2,3 stages in contrast to the S-ROCK method using s=2,…,6 stages per timesteps. 
	
	
	

 §.§ Stochastic heat equation with multiplicative space-time noise
	
	Although our analysis applies only to finite dimensional systems of SDEs, we consider the following stochastic partial differential equation (SPDE) obtained by adding multiplicative noise to the heat equation,
	
    u(t,x)t   =u(t,x)x + u(t,x)Ẇ(t,x),   
    	(t,x)∈[0,T]×[0,1]
    
    	u(0,x)   =5cos(π x),   x∈[0,1],
    
    	u(t,0)   =5,  u(t,1)x=0,   t∈[0,T],

		where Ẇ(t,x) denotes a space-time white noise that we discretize together with the Laplace operator with a standard finite difference formula <cit.>.
	We obtain the following stiff system of SDEs where
u(x_i, t) ≈ u_i(t), with x_i=i Δ x, Δ x = 1/N,

    du_i    =   u_i+1-2u_i+u_i-1/Δ x^2dt+ u_i/√(Δ x)  dw_i,    i=1,…,N,

where the Dirichlet and the Neumann conditions impose u_0=5 and u_N+1=u_N-1, respectively. 
Here, w_1,…, w_N are independent standard Wiener processes and dw_i indicates Itô noise. 
	
		
		
		
		
	
	In Figure <ref>(a), we plot one realization of the SPDE using space stepsize Δ x=1/100 and timestep size Δ t=1/50. Note that the Lipchitz constant associated to the space-discretization of (<ref>) has size ρ=4Δ x^-2, and the stability condition is fulfilled for s=22 stages. For comparison, the standard S-ROCK method would require s=46 stages, while applying the standard Euler-Maruyama with a smaller stable timestep Δ t/s would require s ≥Δ t ρ/2=400 intermediate steps.
	Notice that the initial condition in (<ref>) satisfies the boundary conditions, which permits a smooth solution close to time t=0. Taking alternatively an initial condition that does not satisfy the boundary conditions (for instance u(x,0)=1) yields an inaccurate numerical solution with large oscillations close to the boundary x=0. A simple remedy in such a case is to consider a larger damping parameter η, as described in Remark <ref>.
	
		
		
		
		
		
	
	
		
	
	
		
		
			
			
			
			
		
		
			
			
			
			
			
		
		
			
			
			
		
	
	In Figure <ref>(b), we compare the number of vector field evaluations of the standard S-ROCK and new SK-ROCK methods when applied to the SPDE (<ref>) with finite difference discretization with parameter Δ x=1/100. The better performance of SK-ROCK with damping η=0.05 is due to its larger stability domain with size ≃ 1.94· s^2 compared to the size ≃ 0.33 · s^2 for S-ROCK. Observing the ratio of the two costs in Figure <ref>(b), we see that the new SK-ROCK methods has a reduced cost for stabilization by an asymptotic factor of about √(1.94/0.33)≃ 2.4 for large s and large stepsizes, which confirms the stability analysis of Section <ref>.
	The convergence analysis of the SK-ROCK method for the stochastic heat equation is the topic of future work.
	
Notice that SK-ROCK with s=1 stage has the optimal mean-square stability length (L=2 for η=0) as defined in (<ref>). In contrast, the S-ROCK method with s=1 has the smaller stability length L=3/2, while the standard Euler-Maruyama has L=0.
This explains why for the smallest considered stepsize Δ t=2^-15 in Figure <ref>(b), we have s=1 for SK-ROCK while 
S-ROCK uses s=2 stages.
	
	
	In Figure <ref> we consider again one realization with SK-ROCK of the SPDE problem (<ref>) but with a different initial condition
	u(0,x)=1 not fulfilling the boundary conditions, i.e. that is outside the domain of the Laplace operator, as considered in <cit.>. 
	We compare the result for the same sets of random numbers but for 
	for different values of the damping parameter η. We observe numerically high oscillations in time and space for the small damping value η=0.05  in Figure (<ref>) while the larger damping η=10 yields a smoother solution in Figure (<ref>).
	This illustrates again Remark <ref> showing that the damping parameter η can be increased in the case of
	severely stiff problems, adjusting the stage parameter accordingly with (<ref>).
	
	


Acknowledgements.The work of the first author was partially supported by the Swiss National Foundation, grant
200020_172710. The work of the second and third authors was partially supported by the Swiss National Science Foundation, grants 200020_144313/1 and 200021_162404. 
The computations were performed at University of Geneva on the Baobab cluster.
	

	
	
	
	
	
		
		
		
		
			
			
			
			
			
						
			
		
		


		
		
			
			
			
			
			
			
			
		
		
		
		
		
		
			
		
		
	
	
	
	
	
abbrv

	
	
	
