Symplectic Integrators to Stochastic Hamiltonian Dynamical Systems Derived from Composition Methods

Abstract

“Symplectic” schemes for stochastic Hamiltonian dynamical systems are formulated through “composition methods (or operator splitting methods)” proposed by Misawa (2001). In the proposed methods, a symplectic map, which is given by the solution of a stochastic Hamiltonian system, is approximated by composition of the stochastic flows derived from simpler Hamiltonian vector fields. The global error orders of the numerical schemes derived from the stochastic composition methods are provided. To examine the superiority of the new schemes, some illustrative numerical simulations on the basis of the proposed schemes are carried out for a stochastic harmonic oscillator system.

1. Introduction

The theory of Hamiltonian dynamical systems is understood as one of the fundamental tools for the description of dynamical phenomena treated in physics, engineering, and economics. To investigate the behavior of such a system, numerical simulation is often carried out. Then, as one of the important subjects in the numerical approach, attention has been paid to how to construct numerical schemes preserving Hamiltonian dynamical structures, that is, “symplectic integration schemes" (Hairer et al. [1] and references therein). The symplectic schemes are superior to numerical realization of the phase flows of Hamiltonian dynamical systems on long time intervals, and hence one can apply them to various dynamical systems described in the framework of classical Hamiltonian mechanics.

As an extension of this topic, Milstein et al. [2] have proposed symplectic integrators for “stochastic Hamiltonian dynamical systems (e.g., Misawa [3])" which are governed by stochastic differential equations (SDEs) (Ikeda and Watanabe [4]). The stochastic Hamiltonian dynamical systems are viewed as “open" Hamiltonian systems perturbed by random fluctuations from an external world; from the viewpoint of finance, such systems may also be regarded as dynamical ones describing some risky assets. For the stochastic systems, they have formulated “symplectic schemes of Runge-Kutta type" which are analogous to those in deterministic cases. On the other hand, symplectic integrators for deterministic Hamiltonian systems are often produced through the so-called “composition methods (operator splitting methods)", which are formulated on the basis of Lie algebraic theory and Lie group theory (McLachlan and Quispel [5]). The method is outlined as follows: Let 𝑋 denote vector fields on some space with coordinates 𝑥, with flows exp(𝑡𝑋), that is, the solutions of differential equations of the form 𝑥̇(𝑡)=𝑋(𝑥) are given in the form 𝑥(𝑡)=exp(𝑡𝑋)(𝑥(0)). Assume that one can write 𝑋=𝐴+𝐵 in such a way that exp(𝑡𝐴) and exp(𝑡𝐵) can both be calculated explicitly. (More generally, this can be relaxed by approximations of the exponential maps.) In the most elementary case, the method gives the approximation for 𝑥(𝑡) through𝑥(𝑡)=exp(𝑡𝐴)exp(𝑡𝐵)(𝑥)=𝑥(𝑡)+𝑂(𝑡²); (1.1) the last equality is obtained by using the Baker-Campbell-Hausdorff (BCH) formula, which is well known in Lie algebraic theory. The advantage of this method is that geometric properties of the true flow exp(𝑡𝑋) are often preserved. If 𝑋, 𝐴, and 𝐵 are Hamiltonian vector fields, then both true flow and the approximating flow are symplectic. Therefore, it seems to be quite natural to derive symplectic integrators for stochastic Hamiltonian systems by using these methods.

The purpose of the present paper is to propose some symplectic numerical integration schemes for stochastic Hamiltonian dynamical systems. Then, on the basis of the author's results on composition methods for SDEs (Misawa [6]), the numerical approximation errors of the new stochastic schemes are estimated.

This paper is organized as follows. In Section 2, we first review the composition methods for SDEs and the theorem on error estimation to the schemes. In Section 3, some symplectic schemes by composition methods are proposed for stochastic Hamiltonian dynamical systems described by SDE with a one-dimensional Wiener process. The approximation error of numerical solutions through the schemes is also estimated. In Section 4, in order to examine the superiority of the new schemes, we compare the numerical solutions from the proposed schemes with those from the stochastic Euler-Maruyama scheme through an illustrative example. Some concluding remarks are given in the end of the section.

Finally, we would like to quote here another approach to composition methods for SDEs by K. Burrage and P. M. Burrage [7]. Furthermore, the other contributions related to our topic can also be found in the recent articles by Malham and Wiese [8] and Bou-Rabee and Owhadi [9].

2. Composition Methods for Numerical Integration of SDEs

Now, we start with a review of composition methods for numerical integration of SDEs proposed by Misawa [6]. Let (Ω,ℱ,𝑃) be a probability space. We consider an 𝑛-dimensional stochastic differential equation of the Stratonovich type associated with a one-dimensional Wiener process 𝑊 as follows:𝑑𝑆(𝑡)=𝑏(𝑆(𝑡))𝑑𝑡+𝑔(𝑆(𝑡))∘𝑑𝑊(𝑡), (2.1) where 𝑏 and 𝑔 are 𝑛-dimensional 𝐶^∞ vector functions, respectively. By using these functions, we define the differential operators 𝑋₀ and 𝑋₁ as𝑋₀=^𝑛_𝑖=1𝑏^𝑖𝜕_𝑖,𝑋₁=^𝑛_𝑖=1𝑔^𝑖𝜕_𝑖, (2.2) where 𝜕_𝑖=𝜕/𝜕𝑆^𝑖. According to Kunita [10], a solution of the SDE under an initial value 𝑆(0)=𝑠 can be formally represented in terms of an exponential map derived from a “stochastic" vector field 𝑌_𝑡 as follows:𝑆(𝑡)=(exp𝑌_𝑡)(𝑠), (2.3) where 𝑌_𝑡(𝜔) is the vector field for each 𝑡 and almost surely (a.s.) 𝜔∈Ω given by𝑌_𝑡=𝑡𝑋₀+𝑊(𝑡)𝑋₁+12(^𝑡₀𝑢𝑑𝑊(𝑢)−^𝑡₀𝑊(𝑢)𝑑𝑢)[𝑋₀,𝑋₁]+_{𝐽;3≤|𝐽|}{^∗_Δ𝐽𝑐_Δ𝐽𝑊^Δ𝐽(𝑡)}𝑋^𝐽, (2.4) and 𝑋^𝐽=[⋯[𝑋_𝑗1,𝑋_𝑗2]⋯]𝑋_𝑗𝑚](𝐽=(𝑗₁,…,𝑗_𝑚)), where [𝑋₀,𝑋₁] is a Lie bracket defined as 𝑋₀𝑋₁−𝑋₁𝑋₀. Here, ∑^∗_Δ𝐽 is the sum for all single and double divided indices of a multi-index 𝐽, and 𝑊^Δ𝐽(𝑡) and 𝑐_Δ𝐽 are multiple Wiener-Stratonovich integrals on time interval [0,𝑡] and the coefficients, respectively, which are determined through a rule for the divided indices of 𝐽. Moreover, |𝐽| denotes the length of 𝐽. We note that (2.2) with (2.4) means that the solution 𝑆(𝑡,𝜔) equals 𝜙(1,𝑠,𝜔) a.s., where 𝜙(𝜏,𝑠,𝜔) is the solution of the ordinary differential equation𝑑𝜙(𝜏)𝑑𝜏=𝑌_𝑡(𝜔)(𝜙(𝜏)),𝜙(0)=𝑠, (2.5) regarding 𝑡 and 𝜔 as parameters.

Now, we proceed to a numerical integration of the stochastic equation (2.1) on the discretized time series by using the representation of solutions to SDEs mentioned above. It adopts a uniform discretization of the time interval [0,𝑇] with step size 𝑇Δ𝑡=𝑁 (2.6) for a fixed natural number 𝑁. Let 𝑡_𝑛=𝑛Δ𝑡(𝑛=0,1,2,…,𝑁) be the 𝑛th step point. Then, for all 𝑛∈{0,…,𝑁}, we abbreviate 𝑆_𝑛=𝑆(𝑡_𝑛) and set 𝑆₀=𝑆(0)=𝑠₀. Moreover, we use Δ𝑊_𝑛 for 𝑛=0,1,…,𝑁 to denote the increments 𝑊(𝑡_𝑛+1)−𝑊(𝑡_𝑛); they are independent Gaussian random variables with mean 0 and variance Δ𝑡, that is, 𝑁(0,Δ𝑡)-distributed random variables.

On account of (2.3), we find a series of numerical solutions 𝑆_𝑛(𝑛=0,1,…,𝑁) to SDE (2.1) by using𝑆_𝑛+1=exp(𝑌_𝑛Δ𝑡)(𝑆_𝑛)(𝑛=0,1,2,…,𝑁−1), (2.7) formally, where 𝑌_𝑛Δ𝑡 is a vector field derived by replacing all the multiple Wiener integrals on time [0,𝑡] by those on time interval Δ𝑡 in (2.4). In the theory of ordinary differential equations, exp(𝑌_𝑛Δ𝑡)(⋅) is often called the time-Δ𝑡 map or exponential map. However, it is usually difficult to find the explicit form of the exponential map, and hence, we need to build an approximation for (2.7). Hence, we formulate a new stochastic numerical scheme as the following two procedures, which are composed of the truncation of the vector field (2.4) and a composition method (or operator-splitting method) applied to the exponential map derived from the truncated vector field.

Procedure. We replace the vector field 𝑌_𝑛Δ𝑡 in (2.7) by a “truncated” vector field 𝑌_𝑛Δt with a truncation order 𝛾 which is defined as 𝑌_𝑛Δ𝑡=_{𝐽;1≤|𝐽|≤𝛾}𝐻^𝐽(Δ𝑡)𝑋^𝐽, (2.8) where 𝐻^𝐽(Δ𝑡) is a polynomial function of multiple Wiener integrals for the vector field 𝑋^𝐽 with a multi-index 𝐽. Then, we define a numerical sequence (𝑆_𝑛)^𝑁_𝑛=0 through 𝑆_𝑛+1=exp(𝑌_𝑛Δ𝑡)(𝑆_𝑛)(𝑛=0,1,…,𝑁−1), (2.9) where 𝑆₀=𝑆(0)=𝑠₀.

Procedure. For 𝑆_𝑛+1=exp(𝑌_𝑛Δ𝑡)(𝑆_𝑛), we apply a “composition method" in a way analogous to that in the theory of ordinary differential equations. Suppose that the vector field 𝑌_𝑛Δ𝑡 is of the form 𝑌_𝑛Δ𝑡=𝐴_𝑛Δ𝑡+𝐵_𝑛Δ𝑡, (2.10) where exp(𝐴_𝑛Δ𝑡) and exp(𝐵_𝑛Δ𝑡) can both be explicitly calculated through (2.15). Then an approximation to the exponential map 𝑌_𝑛Δ𝑡 is given by exp(𝐴_𝑛Δ𝑡)exp(𝐵_𝑛Δ𝑡). Hence, the sequence of (𝑆_𝑛)^𝑁_𝑛=0 in Procedure 1 is approximated by 𝑆_𝑛+1=exp(𝐴_𝑛Δ𝑡)exp(𝐵_𝑛Δ𝑡)(𝑆_𝑛)(𝑛=0,1,…,𝑁−1), (2.11) where 𝑆₀=𝑆(0)=𝑠₀.
We regard (𝑆_𝑛)^𝑁_𝑛=0 as a numerical approximation to the exact discretized solutions (𝑆_𝑛)^𝑁_𝑛=0.

We turn to estimate the approximation error of the numerical scheme described previously. In stochastic numerics, it is usually measured through “global errors in the mean-square sense" for the scheme, and hence we first review the following definition (Gard [11], Kloeden and Platen [12]).

Definition 2.1. Suppose that 𝑆(𝑡) and (𝑆_𝑛)^𝑁_𝑛=0 are an exact solution and the numerical approximation solutions to SDE (2.1), respectively. Moreover, let 𝐸_0,𝜉 be the expectation conditioned on starting at 𝜉 at “initial time” 𝜏. Then the global error order 𝜆 is defined by 𝐸_0,𝑠[|𝑆_𝑁−𝑆_𝑁|²]=𝑂((Δ𝑡)^2𝜆)(Δ𝑡↓0), (2.12) where 𝑆_𝑁=𝑆(𝑁Δ𝑡), 𝑁Δ𝑡=𝑇, and |⋅| denotes the Euclidean norm on the space 𝑅^𝑑. Here, the condition in the expectation (2.12) means 𝑆₀=𝑆₀=𝑠.

It is obvious that the accuracy of the numerical scheme improves with increasing global order.

We remark that in the framework of “global error order" defined by Saito and Mitsui [13], the global order for the numerical solutions 𝑆_𝑛 satisfying (2.12) is given by 2𝜆.

It is evident that such an error estimation order depends on a truncated order of 𝑌_Δ𝑡 and the way of decomposition of the truncated vector field 𝑌_Δ𝑡. In [6], Misawa gives a way to calculate the global error order of stochastic numerical schemes determined through the two procedures mentioned above.

Theorem 2.2 (see [6, Theorem 3.2]). Let 𝛾, 𝛼, and 𝛽 be a truncation order of 𝑌_𝑛Δ𝑡 described in Procedure 1, the least values of “mean square orders" of multiple Wiener integrals which are included in coefficients of the operators 𝐴_𝑛Δ𝑡 and 𝐵_𝑛Δ𝑡, respectively. The mean square order 𝑠 of a multiple Wiener integral 𝐼(Δ𝑡) on time interval Δ𝑡 mentioned above is defined as 𝐸[|𝐼(Δ𝑡)|²]=𝑂((Δ𝑡)^𝑠)(Δ𝑡↓0). Then 𝐸[|𝑆_𝑁−𝑆_𝑁|²∣𝑋₀=𝑋₀=𝑥]=𝑂((Δ𝑡)^𝛿−1)(Δ𝑡↓0) (2.13) holds, where 𝛿=min(𝛾+2,𝛼+𝛽). Therefore, the global error order of stochastic numerical schemes provided by Procedures 1 and 2 is equal to (𝛿−1)/2.

As a preparation for later discussions, we review some examples of our schemes and the global error orders given in [6].

Example 2.3. Suppose that a truncated vector field 𝑌_𝑛Δ𝑡 in Procedure 1 is given by 𝑌_𝑛Δ𝑡=Δ𝑡𝑋₀+Δ𝑊_𝑛𝑋₁. (2.14) We further set 𝐴_𝑛Δ𝑡=Δ𝑡𝑋₀ and 𝐵_𝑛Δ𝑡=Δ𝑊_𝑛𝑋₁ in the decomposition in Procedure 2 and assume that the explicit forms of both exponential maps for them are obtained. Then, the numerical scheme (2.11) can be put into the following form.

Scheme 1. One has𝑆_𝑛+1=exp(Δ𝑡𝑋₀)exp(Δ𝑊_𝑛𝑋₁)(𝑆_𝑛). (2.15) In this case, 𝛾, 𝛼 and 𝛽 are equal to 1, 2 and 1, respectively, and therefore the above theorem provides the global order of Scheme 1 as 1.

Example 2.4. For 𝑌_𝑛(Δ𝑡) in Example 2.3, we set 𝐴_𝑛Δ𝑡 =Δ𝑡𝑋^𝐴₀+Δ𝑊_𝑛𝑋^𝐴₁ and 𝐵_𝑛Δ𝑡=Δ𝑡𝑋^𝐵₀+Δ𝑊_𝑛𝑋^𝐵₁, where 𝑋₀=𝑋^𝐴₀+𝑋^𝐵₀ and 𝑋₁=𝑋^𝐴₁+𝑋^𝐵₁. We assume that [𝑋^𝐴₁,𝑋^𝐵₁]≠0 and that the explicit forms of both exponential maps for them are obtained. Then, the numerical scheme (2.11) can be put into the following form.

Scheme 2. One has 𝑆_𝑛+1=exp(Δ𝑡𝑋^𝐴₀+Δ𝑊_𝑛𝑋^𝐴₁)exp(Δ𝑡𝑋^𝐵₀+Δ𝑊_𝑛𝑋^𝐵₁)(𝑆_𝑛). (2.16) In this case, 𝛼 and 𝛽 are both equal to 1; hence the global order of the scheme (2.10) for the above 𝐴_𝑛Δ𝑡 and 𝐵_𝑛Δ𝑡 becomes 0.5. Note that this order is equal to that of the Euler-Maruyama scheme for SDEs which is the most famous scheme in stochastic numerics.

Remark 2.5. By further manipulating the BCH formula to eliminate higher-order terms, we can obtain various schemes which give higher-order approximations to the exponential map. For example, the scheme corresponding to “leapfrog", which is well known in deterministic numerical analysis, is given by exp((Δ𝑡)(𝑋+𝑌))=exp(Δ𝑡𝑌2)exp(Δ𝑡𝑋)exp(Δ𝑡𝑌2)+𝑂((Δ𝑡)³). (2.17) In a way analogous to that in (2.11), we define a stochastic leapfrog scheme as follows: 𝑆_𝑛+1𝐵=exp(^𝑛Δ𝑡2)exp(𝐴_𝑛Δ𝑡𝐵)exp(^𝑛Δ𝑡2)(𝑆_𝑛)(𝑛=0,1,…,𝑁−1). (2.18) Then, using the BCH formula, we can produce numerical schemes having a better error order than that of scheme (2.11). In detail, see Remark 3.3 and Example 3.2 in [6].

3. Symplectic Integrators of Stochastic Hamiltonian DynamicalSystems

Now, we proceed to derive symplectic integrators for stochastic Hamiltonian dynamical systems in terms of composition methods mentioned in the previous section. We start with the definition of stochastic Hamiltonian dynamical systems defined in [3]. Consider the following 2ℓ-dimensional stochastic dynamical system with one Wiener process:𝑑𝑥^𝑖𝑥(𝑡)_ℓ+𝑖=𝜕(𝑡)^ℓ+𝑖𝐻₀(𝑥(𝑡))−𝜕_𝑖𝐻₀𝜕(𝑥(𝑡))𝑑𝑡+^ℓ+𝑖𝐻₁(𝑥(𝑡))−𝜕_𝑖𝐻₁(𝑥(𝑡))∘𝑑𝑊(𝑡)(𝑖=1,…,ℓ), (3.1) where 𝑥=(𝑥^𝑘)^2ℓ_𝑘=1 and 𝜕_𝑗=𝜕/𝜕𝑥^𝑗(𝑗=1,2,…,2ℓ), respectively. In (3.1), 𝐻_𝛼(𝑥)(𝛼=0,1) are smooth scalar functions on 𝑅^2ℓ. Formally, one may regard this as a Hamiltonian dynamical system 𝑑𝑥𝑑𝑡𝑥^𝑖ℓ+𝑖=𝜕_ℓ+𝑖𝐻(𝑥)−𝜕_𝑖𝐻(𝑥)(𝑖=1,…,ℓ) (3.2) with a “randomized" Hamiltonian 𝐻 given by 𝐻=𝐻₀+𝐻₁𝛾_𝑡, (3.3) where 𝛾_𝑡 is a one-dimensional Gaussian white noise. Hence, we call (3.1) and 𝐻_𝛼(𝛼=0,1) an (ℓ-dimensional) stochastic Hamiltonian dynamical system and the Hamiltonian, respectively. In the following, for simplicity, we set ℓ=1 and denote 𝑥¹(𝑡) and 𝑥²(𝑡) by 𝑞(𝑡) and 𝑝(𝑡), respectively, that is, we treat the following stochastic Hamiltonian dynamical system:𝑑=𝜕𝑞(𝑡)𝑝(𝑡)𝐻𝜕𝑝⁰−𝜕(𝑝(𝑡),𝑞(𝑡))𝐻𝜕𝑞⁰𝜕(𝑝(𝑡),𝑞(𝑡))𝑑𝑡+𝐻𝜕𝑝¹−𝜕(𝑝(𝑡),𝑞(𝑡))𝐻𝜕𝑞¹(𝑝(𝑡),𝑞(𝑡))∘𝑑𝑊(𝑡). (3.4)

Let 𝑡_𝑛=𝑛Δ𝑡(𝑛=0,1,2,…,𝑁) be the 𝑛th step-point in a uniform discretization of the time interval [0,𝑇] with step size Δ𝑡=𝑇/𝑁 for fixed natural number 𝑁. Suppose that {(𝑝_𝑛,𝑞_𝑛)}^𝑁_𝑛=0 is a series of the numerical solutions for the above Hamiltonian system derived by some numerical scheme on the above time discretization. In the same manner as that in the case of the numerical scheme on deterministic Hamiltonian systems, we call the numerical scheme for stochastic Hamiltonian dynamical systems “symplectic", which produces numerical solutions satisfying det(𝜕(𝑝_𝑛+1,𝑞_𝑛+1)𝜕(𝑝_𝑛,𝑞_𝑛))=1 (3.5) on any step point 𝑛. In consideration of the advantages of using symplectic schemes for usual deterministic systems, we also may expect to obtain more stable numerical orbits by using these schemes than by using ordinary numerical schemes, for example, the Euler scheme.

Now, we produce the symplectic schemes by using the composition methods mentioned in Section 2. First, as an example of symplectic scheme on the basis of Scheme 1, we can provide the following scheme.

Scheme 3. One has 𝑞_𝑛+1𝑝_𝑛+1=exp(𝐴_𝑛Δ𝑡)exp(𝐵_𝑛Δ𝑡)𝑞_𝑛𝑝_𝑛, (3.6) where a truncated vector field 𝑌_𝑛Δ𝑡 shown by (2.14) and the splitting vecter fields 𝐴_𝑛Δ𝑡 and 𝐵_𝑛Δ𝑡 are given by 𝑌_𝑛Δ𝑡=Δ𝑡𝑋₀+Δ𝑊_𝑛𝑋₁=(𝜕_𝑝𝐻₀𝜕_𝑞−𝜕_𝑞𝐻₀𝜕_𝑝)Δ𝑡+(𝜕_𝑝𝐻₁𝜕_𝑞−𝜕_𝑞𝐻₁𝜕_𝑝)Δ𝑊_𝑛,𝐴_𝑛Δ𝑡=Δ𝑡𝑋₀=(𝜕_𝑝𝐻₀𝜕_𝑞−𝜕_𝑞𝐻₀𝜕_𝑝𝐵)Δ𝑡,_𝑛Δ𝑡=Δ𝑊_𝑛𝑋₁=(𝜕_𝑝𝐻₁𝜕_𝑞−𝜕_𝑞𝐻₁𝜕_𝑝)Δ𝑊_𝑛, (3.7)(3.8)(3.9) respectively.

It is easy to see that the above composition scheme is symplectic, since for each fixed stochastic parameter 𝜔, 𝐴_𝑛Δ𝑡, and 𝐵_𝑛Δ𝑡 are Hamiltonian vector fields with the Hamiltonians 𝐻₀Δ𝑡 and 𝐻₁Δ𝑊_𝑛 on time interval [𝑡_𝑛,𝑡_𝑛+1], respectively, and hence the composition of their exponential maps is symplectic (McLachlan and Quispel [5]). Moreover, as mentioned in Section 2, the global error order of this scheme is 1.

Next, we give an example of the symplectic scheme on the basis of Scheme 2.

Scheme 4. One has 𝑞_𝑛+1𝑝_𝑛+1=exp(𝐴_𝑛Δ𝑡)exp(𝐵_𝑛Δ𝑡)𝑞_𝑛𝑝_𝑛, (3.10) where for the same truncated vector field 𝑌_𝑛Δ𝑡 as Scheme 3, we choose 𝐴_𝑛Δ𝑡 and 𝐵_𝑛Δ𝑡 as follows: 𝐴_𝑛Δ𝑡=Δ𝑡𝑋^𝐴₀+Δ𝑊_𝑛𝑋^𝐴₁=Δ𝑡(𝜕_𝑝𝐻^𝐴₀𝜕_𝑞−𝜕_𝑞𝐻^𝐴₀𝜕_𝑝)+Δ𝑊_𝑛(𝜕_𝑝𝐻^𝐴₁𝜕_𝑞−𝜕_𝑞𝐻^𝑎₁𝜕_𝑝𝐵),_𝑛Δ𝑡=Δ𝑡𝑋^𝐵₀+Δ𝑊_𝑛𝑋^𝐵₁=Δ𝑡(𝜕_𝑝𝐻^𝐵₀𝜕_𝑞−𝜕_𝑞𝐻^𝐵₀𝜕_𝑝)+Δ𝑊_𝑛(𝜕_𝑝𝐻^𝐵₁𝜕_𝑞−𝜕_𝑞𝐻^𝐵₁𝜕_𝑝), (3.11)(3.12) under the conditions 𝐻^𝐴₀+𝐻^𝐵₀=𝐻₀ and 𝐻^𝐴₁+𝐻^𝐵₁=𝐻₁.
We can also prove that this scheme is symplectic, since for each fixed stochastic parameter 𝜔, 𝐴_𝑛Δ𝑡, and 𝐵_𝑛Δ𝑡 are also Hamiltonian vector fields with the Hamiltonians Δ𝑡𝐻^𝐴₀+Δ𝑊_𝑛𝐻^𝐴₁ and Δ𝑡𝐻^𝐵₀+Δ𝑊_𝑛𝐻^𝐵₁ on time interval [𝑡_𝑛,𝑡_𝑛+1]. We note that the global error order of this scheme is 0.5 as mentioned in Section 2.

Remark 3.1. If concrete Hamiltonians 𝐻₀ and 𝐻₁ are given, one can also check the condition (3.5) for each scheme mentioned above by direct calculation.

4. Illustrative Examples

In this section, as an illustrative example of Schemes 3 and 4, we will show the numerical solutions to a stochastic harmonic oscillator system by the symplectic schemes and compare them with the numerical results by the Euler scheme.

Let us consider the following stochastic Hamiltonians dynamical system (3.4) with the Hamiltonian 𝐻₀=(1/2)𝑝²+(1/2)𝑞² and 𝐻₁=𝜎𝐻₀:𝑑=𝑞(𝑡)𝑝(𝑡)𝑝(𝑡)−𝑞(𝑡)𝑑𝑡+𝜎𝑝(𝑡)−𝑞(𝑡)∘𝑑𝑊(𝑡), (4.1) where 𝜎 is a constant. It is evident that this is a stochastic system with the “conserved quantity" 𝐻₀ [3, 6]; that is, the solution (𝑞(𝑡),𝑝(𝑡)) randomly moves on a circle determined by 𝐻₀(𝑞(𝑡),𝑝(𝑡))=𝐻₀(𝑞(0),𝑝(0))=(1/2)(𝑞(0)²+𝑝(0)²) on 𝑝-𝑞 plane. We call such a system a stochastic harmonic oscillator system.

Now, in order to investigate the advantage of our symplectic schemes, we will produce numerical solutions to this system by the new schemes and the “Euler scheme" [11–13]. Then, by straightforward calculation of the definitions of the schemes together with the exponential maps (2.5), the Euler scheme, Schemes 3 and 4 under setting 𝐻^𝐴₀, 𝐻^𝐴₁, 𝐻^𝐵₀, and 𝐻^𝐵₁ as 𝐻^𝐴₀=(1/𝜎)𝐻^𝐴₁=(1/2)𝑝² and 𝐻^𝐵₀=(1/𝜎)𝐻^𝐵₁=(1/2)𝑞² for this system are given as follows.

Euler Scheme:
𝑞_𝑛+1𝑝_𝑛+1=𝑞_𝑛𝑝_𝑛+𝑝_𝑛−12𝑞_𝑛−𝑞_𝑛−12𝑝_𝑛𝑝Δ𝑡+^𝑛−𝑞_𝑛𝜎Δ𝑊_𝑛. (4.2)Scheme 5 (Scheme 3 for (4.1)). One has 𝑞_𝑛+1𝑝_𝑛+1=cos(Δ𝑡)−sin(Δ𝑡)sin(Δ𝑡)cos(Δ𝑡)cos(𝜎Δ𝑊_𝑛)−sin(𝜎Δ𝑊_𝑛)sin(𝜎Δ𝑊_𝑛)cos(𝜎Δ𝑊_𝑛)𝑞_𝑛𝑝_𝑛. (4.3)Scheme 6 (Scheme 4 for (4.1)). One has 𝑞_𝑛+1𝑝_𝑛+1=𝑝_𝑛(Δ𝑡+𝜎Δ𝑊_𝑛)+𝑞_𝑛𝑝_𝑛, (4.4) where 𝑞_𝑛𝑝_𝑛=𝑞_𝑛−𝑞_𝑛(Δ𝑡+𝜎Δ𝑊_𝑛)+𝑝_𝑛. (4.5) Here we note that we can directly examine the symplectic condition (3.5) which holds for Schemes 5 and 6mentioned above.

In these schemes, Δ𝑊_𝑛 are numerically realized in terms of the independent 𝑁(0,1) random numbers 𝛾_𝑛 as [12]Δ𝑊_𝑛=𝛾_𝑛√Δ𝑡. (4.6)Moreover, we set 𝜎, Δ𝑡, and 𝑁 as 𝜎=0.2, Δ𝑡=0.05, and 𝑁=1000, respectively, and we choose the initial values of 𝑝₀ and 𝑞₀ as 1 and 0, respectively.

Figure 1 denotes the numerical time series of “𝑝" by these schemes, and Figures 2(a)–2(c) display sample paths of numerical solutions from Euler scheme andSchemes 6 and 5, respectively. Here, as mentioned, it should be remembered that the exact solution to the original stochastic system should run on the circle 𝐻₀(𝑞(𝑡),𝑝(𝑡))=(1/2)(𝑞(0)²+𝑝(0)²)=1, that is, a unit circle on the 𝑝-𝑞 plane. Figures 1 and 2(a), however, indicate that a series of the numerical solutions by the Euler scheme gradually goes far from the circle. In contrast, Figures 1, 2(b), and 2(c) show that the numerical solutions by our symplecticSchemes 5 and 6move around the circle stably, and especially, a series of the solution byScheme 5runs completely on the circle. Thus, as in the deterministic case, the new stochastic symplectic schemes also numerically preserve the Hamiltonian dynamical structure; this is the superiority of the new stochastic schemes over the ordinary stochastic schemes.

Figure 1 

Numerical sample paths of “$p$” of a stochastic harmonic oscillator system (SHOS) via Euler scheme, Composition 1 (Scheme 6) and Composition 2 (Scheme 5) ($dt$ = 0.05, sigma = 0.2).

(a)

(b)

(c)

Figure 2 

(a) A numerical sample path of SHOS in phase space via Euler scheme. (b) A numerical sample path of SHOS in phase space via Composition 1 (Scheme 6). (c) A numerical sample path of SHOS in phase space via Composition 2 (Scheme 5).

Finally, we give some concluding remarks on our method.(1)In this paper, we have treated the stochastic Hamiltonian dynamical systems with a one-dimensional Wiener process. However, it often happens that the error order of a numerical method collapses, if there is more than one Wiener process. Hence, it would be important to investigate the error orders of our schemes in the case of stochastic Hamiltonian dynamical systems with a “multidimensional Wiener process". Therefore, we should be able to improve our composition methods and offer new schemes with high order for SDEs with a multi-dimensional Wiener process.(2)Moreover, we have referred only to a single composition method to SDEs so far. It seems that one would readily be able to extend the method to a multi-composition one, since such compositions are obtained by iterations of single compositions. Using the multicompositions, we may easily produce a symplectic scheme for the more complicated stochastic Hamiltonian systems.

In future papers, the topics mentioned above will be reported.

Acknowledgment

This work was supported by Grants-in-Aid for Scientific Research nos. 18540134 and 21540140 from the Japan Society for the Promotion of Science.