Mathematical Problems in Engineering

Volume 2010, Article ID 384937, 12 pages

http://dx.doi.org/10.1155/2010/384937

## Symplectic Integrators to Stochastic Hamiltonian Dynamical Systems Derived from Composition Methods

Faculty of Economics, Nagoya City University, Nagoya 467-8501, Japan

Received 25 December 2009; Accepted 20 June 2010

Academic Editor: Katica R. (Stevanovic) Hedrih

Copyright © 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.

#### 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 , that is, the solutions of differential equations of the form are given in the form . Assume that one can write in such a way that and 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
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 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:
where and are -dimensional vector functions, respectively. By using these functions, we define the differential operators and as
where . According to Kunita [10], a solution of the SDE under an initial value can be formally represented in terms of an exponential map derived from a “stochastic" vector field as follows:
where is the vector field for each and almost surely (a.s.) given by
and , 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 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 a.s., where is the solution of the * ordinary* differential equation
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 with step size for a fixed natural number . Let be the th step point. Then, for all , we abbreviate and set . Moreover, we use for to denote the increments ; they are independent Gaussian random variables with mean 0 and variance , that is, -distributed random variables.

On account of (2.3), we find a series of numerical solutions to SDE (2.1) by using formally, where is a vector field derived by replacing all the multiple Wiener integrals on time by those on time interval in (2.4). In the theory of ordinary differential equations, 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 with a truncation order which is defined as
where is a polynomial function of multiple Wiener integrals for the vector field with a multi-index . Then, we define a numerical sequence through
where .

*Procedure. *For , 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
where and can both be explicitly calculated through (2.15). Then an approximation to the exponential map is given by . Hence, the sequence of in Procedure 1 is approximated by
where .

We regard as a numerical approximation to the exact discretized solutions .

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 are an exact solution and the numerical approximation solutions to SDE (2.1), respectively. Moreover, let be the expectation conditioned on starting at at * “initial time” *. Then * the global error order * is defined by
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 .

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 ]). *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 . Then
**
holds, where . Therefore, the global error order of stochastic numerical schemes provided by Procedures 1 and 2 is equal to .*

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
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
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 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
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
In a way analogous to that in (2.11), we define a stochastic leapfrog scheme as follows:
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 and Example 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:
where and , respectively. In (3.1), are smooth scalar functions on . Formally, one may regard this as a Hamiltonian dynamical system
with a “randomized" Hamiltonian given by
where is a one-dimensional Gaussian white noise. Hence, we call (3.1) and an (-dimensional) * stochastic Hamiltonian dynamical system* and the * Hamiltonian*, respectively. In the following, for simplicity, we set and denote and by and , respectively, that is, we treat the following stochastic Hamiltonian dynamical system:

Let be the th step-point in a uniform discretization of the time interval with step size for fixed natural number . Suppose that 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 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
where a truncated vector field shown by (2.14) and the splitting vecter fields and are given by
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 , 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
where for the same truncated vector field as Scheme 3, we choose and as follows:
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 . 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 and : 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 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 and for this system are given as follows.

*Euler Scheme:**Scheme 5 (Scheme 3 for (4.1)). *One has
*Scheme 6 (Scheme 4 for (4.1)). *One has
where
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 random numbers as [12]Moreover, we set , , and as , , and , respectively, and we choose the initial values of and as and , 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 , 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.

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.

#### References

- E. Hairer, C. Lubich, and G. Wanner,
*Geometric Numerical Integration*, vol. 31 of*Springer Series in Computational Mathematics*, Springer, Berlin, Germany, 2nd edition, 2006, Structure-Preserving Algorithms for Ordinary Differential Equation. View at MathSciNet - G. N. Milstein, Yu. M. Repin, and M. V. Tretyakov, “Numerical methods for stochastic systems preserving symplectic structure,”
*SIAM Journal on Numerical Analysis*, vol. 40, no. 4, pp. 1583–1604, 2002. View at Publisher · View at Google Scholar · View at MathSciNet - T. Misawa, “Conserved quantities and symmetries related to stochastic dynamical systems,”
*Annals of the Institute of Statistical Mathematics*, vol. 51, no. 4, pp. 779–802, 1999. View at Publisher · View at Google Scholar · View at MathSciNet - N. Ikeda and S. Watanabe,
*Stochastic Differential Equations and Diffusion Processes*, vol. 24 of*North-Holland Mathematical Library*, North-Holland, Amsterdam, The Netherlands, 2nd edition, 1989. View at MathSciNet - R. I. McLachlan and G. R. W. Quispel, “Splitting methods,”
*Acta Numerica*, vol. 11, pp. 341–434, 2002. View at Publisher · View at Google Scholar · View at MathSciNet - T. Misawa, “A Lie algebraic approach to numerical integration of stochastic differential equations,”
*SIAM Journal on Scientific Computing*, vol. 23, no. 3, pp. 866–890, 2001. View at Publisher · View at Google Scholar · View at MathSciNet - K. Burrage and P. M. Burrage, “High strong order methods for non-commutative stochastic ordinary differential equation systems and the Magnus formula,”
*Physica D*, vol. 133, no. 1–4, pp. 34–48, 1999. View at Publisher · View at Google Scholar · View at MathSciNet - S. J. A. Malham and A. Wiese, “Stochastic Lie group integrators,”
*SIAM Journal on Scientific Computing*, vol. 30, no. 2, pp. 597–617, 2008. View at Publisher · View at Google Scholar · View at MathSciNet - N. Bou-Rabee and H. Owhadi, “Stochastic variational integrators,”
*IMA Journal of Numerical Analysis*, vol. 29, no. 2, pp. 421–443, 2009. View at Publisher · View at Google Scholar · View at MathSciNet - H. Kunita, “On the representation of solutions of stochastic differential equations,” in
*Séminaire de Probabilités XIV 1978/79*, vol. 784 of*Lecture Notes in Mathematics*, pp. 282–304, Springer, Berlin, Germany, 1980. View at Publisher · View at Google Scholar · View at MathSciNet - T. C. Gard,
*Introduction to Stochastic Differential Equations*, vol. 114 of*Monographs and Textbooks in Pure and Applied Mathematics*, Marcel Dekker, New York, NY, USA, 1988. View at MathSciNet - P. E. Kloeden and E. Platen,
*Numerical Solution of Stochastic Differential Equations*, vol. 23 of*Applications of Mathematics (New York)*, Springer, Berlin, Germany, 1992. View at MathSciNet - Y. Saito and T. Mitsui, “Simulation of stochastic differential equations,”
*Annals of the Institute of Statistical Mathematics*, vol. 45, no. 3, pp. 419–432, 1993. View at Publisher · View at Google Scholar · View at MathSciNet