High order numerical integrators for single integrand Stratonovich SDEs

We show that applying any deterministic B-series method of order p d with a random step size to single integrand SDEs gives a numerical method converging in the mean-square and weak sense with order (cid:3) p d / 2 (cid:4) . As an application, we derive high order energy-preserving methods for stochastic Poisson systems as well as further geometric numerical schemes for this wide class of Stratonovich SDEs. © 2020 The Author(s). Published by Elsevier B.V. on behalf of IMACS. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

Stochastic differential equations (SDEs) of this form arise, e.g., when the right-hand side of an ordinary differential equation (ODE) model is randomly perturbed. Interesting applications are: Example 1 (Fatigue cracking [34]).
Example 2 (Stochastic Hamiltonian systems [25,35]). dX(t) = J −1 ∇ H(X(t)) ( dt + c • dW (t)) , Example 3 (Stochastic perturbations of Poisson systems [9]). dX(t) = B(X(t))∇ H(X(t)) ( dt + c • dW (t)) , In the present communication, we generalize the main results in [14] from Runge-Kutta methods to B-series methods. In particular, this permits to show that when applying any (deterministic) B-series numerical integrators of order p d to (1) by replacing the deterministic step size h in the numerical scheme by λh + σ (W (t + h) − W (t)) when calculating one step of the approximation on the time interval [t, t + h] respectively, one then automatically obtains a time integrator with step size h of mean-square as well as weak order p d /2 for the SDE (1). We recall that the notation x denotes the floor function of a real number x.
On top of that, it can easily be observed that, in general, when a (deterministic) numerical integrator possesses some geometric properties, then the same geometric properties hold for the time integrator applied to single integrand SDEs (1). This observation thus allows to carry forward various results already obtained in the literature as well as new ones, see for instance Section 3 below.

Convergence of B-series methods applied to single integrand SDEs
B-series for solutions to SDEs and their numerical solution by stochastic Runge-Kutta methods have been developed in [4,5] to study strong convergence in the Stratonovich case, in [32,33] for strong converge in the Itô as well as the Stratonovich case, in [20] and [21] to study weak convergence in the Stratonovich case and in [30,31] to study weak convergence in both the Itô and the Stratonovich case. A uniform and self-contained theory for the construction of stochastic B-series for the exact solution of SDEs and its numerical approximation by stochastic Runge-Kutta methods is given in [11]. Based on the notation used there, in [14] the B-series for the exact solution and Runge-Kutta approximations of singleintegrand SDEs were derived. For convenience, we summarize the results we will need in the following. Due to the single integrand we do, similar to the ODE case [6], only need non-colored trees in the expansion of the solution.

Definition 1 (Trees).
The set of rooted trees T related to single-integrand SDEs is recursively defined as follows: a) The empty tree ∅ and the graph • = [∅] with only one vertex belong to T .
We are now able to state the B-series of the exact solution to the SDE (1), see also [11,32]. In the following, ρ(τ ) denotes the number of nodes in a tree τ . Theorem 1 ([14]). Let γ : T → N be given by Then the solution X(t+h) of (1) starting at the point (t, x) can be written as a B-series That is, the B-series of the exact solution of (1) coincides with the one of the exact solution of an ODE (i.e. when λ = 1, Next, we derive the B-series for the numerical solution of the single integrand SDE (1). To do this, we consider B-series methods for the ODE d X(t) = f (X(t)) dt. Such methods can be written as [15] consider therefore the following B-series method for solving (1):

Definition 4.
Given an ODE-B-series method (4) for the ODE d X(t) = f (X(t)) dt, the corresponding ODE-B-series method for (1) is given by where Now we give the definitions of both weak and strong convergence used in this note. Rd) fulfilling a polynomial growth condition [19] and Ih be the discretized time interval on which the numerical approximations are calculated.

Definition 5.
A time discrete approximation Y = (Y (t)) t∈Ih converges weakly with order p to X at time t ∈ Ih as the maximum step size h → 0 if for each g ∈ C 2(p+1) P (R d , R) there exist a constant C g and a finite δ 0 > 0 such that | E(g(Y (t))) − E(g(X(t)))| ≤ C gh p holds for each h ∈ ]0, δ 0 [. Whereas weak approximation methods are used to estimate the expectation of functionals of the solution, strong approximation methods approach the solution path-wise. In this article, next to weak convergence, we will consider mean-square instead of strong convergence. Definition 6. A time discrete approximation Y = (Y (t)) t∈Ih converges in the mean-square sense with order p to X at time t ∈ Ih as the maximum step size h → 0 if there exist a constant C and a finite δ 0 > 0 such that E( Y (t) − X(t) 2 ) ≤ Ch p holds for each h ∈ ]0, δ 0 [ . Observe that, by Jensen's inequality, mean-square convergence implies strong convergence of the same order. We are now in position to state the main result of the present publication.

Theorem 2. Assume that the B-series method (5) is of deterministic order p d and let p
and f and f f fulfill a Lipschitz condition. Finally, assume • for mean-square convergence, that all elementary differentials F (τ ) fulfill a linear growth condition, • respectively for weak convergence, that f ∈ C Then this very same B-series method is of mean-square as well as weak order p μ when applied to the single integrand SDE (1) with step size t,t+h μ. For weak convergence, it suffices that t,t+h μ is chosen such that at least the first 2p μ + 1 moments coincide with those of μ t (h), and all the others are in O(h p μ +1 ).
Proof. To prove this result, one first observes that the corresponding Theorem for Runge-Kutta methods [14, Theorem 1.1] only uses that the B-series of the exact and numerical solutions fulfil the properties summarized in Definition 4. One can then directly extend the proof to the present situation of B-series methods.

Applications to geometric numerical integration of single integrand SDEs
Geometric properties of numerical schemes for single integrand SDEs (1) are closely related to those of numerical schemes for their corresponding ODEs. This relation, which holds for B-series as well as non B-series (e.g. volume preserving) methods, can be specified furthermore: In principle, when applying any deterministic geometric numerical integrators to single integrand SDEs (1) with random step size t,t+h μ, one then obtains the same geometric property as by the corresponding deterministic numerical scheme. This is because the random perturbation in some sense respects the geometric structure of the phase space. Hence, proofs of conservation properties in the deterministic setting can be adapted to single integrand Stratonovich SDEs (1). With the main idea used in this paper, one easily derives several results already obtained in the literature as well as new ones for geometric numerical integrators of single integrand SDEs (1): 1. If the deterministic numerical method preserves linear or quadratic invariants, so does it for single integrand SDEs. See the example below and e.g. [18]. 2. If the deterministic scheme is energy or invariant-preserving, then it is also energy or invariant-preserving for single integrand SDEs. See the example below and e.g. [9,22,23,26]. 3. If the deterministic numerical method is symmetric, then it is also symmetric for the SDE (1).

If the deterministic time integrator is symplectic and f
, then it is symplectic for single integrand SDEs, e.g. [25,27,35]. 5. Deterministic (symmetric) projection methods can be directly applied to single integrand SDEs with constraints, e.g.
[36]. 6. If the deterministic numerical method is volume preserving, then it is also volume preserving for single integrand SDEs.

If the deterministic numerical scheme is a Poisson integrator and f (x) = B(x)∇ H(x), where B(x) represents a Poisson
bracket, then it is also a Poisson integrator for single integrand SDEs.

The Hamiltonian
Note that the right hand side f of (6) is not globally Lipschitz continuous and does thus a priori not fulfill the conditions of Theorem 2. However, for methods that conserve the Casimir, we can replace f by a function being zero outside a suitable ball and fulfilling the conditions of Theorem 2, e.g. replacing f bỹ Note that f (X) = f (X) for all X 2 ≤ √ 2 X(0) and f (X) = 0 for all X 2 ≥ 2 X(0) (see e.g. [28]), and all derivatives of f remain bounded. A similar modification can be made for numerical methods that preserve the Hamiltonian of the deterministic rigid body.
As a starting deterministic method, we choose the linear energy-preserving integrators of orders 2, 4 and 6 from [10]. By Theorem 2, we thus expect a first, resp. second, resp. third strong and weak order energy-preserving numerical integrator for the stochastic rigid body (6), denoted by EPs1, EPs2, EPs3. Recall that these novel numerical methods are simply given by replacing the step size h by λh + σ W n in the deterministic scheme. The orders of convergence can be seen in Figs  Furthermore, as it can be seen in Fig. 3, for the example of the third order method, these numerical schemes not only preserve the energy but also the quadratic Casimir C (X). This follows from the discussion in Section 3, as the numerical integrators from [10] are known to preserve quadratic Casimirs.