Pergamon
Nonlinear Analysis, Theory, Methods & Applications, Vol. 30, No. 3, pp. 1887-1892, 1997 Proc. 2nd World Congress of Nonlinear Analysts
PII: S 0 3 6 2 - 5 4 6 X ( 9 7 ) O O 1 5 1 - X
© 1997 Elsevier Science Ltd Printed in Great Britain. All fights reserved 0362-546X/97 $17.00 + 0.00
HAMILTONIAN SYSTEMS AND SYMPLECTIC INTEGRATORS PETER GC)RTZ and RUDOLF SCHERER Institut fiir Praktische Mathematik, Universit£t Karlsruhe, 76128 Karlsruhe, Germany 1
Key words and phrases: Hamiltonian systems, symplectic integrators, Runge-Kutta type methods, symplectic stability, conservation of energy.
1. INTRODUCTION Hamiltonian systems arise in many areas of physics, mechanics and engineering sciences, as well as in pure and applied mathematics. In the last years the interest in integrating these systems numerically has increased strongly. Numerical integrators have to conserve the basic and important properties of Hamiltonian systems, and up to now many symplectic integrators are known. But there are further properties to be studied such as stability and conservation of energy. When integrating dynamical systems it is important to study the properties of the method in comparison with those of the dynamical systems. The numerical methods represent discrete dynamical systems, and hence first, common properties of continuous and discrete dynamical systems in the linear and nonlinear case have to be analyzed. Then mostly a very difficult task is to translate the desired properties to algebraic conditions for the coefficients of the numerical methods. When solving Hamiltonian systems the property of symplecticness is basic and symplectic integrators of Runge-Kutta type have been studied in detail in the last years (cf. Sanz-Serna and Calvo [12], GSrtz and Scherer [4]). Explicit Runge-Kutta methods do not have this property but there exist symplectic Runge-Kutta type methods, which are explicit. In this field one approach was derived from Runge-Kutta methods and the other from NystrSm methods applied to second order systems, and the results are closely related to each other. Stability properties have been studied for second order systems (cf. Hairer [5], van der Houwen and S0mmeijer [6], [7]) but not in detail for Hamiltonian systems (cf. McLachlan [8], Qin and Zhang [10]). Conservation of energy means stability, too. It is known that Runge-Kutta methods with stability matrix M = 0 or with symmetric A-stability function conserve energy when applied to linear Hamiltonian systems, but for the other Runge-Kutta type methods the situation is more complicated (cf. [12]). In this paper some new ideas about symplectic stability and the relation to energy conservation for Runge-Kutta type methods are discussed. As a test equation the harmonic oscillator 1e-mail:
[email protected]
1887
1888
Second World Congress of Nonlinear Analysts
is used, which is stable but not asymptotically stable. A symplectic Runge-Kutta type method holds this stability property at least for sufficiently small stepsize h. Especially, Runge-Kutta methods are emphasized but the results partially hold for other types, too (el. G6rtz [3]). Consider the autonomous Hamiltonian system (cf. Arnold [1]) .
.
OH(p, q) . . Oq ,
4 =
OH(p, q) Op
with Hamiltonian H = H(p, q) being sufficiently smooth. A Runge-Kutta method (c, A, b) will be applied to the general system, a partitioned RungeKutta method (c (1), A (1), b(D; c (2), A (2), b(2)) to the separable system with H(p, q) = T(p)+V(q), a Runge-Kutta-NystrSm method (c, A, b, b) to the separable and partially linear system with T(p) = ½pTKp and a Runge-Kutta method with symmetric A-stability function to the linear system with H(u) = ½uTSu. Hence, if we refer to symplectic Runge-Kutta type methods, keep in mind the four different cases of Hamiltonian systems. The greatest common measure is the system with 1 Tiq, H(p,q) = lpTKp + gq
which in the case K and N positive definite will be reduced to the test equation of the harmonic oscillator
°1)0
with Hamiltonian representing the total energy
H(u) = ~uT u,
(1.2)
which is a conserved quantity, i.e., H(u(t)) = H(u(O)) for all t > 0. Naturally (1.1) is closely related to ~ = -w2q ([5], [6], [7]). A Runge-Kutta type method, which is always assumed to be consistent, applied to the test equation (1.1) with initial value u(0) = u0 yields the discrete system Urn+ 1 =
G(hu)um, m = O, 1,2,...,
(1.3)
where with x = hu the matrix
g l(x) consists of rational functions go(x), in case of explicit methods only of polynomials ([3]). In the limiting case it holds G(0) = I. When the solution u(t) = (p(t), q(t)) T of (1.1) with u(0) = u0 is plotted in the phase (p, q)plane, the parametric curve describes the circle H(p,q) = H(uo). The solution {Urn}m>_0 of (1.3) should have a similar behaviour in the phase (p, q)-plane. But conservation of energy for the discrete system, which means
H(u,,) = H(uo), m = 1,2,... cannot always be expected.
for all initial values u(O) = Uo,
(1.4)
Second World Congress of Nonlinear Analysts
1889
For R u n g e - K u t t a methods it is G(x) = R ( - x J ) , where R(z) is the A-stability function, and especially, the matrix G(x) has the simple configuration gu(x) =
g22(x)
g21(X) =
--g12(X)
= 1 - x2bT(I + x2A2)-lc,
(1.5)
= xbT([ + x2A2)-le,
(1.6)
where e = ( 1 , . . . , 1) T, and the determinant and trace of G(z) satisfy det G(x) = 1 + ((bT e) 2 -- 2bTc) x 2 + O(x 4) for z ---, 0
(1.7)
tr a ( z ) = 2 - 2b~'c x ~ + O(x 4) for ~ -~ 0,
(1.S)
and respectively. The transformation (1.3) is symplectic if and only if det G(x) = 1 for all permissible x
(1.9)
(poles of glj(x) have to be taken in consideration), which immediately follows from condition G T j G = J. Hence, symplectic Runge-Kutta type methods satisfy (1.9). Without property (1.9) there exist at most a finite number of x such that det G(x) = 1, i.e., symplecticness does not hold for x = hv -.-+ O. For a consistent and symplectic Runge-Kutta method the relation (1.7) implies 2bTc = 1, i.e., the quadrature formula (c, b) has at least order two.
2. STABILITY The test equation (1.1) having the eigenvalues #1/2 = +iv is stable but not asymptotically stable for all v. A linear Hamiltonian system never can be asymptotically stable in the origin. The discrete system (1.3) must have the same property as the test equation. It is stable but not asymptotically stable if and only if the eigenvalues A1/2(x) of G(x) satisfy IA1/~(x)l = 1, and .~l(X) : /~2(X) = "4-1 implies G(x) = +I.
(2.1)
This request is studied under the assumption that the method is symplectic, which means that at least (1.9) is satisfied. The eigenvalues A1/2(x) are expressed by trace of G(x). LEMMA 2.1. For a symplectic Runge-Kutta type method the discrete system (1.3) with x = hv is stable but not asymptotically stable if and only if I tr G(x)l <_ 2, and ] tr G(x)l = 2 implies G(x) = +I. PROOF. With (1.9) the eigenvalues of G(x) are given by A,/2(x) = ~1 tr G(x) + i ~ / 1 - ' ~( tr C(x))~.
(2.2)
1890
Second World Congress of Nonlinear Analysts
Hence, condition I tr G(x)l < 2 implies
and I tr G(x)l = 2 implies )h(x) = ~2(x) = -t-1. The reversal is given in a similar manner. The matrix G depends on x = hv, where v is given in (1.1) and the stepsize h of the method has to be chosen. We are interested in determining the size of h, where stability is given. THEOREM 2.2. For a symplectic Runge-Kutta type method the discrete system (1.3) satisfies
Itr G(x)I < 2 for all sufficiently small x = hv > O.
(2.3)
PROOF. From (1.8) the desired statement follows, since a consistent and symplectic RungeK u t t a method satisfies 2brc = 1. For partitioned Runge-Kutta methods and R u n g e - K u t t a Nystrhm methods the proof is given in a similar way ([3]). [q DEFINITION 2.3. For a symplectic Runge-Kutta type method F = {x > 0: Itr G(~)I _< 2, and Itr G(x)l = 2 implies G(x) = : k I }
(2.4)
is called the region of symplectic stability. Treating the test equation (1.1) then the stepsize h of the method has to be chosen such that x = hu E P. The region of symplectic stability r may be a finite or infinite intervall or a junction of intervalls. Since tr G(x) is a polynomial, an explicit Runge-Kutta type method never has an infinite region I'. COROLLARY 2.4. For a symplectic Runge--Kutta method there always holds I" = [0, co) and conservation of energy (1.4) is satisfied. PROOF. Notice the simple configuration of a(x) (see (1.5) and (1.6)) then it follows )~l/2(X) = g l l ( x ) q" ilg2
(x)l,
-- g a(X) + g l(X) = det
and hence,
I tr G(x)l
< I~l(x)l + 1~2(x)l = 2
for all x > 0; further A,(x) = A2(x) = +1 implies G(x) = :i:I. Since
aT(x)a(x) = det G(x)I = I it follows
H(um+l) = ~(G(x)um)TG(x)u,~ --- H(um), which means conservation of energy.
G(x)
= 1,
Second World Congress of Nonlinear Analysts
1891
This result cannot be completely extended to partitioned Runge-Kutta and R u n g e - K u t t a NystrSm methods, namely there is no conservation of energy. Examples show that the discrete solution {urn = (Pm,qm)T}~>_O plotted in the phase (p,q)-plane corresponds to an ellipse if and only if hu E F, in the other case {llumll},~>o grows up very rapidly. EXAMPLES: i) The Gauss-Runge-Kutta methods satisfy i = 0 (cf. Butcher [2]) and conserve energy when applied to linear Hamiltonian systems, hence, according to Corollary 2.4 the stability region is given by F = [0, oc). Consider the one and two stage methods with 1 - ~ x 2 + l-~x tr G(x) = 21 1 ~x 2 1 4 1 T l x 2 and t r G ( x ) = 2 1 + i x 2+1x4, respectively, where tr G ( 2 v ~ ) = - 2
and G ( 2 v ~ ) = - I
for the two stage method.
ii) The trapezoidal rule having a symmetric A-stability function possesses the same region F = [0, oo) as the one stage Gauss method. iii) The explicit partitioned Runge-Kutta method matrix G(x) with
(111111]0{011) (Ruth [11]) yields the
gll(X) ~-- 1, g 2 2 ( X ) ~--- 1 - x 2, g21(x) = -g12(x) = x, and det G ( x ) = l , t r G ( x ) = 2 - x 2, G ( 2 ) - ~ + I , and hence, F = [ 0 , 2 ) . For example, if x = 1 then the points um = Gm(1)Uo, m = 0 , 1 , 2 , . . . , with twisted axes, e.g.,
lie on an ellipse
?A0 ~-- ( 1 , 0 ) T , It1 : (1, 1)T,u2 = (0, 1)T, u3 = (--1,0)T, u4 = (--1,--1)T, us = (0,--1)T, u6 = Uo lie on p2 _
p q + q2 =
1.
iv) The explicit partitioned Runge-Kutta method (McLachlan and Atela [9]) 1-~11-~ 1 1
a
0 1-a a a
0]0 a a
0 a 0 1-a
)
lx/~ , a=
yields the matrix G(x) with
g~l(X) = 1 - 1 ~X 2 , g 2 2 ( X ) --~ g l l ( X ) -~- 1 ( 3 -- 4 a ) X 4, g21(-T) = - - g l 2 ( X ) = X -- 1 ( 1 -- O~)X 3, and hence, det G(x) = 1 and tr G(x) --- 2 - x 2 + 1(3 - 4a)x 4. With the abbreviations /3 = 3 - 4a, 5 = ~ the
+ 1, pl
----
20-~ B ),
P2 =
2(1+~) and p3 ~
--
stability region is given by F = [0, v @ ) t2 ( v ~ ,
v ~ ) = [0, 2.26...) U (4.26..., 4.S2...).
4
1892
SecondWoddCongressof NonlinearAnalysts 3. OUTVIEW
The behaviour of the discrete solutions obtained by partitioned Runge--Kutta and RungeKutta-Nystrbm methods with hu E F should be studied in detail in context with backward error analysis. The elliptical form of the approximate solution in the phase (p, q)-plane indicates that the discrete system conserves energy of a perturbed Hamiltonian problem. Furthermore, the stability investigations should be extended to nonlinear Hamiltonian systems with equilibrium points, for example with Hamiltonian H(p, q) = ~] v~(p~ + q~) + V(q), where P'(q) contains higher order terms of the qi, such as the H6non-Heiles problem or a simple type of the Duffing equation. These topics are investigated in a forthcoming paper. REFERENCES
1. ARNOLD, V. I., Mathematical Methods of Classical Mechanics, Second ed., Springer, New York (1989). 2. BUTCHER, J. C., The Numerical Analysis of Ordinary Differential Equations, John Wiley & Sons, Chichester (1987). 3. G()RTZ, P., Symplektische Stabilit£tstheorie zur numerischen Integration von HamiltonSystemen, Dissertation, Karlsruhe (1995). 4. G()RTZ, P. & SCHERER, R., Reducibility and characterization of symplectic Runge-Kutta methods, ETNA 2, 194-204 (1994). 5. ItAIRER, E., Unconditionally stable methods for second order differential equations, Numer. Math. 32, 373-379 (1979). 6. VAN DER HOUWEN, P. & SOMMEIJER, B. P., Explicit Runge-Kutta(-Nystrbm) methods with reduced phase errors for computing oscillating solutions, SIAM J. Numer. Anal. 24, 595-617 (1987). 7. VAN DER HOUWEN, P. & SOMMEIJER, B. P., Diagonally implicit Runge-Kutta-Nystrbm methods for oscillatory problems, SIAM J. Numer. Anal. 26, 414-429 (1989). 8. MCLACHLAN, R. I., Symplectic integration of Hamiltonian wave equations, Numer. Math. 66, 465-492 (1994). 9. McLACHLAN, R. I. & ATELA, P., The accuracy of symplectic integrators, Nonlinearity 5, 541-562 (1992). 10. QIN, M. Z. & ZHANG, M. Q., Multi-stage symplectic schemes of two kinds of Hamiltonian systems for wave equations, Computers Math. Applic. 9, 51-62 (1990). 11. RUTH, R. D., A canonical integration technique, IEEE Trans. Nucl. Sci., NS-30, 2669-2671 (1983). 12. SANZ-SERNA, J. M. & CALVO, M. P., Numerical ttamiltonian Problems, Chapman & Hall, London (1994).