Continuous stage stochastic Runge–Kutta methods

In this work, a version of continuous stage stochastic Runge–Kutta (CSSRK) methods is developed for stochastic differential equations (SDEs). First, a general order theory of these methods is established by the theory of stochastic B-series and multicolored rooted tree. Then the proposed CSSRK methods are applied to three special kinds of SDEs and the corresponding order conditions are derived. In particular, for the single integrand SDEs and SDEs with additive noise, we construct some specific CSSRK methods of high order. Moreover, it is proved that with the help of different numerical quadrature formulas, CSSRK methods can generate corresponding stochastic Runge–Kutta (SRK) methods which have the same order. Thus, some efficient SRK methods are induced. Finally, some numerical experiments are presented to demonstrate those theoretical results.


Introduction
Stochastic differential equations (SDEs) have wide applications in many disciplines like biology, economics, medicine, engineering and finance (see, e.g., [1][2][3]). However, most SDEs arising in practice are nonlinear, and cannot be solved explicitly. There has been tremendous interests in developing effective and reliable numerical methods for SDEs during the last few decades, for example see [4][5][6][7][8][9][10][11][12][13][14]. Runge-Kutta (RK) methods with continuous stage were firstly presented by Butcher in 1970s [15], and they have been investigated and discussed by several authors recently because of the great advantages in conserving symplecticity [16], preserving energy [17] and so on. So constructing continuous stage stochastic Runge-Kutta (CSSRK) methods for SDEs is a valuable task. This paper mainly aims to construct CSSRK methods for SDEs which can be written in integral form as with a d-dimensional Winer process (W (t)) t≥0 = ((W 1 (t), . . . , W d (t)) T ) t≥0 . Here g k : R m → R m , k = 1, . . . , d, are Borel measurable, satisfying Lipschitz condition and linear growth condition, x(t) is the stochastic process which denotes the solution of (1.1).
In engineering, the convergence order of applied numerical methods plays an important role in efficient operation. And B-series appear as a fundamental tool to do error analysis on a wide range of problems. Therefore, a general order theory for CSSRK methods based on stochastic B-series and multicolored rooted tree is developed. Then, for the above SDE (1.1) with d = 1, we discuss the order conditions, and construct relevant CSSRK methods of order 1.0. Furthermore, CSSRK methods are also constructed for problems with some special structures, such as single integrand SDEs, SDEs with additive noise.
Single integrand SDEs are given by where λ ∈ {0, 1} and σ k ∈ R, k = 1, 2, . . . , d, are given constants. Some well known examples for single integrand SDEs are the SDE describing fatigue cracking [18], the stochastic Van der Pol equation [19]. A theory for the derivation of B-series of the numerical approximation by CSSRK methods is given. Due to the single integrant, it is similar to the ordinary differential equations (ODEs) case [20]. The deterministic continuous stage Runge-Kutta (DCSRK) methods of deterministic order p d are corresponding with the CSSRK methods of mean square p μ = p d /2 , which can be viewed as the stochastic generalization of the DCSRK methods.
Concerning the specific problems with additive noise, thus where σ k , i = 1, 2, . . . , d, are constants. Based on stochastic B-series, we propose a class of CSSRK methods, totally derivative-free for (1.3), which are able to reach mean-square order 1.5. Then we focus on the special second-order systems with additive noise of the following form: the corresponding CSSRK methods are able to reach mean-square order 2.0. Interestingly, it is verified that a great deal of SRK methods could be derived from the CSSRK methods by means of the adoption of different numerical quadrature formulas. Furthermore, the conditions on algebraic precision of the numerical quadrature formulas and the degree of the coefficient polynomials of CSSRK methods are analyzed to guarantee that SRK methods share the same order with the CSSRK one. Thus, some efficient SRK methods are obtained.
The outline of this paper is organized as follows. In the next section, a family of CSSRK methods is proposed. Based on the multicolored rooted tree and stochastic B-series, the general order conditions are calculated. Up to strong global order 1.0 of SRK methods are proposed in Sect. 3. In Sect. 4, we consider a particular class of Stratonovich SDEs, single integrand SDEs. Sections 5 is devoted to discuss the CSSRK methods solving SDEs with additive noise. Finally, numerical experiments are reported in Sect. 6.

CSSRK methods and their order conditions
In this section, we first introduce the definition of CSSRK methods, then show some results which are useful in receiving the order theory.

CSSRK methods
We consider the following methods given by y 0 = x 0 and the step size h > 0, let A (i,k) τ ,ξ be a polynomial in τ and ξ , B (i,k) τ be a polynomial in τ . The one-step method h : y 0 → y 1 given by τ , then another presentation can be obtained, the coefficients Z (k) τ ,ξ and z (k) τ include random variables that depend on the stepsize h. Then it defines a m-dimensional approximation process y n with y n ≈ x n which called CSSRK methods.

Order theory
The basic tool of constructing our numerical methods is the multicolored rooted tree theory in [21]. Thus, we briefly list some definitions and theorems used in constructing numerical schemes later. Our first goal is to find B-series representations of (2.2a)-(2.2b) which can be written as a B-series, where T is the set of multicolored rooted tree, α(t) are combination terms, the elementary weight functions φ(t)(h) are stochastic integrals or random variables, and F(t)(x 0 ) are the elementary differentials.

Definition 2.1 (Trees and combinatorial coefficients) The set of multicolored rooted trees
is recursively defined as follows: (a) The graph • k = [∅] k with only one vertex of color k belongs to T k . Let t = [t 1 , t 2 , . . . , t l ] k be the tree formed by joining the subtrees t 1 , t 2 , . . . , t l each by a single branch to a common root of color k.
Thus, T k is the set of trees with a k-colored root, and T is the union of these sets. Further, we define α(t) as where r 1 , r 2 , . . . , r q count equal trees among t 1 , t 2 , . . . , t l .

Definition 2.2 (Elementary differentials)
For a tree t ∈ T, the elementary differential is a mapping The next lemma proves that, if Y (h) can be written as a B-series, then f (Y (h)) can be written as a similar series. The lemma is fundamental for deriving B-series for the exact and the numerical solution.
where U f are the collection derived from T.

Lemma 2.2 ([21])
The exact solution x(h) of (1.1) can be written as a B-series A similar result can be found for the numerical solution of (1.1) by the CSSRK methods (2.2a)-(2.2b).

Theorem 2.1 The numerical solution y 1 as well as the continuous stage values can be written in terms of B-series
Proof Write Y τ as a B-series, that is, for τ ∈ [0, 1], use the definition of the method (2.2a) and the term by term comparison yielding (2.2a) to obtain where t = [t 1 , t 2 , . . . , t l ] k ∈ T k . The proof of (2.5) is similar. Now, the local order of accuracy of the CSSRK methods can be decided by comparing the B-series of the exact and the numerical solution. First, we need to define the tree order.

Definition 2.3
The order of a tree t ∈ T is defined by The following lemma relates the global order of accuracy to the local order. Here, we assume that (2.5) is constructed such that ϕ(t)(h) = O(h ρ(t) ) for all t ∈ T for mean square convergence.

Lemma 2.3 ([21]) The method has mean square global order p if
Here, the O(·)-notation refers to h → 0 and, especially in (2.6), to the L 2 -norm.

Construction of CSSRK methods for SDEs with one random variable
Without loss of generality, we restrict consideration to autonomous systems (the coefficients of SDEs do not depend on t explicitly) in this part. We write then the CSSRK methods for system (3.1) are given by Based on the multicolored rooted tree theory in Sect. 2, we are able to obtain a set of order conditions guaranteeing that CSSRK methods (3.2a)-(3.2b) of mean-square order 1.0 as detailed below. A list of all trees with ρ(t) ≤ 1.5 and their corresponding functions are given in Table 1. Above all, the following conditions need to be satisfied: . Through calculation, if the coefficients of CSSRK methods (3.2a)-(3.2b) satisfy the conditions 3) Table 1 Multicolored rooted trees for (3.2a)-(3.2b) with order less than or equal to 1.5 then it is of mean-square order 1.0. We can choose some polynomials satisfying (3.3), for example It is almost mandatory that the practical implementation of the CSSRK methods (3.2a)-(3.2b) uses a numerical quadrature formula. By applying the quadrature formula where which can be formulated by the following Butcher tableau: The following result implies that we can construct SRK methods with the same order of CSSRK methods with the help of suitable numerical quadrature formulas.  3 -order polynomial, based on algebra analysis, the numerical quadrature formulas (b i , c i ) s i=1 with accuracy up to degree max{m 2 , m 1 + m 3 } are exactly established for (3.3). Then we have thus, the associated SRK methods (3.9) are of order 1.0.
As an application, for the constructed CSSRK methods of mean-square order 1.0, if we use more quadrature points then we can get lots of new SRK methods which are not found in the literature. For example, it gives the SRK schemes shown in the following.

High order CSSRK methods for single integrand SDEs
In this section, we consider a particular class of Stratonovich SDEs, single integrand SDEs (1.2). All results in this section can straightforwardly be considered as the case of onedimensional Wiener process, Then (4.1) is solved by (3.2a)-(3.2b), we can obtain Due to the following lemmas, the main result is that the B-series of the exact and the numerical solution are exactly as in the ODEs case, with the exception that integration is now performed with respect to μ instead of h.
Proof Write Y τ as a B-series, that is,
In the following, the main results of this section are presented.
and all elementary differentials F(t) fulfill the linear growth condition. Assume that μ is either sampled from the exact distribution or chosen such that at least the first 2p μ + 1 moments coincided with the ones of the exact distribution.

Due to E(μ(h) 2ρ(t) ) = O(hρ (t) ) by Lemma 4.3, (4.3) is then by Lemma 4.2 and Theorem 4.1 automatically fulfilled for all t ∈ T withρ (t) ≥ 2p μ + 1, and satisfying the remaining trees if and only if
, t ∈ T,ρ(t) ≤ 2p μ . (4.5) Note that (4.5) is just the condition that for the order p d of the DCSRK methods applied to a deterministic system (σ = 0) with p d = 2p μ . Similarly, (4.4) is automatically fulfilled for all t ∈ T witĥ satisfying the remaining trees if and only if Thus, the methods will be mean-square consistent of order p μ if the deterministic order is or, vice versa, the methods of deterministic order p d will converge with mean-square order of p d 2 .
As is well known, when solving the deterministic differential equation, the DCSRK method with the polynomial is of mean-square order 4.0 [22]. Hence, the corresponding CSSRK method is of order 2.0. We can also choose multiple numerical quadrature formulas to get different SRK methods of high order.

CSSRK methods for SDEs with additive noise
Concerning the specific problems with additive noise, thus (1.3) is solved by (2.1a)-(2.1b) defined by Because the noise terms are additive, we find that a elementary differential vanishes if its multicolored rooted tree contains a node following a stochastic node directly except if the deterministic node t 0 is the only succeeding end node. Table 2 consists of multicolored rooted trees whose elementary differentials are nonzero. We are able to obtain a set of order conditions guaranteeing that the CSSRK methods (5.1a)-(5.1b) obtained mean-square order 1.5 as detailed below. Table 2 Multicolored rooted trees for (5.1a)-(5.1b) with order less than or equal to 2.0   (2) τ ,ξ of the CSSRK methods (5.1a)-(5.1b) satisfy conditions then they are of mean-square order 1.5.
The condition (f ) in Theorem 5.1 is derived from the tree 6 in Table 2, which is essential for general system. However, for system (5.8), the corresponding elemental differential of tree 6 is because the nth component of it is where (f (x)) n means the nth element of the vector function f (x). Thus tree 6 is unnecessary in this situation. With the same analysis, the elementary differential of tree 5 in Table 2 vanishes as well. Moreover, for the additional trees 7, 8 in Table 3, we can check that

Numerical experiments
In this section, we perform numerical tests to verify the mean-square convergence order proposed in Sects. 3, 4 and 5. The first group consists of method (3.4)-(3.7) and (3.11), (3.14) of different free parameters. The second group of methods consists of the methods (4.6) of deterministic order p d = 4.0 for solving deterministic differential equation, in which case the predicted stochastic order p μ = 2. Finally, we consider methods (5.9) for SDEs with additive noise. In each example, the solution is approximated with step sizes 2 -5 , . . . , 2 -9 and the sample average of M = 2000 independent simulated realizations of the absolute error at the terminal time T = 1 is calculated in order to estimate the expectation, accordingly.

Kubo stochastic oscillator
The Kubo stochastic oscillator in the sense of Stratonovich SDE is described by (6.8) Here a and b are constants. The general solution of (6.8) is given by We choose the coefficients of (6.8) as a = 0, b = 1 with initial value p 0 = 0.5, q 0 = 0. By using appropriate numerical quadrature formula, we can get the classical SRK methods of high order. For (4.6), we can choose Gaussian quadrature with 2 nodes to get the following method of order 2.0: . (6.9) The average sample errors at terminal time T = 1 are expressed by  Table 5. Figure 2 shows the results of Table 5 in a log-log plot.

Linear stochastic oscillator with additive noise
Then we consider the following linear stochastic oscillator with additive noise: where x is the position and y is the velocity of a particle under the simple harmonic restoring force. Since (6.10) is a SDE with additive noise, its Itô and Stratonovich form are identical. We consider it as a SDE of Stratonovich type. Choose the free parameters of (5.3) Table 5 The endpoint average sample errors of method (4.6) and (6.9) for solving (6.8)   To check the order of method (6.11), (6.12), the average sample errors at terminal time T = 1 (i.e., 2000 i=1 |x(1, ω i )x N (ω i )| 2 + |y(1, ω i )y N (ω i )| 2 /2000) for method (6.11), (6.12) in Table 6. Figure 3 shows the results of Table 6 in a log-log plot and the slope of the reference line is 2.0.

Conclusions
This paper presents the extension work of DCSRK methods which are upgraded to the stochastic counterpart. For general autonomous SDEs in the Stratonovich sense, a class of CSSRK methods are proposed. The general order conditions are obtained by the use of colored rooted tree and stochastic B-series. The order conditions up to strong global order 2.0 of the CSSRK methods for solving single integrand SDEs are simplified. And for solving SDEs with additive noise, the CSSRK methods of strong global order 1.5 and 2.0 are given. It can be proved that a series of SRK methods of the same order with one CSSRK method can be obtained by the rational application of numerical quadrature formulas.
The numerical experiments are examined to verify the results of our theoretical analysis and show that the schemes are useful in very long-time numerical simulation of SDEs. The properties of more numerical SRK methods can be obtained by the CSSRK methods with all of these numerical experiments demonstrating significant order convergence of the numerical methods.