Symplectic Integrators to Stochastic Hamiltonian Dynamical Systems Derived from Composition Methods

Copyright q 2010 Tetsuya Misawa. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. " 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.


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

Mathematical Problems in Engineering
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 X denote vector fields on some space with coordinates x, with flows exp tX , that is, the solutions of differential equations of the form ẋ t X x are given in the form x t exp tX x 0 .Assume that one can write X A B in such a way that exp tA and exp tB 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 x t through x t exp tA exp tB x x t O t 2 ; 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 tX are often preserved.If X, A, and B 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 .

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 Ω, F, P be a probability space.We consider an n-dimensional stochastic differential equation of the Stratonovich type associated with a one-dimensional Wiener process W as follows: where b and g are n-dimensional C ∞ vector functions, respectively.By using these functions, we define the differential operators X 0 and X 1 as where ∂ i ∂/∂S i .According to Kunita 10 , a solution of the SDE under an initial value S 0 s can be formally represented in terms of an exponential map derived from a "stochastic" vector field Y t as follows: where Y t ω is the vector field for each t and almost surely a.s.ω ∈ Ω given by Here, * ΔJ is the sum for all single and double divided indices of a multi-index J, and W ΔJ t and c ΔJ are multiple Wiener-Stratonovich integrals on time interval 0, t and the coefficients, respectively, which are determined through a rule for the divided indices of J.Moreover, |J| denotes the length of J.We note that 2.2 with 2.4 means that the solution S t, ω equals φ 1, s, ω a.s., where φ τ, s, ω is the solution of the ordinary differential equation regarding t 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, T with step size for a fixed natural number N. Let t n nΔt n 0, 1, 2, . . ., N be the nth step point.Then, for all n ∈ {0, . . ., N}, we abbreviate S n S t n and set S 0 S 0 s 0 .Moreover, we use ΔW n for n 0, 1, . . ., N to denote the increments W t n 1 − W t n ; they are independent Gaussian random variables with mean 0 and variance Δt, that is, N 0, Δt -distributed random variables.On account of 2.3 , we find a series of numerical solutions S n n 0, 1, . . ., N to SDE 2.1 by using formally, where Y nΔt is a vector field derived by replacing all the multiple Wiener integrals on time 0, t by those on time interval Δt in 2.4 .In the theory of ordinary differential equations, exp Y nΔt • is often called the time-Δt 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 1.We replace the vector field Y nΔt in 2.7 by a "truncated" vector field Y nΔt with a truncation order γ which is defined as where H J Δt is a polynomial function of multiple Wiener integrals for the vector field X J with a multi-index J.Then, we define a numerical sequence S n N n 0 through where S 0 S 0 s 0 .
Procedure 2. For S n 1 exp Y nΔt S n , we apply a "composition method" in a way analogous to that in the theory of ordinary differential equations.Suppose that the vector field Y nΔt is of the form where exp A nΔt and exp B nΔt can both be explicitly calculated through 2.15 .Then an approximation to the exponential map Y nΔt is given by exp A nΔt exp B nΔt .Hence, the sequence of S n N n 0 in Procedure 1 is approximated by where S 0 S 0 s 0 .
We regard S n N n 0 as a numerical approximation to the exact discretized solutions S n N n 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 meansquare sense" for the scheme, and hence we first review the following definition Gard 11 , Kloeden and Platen 12 .
Definition 2.1.Suppose that S t and S n N n 0 are an exact solution and the numerical approximation solutions to SDE 2.1 , respectively.Moreover, let E 0,ξ be the expectation conditioned on starting at ξ at "initial time" τ.Then the global error order λ is defined by where S N S NΔt , NΔt T , and | • | denotes the Euclidean norm on the space R d .Here, the condition in the expectation 2.12 means S 0 S 0 s.
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 S n satisfying 2.12 is given by 2λ.
It is evident that such an error estimation order depends on a truncated order of Y Δt and the way of decomposition of the truncated vector field Y Δt .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 Y nΔt described in Procedure 1, the least values of "mean square orders" of multiple Wiener integrals which are included in coefficients of the operators A nΔt and B nΔt , respectively.The mean square order s of a multiple Wiener integral I Δt on time interval Δt mentioned above is defined as 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 Y nΔt in Procedure 1 is given by

2.14
We further set A nΔt ΔtX 0 and B nΔt ΔW n X 1 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
S n 1 exp ΔtX 0 exp ΔW n X 1 S n .

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.
We assume that X A 1 , X B 1 / 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.

2.16
In this case, α and β are both equal to 1; hence the global order of the scheme 2.10 for the above A nΔt and B nΔt 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

2.17
In a way analogous to that in 2.11 , we define a stochastic leapfrog scheme as follows:

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 .

Symplectic Integrators of Stochastic Hamiltonian Dynamical Systems
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: where x x k 2 k 1 and ∂ j ∂/∂x j j 1, 2, . . ., 2 , respectively.In 3.1 , H α x α 0, 1 are smooth scalar functions on R 2 .Formally, one may regard this as a Hamiltonian dynamical system with a "randomized" Hamiltonian H given by where γ t is a one-dimensional Gaussian white noise.Hence, we call 3.1 and H α α 0, 1 an -dimensional stochastic Hamiltonian dynamical system and the Hamiltonian, respectively.In the following, for simplicity, we set 1 and denote x 1 t and x 2 t by q t and p t , respectively, that is, we treat the following stochastic Hamiltonian dynamical system: Let t n nΔt n 0, 1, 2, . . ., N be the nth step-point in a uniform discretization of the time interval 0, T with step size Δt T/N for fixed natural number N. Suppose that { p n , q n } N n 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 ∂ p n 1 , q n 1 ∂ p n , q n 1 3.5 on any step point n.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.
It is easy to see that the above composition scheme is symplectic, since for each fixed stochastic parameter ω, A nΔt , and B nΔt are Hamiltonian vector fields with the Hamiltonians H 0 Δt and H 1 ΔW n on time interval t n , t n 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
q n 1 p n 1 exp A nΔt exp B nΔt q n p n , 3.10 where for the same truncated vector field Y nΔt as Scheme 3, we choose A nΔt and B nΔt as follows: ΔW n H B 1 on time interval t n , t n 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 H 0 and H 1 are given, one can also check the condition 3.5 for each scheme mentioned above by direct calculation.Here we note that we can directly examine the symplectic condition 3.5 which holds for Schemes 5 and 6 mentioned above.
In these schemes, ΔW n are numerically realized in terms of the independent N 0, 1 random numbers γ n as 12 Moreover, we set σ, Δt, and N as σ 0.2, Δt 0.05, and N 1000, respectively, and we choose the initial values of p 0 and q 0 as 1 and 0, respectively.Figure 1 denotes the numerical time series of "p" 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 H 0 q t , p t 1/2 q 0 2 p 0 2 1, that is, a unit circle on the p-q 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 6 move around the circle stably, and especially, a series of the solution byScheme 5 runs 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.
Finally, we give some concluding remarks on our method.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.
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.

Figure 1 :
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 .

Figure 2 :
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 .