Advanced
Stability of Explicit Symplectic Partitioned Runge-Kutta Methods
Stability of Explicit Symplectic Partitioned Runge-Kutta Methods
Journal of Information and Communication Convergence Engineering. 2014. Mar, 12(1): 39-45
Copyright © 2014, The Korea Institute of Information and Commucation Engineering
This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.
  • Received : October 09, 2013
  • Accepted : December 02, 2013
  • Published : March 31, 2014
Download
PDF
e-PUB
PubReader
PPT
Export by style
Share
Article
Author
Metrics
Cited by
TagCloud
About the Authors
Toshiyuki Koto
Department of Information Systems and Mathematical Sciences, Nanzan University, Seto 489-0863, Japan
koto@ms.nanzan-u.ac.jp
Eunjee Song
Department of Computer Science, Namseoul University, Cheonan 331-707, Korea

Abstract
A numerical method for solving Hamiltonian equations is said to be symplectic if it preserves the symplectic structure associated with the equations. Various symplectic methods are widely used in many fields of science and technology. A symplectic method preserves an approximate Hamiltonian perturbed from the original Hamiltonian. It theoretically supports the effectiveness of symplectic methods for long-term integration. Although it is also related to long-term integration, numerical stability of symplectic methods have received little attention. In this paper, we consider explicit symplectic methods defined for Hamiltonian equations with Hamiltonians of the special form, and study their numerical stability using the harmonic oscillator as a test equation. We propose a new stability criterion and clarify the stability of some existing methods that are visually based on the criterion. We also derive a new method that is better than the existing methods with respect to a Courant-Friedrichs-Lewy condition for hyperbolic equations; this new method is tested through a numerical experiment with a nonlinear wave equation.
Keywords
I. INTRODUCTION
A symplectic integration method is a numerical method for solving Hamiltonian equations,
PPT Slide
Lager Image
a special class of differential equations related to classical mechanics and symplectic geometry. Various symplectic methods are designed and widely used in celestial mechanics, molecular dynamics, electromagnetic field analysis, etc., particularly for the longterm integration of Hamiltonian equations.
The time evolution of Hamiltonian equations preserves a special differential 2-form dp dq called the symplectic form. A numerical method is said to be symplectic if it also preserves the symplectic form. Since the concept of symplectic integration methods was proposed in the mid- 1980s [1] , many mathematical researches have been carried out [2 - 4] . In particular, it has been revealed that a symplectic method preserves an approximate Hamiltonian perturbed from the original Hamiltonian [5 , 6] . It theoretically supports the effectiveness of symplectic methods for long-term integration.
On the other hand, the numerical stability of symplectic methods has received little attention, although it is also related to long-term integration; only a few papers [7 , 8] deal with this subject. It is certain that many outstanding symplectic methods are implicit and possess originally superior stability. However, in a large-scale computation, e.g., in the solution of partial differential equations, explicit methods are still effective tools. A study of their stability has significance for practical computation because the stability of numerical methods is closely related to step size restrictions, such as a Courant-Friedrichs-Lewy (CFL) condition for hyperbolic equations.
In this paper, we study the stability of an explicit symplectic method by using the harmonic oscillator as a test equation, following [8] . An outline of this paper is as follows: In Section II, we describe the fundamental concept and notation concerning explicit symplectic methods and their numerical stability. In Section III, we propose a new stability criterion for the symplectic methods and discuss the stability of the basic methods on the basis of this criterion. In Section IV, we continue to analyze more advanced methods and derive a new method, which is tested through a numerical experiment with the sine-Gordon equation, a nonlinear wave equa-tion in Section V.
II. PRELIMINARIES
- A. Explicit Symplectic Methods
We consider a Hamiltonian of the special form
PPT Slide
Lager Image
and the initial value problem
PPT Slide
Lager Image
for the corresponding Hamiltonian equation, where
PPT Slide
Lager Image
In mechanics, T and U represent kinetic energy and potential energy, respectively.
In general, symplectic methods are implicit; i.e., it is necessary to solve nonlinear equations for the implementation of these methods. For problem (3), there are explicit symplectic methods by virtue of the special form (2). A well-known instance is a symplectic partitioned Runge-Kutta method, whose general form is as follows (see, e.g., [2 , 4] ):
PPT Slide
Lager Image
Here, Δ t > 0 is the time step size, tn = nΔt ( n = 0,1,… ), and qn pn are approximate values for q(tn) and p(tn) , respectively. Further, b 1 , b 2 , … b s ,
PPT Slide
Lager Image
and
PPT Slide
Lager Image
are parameters of the method, and Qi and Pi are intermediate variables for computation. The parameters of the method, determined from order conditions [2 , 4] , are often written as
PPT Slide
Lager Image
- B. Test Equation for Stability Analysis
To study the stability of the symplectic method (5), we adopt the harmonic oscillator
PPT Slide
Lager Image
as a test equation ( [8] ; see also [7] for another test equation). This is a Hamiltonian equation with the Hamiltonian H(p,q) = (ω /2 )( p 2 + q 2 ), ω ≥ 0. We also adopt the scaled step size
PPT Slide
Lager Image
as a basic parameter for the stability analysis. Upon the restriction of the frequency ω ≥ 0, the range of the parameter is θ ≥ 0.
It should be noted that exact solutions to (6) satisfy
PPT Slide
Lager Image
The matrix M ( ) is an orthogonal matrix, and its eigenvalues are
PPT Slide
Lager Image
, both of which have unit modulus.
In the case f ( p ) = ωp and g ( t,q ) = −ω q , the equations for the intermediate variables in (5) becomes
PPT Slide
Lager Image
The substitution of the first equation into the second equation gives
PPT Slide
Lager Image
Hence, (9) is rewritten as
PPT Slide
Lager Image
and application of method (5) to test equation (6) yields an analogue to (8),
PPT Slide
Lager Image
It is clear that det Mi ( θ ) = 1. Hence M * ( θ ) = 1 holds for any method of the form(5). The Characteristic equation M * ( θ ) is written as
PPT Slide
Lager Image
and the eigen values are given by
PPT Slide
Lager Image
where tr M * ( θ ) denotes the trace of the matrix M * ( θ ). If │ tr M * ( θ )│ < 2, the eigenvalues are complex numbers with |λ|=1. If tr M * ( θ ) = 2, then λ = 1, and if tr M * ( θ ) = − 2, then λ = − 1. If │tr M * ( θ )│ > 2, the eigenvalues are real, and one of them satisfies |λ| > 1. The set { θ ≥ 0 : │ tr M * ( θ )│ ≤ 2} is a union of closed intervals. The connected component of the set that contains the origin is called the stability interval of method (5), which has been used for comparing the stability of numerical methods [8] .
III. STABILITY CRITERION
If │tr M * ( θ )│ < 2, M * ( θ ) has complex conjugate eigenvalues λ,
PPT Slide
Lager Image
which satisfy │λ│=│
PPT Slide
Lager Image
│= 1 and λ ≠
PPT Slide
Lager Image
Hence, M * ( θ ) is represented in form
PPT Slide
Lager Image
with some nonsingular matrix T. Since
PPT Slide
Lager Image
and │λ│=│
PPT Slide
Lager Image
│ = 1, we have ║ M * ( θ ) n ║ ≤║ T ║ ║ T -1 ║ for any integer n ≥ 0, where ║•║ denotes the matrix norm induced from the Euclidean norm. The upper bound ║ T ║║ T -1 ║ is represented as follows.
Theorem 1. Let a,b,c,d be real numbers, Assume that
PPT Slide
Lager Image
satisfies det M = 1 and │ tr M │ < 2. Then, we have
PPT Slide
Lager Image
for any integer n ≥ 0, where
PPT Slide
Lager Image
The proof of the theorem is obtained by a simple but tiresome computation. We omit the proof (cf. the proof of Theorem 3.1 in [9] ). As shown below, ϕ in Theorem 1 is used as a criterion for the stability of the numerical methods.
In the case s = 1 and
PPT Slide
Lager Image
(5) is reducde to
PPT Slide
Lager Image
This is called the symplectic Euler method and is of the order 1 in accuracy. In the case of the symplectic Euler method, we have
PPT Slide
Lager Image
Since tr M ( θ ) = 2 - θ 2 , the stability interval of the method is [0, 2]. For 0 < θ < 2, ϕ in Theorem 1 is computed as
PPT Slide
Lager Image
In the case s = 2, method (5) is rewritten as
PPT Slide
Lager Image
which is of the order 2 if the parameters satisfy
PPT Slide
Lager Image
In particular, the parameter values
PPT Slide
Lager Image
satisfy the condition, and the corresponding method is known as the Störmer - Verlet method [4 , 8] .
For this method, we have
PPT Slide
Lager Image
Since tr M ( θ ) = 2 − θ 2 , the stability interval of the Stormer - Verlet method is [0, 2], which is the same as that of the symplectic Euler method. However, since (2 θ - θ 3 /4) 2 - {4 - (2 - θ 2 ) 2 } = θ 6 /16 , we have, for 0 < θ < 2,
PPT Slide
Lager Image
Fig. 1 shows the functions ϕ for the two methods. Function (25) for the Störmer-Verlet method is closer to the line ϕ = 1 than (20) for the symplectic Euler method. The matrix M ( θ ) in (8) is an orthogonal matrix and satisfies ║ M ( θ ) n ║ = 1 for any θ ≥ 0 and any integer n ≥ 0. Since (25) reflects this property more appropriately than (20), we can consider the Störmer-Verlet method has a better stability property than the symplectic Euler method although the two methods have the same stability intervals.
Table 1 presents ϕ and ϕ 100 = max 0≤n ≤100 M * ( θ ) n ║, computed numerically, for several values of θ . This shows that ϕ gives an appropriate approximation to sup n≥0 M * ( θ ) n ║ except θ = 1
PPT Slide
Lager Image
Functions ϕ for the symplectic Euler and Störmer-Verlet methods.
Comparison betweenϕandϕ= max0≤n≤100║M*(θ)n║
PPT Slide
Lager Image
Comparison between ϕ and ϕ = max 0≤n≤100M*(θ)n
IV. STABILITY OF METHODS OF ORDER 3 AND ORDER 4
Method (5) for s = 3 corresponding to the parameter values
PPT Slide
Lager Image
is called Ruth’s method, which is of the order 3 in accuracy. For Ruth’s method, we have
PPT Slide
Lager Image
PPT Slide
Lager Image
The stability interval is
PPT Slide
Lager Image
PPT Slide
Lager Image
≈2.50748 where
PPT Slide
Lager Image
denotes a root of tr M * ( θ ) = −2.
To try to improve Ruth’s method with respect to stability, we consider (5) for s = 4 with
PPT Slide
Lager Image
, which is reduced to
PPT Slide
Lager Image
At first glance, it appears that (29) needs more evaluation of f than (5) with s = 3, but f ( p n+1 ) for the computation of q n+1 is again used for the computation of Q 1 at the next step t = t n+1 . Hence, from the perspective of function evaluation, the work needed for (29) is the same as that for (5) with s = 3 e. g., Ruth’s method. This idea is called first same as last and is often utilized in the numerical analysis of differential equations [2] .
Method (29) is of the order 3 if the parameters satisfy
PPT Slide
Lager Image
These are too complicated to treat. We thus introduce the simplifying condition
PPT Slide
Lager Image
By virtue of this condition, the coefficient of θ 6 in tr M * ( θ ) becomes 0, and the trace is reduced to
PPT Slide
Lager Image
The stability interval becomes
PPT Slide
Lager Image
PPT Slide
Lager Image
which is larger than that of Ruth’s method.
Eqs. (30) and (31) form a system of 6 equations with 7 unknown variables, which has solutions with a free parameter, e.g., b 1 . Letting b 1 = 1/3, we obtain the following :
PPT Slide
Lager Image
We refer to the corresponding method as the stabilized 3rd-order method. In Fig. 2 , the functions ϕ for Ruth’s method and the stabilized 3rd-order method are presented. For θ ≤ 2.37, ϕ for Ruth’s method is smaller than ϕ for the stabilized 3rd-order method, but the latter has finite values up to
PPT Slide
Lager Image
Several symplectic methods of the order 4 are known. Among them, a method of the form (29) corresponding to the parameter values (see, e.g., [4] , p. 109)
PPT Slide
Lager Image
For this method, we have the following:
PPT Slide
Lager Image
PPT Slide
Lager Image
The stability interval is
PPT Slide
Lager Image
PPT Slide
Lager Image
≈ 1.57340, where
PPT Slide
Lager Image
is a root of tr M * ( θ ) = 2 . The stability interval is smaller than that of the symplectic Euler method ( Fig. 2 ).
PPT Slide
Lager Image
Functions ϕ for the three symplectic methods
V. NUMERICAL ILLUSTRATION
To test our numerical method, we consider the sine- Gordon equation
PPT Slide
Lager Image
This equation has the solitary wave solution (see, e.g., [10] , chapter 17).
PPT Slide
Lager Image
By introducing a new variable v = 𝜕u/𝜕t and restricting the space variable x to -5 ≤ x ≤ 5, we get the problem
PPT Slide
Lager Image
where φ 0 ( t ) and φ 1 ( t ) are given so that (38) satisfies (39). Moreover, we apply the method of lines approximation to problem (39) by using a mesh of the form xj = -5 + jΔ x , j = 0,1 …, M , Δx = 10/ M , enotes a positive integer. As usual, we denote approximate functions to u ( t,xj ) and v ( t,xj ) by u j ( t ) and v j ( t ), respectively. By approximating 𝜕 2 u /𝜕 x 2 with the standard central difference scheme, we get a Hamiltonian equation
PPT Slide
Lager Image
where q ( t ) = [ u 1 ( t ), u 2 ( t ), …, u M-1 ( t )] T , p ( t ) = [ v 1 ( t ), v 2 ( t ), …, v M-1 ( t )] T
PPT Slide
Lager Image
PPT Slide
Lager Image
The matrix L Δx has eigenvalues represented as
PPT Slide
Lager Image
By using a linear transform, we change the linear part of (40) into equations of the form
PPT Slide
Lager Image
Since ω M-1 is the largest among ωj ’s, a symplectic method is stable for the linear part of (40) if ω M-1 Δ t is included in the interior of the stability interval. Denoting the stability interval by [0, θ 0 ] we express this condition as
PPT Slide
Lager Image
which gives, M → ∞ a CFL condition
PPT Slide
Lager Image
We now consider time step sizes of the form Δ t = 10/ N , where N is a positive integer, and assume 3 N = 2 M for M and N . Then, since Δ t / Δ x = 3/2, among the specific symplectic methods in Sections 2 and 3, only the stabilized 3rd-order method satisfies the CFL condition (46).
Table 2 shows the errors
PPT Slide
Lager Image
for M = 150, 300, 600, 1200, in the case γ = 1 ⁄2 . It is observed that the numerical solution converges to the exact solution (38) with O x 2 ). For this selection of Δ x and Δ t , the other methods bring no significant numerical results because of overflow.
Numerical results by the stabilized 3rd-order method
PPT Slide
Lager Image
Numerical results by the stabilized 3rd-order method
Acknowledgements
The authors would like to thank Ms. Wakana Tamaru, a student of Nanzan University for her help with checking the mathematical expressions and numerical results presented in the paper.
BIO
Toshiyuki Koto received his B.S. degree in 1984 and M.S. degree in 1986 from Department of Mathematics, the University of Tokyo. In 1992, he received his Ph.D. in Engineering from Nagoya University. From April 1986 to March 1991, he was with Fujitsu Limited; from April 1991 to March 2004, with the University of Electro- Communications; and from April 2004 to March 2009, with Nagoya University. Since April 2009, he has been a professor at Nanzan University. His research interests include numerical analysis and applied mathematics. He is a member of the Mathematical Society of Japan, the Society of Information Processing, and the Japan Society of Industrial and Applied Mathematics.
Eun-Jee Song received her B.S. degree from the Department of Mathematics, Sookmyung Women’s University, in 1984. She earned her M.S. and Ph.D. degrees from the Department of Information Engineering, Nagoya University, Japan in 1988 and 1991, respectively. She was an exchange professor at the Department of Computer Science, the University of Auckland, New Zealand, in 2007. She is currently a full and tenured professor of the Department of Computer Science, Namseoul University, Cheonan, Korea.
References
Feng K. 1985 “On the difference schemes and symplectic geometry,” in Proceedings of the 1984 Beijing Symposium on Differential Geometry and Differential Equations 42 - 58
Hairer E. , Norsett S. P. , Wanner G. 1993 Solving Ordinary Differential Equations I: Nonstiff Problems 2nd ed Springer-Verlag Heidelberg, Germany
Leimkuhler B. , Reich S. 2004 Simulating Hamiltonian Dynamics Cambridge University Press Cambridge, MA
Sanz-Serna J. M. , Calvo M. P. 1994 Numerical Hamiltonian Problems Chapman & Hall London
Hairer E. 1994 “Backward analysis of numerical integrators and symplectic methods,” Annals of Numerical Mathematics 1 107 - 132
Yoshida H. 1993 “Recent progress in the theory and application of symplectic integrators,” Celestial Mechanics and Dynamical Astronomy http://dx.doi.org/10.1007/BF00699717 56 (1-2) 27 - 4
Liu F. Y. , Wu X. , Lu B. K. 2007 “On the numerical stability of some symplectic integrators,” Chinese Astronomy and Astrophysics http://dx.doi.org/10.1016/j.chinastron.2007.04.006 31 (2) 172 - 186
Lopez-Marcos M. A. , anz-Serna J. M. , Skeel R. D. 1996 “An explicit symplectic integrator with maximal stability interval,” in Numerical Analysis: A. R. Mitchell 75th Birthday Volume World Scientific Singapore 163 - 175
Murai D. , Koto T. 2011 “Stability and convergence of staggered Runge-Kutta schemes for semilinear wave equations,” Journal of Computational and Applied Mathematics http://dx.doi.org/10.1016/j.cam.2011.03.020 235 (14) 4251 - 4264
Whitham G. B. 1974 Linear and Nonlinear Wave John Wiley & Sons New York, NY