- Email: [email protected]
Symmetric and symplectic exponentially fitted Runge–Kutta–Nyström methods for Hamiltonian problems
Available online at www.sciencedirect.com
Mathematics and Computers in Simulation 94 (2013) 76–95
Original article
Symmetric and symplectic exponentially fitted Runge–Kutta–Nyström methods for Hamiltonian problems Xiong You a,c,∗ , Bingzhen Chen b a Department of Applied Mathematics, Nanjing Agricultural University, Nanjing 210095, PR China Department of Mathematics, Beijing Jiaotong University Haibin College, Huanghua 061100, PR China c State Key Laboratory for Novel Software Technology at Nanjing University, Nanjing 210093, PR China
b
Received 19 April 2012; received in revised form 11 December 2012; accepted 30 May 2013 Available online 19 June 2013
Abstract The construction of symmetric and symplectic exponentially fitted modified Runge–Kutta–Nyström (SSEFRKN) methods is considered. Based on the symmetry, symplecticity, and exponentially fitted conditions, new explicit modified RKN integrators with FSAL property are obtained. The new integrators integrate exactly differential systems whose solutions can be expressed as linear combinations of functions from the set { exp(± iωt)}, ω > 0, i2 = −1, or equivalently from the set { cos(ωt), sin(ωt)}. The phase properties of the new integrators are examined and their periodicity regions are obtained. Numerical experiments are accompanied to show the high efficiency and competence of the new SSEFRKN methods compared with some highly efficient nonsymmetric symplecti EFRKN methods in the literature. © 2013 IMACS. Published by Elsevier B.V. All rights reserved. MSC: 65L05; 65L06; 65M20; 65N40 Keywords: Runge–Kutta–Nyström method; Symmetry; Symplecticity; Exponential fitting; Hamiltonian system
1. Introduction In this paper we focus on the initial-value problem (IVP) of a system of second-order ODEs in the form y¨ = f (t, y),
y(t0 ) = y0 ,
y˙ (t0 ) = y˙ 0 ,
(1)
whose solution exhibits an oscillatory character, where y ∈ Rd , f : R × Rd → Rd is a smooth function. Problems of this type arise naturally in a variety of applied sciences such as molecular dynamics, orbital mechanics and electronics, chemistry, biology and engineering. High accuracy of integration is often required in these areas. Until now there have been broadly two categories of approaches to numerical integration of the IVP (1): indirect and direct. On one hand,
∗ Corresponding author at: Department of Applied Mathematics, Nanjing Agricultural University, Nanjing 210095, PR China. Tel.: +86 13605184263. E-mail addresses: [email protected] (X. You), [email protected] (B. Chen).
0378-4754/$36.00 © 2013 IMACS. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.matcom.2013.05.010
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
77
by introducing a new variable w to represent the first derivative y˙ , the IVP (1) is turned into a partitioned system of first-order equations y˙ = w,
˙ = f (t, y), w
y(t0 ) = y0 ,
w(t0 ) = y˙ 0 ,
(2)
and the problem can be solved by the general-purpose Runge–Kutta (RK) methods or partitioned Runge–Kutta (PRK) methods (see [3,4,6,10,13,18,19]). On the other hand, Runge–Kutta–Nyström (RKN) methods, initiated by Nyström in 1925, are designed for attacking the second-order problem (1) directly (see [9]). Some variations of the standard RKN methods are given by Butcher [5], Simos [21], Franco [14] and so forth. Recently, Yang et al. [26] proposed a new family of extend Runge–Kutta–Kutta–Nyström (ERKN) methods which have been shown to be very effective and efficient when applied to perturbed oscillators. On many occasions, the problem under consideration takes the form of a Hamiltonian system p˙ = −
∂ U(t, q), ∂q
q˙ = M −1 p
(3)
with the Hamiltonian H(t, p, q) =
1 T −1 p M p + U(t, q), 2
(4)
where M is a symmetric positive definite constant matrix. This system is equivalent to the second-order equation (1) with ∂ f (t, q) = −M −1 ∂q U(t, q). Note that if the system (3) is autonomous, its Hamiltonian is invariant along any solution. Symplecticity is one of the characteristic properties of Hamiltonian systems and it has been widely acknowledged that the so-called symplectic integrators have the advantage of preserving qualitative properties such as the symplecticity of the exact flow over long time when they are applied to Hamiltonian systems [15,19,23]. Among pioneers in the investigation of symplectic integrators are Vogelaere, Ruth and Feng (see Hairer et al. [15]). From then on, identification and analysis of symplectic methods for solving IVPs (1) have received more and more attention. As is known, the symplectic Runge–Kutta–Nyström methods share the pronounced property of being zerodissipative, which is an important requirement for solving oscillatory and Hamiltonian problems. Simos and Vigo-Aguiar [21] considered symplectic methods of RKN type which are adapted to certain types of oscillators. In particular, they have constructed a symplectic modified RKN method that integrates exactly the linear differential equation y¨ = −ω2y , and whose coefficients are frequency dependent. This method has algebraic order two, and it performs much more efficiently than the classical symplectic RKN second-order method introduced by Calvo and Sanz-Serna [6] when it is applied to oscillatory problems (see the numerical results presented in [21]). On the other hand, exponentially fitted methods, which intend to integrate exactly differential equations whose solutions can be expressed as linear combinations of functions { exp(± iωt)}, ω > 0, i2 = −1, or equivalently from the set { cos(ωt), sin(ωt)}, share excellent behavior when applied to oscillatory problems. The construction of exponentially fitted RK(N) methods is originally due to Paternoster [17] and a detailed exposition of exponentially fitted methods with an extensive bibliography on this subject can be found in Ixaru and Vanden Berghe [16]. Recently, Van de Vyver [24] and Franco [14] considered symplectic exponentially fitted modified RKN method. Unfortunately they did not take into account the symmetry conditions. As is pointed out in Chapters V and XI of [15], symmetric methods show a better long time behavior than non symmetric ones when applied to reversible differential systems, as it is the case for conservative mechanical systems. Very recently, Chen et al. [8] derived the symmetry and symplecticity conditions for ERKN (extended RKN) methods solving oscillatory Hamiltonian systems and constructed two examples of symmetric and symplectic ERKN (SSERKN) methods of orders two and four respectively. Motivated by this work, we derive a family of symmetric and symplectic exponentially fitted modified RKN methods in this paper. The rest of this paper is organized as follows: in Section 2 we present the symmetry conditions, symplecticity conditions, exponential fitting conditions and algebraic order (up to third-order) conditions for modified RKN methods. In Section 3 we derive exponentially fitted symmetric and symplectic integrators with two, three and four stages per step. We will also reveal relationship between the symmetry and symplecticity. In Section 4 we analyze the phase and stability properties of the new methods, obtaining generalized periodicity regions for the classical second-order linear test model. In Section 5 we carry out five numerical experiments to show the efficiency of the new methods compared with the symplectic integrators given in [6,14,21,24]. Section 6 is devoted to conclusive comments.
78
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
2. Conditions for symmetry, symplecticity, exponential fitting and consistency of modified RKN methods In this paper, we consider the method (5) whose coefficients are z-dependent, like exponentially fitted (EF)-type or other adapted-type methods (see, for example, [3,4,7,20,21]). We consider a class of RKN-type methods studied by Butcher [5], Simos and Vigo-Aguiar [21], and Van de Vyver [24]. Definition 2.1. An s-stage modified RKN method for the second-order IVP (1) is a scheme of the form Yi = y0 + ci (z)γi (z)h˙y0 + h2 y1 =
s
aij (z)f (t0 + cj (z)h, Yj ), j=1 s g1 (z)y0 + g2 (z)h˙y0 + h2 bi (z)f (t0 + ci (z)h, Yi ), i=1 s
y˙ 1 = g3 (z)˙y0 + h
i = 1, . . . , s, (5)
bi (z)f (t0 + ci (z)h, Yi ),
i=1
where the nodes ci (z), i = 1, . . ., s and the coefficients g1 (z), g2 (z), g3 (z) and γi (z), bi (z), bi (z), aij (z), i, j = 1, . . ., s are assumed to be even functions of z = ihω with ω > 0 the principal frequency of the problem. Scheme (5) can be expressed by the following Butcher tableau c(z)
e g1 (z)
γ(z) g2 (z) g3 (z)
c1 (z) .. A(z) . ¯bT (z) = cs (z) bT (z)
1 .. .
1 g1 (z)
γ1 (z) .. . γs (z) g2 (z) g3 (z)
a11 (z) .. . as1 (z) ¯b1 (z) b1 (z)
... .. . ... ... ...
a1s (z) .. . ass (z) ¯bs (z) bs (z)
If the matrix A(z) is strictly lower triangular, i.e., aij (z) = 0, 1 ≤ i < j ≤ s, the method (5) is explicit, otherwise it is implicit. If g1 (z) = 1, cs (z)γ s (z) = g2 (z) and bj (z) = asi (z), j = 1, . . . , s, the method is said to have the FSAL-property (the last evaluation at any step is the same as the first evaluation at the next step). A s-stage modified RKN method having the FSAL-property actually requires only s − 1 function evaluations at each step. In convention, we assume that as z → 0, the coefficients g1 (z), g2 (z), g3 (z), γ i (z), i = 1, . . ., s all approach 1, so that the limit method of Scheme (5) coincides with a classical s-stage RKN formula. The objective of this section is to find out when the modified RKN method (5) is symmetric, symplectic and exponentially fitted and to construct practical methods with these favorable properties. 2.1. Symmetry conditions The key to understanding symmetry is the concept of the adjoint method. Definition 2.2. The adjoint method h * of a one-step method h is the inverse map of the original method with reversed time step −h, i.e., h * : = −h −1 . A method for which h * = h is said to be symmetric. (See Hairer [15].) Theorem 2.1. The modified RKN method (5) is symmetric if its coefficients satisfy the following conditions ci (z) + cs+1−i (z) = 1, g1 (z) = g3 (z) = 1, ci γi (z) + cs+1−i (z)γs+1−i (z) = g2 (z), bi (z) = bs+1−i (z), bi (z) + bs+1−i (z) = g2 (z)bs+1−i (z), aij (z) − ci (z)γi (z)bj (z) = as+1−i,s+1−j (z) − bs+1−j (z), for all i, j = 1, . . ., s.
(6)
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
79
Proof. Following the idea of Hairer et al. [15], we interchange h ↔ − h, z = iωh ↔ − z = − iωh, t0 ↔ t0 + h, y0 ↔ y1 and y˙ 0 ↔ y˙ 1 , and obtain the following equations which define the adjoint method of (5)
Yi = y1 − ci (z)γi (z)h˙y1 + h2 y0 =
s
aij (z)f (t0 + (1 − cj (z))h, Yj ), j=1 s g1 (z)y1 − g2 (z)h˙y1 + h2 bi (z)f (t0 + (1 − ci (z))h, Yi ), i=1 s
y˙ 0 = g3 (z)˙y1 − h
i = 1, . . . , s, (7)
bi (z)f (t0 + (1 − ci (z))h, Yi ).
i=1
From Eqs. (7), we have bi (z) 1 y˙ 0 + h f (t0 + (1 − ci (z))h, Yi ) g3 (z) g (z) i=1 3 s 1 bs+1−i (z) = y˙ 0 + h f (t0 + (1 − cs+1−i (z))h, Ys+1−i ), g3 (z) g3 (z) s
y˙ 1 =
(8)
i=1
bi (z) 1 g2 (z) y0 + h˙y1 − h2 f (t0 + (1 − ci (z))h, Yi ) g1 (z) g1 (z) g (z) i=1 1 s 1 bs+1−i (z) g2 (z) = y0 + h˙y1 − h2 f (t0 + (1 − cs+1−i (z))h, Ys+1−i ) g1 (z) g1 (z) g1 (z) i=1 s 1 bs+1−i (z) g2 (z) 1 h y˙ 0 + h f (t0 + (1 − cs+1−i (z))h, Ys+1−i ) y0 + = g1 (z) g1 (z) g3 (z) g3 (z) s
y1 =
i=1
s 1 bs+1−j (z)f (t0 + (1 − cs+1−i (z))h, Ys+1−i ) −h2 g (z) j=1 1 1 g2 (z) = y0 + h˙y0 g1 (z) g1 (z)g3 (z) s 1 g2 (z) +h2 ( bs+1−i (z) − bs+1−i (z))f (t0 + (1 − cs+1−i (z))h, Ys+1−i ) g1 (z) g3 (z)
(9)
i=1
and
Ys+1−i = y1 − cs+1−i (z)γs+1−i (z)h˙y1 + h2
s
as+1−i,s+1−j (z)f (t0 + (1 − cj (z))h, Yj ),
j=1
=
1 1 g2 (z) y0 + ( − cs+1−i (z)γs+1−i (z))hy0 g1 (z) g3 (z) g1 (z) s 1 g2 (z) 1 +h2 ( bs+1−j (z) ( − cs+1−i (z)γs+1−i (z))bs+1−j (z) − g3 (z) g1 (z) g1 (z) j=1
+as+1−i,s+1−j (z)) · f (t0 + (1 − cs+1−j (z))h, Ys+1−j ).
(10)
80
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
Denoting Ys+1−i = Y˜ i , we know that the scheme defined by Eqs. (8)–(10) coincides with Scheme (5) if the following conditions are satisfied 1 − cs+1−i (z) = ci (z), 1 = 1, g1 (z) 1 g2 (z) − cs+1−i (z)γs+1−i (z) = ci (z)γi (z), g3 (z) g1 (z) 1 g2 (z) 1 bs+1−j (z) − cs+1−i (z)γs+1−i (z) bs+1−j (z) − g3 (z) g1 (z) g1 (z) +as+1−i,s+1−j (z) = aij (z), 1 = g1 (z), g1 (z) g2 (z) = g2 (z), g1 (z)g3 (z) 1 g2 (z) ( bs+1−i (z) − bs+1−i (z)) = bi (z), g1 (z) g3 (z) 1 = g3 (z), g3 (z) bs+1−i (z) = bi (z). g3 (z)
(11)
It is easy to verify that Eqs. (11) are equivalent to Eqs. (6). Corollary 2.1. As z → 0, the conditions (6) reduce to the symmetry conditions for traditional RKN methods with constant coefficients ci + cs+1−i = 1, bi = bs+1−i , bi + bs+1−i = bi ,
(12)
aij − ci bj = as+1−i,s+1−j − bs+1−j , for i, j = 1, . . ., s . 2.2. Symplecticity conditions Now we turn to the symplecticity conditions for Scheme (5). Recall that a one-step method h : (q0 , p0 )T → (q1 , p1 )T is symplectic if for every smooth Hamiltonian function H and for every step size h, the numerical flow preserves the differential 2-form dp ∧ dq = di=1 dpi ∧ dqi , i.e., dp1 ∧ dq1 = dp0 ∧ dq0 . In analog, for the second-order system (1) as a Hamiltonian system, the RKN scheme (5) h : is said to be symplectic if dy1 ∧ d˙y1 = dy0 ∧ d˙y0 . From Tocino and Vigo-Aguiar [22], we have the following theorem.
(y0 , y˙ 0 )T → (y1 , y˙ 1 )T (13)
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
81
Theorem 2.2. The modified RKN method (5) is symplectic if its coefficients satisfy the following conditions g1 (z)g3 (z) = 1, g3 (z)bi (z) + (g1 (z)ci (z)γi (z) − g2 (z))bi (z) = 0,
i = 1, . . . , s,
bi (z)(bj (z) − g1 (z)aij (z)) = bj (z)(bi (z) − g1 (z)aji (z)),
(14)
1 ≤ j < i ≤ s.
When g1 (z) = g3 (z) = 1, Eqs. (14) become bi (z) = (g2 (z) − ci (z)γi (z))bi (z),
i = 1, . . . , s, (15)
bi (z)(bj (z) − aij (z)) = bj (z)(bi (z) − aji (z)),
1 ≤ j < i ≤ s.
Corollary 2.2. As z → 0, the conditions (14) or (15) reduce to the symplectic conditions for traditional RKN methods with constant coefficients bi = (1 − ci )bi ,
i = 1, . . . , s, (16)
bi aij − bj aji = bi bj (ci − cj ),
1 ≤ j < i ≤ s.
We refer to Eqs. (6) and (15) as SS conditions. 2.3. Relationship between symmetry, symplecticity and FSAL-property of explicit methods
Theorem 2.3. If an s-stage explicit modified RKN method (5) is symmetric and the sixth equation in (6) is also true for i = j, then the method is symplectic. Proof. The first condition in (14) is satisfied immediately. We only need to check Eqs. (15) for i > j. The first equation in (15) follows from the third and sixth equations in (14) with j = i. On the other hand, with the third and fourth equations in (14) and the first equation in (15), for i > j, the sixth equation in (6) becomes aij (z) = ci (z)γi (z)bj (z) − bs+1−j (z) = ci (z)γi (z)bs+1−j (z) − cj (z)γj (z)bs+1−j (z) = (ci (z)γi (z) − cj (z)γj (z))bs+1−j (z) = (cs+1−j (z)γs+1−j (z) − cs+1−i (z)γs+1−i (z))bj (z). Multiplying both sides of the above equation by bi (z), we obtain bi (z)aij (z)
= bi (z)cs+1−j (z)γs+1−j (z)bj (z) − bj (z)cs+1−i (z)γs+1−i (z)bi (z) = bi (z)bj (z) − bj (z)bi (z).
This is exactly the second equation in (15) for explicit methods with i > j. According to Theorem 2.3, Eqs. (6) for 1 ≤ j ≤ i ≤ s become the SS conditions for the explicit modified RKN method (5). The following theorem shows the connection between symmetry and the FSAL-property.
82
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
Theorem 2.4. Suppose that an s-stage explicit modified RKN method (5) is symmetric and sixth equation in (6) is also true for j = i, i = 1, . . ., s. Suppose c1 (z) = 0 and bi (z) = / 0, i = 1, . . ., s. Then the method has the FSAL-property. Proof. By Theorem 2.3, Eqs. (6) for 1 ≤ j ≤ i ≤ s imply Eqs. (15). From the third equation in (6) and the first equation in (15) we have bs (z) = c1 (z)γ1 (z)bs (z) = 0. Then in the case of i = s, the second equation in (15) gives asj (z) = bj (z), j = 1, . . . , s − 1. By imposing i = j = 1 in the sixth equation in (6), we obtain c1 (z)γ 1 (z) = 0 which implies that cs (z)γ s (z) = g2 (z) by the third equation in (6). Hence, the method satisfies the FSAL-property. 2.4. Exponential fitting conditions Following Albrecht’s lines in [1,2], each stage of the explicit scheme (5) can be viewed as a linear multistep method on a non-equidistant grid. Exponential fitting of the modified RKN method (5) can be presented in terms of a particular type of linear operators on C2 ([t0 , T]), the space of second-order continuously differentiable functions on the interval [t0 , T]. With each stage we associate a linear operator as follows: • for the internal stages,
Li [A(z), h]y(x) = y(x + ci h) − y(x) − ci (z)γi h˙y(x) − h2
s
aij (z)¨y(x + cj (z)h),
i = 1, . . . , s;
j=1
• for the final stages,
L[b(z), h]y(x) = y(x + h) − g1 (z)y(x) − hg2 (z)˙y(x) − h2 L[b(z), h]y(x) = y˙ (x + h) − g3 (z)˙y(x) − h
s
s
bi (z)¨y(x + ci (z)h),
i=1
bi (z)¨y(x + ci (z)h).
i=1
Definition 2.3. The modified RKN method (5) is said to be exponentially fitted if the space span{ exp(± iωt)} is a subspace of the kernels kerL and kerLi , i = 1, . . . , s. By definition, requiring the modified RKN method (5) to be exponentially fitted is equivalent to requiring the operators L and Li , i = 1, . . . , s to vanish for the functions from the set { exp(± iωt)} leading to the following equations exp(±ci (z)z) = 1 ± ci (z)γi (z)z + z2 exp(±z) = g1 (z) ± g2 (z)z + z2 exp(±z) = g3 (z) ± z
s i=1
s
s
aij (z) exp(±cj (z)z),
i = 1, . . . , s,
j=1
bi (z) exp(±ci (z)z),
i=1
bi (z) exp(±ci (z)z),
z = iωh.
(17)
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
83
Using the definitions cosh(z) = (exp(z) + exp(− z))/2 and sinh(z) = (exp(z) − exp(− z))/2, the equations (17) imply that s
aij (z) cosh(cj (z)z) =
cosh(ci (z)z) − 1 , z2
aij (z) sinh(cj (z)z) =
sinh(ci (z)z) − ci (z)γi (z)z , z2
bi (z) cosh(ci (z)z) =
cosh(z) − g1 (z) , z2
j=1
s j=1 s i=1
s i=1
s
(18)
sinh(z) − zg2 (z) bi (z) sinh(ci (z)z) = , z2 bi (z) sinh(ci (z)z) =
cosh(z) − g3 (z) , z
bi (z) cosh(ci (z)z) =
sinh(z) . z
i=1
s
i = 1, . . . , s,
i=1
Eqs. (18) are referred to as the EF conditions. 2.5. Algebraic order conditions In this section, we will present algebraic order conditions for the modified RKN method (5). Like the case of a classical RKN method, for an modified RKN method, the local truncation errors in the approximations of the solution and its derivative can be expressed as ⎛ ⎞ kj p−1 e1 = y(x0 + h) − y1 = hj+1 ⎝ di (j+1) F (j) (τji )(y0 )⎠ + O(hp+1 ), i=1
j=1
e˙ 1 = y˙ (x0 + h) − y˙ 1 =
p
kj (j) j h ( d˙ i F (j) (τji )(y0 )) + O(hp+1 ),
j=1
i=1
where F(j) (τ ji )(y0 ) denotes the elementary differential associated to the Nyström tree τ ji of order j, kj is the number of (j) Nyström trees of order j, and the terms di (j+1) and d˙ i depend on the coefficients of the EFRKN method. A method (5) has order p if, for every sufficiently smooth IVP (1) and for every small step size h, the local truncation errors of the numerical solutions satisfy e1 = y(t0 + h) − y1 = O(hp+1 ), e˙ 1 = y˙ (t0 + h) − y˙ 1 = O(hp+1 ), or equivalently, di (j+1) = 0, (j) d˙ i = 0,
i = 1, . . . , kj ,
i = 1, . . . , kj ,
j = 1, . . . , p − 1,
j = 1, . . . , p.
Following the idea of Franco [14], we assume the following expansions bi (z) = bi
(0)
(2)
(4)
+ bi h2 + bi h4 + · · ·,
γi (z) = 1 + γi (2) h2 + γi (4) h4 + · · ·, g1 (z) = 1 + g1 (2) h2 + g1 (4) h4 + · · ·, g3 (z) = 1 + g3 (2) h2 + g3 (4) h4 + · · ·.
bi (z) = bi (0) + bi (2) h2 + bi (4) h4 + · · ·,
aij (z) = aij (0) + aij (2) h2 + aij (4) h4 + · · ·, g2 (z) = 1 + g2 (2) h2 + g2 (4) h4 + · · ·,
84
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
In the sequel, we only consider methods with constant nodes. Then the order conditions up to order three for the modified RKN method (5) can be listed as follows: Order 1 requires:
(1) d˙ 1 :=
bi (0) − 1 = 0.
(19)
i
Order 2 requires in addition:
d1 (2)
1 (2) = 0, d˙ 2 := g3 (2) = 0, 2 i (0) 1 := bi − = 0, d2 (2) := g1 (2) = 0. 2
(2) d˙ 1 :=
bi (0) ci −
(20)
i
Order 3 requires in addition:
(3) d˙ 1 (3) d˙ 3
1 (0) 2 1 1 (3) := bi ci − bi (0) aik (0) − = 0, = 0, d˙ 2 := 2 3 6 i i k (0) 1 := bi (2) = 0, d1 (3) := bi ci − = 0, 6 i
(21)
i
d2 (3) := g2 (2) = 0.
From Theorem 2.1 in [14], we know that the modified RKN method (5) has algebraic order at least 2.
3. Construction of symmetric and symplectic explicit EFRKN methods In this section we proceed to construct practical explicit EFRKN methods with constant nodes under the symmetry, symplecticity, exponential fitting and algebraic order conditions obtained in the previous section. We refer to these methods as SSEFRKN methods.
3.1. The case s = 2 For s = 2, the SS conditions (6) for 1 ≤ j ≤ i ≤ 2 become c1 + c2 = 1,
g1 (z) = g3 (z) = 1,
b1 (z) + b2 (z) = g2 (z)b1 (z),
c1 γ1 (z) + c2 γ2 (z) = g2 (z),
b1 (z) = c2 γ2 (z)b1 (z),
a21 (z) = c2 γ2 (z)b1 (z) − b2 (z),
b1 (z) = b2 (z),
b2 (z) = c1 γ1 (z)b2 (z),
0 = a21 (z) + c1 γ1 (z)b2 (z) − b1 (z).
(22)
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
85
On the other hand, for s = 2 the EF conditions (18) give cosh(c1 z) − 1 sinh(c1 z) − c1 γ1 (z)z = 0, = 0, 2 z z2 cosh(z) − 1 a21 (z) cosh(c1 z) = , z2 sinh(z) − c2 γ2 (z)z a21 (z) sinh(c1 z) = , z2 cosh(z) − g1 (z) b1 (z) cosh(c1 z) + b2 (z) cosh(c2 z) = , z2 sinh(z) − g2 (z) b1 (z) sinh(c1 z) + b2 (z) sinh(c2 z) = , z2 cosh(z) − g3 (z) b1 (z) cosh(c1 z) + b2 (z) cosh(c2 z) = , z sinh(z) b1 (z) cosh(c1 z) + b2 (z) cosh(c2 z) = . z
(23)
Solving (22) and (23), we obtain g1 (z) = g3 (z) = 1,
c1 = 0,
c2 = 1,
γ2 (z) =
sinh(z) , z
g2 (z) = γ2 (z),
cosh(z) − 1 , b1 (z) = a21 (z), b2 (z) = 0, z2 sinh(z) b1 (z) = , b2 (z) = b1 (z). z(cosh(z) + 1)
a21 (z) =
(24)
The two-stage method of order two is denoted by SSEFRKN2. It is exactly the symplectic exponentially fitted integrator EFRKN2 presented in [24]. However, unlike [24], we do not artificially assume c2 = 1 to find the coefficients.
3.2. The case s = 3 In a similar way we can obtain the coefficients of a three-stage SSERKN method of order two c1 = 0,
c2 =
1 , 2
g1 (z) = g3 (z) = 1,
c3 = 1, g2 (z) = γ3 (z) = γ2 (z) =
2 sinh(z/2) , z
cosh(z/2) − 1 , z2 2(cosh(z/2) − 1) a31 (z) = , a32 (z) = a31 (z), z2 b1 (z) = a31 (z), b2 (z) = a32 (z), b3 (z) = 0, a21 (z) =
b1 (z) =
cosh(z/2) − 1 , z sinh(z/2)
b2 (z) = 2b1 (z),
We denote this method as SSEFRKN3.
b3 (z) = b1 (z).
(25)
86
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
3.3. The case s = 4 Our main interest is in the case s = 4. The SS conditions (6) for 1 ≤ j ≤ i ≤ 4 become
g1 (z) = g3 (z) = 1,
c1 + c4 = 1,
c1 γ1 (z) + c4 γ4 (z) = g2 (z), b1 (z) = b4 (z),
c2 + c3 = 1,
c2 γ2 (z) + c3 γ3 (z) = g2 (z),
b2 (z) = b3 (z),
b1 (z) + b4 (z) = g2 (z)b1 (z), b1 = c4 γ4 (z)b4 (z),
b2 (z) + b3 (z) = g2 (z)b2 (z),
b2 = c3 γ3 (z)b3 (z),
b3 = c2 γ2 (z)b2 (z),
(26) b4 = c1 γ1 (z)b1 (z),
a21 (z) = c2 γ2 (z)b1 (z) − b4 (z),
a31 (z) = c3 γ3 (z)b1 (z) − b4 (z),
a32 (z) = c3 γ3 (z)b2 (z) − b3 (z),
a41 (z) = c4 γ4 (z)b1 (z) − b4 (z),
a42 (z) = c4 γ4 (z)b2 (z) − b3 (z),
a43 (z) = c3 γ4 (z)b3 (z) − b2 (z).
On the other hand, the EF conditions (18) for s = 4 can be become
cosh(c1 z) − 1 = 0, z2
sinh(c1 z) − c1 γ1 (z)z = 0, z2
a21 (z) cosh(c1 z) =
cosh(c2 z) − 1 , z2
a21 (z) sinh(c1 z) =
sinh(c2 z) − c2 γ2 (z)z , z2
a31 (z) cosh(c1 z) + a32 (z) cosh(c2 z) = a31 (z) sinh(c1 z) + a32 (z) sinh(c2 z) =
cosh(c3 (z)z) − 1 , z2 sinh(c3 z) − c3 γ3 (z)z , z2
a41 (z) cosh(c1 z) + a42 (z) cosh(c2 z) + a43 (z) cosh(c3 z) = a41 (z) sinh(c1 z) + a42 (z) sinh(c2 z) + a43 (z) sinh(c3 z) =
cosh(c4 z) − 1 , z2
sinh(c4 z) − c4 γ4 (z)z . z2
cosh(z) − 1 , z2 sinh(z) − g2 (z)z b1 (z) sinh(c1 z) + b2 (z) sinh(c2 z) + b3 (z) sinh(c3 z) = , z2 sinh(z) b1 (z) cosh(c1 z) + b2 (z) cosh(c2 z) + b3 cosh(c3 z) + b4 (z) cosh(c4 z) = , z cosh(z) − 1 b1 (z) sinh(c1 z) + b2 sinh(c2 z) + b3 (z) sinh(c3 z) + b4 (z) sinh(c4 z) = . z b1 (z) cosh(c1 z) + b2 (z) cosh(c2 z) + b3 (z) cosh(c3 z) =
(27)
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
87
We solve the system of Eqs. (26) and (27) for the coefficients with c2 as a parameter as follows c1 = 0,
c4 = 1,
c3 = 1 − c2 ,
γ2 (z) =
sinh(c2 z) , c2 z
a21 (z) =
g1 (z) = g3 (z) = 1,
cosh(c2 z) − 1 , z2
z sinh(z/2) − b1 (z)z cosh( ) a21 (z) 2 b1 (z) = , , b2 (z) = c2 γ2 (z) z cosh(z(c3 − c2 )/2) b3 (z) = c2 γ2 (z)b3 (z), N0 b1 (z)z − N1 , b1 (z)b2 (z)z sinh(c2 z) − b2 (z) cosh(c2 z) N1 − g2 (z)b2 (z) cosh(c2 z) γ3 (z) = c3 b1 (z) a31 (z) = c3 γ3 (z)b1 (z),
g2 (z) =
b3 (z) = b2 (z),
b4 (z) = b1 (z),
γ4 (z) = g2 (z),
(28)
a32 (z) = g2 (z)b2 (z) − 2c2 γ2 (z)b2 (z), b1 (z) = g2 (z)b1 (z), a41 (z) = b1 (z),
b2 (z) = c3 γ3 (z)b2 (z),
a42 (z) = b2 (z),
a43 (z) = b3 (z).
where N0 =
sinh(c3 z) + 2b3 (z) sinh(c2 z), z2
N1 =
cosh(c3 z) − 1 + 2b3 (z) cosh(c2 z). z2 √ 3
√ 3
Substituting ci , bi (z), i = 1, . . ., 4 into the first equation in (21) are satisfied, we obtain c2 = 32 + 64 + 23 ≈ 1.351207191959658. With this value of c2 , it is easy to check that (19), (20) and the other equations in (21) are satisfied. Thus the method defined by (28) is of order three. Since a symmetric method must be of an even order, then the method is actually of order 4 (see Chapter II of [15]). We denote this method as SSEFRKN4. It should be noted that the new methods SSEFRKN2 and SSEFRKN3, SSEFRKN4 derived in this section all have the FSAL-property.
4. Phase-lags and periodicity regions of the new methods Now we set out to analyze the phase and stability properties of our new methods. Stability means that the numerical solutions remain bounded as we move further away from the starting point. For classical RKN methods, the phase and stability properties are checked using the second-order linear test model y¨ (t) = −μ2 y(t),
μ > 0.
(29)
Recall that the new symmetric and symplectic exponentially fitted modified RKN methods (SEFRKN) derived in the previous section are dependent on the complex number λ = iω, where ω > 0 is an estimate of the principal frequency. Applying an s-stage SEFRKN method (5) to the test model (29) yields
y1 h˙y1
= M(η , ν ) 2
2
y0 h˙y0
,
(30)
88
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
where
M(η2 , ν2 ) =
g1 (z) − η2 b(z)T N(z)−1 e
g2 (z) − η2 b(z)T N(z)−1 c(z)γ(z)
, −η2 b(z)T N(z)−1 e g3 (z) − η2 b(z)T N(z)−1 c(z)γ(z) ν = ωh, η = μh, N(z) = I + η2 A(z), b(z)T = (b1 (z), . . . , bs (z)), bT (z) = (b1 (z), . . . , bs (z)), ⎛ ⎞ 0 0 ... 0 ⎛ ⎞ ⎛ ⎞ c (z)γ (z) 1 1 1 ⎜ a21 (z) 0 ... 0⎟ ⎜ ⎟ ⎜ ⎟ ⎜ .. ⎟ .. ⎟ , A(z) = ⎜ ⎟. e = ⎝ . ⎠ , c(z)γ(z) = ⎜ ⎜ ⎝ ⎠ . . . .. .. .. ⎟ ⎝ .. ⎠ . . 1 cs (z)γs (z) as1 (z) . . . as,s−1 (z) 0
The phase and stability behavior of the numerical solution depends on the eigenvalues or the spectrum of the stability matrix M = M(η2 , ν2 ). Eliminating y˙ 0 and y˙ 1 from (30) and the equation that is obtained from (30) by replacing the subscript 0 by 1 gives the difference equation y2 − tr(M(η2 , ν2 ))y1 + det(M(η2 , ν2 ))y0 = 0. Accordingly, the characteristic equation is given by ξ 2 − tr(M)ξ + det(M) = 0,
(31)
where tr(M) and det(M) are the trace and the determinant of M, respectively. Definition 4.1. For the characteristic Eq. (31) obtained by applying an exponentially fitted method to Eq. (29), set the ratio r = ω/μ. The two quantities tr(M) √ φ(η, r) = η − arccos , 2 det(M) √ d(η, r) = 1 − det(M) are called the dispersion error and the dissipation error of the numerical method under consideration, respectively. The method is said to be dispersive of order p and dissipative of order q, if φ(η, r) = O(ηp+1 ) and d(η, r) = O(ηq+1 ), respectively. The method is called zero-dissipative if d(η, r) = 0. For the method SSEFRKN2 method, we have r2 − 1 3 η + O(η5 ), 24 which shows that SSEFRKN2 is dispersive of order 2. For the method SSEFRKN3, we have φ(η, r) =
r2 − 1 3 η + O(η5 ), 96 which shows that SSEFRKN3 is dispersive of order 2. For the method SSEFRKN4 method, we have φ(η, r) =
12c2 4 − 18c2 3 + 6c2 − 1 2 3 r η + O(η5 ), 24 which shows that SSEFRKN4 is dispersive of order 2. It is noted that for the methods SSEFRKN2, SERRKN3 and SSEFRKN4, det(M) ≡ 1, since they are all symplectic. Therefore these methods are all zero-dissipative. The stability properties of an EFRKN method are characterized by the spectral radius ρ(M). Having in mind that the matrix M depends on the variables η and ν, the periodicity and stability intervals typical in the classical RKN methods now become two-dimensional regions in the (H, ν)-plane. So, for the new methods we use an adaptation of the terminology introduced by Coleman and Ixaru [10]: φ(η, r) =
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95 SSEFRKN2
30
SSEFRKN4
8
SSEFRKN3
30
89
25
20
20
5
15
15
4
6
ν
ν
ν
7 25
3 10
10
5
5
2
0 0
5
10
15 η
20
25
30
0 0
1 5
(a)
10
15 η
(b)
20
25
30
0 0
1
2
3
4
η
5
6
7
8
(c)
Fig. 1. Periodicity regions for the methods (a) SSEFRKN2, (b) SSEFRKN3 and (c) SSEFRKN4.
(i) (ii) (iii) (iv)
Rs = {η > 0, ν > 0|ρ(M) < 1} is the region of stability; Rp = {η > 0, ν > 0|ρ(M) = 1andtr(M)2 < 4det(M)} is the region of periodicity; If Rs = (0, ∞) × (0, ∞) except possibly for a discrete set of curves, the method is A-stable; If Rp = (0, ∞) × (0, ∞) except possibly for a discrete set of curves, the method is P-stable.
Note again that the symplecticity implies that det(M) ≡ 1, and hence symplectic methods have nonempty regions of periodicity. The regions of periodicity of the methods SSEFRKN2, SSEFRKN3 and SSEFRKN4 are depicted in Fig. 1.
5. Numerical experiments To test the numerical performance of the new modified RKN methods, we carry out experiments on five model problems to illustrate the effectiveness and efficiency of the new methods SSEFRKN3 and SSEFRKN4. The codes used in the comparisons are: • SIMOS2: The symplectic modified two-stage RKN method of order two proposed by Simos and Vigo-Aguiar [21]. • SSEFRKN2: The symplectic exponentially fitted modified two-stage RKN method of order two with FSAL-property proposed by Van de Vyver [24], that is, the method given by (24). • EFRKNS23F: The symplectic exponentially fitted modified three-stage RKN method of order two with FSALproperty proposed by Franco [14]. • EFRKNS4: The symplectic exponentially fitted modified three-stage RKN method of order four proposed by Franco [14]. • SSEFRKN3: The symmetric and symplectic exponentially fitted modified three-stage RKN method of order two with FSAL-property proposed in this paper. • SSEFRKN4: The symmetric and symplectic exponentially fitted modified four-stage RKN method of order four with FSAL-property proposed in this paper. The criterion used in the numerical comparisons is efficiency in terms of the decimal logarithm of the maximum global error (log10 (MGE)) versus the computational effort measured in the number of function evaluations required by each method. Problem 1. Consider the “almost periodic” orbit problem (studied by Franco and Palacios [12]) z¨ + z = ε exp(iψx),
z(0) = 1,
z˙ (0) = i,
90
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
Fig. 2. Efficiency curves for solution in (a) Problem 1, (b) Problem 2 and (c) Problem 3.
where ε = 0.001 and ψ = 0.01. The analytical solution z(x) = u(x) + iv(x) is given by ε 1 − ε − ψ2 cos(x) + cos(ψx), 2 1−ψ 1 − ψ2 1 − εψ − ψ2 ε v(x) = sin(x) + sin(ψx). 2 1−ψ 1 − ψ2 u(x) =
The problem is a nonhomogeneous Hamiltonian system H(p, q) =
1 1 (p1 2 + p2 2 ) + (q1 2 + q2 2 ) − ε(cos(ψt)q1 + sin(ψt)q2 ), 2 2
where q = (q1 , q2 )T = (u, v)T , p = (p1 , p2 )T = (˙q1 , q˙ 2 )T . The initial values are q(0) = (1, 0)T , q˙ (0) = (0, 1)T . We integrate the problem on the interval [0, 1000] with fitting frequency ω = 1. The stepsizes are taken as h = 1/2j , j = 1, . . ., 5 for the methods SIMOS2, SSEFRKN2, SSEFRKN3, and 3h/2 for the other three methods EFRKNS23F, EFRKNS4 and SSERKN4. We will follow this convention in the following experiments. The numerical results are displayed in Fig. 2(a).
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
91
Problem 2. Consider the nonhomogeneous Duffing equation (studied in [21]) q + q¨ + q3 = 0.002 cos(1.01t),
q(0) = 0.200426728069666,
q˙ (0) = 0.
The Hamiltonian of this system is 1 2 1 2 1 4 p + q + q − 0.002 cos(1.01t)q. 2 2 4 An “exact” solution was first given by Van Dooren [11], but a more accurate solution is given by Zhao et al. [27] as H(p, q) =
q(t) =
4
A2i+1 cos[(2i + 1)1.01t],
i=0
with (A1 , A3 , A5 , A7 , A9 ) = (0.20017947753661502, 2.46946143255559 × 10−4 , 3.0401498519692437 × 10−7 , 3.743490701609247 × 10−10 , 4.609682949622697 × 10−13 ). The equation is integrated on the interval [0, 1000] with ω = 1, h = 1/2j , j = 0, . . ., 4. The numerical results are displayed in Fig. 2(b). Problem 3. Consider the perturbed orbit problem (studied in [14]) q¨ 1 = − q¨ 1 = −
q1 (q1 2 + q2 q1
2 )3/2
(q1 + q2
2 )3/2
2
− −
(2ε + ε2 )q1 5/2
(q1 2 + q2 2 ) (2ε + ε2 )q1 (q1 + q2 2
2 )5/2
,
q1 (0) = 1,
q1 (0) = 0,
,
q2 (0) = 0,
q2 (0) = 1 + ε.
with ε = 0.001 . The Hamiltonian of this problem is given by H(p, q) =
1 1 2ε + ε2 (p1 2 + p2 2 ) − − . 1/2 3/2 2 (q1 2 + q2 2 ) 3(q1 2 + q2 2 )
The analytic solution of the problem is given by q1 (t) = cos((1 + ε)t),
q2 (t) = sin((1 + ε)t).
representing a periodic motion with two dominant frequencies and a small perturbation of low frequency. We take the fitting frequency ω = 1. The system is integrated on the interval [0, 100] with the step sizes h = 1/2j , j = 4, . . ., 7. The numerical results are displayed in Fig. 2(c). In the following two experiments, since the exact solutions are not available in a closed form, we compare the errors of the Hamiltonian produced by the six methods. Problem 4. Consider the pendulum with the Hamiltonian (studied in [19]) given by p2 − a cos(q), a > 0. 2 The Hamiltonian system is equivalent to the following second-order equation H(p, q) =
q¨ = −a sin(q). The initial values are q(0) = 0, q˙ (0) = 1. For these initial conditions the constant value of the Hamiltonian is 3 5 H0 = (1/2) − a. Since for |q| < 1√the equation can be expanded as q¨ + aq = a( q3! − q5! + · · ·), a reasonable choice for the fitted frequency is ω = a. We take the parameter a = 0.2. The problem is integrated on the interval [0, 5000] and the step size h = 1/4. Fig. 3 presents the time evolution of numerical Hamiltonian for each methods on the interval [4930, 5000]. From this figure we see that the symmetric methods SSEFRKN2, SSEFRKN3 and SSEFRKN4 are much more precise in conserving the Hamiltonian than the non-symmetric methods SIMOS2, EFRKNS23F and EFRKNS4. Then we integrate the system on the interval [0, 100] with the step h = 1/2j , j = 2, . . ., 5. Fig. 4 presents the efficiency of each method in the preservation of the Hamiltonian H(p, q).
92
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95 Problem 4: Time−evolution of Hamiltonian for SSEFRKN2 0.55
0.5
0.5
0.45
0.45
0.4
0.4
0.35
0.35
Hamiltonian
Hamiltonian
Problem 4: Time−evolution of Hamiltonian for SIMOS2 0.55
0.3 0.25
0.3 0.25
0.2
0.2
0.15
0.15
0.1 0.05 4930
0.1 4940
4950
4960
4970
4980
4990
0.05 4930
5000
4940
4950
4960
Time
(a) 0.5
0.5
0.45
0.45
0.4
0.4
0.35
0.35
0.3 0.25
0.25
0.2
0.2 0.15
0.1
0.1 4960
4970
4980
4990
0.05 4930
5000
4940
4950
4960
Time
(c)
4980
4990
5000
(d) Problem 4: Time−evolution of Hamiltonian for SSEFRKN4
0.55
0.55
0.5
0.5
0.45
0.45
0.4
0.4
0.35
0.35
Hamiltonian
Hamiltonian
4970 Time
Problem 4: Time−evolution of Hamiltonian for EFRKNS4
0.3 0.25
0.3 0.25
0.2
0.2
0.15
0.15
0.1 0.05 4930
5000
0.3
0.15
4950
4990
Problem 4: Time−evolution of Hamiltonian for EFRKNS23F 0.55
Hamiltonian
Hamiltonian
Problem 4: Time−evolution of Hamiltonian for SSEFRKN3
4940
4980
(b)
0.55
0.05 4930
4970 Time
0.1 4940
4950
4960
4970
4980
4990
5000
0.05 4930
4940
4950
Time
4960
4970
4980
4990
5000
Time
(e)
(f)
Fig. 3. Numerical Hamiltonian for the methods (a) SIMOS2, (b) SSEFRKN2, (c) SSEFRKN3, (d) EFRKNS23F, (e) EFRKNS4 and (f) SSEFRKN4 in Problem 4.
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
93
Problem 4: Efficiency curves for Hamiltonian −2 SIMOS2 SSERKN2 EFRKN23F EFRKNS4 SSEFRKN3 SSEFRKN4
log10(Error of Hamiltonian)
−3
−4
−5
−6
−7
−8
−9
0
1000
2000 3000 4000 5000 Number of function evaluations
6000
7000
Fig. 4. Efficiency curve for Hamiltonian in Problem 4.
Problem 5. The two-body problem (studied in Hairer et al. [15]) consists of finding the positions and velocities of two massive bodies that attract each other gravitationally, given their masses, positions, and velocities at some initial time. The first body is located in the origin while the second body is located in the plane with Cartesian coordinates (q1 , q2 ). The velocity of the second body is given by (˙q1 , q˙ 2 ) = (p1 , p2 ). The Hamiltonian is given by H = T + V,
T =
1 (p1 2 + p2 2 ), 2
1 V = − . q1 2 + q 2 2
The Hamiltonian system is equivalent to the following second-order problem q¨ 1 + q¨ 2 +
q1 3/2
= 0,
3/2
= 0.
(q1 + q2 2 ) q2 2
(q1 2 + q2 2 )
For 0 < e < 1, the initial values are given as follows: q1 (0) = 1 − e,
q2 (0) = 0,
q˙ 1 (0) = 0,
q˙ 2 (0) =
1+e . 1−e
(32)
The exact solution of the problem is q1 = cos(E) − e, q2 = 1 − e2 sin(E). where e is the eccentricity of the orbit and the eccentric anomaly E is expressed as an implicit function of the independent variable t by Kepler’s equation t = E − e sin(E). The system has the conservative energy H and the conservative angular √ momentum M = q1 p2 − q2 p1 . The initial energy H0 = − (1/2(1 − e)) and the initial momentum M0 = 1 − e2 . In this experiment, we choose e = 10−2 . Fig. 5 presents the efficiency curves for the Hamiltonian and for the momentum on the interval [0, 100] with the stepsizes h = 1/2i , i = 3, . . ., 6. From Fig. 2, we can see that the symmetric methods SSEFRKN2, SSSEFRKN3 and SSEFRKN4 are more precise than the non-symmetric methods of the same algebraic order when they are applied to Problems 1 to 3. In Figs. 4 and 5(a), we see the same result in energy-preservation for Problems 4 and 5. Among all the methods, the SSEFRKN4 is the most efficient. It is noted that, in Problems 1 and 4, the second-order methods SSEFRKN2 and SSEFRKN3 is even more efficient than the fourth-order method EFRKNS4. Fig. 5(b) shows that, for Problem 5, among all the methods we select, the new second-order method SSEFRKN3 is the most efficient in preserving the momentum.
94
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
Fig. 5. Efficiency curve for (a) Hamiltonian and (b) momentum in Problem 5.
6. Conclusions A new approach to constructing RKN-type methods for oscillatory Hamiltonian systems is presented in this paper. This approach is based on a combination of the symmetry, symplecticity and exponential fitting conditions. Like the existing RKN integrators adapted to oscillators (see [14], for example), the coefficients of the new SSEFRKN methods depend on the product of the dominant frequency ω and the step size h. When the parameter z approaches zero, the SSEFRKN methods reduce to classical RKN methods, and symmetric RKN methods and the symmetry and symplecticity conditions for SSERKN methods reduce to those for RKN methods. The numerical experiments carried out on several practical oscillatory Hamiltonian systems show that the new SSEFRKN methods are more efficient than non-symmetric symplectic integrators of the same order. In every experiment carried out, the method SSEFRKN4 is shown to be the most efficient among the methods we use for comparison. The FSAL property also contributes to the high efficiency of the new methods SSEFRKN3 and SSEFRKN4. Following the approach of this paper, higher order and multidimensional symmetric and symplectic exponentially fitted RKN methods can be constructed (see Section 3.2 of the monograph by Wu et al. [25]). On the other hand, the development of SSEFRKN methods with variable nodes and embedded pairs will be interesting topics for future work. Acknowledgements The authors are grateful to the anonymous referees for their careful reading of the manuscript and for their invaluable comments and suggestions, which largely help to improve this paper. This research was jointly supported by National Natural Science Foundation of China under Grant No. 11171155, the Fundamental Research Fund for the Central Universities under Grant No. Y0201100265. References [1] P. Albrecht, The extension of the theory of A-methods to RK methods., in: K. Strehmel (Ed.), Numerical Treatment of Differential Equations, Proceedings of the 4th Seminar NUMDIFF-4, Tuebner-Texte Zur Mathematik, Tuebner, Leipzig, 1987, pp. 8–18. [2] P. Albrecht, A new theoretical approach to Runge Kutta methods, SIAM Journal on Numerical Analysis 24 (1987) 391–406. [3] G. Vanden Berghe, M. Van Daele, H. Van de Vyver, Exponential fitted Runge–Kutta methods of collocation type: fixed or variable knot points? Journal of Computational and Applied Mathematics 159 (2003) 217–239. [4] G. Vanden Berghe, H. De Meyer, M. Van Daele, T. Van Hecke, Exponentially fitted Runge–Kutta methods, Journal of Computational and Applied Mathematics 125 (2000) 107–115. [5] J.C. Butcher, Numerical Methods for Ordinary Differential Equations, Second Edition, John Wiley & Sons Ltd, Chichester, 2008. [6] M.P. Calvo, J.M. Sanz-Serna, High-order symplectic Runge–Kutta–Nyström methods, SIAM Journal on Scientific Computing 14 (1993) 1237–1252.
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
95
[7] J. Candy, W. Rozmus, A symplectic integration algorithm for separable Hamiltonian functions, Journal of Computational Physics 92 (1991) 230–256. [8] Z. Chen, X. You, W. Shi, Z. Liu, Symmetric and symplectic ERKN methods for oscillatory Hamiltonian systems, Computer Physics Communications 183 (2012) 86–98. [9] C. Burnton, R. Scherer, Gauss–Runge–Kutta–Nyström methods, BIT 38 (1) (1998) 12–21. [10] J.P. Coleman, L.Gr. Ixaru, P-stability and exponential-fitting methods for y = f(x, y), IMA Journal of Numerical Analysis 16 (1996) 179–199. [11] R. Van Dooren, Stabilization of Cowell’s classical finite difference method for numerical integration, Journal of Computational Physics 16 (1974) 186–192. [12] J.M. Franco, M. Palacios, High-order P-stable multistep methods, Journal of Computational and Applied Mathematics 30 (1990) 1–10. [13] J.M. Franco, Exponentially fitted explicit Runge–Kutta–Nyström methods, Journal of Computational and Applied Mathematics 167 (2004) 1–19. [14] J.M. Franco, Exponentially fitted symplectic integrators of RKN type for solving oscillatory problems, Computer Physics Communications 177 (2007) 479–492. [15] E. Hairer, S.P. Nørsett, G. Wanner, Solving Ordinary Differential Equations I, Nonstiff Problems, Second revised edition, Springer-Verlag, Berlin, 1993. [16] L.Gr. Ixaru, G. Vanden Berghe, Exponential Fitting, Kluwer, Boston, Dordrecht, London, 2004. [17] B. Paternoster, Runge–Kutta(–Nyström) methods for ODEs with periodic solutions based on trigonometric polynomials, Applied Numerical Mathematics 28 (1998) 401–412. [18] J.M. Sanz-Serna, Symplectic integrators for Hamiltonian problems: an overview, Acta Numerica 1 (1992) 243–286. [19] J.M. Sanz-Serna, M.P. Calvo, Numerical Hamiltonian Problems, Chapman and Hall, London, 1994. [20] T.E. Simos, An exponentially-fitted Runge–Kutta method for the numerical integration of initial-value problems with periodic or oscillating solutions, Computer Physics Communications 115 (1998) 1–8. [21] T.E. Simos, J. Vigo-Aguiar, Exponentially fitted symplectic integrator, Physical Review E 67 (2003) 1–7. [22] A. Tocino, J. Vigo-Aguiar, Symplectic Conditions for Exponential Fitting Runge–Kutta–Nyström Methods, Mathematical and Computer Modelling 42 (2005) 873–876. [23] H. Van de Vyver, A fourth-order symplectic exponentially fitted integrator, Computer Physics Communications 174 (2006) 255–262. [24] H. Van de Vyver, A symplectic exponentially fitted modified Runge–Kutta–Nyström method for the numerical integration of orbital problems, New Astronomy 10 (2005) 261–269. [25] X. Wu, X. You, B. Wang, Structure-Preserving Algorithms for Oscillatory Differential Equations, Springer, Berlin, Heidelberg, 2013. [26] H. Yang, X. Wu, X. You, Y. Fang, Extended RKN-type methods for numerical integration of perturbed oscillators, Computer Physics Communications 180 (2009) 1777–1794. [27] D. Zhao, Z. Wang, Y. Dai, Importance of the first-order derivative formula in the Obrechkoff methods, Computer Physics Communications 167 (2005) 65–75.
Mathematics and Computers in Simulation 94 (2013) 76–95
Original article
Symmetric and symplectic exponentially fitted Runge–Kutta–Nyström methods for Hamiltonian problems Xiong You a,c,∗ , Bingzhen Chen b a Department of Applied Mathematics, Nanjing Agricultural University, Nanjing 210095, PR China Department of Mathematics, Beijing Jiaotong University Haibin College, Huanghua 061100, PR China c State Key Laboratory for Novel Software Technology at Nanjing University, Nanjing 210093, PR China
b
Received 19 April 2012; received in revised form 11 December 2012; accepted 30 May 2013 Available online 19 June 2013
Abstract The construction of symmetric and symplectic exponentially fitted modified Runge–Kutta–Nyström (SSEFRKN) methods is considered. Based on the symmetry, symplecticity, and exponentially fitted conditions, new explicit modified RKN integrators with FSAL property are obtained. The new integrators integrate exactly differential systems whose solutions can be expressed as linear combinations of functions from the set { exp(± iωt)}, ω > 0, i2 = −1, or equivalently from the set { cos(ωt), sin(ωt)}. The phase properties of the new integrators are examined and their periodicity regions are obtained. Numerical experiments are accompanied to show the high efficiency and competence of the new SSEFRKN methods compared with some highly efficient nonsymmetric symplecti EFRKN methods in the literature. © 2013 IMACS. Published by Elsevier B.V. All rights reserved. MSC: 65L05; 65L06; 65M20; 65N40 Keywords: Runge–Kutta–Nyström method; Symmetry; Symplecticity; Exponential fitting; Hamiltonian system
1. Introduction In this paper we focus on the initial-value problem (IVP) of a system of second-order ODEs in the form y¨ = f (t, y),
y(t0 ) = y0 ,
y˙ (t0 ) = y˙ 0 ,
(1)
whose solution exhibits an oscillatory character, where y ∈ Rd , f : R × Rd → Rd is a smooth function. Problems of this type arise naturally in a variety of applied sciences such as molecular dynamics, orbital mechanics and electronics, chemistry, biology and engineering. High accuracy of integration is often required in these areas. Until now there have been broadly two categories of approaches to numerical integration of the IVP (1): indirect and direct. On one hand,
∗ Corresponding author at: Department of Applied Mathematics, Nanjing Agricultural University, Nanjing 210095, PR China. Tel.: +86 13605184263. E-mail addresses: [email protected] (X. You), [email protected] (B. Chen).
0378-4754/$36.00 © 2013 IMACS. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.matcom.2013.05.010
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
77
by introducing a new variable w to represent the first derivative y˙ , the IVP (1) is turned into a partitioned system of first-order equations y˙ = w,
˙ = f (t, y), w
y(t0 ) = y0 ,
w(t0 ) = y˙ 0 ,
(2)
and the problem can be solved by the general-purpose Runge–Kutta (RK) methods or partitioned Runge–Kutta (PRK) methods (see [3,4,6,10,13,18,19]). On the other hand, Runge–Kutta–Nyström (RKN) methods, initiated by Nyström in 1925, are designed for attacking the second-order problem (1) directly (see [9]). Some variations of the standard RKN methods are given by Butcher [5], Simos [21], Franco [14] and so forth. Recently, Yang et al. [26] proposed a new family of extend Runge–Kutta–Kutta–Nyström (ERKN) methods which have been shown to be very effective and efficient when applied to perturbed oscillators. On many occasions, the problem under consideration takes the form of a Hamiltonian system p˙ = −
∂ U(t, q), ∂q
q˙ = M −1 p
(3)
with the Hamiltonian H(t, p, q) =
1 T −1 p M p + U(t, q), 2
(4)
where M is a symmetric positive definite constant matrix. This system is equivalent to the second-order equation (1) with ∂ f (t, q) = −M −1 ∂q U(t, q). Note that if the system (3) is autonomous, its Hamiltonian is invariant along any solution. Symplecticity is one of the characteristic properties of Hamiltonian systems and it has been widely acknowledged that the so-called symplectic integrators have the advantage of preserving qualitative properties such as the symplecticity of the exact flow over long time when they are applied to Hamiltonian systems [15,19,23]. Among pioneers in the investigation of symplectic integrators are Vogelaere, Ruth and Feng (see Hairer et al. [15]). From then on, identification and analysis of symplectic methods for solving IVPs (1) have received more and more attention. As is known, the symplectic Runge–Kutta–Nyström methods share the pronounced property of being zerodissipative, which is an important requirement for solving oscillatory and Hamiltonian problems. Simos and Vigo-Aguiar [21] considered symplectic methods of RKN type which are adapted to certain types of oscillators. In particular, they have constructed a symplectic modified RKN method that integrates exactly the linear differential equation y¨ = −ω2y , and whose coefficients are frequency dependent. This method has algebraic order two, and it performs much more efficiently than the classical symplectic RKN second-order method introduced by Calvo and Sanz-Serna [6] when it is applied to oscillatory problems (see the numerical results presented in [21]). On the other hand, exponentially fitted methods, which intend to integrate exactly differential equations whose solutions can be expressed as linear combinations of functions { exp(± iωt)}, ω > 0, i2 = −1, or equivalently from the set { cos(ωt), sin(ωt)}, share excellent behavior when applied to oscillatory problems. The construction of exponentially fitted RK(N) methods is originally due to Paternoster [17] and a detailed exposition of exponentially fitted methods with an extensive bibliography on this subject can be found in Ixaru and Vanden Berghe [16]. Recently, Van de Vyver [24] and Franco [14] considered symplectic exponentially fitted modified RKN method. Unfortunately they did not take into account the symmetry conditions. As is pointed out in Chapters V and XI of [15], symmetric methods show a better long time behavior than non symmetric ones when applied to reversible differential systems, as it is the case for conservative mechanical systems. Very recently, Chen et al. [8] derived the symmetry and symplecticity conditions for ERKN (extended RKN) methods solving oscillatory Hamiltonian systems and constructed two examples of symmetric and symplectic ERKN (SSERKN) methods of orders two and four respectively. Motivated by this work, we derive a family of symmetric and symplectic exponentially fitted modified RKN methods in this paper. The rest of this paper is organized as follows: in Section 2 we present the symmetry conditions, symplecticity conditions, exponential fitting conditions and algebraic order (up to third-order) conditions for modified RKN methods. In Section 3 we derive exponentially fitted symmetric and symplectic integrators with two, three and four stages per step. We will also reveal relationship between the symmetry and symplecticity. In Section 4 we analyze the phase and stability properties of the new methods, obtaining generalized periodicity regions for the classical second-order linear test model. In Section 5 we carry out five numerical experiments to show the efficiency of the new methods compared with the symplectic integrators given in [6,14,21,24]. Section 6 is devoted to conclusive comments.
78
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
2. Conditions for symmetry, symplecticity, exponential fitting and consistency of modified RKN methods In this paper, we consider the method (5) whose coefficients are z-dependent, like exponentially fitted (EF)-type or other adapted-type methods (see, for example, [3,4,7,20,21]). We consider a class of RKN-type methods studied by Butcher [5], Simos and Vigo-Aguiar [21], and Van de Vyver [24]. Definition 2.1. An s-stage modified RKN method for the second-order IVP (1) is a scheme of the form Yi = y0 + ci (z)γi (z)h˙y0 + h2 y1 =
s
aij (z)f (t0 + cj (z)h, Yj ), j=1 s g1 (z)y0 + g2 (z)h˙y0 + h2 bi (z)f (t0 + ci (z)h, Yi ), i=1 s
y˙ 1 = g3 (z)˙y0 + h
i = 1, . . . , s, (5)
bi (z)f (t0 + ci (z)h, Yi ),
i=1
where the nodes ci (z), i = 1, . . ., s and the coefficients g1 (z), g2 (z), g3 (z) and γi (z), bi (z), bi (z), aij (z), i, j = 1, . . ., s are assumed to be even functions of z = ihω with ω > 0 the principal frequency of the problem. Scheme (5) can be expressed by the following Butcher tableau c(z)
e g1 (z)
γ(z) g2 (z) g3 (z)
c1 (z) .. A(z) . ¯bT (z) = cs (z) bT (z)
1 .. .
1 g1 (z)
γ1 (z) .. . γs (z) g2 (z) g3 (z)
a11 (z) .. . as1 (z) ¯b1 (z) b1 (z)
... .. . ... ... ...
a1s (z) .. . ass (z) ¯bs (z) bs (z)
If the matrix A(z) is strictly lower triangular, i.e., aij (z) = 0, 1 ≤ i < j ≤ s, the method (5) is explicit, otherwise it is implicit. If g1 (z) = 1, cs (z)γ s (z) = g2 (z) and bj (z) = asi (z), j = 1, . . . , s, the method is said to have the FSAL-property (the last evaluation at any step is the same as the first evaluation at the next step). A s-stage modified RKN method having the FSAL-property actually requires only s − 1 function evaluations at each step. In convention, we assume that as z → 0, the coefficients g1 (z), g2 (z), g3 (z), γ i (z), i = 1, . . ., s all approach 1, so that the limit method of Scheme (5) coincides with a classical s-stage RKN formula. The objective of this section is to find out when the modified RKN method (5) is symmetric, symplectic and exponentially fitted and to construct practical methods with these favorable properties. 2.1. Symmetry conditions The key to understanding symmetry is the concept of the adjoint method. Definition 2.2. The adjoint method h * of a one-step method h is the inverse map of the original method with reversed time step −h, i.e., h * : = −h −1 . A method for which h * = h is said to be symmetric. (See Hairer [15].) Theorem 2.1. The modified RKN method (5) is symmetric if its coefficients satisfy the following conditions ci (z) + cs+1−i (z) = 1, g1 (z) = g3 (z) = 1, ci γi (z) + cs+1−i (z)γs+1−i (z) = g2 (z), bi (z) = bs+1−i (z), bi (z) + bs+1−i (z) = g2 (z)bs+1−i (z), aij (z) − ci (z)γi (z)bj (z) = as+1−i,s+1−j (z) − bs+1−j (z), for all i, j = 1, . . ., s.
(6)
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
79
Proof. Following the idea of Hairer et al. [15], we interchange h ↔ − h, z = iωh ↔ − z = − iωh, t0 ↔ t0 + h, y0 ↔ y1 and y˙ 0 ↔ y˙ 1 , and obtain the following equations which define the adjoint method of (5)
Yi = y1 − ci (z)γi (z)h˙y1 + h2 y0 =
s
aij (z)f (t0 + (1 − cj (z))h, Yj ), j=1 s g1 (z)y1 − g2 (z)h˙y1 + h2 bi (z)f (t0 + (1 − ci (z))h, Yi ), i=1 s
y˙ 0 = g3 (z)˙y1 − h
i = 1, . . . , s, (7)
bi (z)f (t0 + (1 − ci (z))h, Yi ).
i=1
From Eqs. (7), we have bi (z) 1 y˙ 0 + h f (t0 + (1 − ci (z))h, Yi ) g3 (z) g (z) i=1 3 s 1 bs+1−i (z) = y˙ 0 + h f (t0 + (1 − cs+1−i (z))h, Ys+1−i ), g3 (z) g3 (z) s
y˙ 1 =
(8)
i=1
bi (z) 1 g2 (z) y0 + h˙y1 − h2 f (t0 + (1 − ci (z))h, Yi ) g1 (z) g1 (z) g (z) i=1 1 s 1 bs+1−i (z) g2 (z) = y0 + h˙y1 − h2 f (t0 + (1 − cs+1−i (z))h, Ys+1−i ) g1 (z) g1 (z) g1 (z) i=1 s 1 bs+1−i (z) g2 (z) 1 h y˙ 0 + h f (t0 + (1 − cs+1−i (z))h, Ys+1−i ) y0 + = g1 (z) g1 (z) g3 (z) g3 (z) s
y1 =
i=1
s 1 bs+1−j (z)f (t0 + (1 − cs+1−i (z))h, Ys+1−i ) −h2 g (z) j=1 1 1 g2 (z) = y0 + h˙y0 g1 (z) g1 (z)g3 (z) s 1 g2 (z) +h2 ( bs+1−i (z) − bs+1−i (z))f (t0 + (1 − cs+1−i (z))h, Ys+1−i ) g1 (z) g3 (z)
(9)
i=1
and
Ys+1−i = y1 − cs+1−i (z)γs+1−i (z)h˙y1 + h2
s
as+1−i,s+1−j (z)f (t0 + (1 − cj (z))h, Yj ),
j=1
=
1 1 g2 (z) y0 + ( − cs+1−i (z)γs+1−i (z))hy0 g1 (z) g3 (z) g1 (z) s 1 g2 (z) 1 +h2 ( bs+1−j (z) ( − cs+1−i (z)γs+1−i (z))bs+1−j (z) − g3 (z) g1 (z) g1 (z) j=1
+as+1−i,s+1−j (z)) · f (t0 + (1 − cs+1−j (z))h, Ys+1−j ).
(10)
80
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
Denoting Ys+1−i = Y˜ i , we know that the scheme defined by Eqs. (8)–(10) coincides with Scheme (5) if the following conditions are satisfied 1 − cs+1−i (z) = ci (z), 1 = 1, g1 (z) 1 g2 (z) − cs+1−i (z)γs+1−i (z) = ci (z)γi (z), g3 (z) g1 (z) 1 g2 (z) 1 bs+1−j (z) − cs+1−i (z)γs+1−i (z) bs+1−j (z) − g3 (z) g1 (z) g1 (z) +as+1−i,s+1−j (z) = aij (z), 1 = g1 (z), g1 (z) g2 (z) = g2 (z), g1 (z)g3 (z) 1 g2 (z) ( bs+1−i (z) − bs+1−i (z)) = bi (z), g1 (z) g3 (z) 1 = g3 (z), g3 (z) bs+1−i (z) = bi (z). g3 (z)
(11)
It is easy to verify that Eqs. (11) are equivalent to Eqs. (6). Corollary 2.1. As z → 0, the conditions (6) reduce to the symmetry conditions for traditional RKN methods with constant coefficients ci + cs+1−i = 1, bi = bs+1−i , bi + bs+1−i = bi ,
(12)
aij − ci bj = as+1−i,s+1−j − bs+1−j , for i, j = 1, . . ., s . 2.2. Symplecticity conditions Now we turn to the symplecticity conditions for Scheme (5). Recall that a one-step method h : (q0 , p0 )T → (q1 , p1 )T is symplectic if for every smooth Hamiltonian function H and for every step size h, the numerical flow preserves the differential 2-form dp ∧ dq = di=1 dpi ∧ dqi , i.e., dp1 ∧ dq1 = dp0 ∧ dq0 . In analog, for the second-order system (1) as a Hamiltonian system, the RKN scheme (5) h : is said to be symplectic if dy1 ∧ d˙y1 = dy0 ∧ d˙y0 . From Tocino and Vigo-Aguiar [22], we have the following theorem.
(y0 , y˙ 0 )T → (y1 , y˙ 1 )T (13)
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
81
Theorem 2.2. The modified RKN method (5) is symplectic if its coefficients satisfy the following conditions g1 (z)g3 (z) = 1, g3 (z)bi (z) + (g1 (z)ci (z)γi (z) − g2 (z))bi (z) = 0,
i = 1, . . . , s,
bi (z)(bj (z) − g1 (z)aij (z)) = bj (z)(bi (z) − g1 (z)aji (z)),
(14)
1 ≤ j < i ≤ s.
When g1 (z) = g3 (z) = 1, Eqs. (14) become bi (z) = (g2 (z) − ci (z)γi (z))bi (z),
i = 1, . . . , s, (15)
bi (z)(bj (z) − aij (z)) = bj (z)(bi (z) − aji (z)),
1 ≤ j < i ≤ s.
Corollary 2.2. As z → 0, the conditions (14) or (15) reduce to the symplectic conditions for traditional RKN methods with constant coefficients bi = (1 − ci )bi ,
i = 1, . . . , s, (16)
bi aij − bj aji = bi bj (ci − cj ),
1 ≤ j < i ≤ s.
We refer to Eqs. (6) and (15) as SS conditions. 2.3. Relationship between symmetry, symplecticity and FSAL-property of explicit methods
Theorem 2.3. If an s-stage explicit modified RKN method (5) is symmetric and the sixth equation in (6) is also true for i = j, then the method is symplectic. Proof. The first condition in (14) is satisfied immediately. We only need to check Eqs. (15) for i > j. The first equation in (15) follows from the third and sixth equations in (14) with j = i. On the other hand, with the third and fourth equations in (14) and the first equation in (15), for i > j, the sixth equation in (6) becomes aij (z) = ci (z)γi (z)bj (z) − bs+1−j (z) = ci (z)γi (z)bs+1−j (z) − cj (z)γj (z)bs+1−j (z) = (ci (z)γi (z) − cj (z)γj (z))bs+1−j (z) = (cs+1−j (z)γs+1−j (z) − cs+1−i (z)γs+1−i (z))bj (z). Multiplying both sides of the above equation by bi (z), we obtain bi (z)aij (z)
= bi (z)cs+1−j (z)γs+1−j (z)bj (z) − bj (z)cs+1−i (z)γs+1−i (z)bi (z) = bi (z)bj (z) − bj (z)bi (z).
This is exactly the second equation in (15) for explicit methods with i > j. According to Theorem 2.3, Eqs. (6) for 1 ≤ j ≤ i ≤ s become the SS conditions for the explicit modified RKN method (5). The following theorem shows the connection between symmetry and the FSAL-property.
82
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
Theorem 2.4. Suppose that an s-stage explicit modified RKN method (5) is symmetric and sixth equation in (6) is also true for j = i, i = 1, . . ., s. Suppose c1 (z) = 0 and bi (z) = / 0, i = 1, . . ., s. Then the method has the FSAL-property. Proof. By Theorem 2.3, Eqs. (6) for 1 ≤ j ≤ i ≤ s imply Eqs. (15). From the third equation in (6) and the first equation in (15) we have bs (z) = c1 (z)γ1 (z)bs (z) = 0. Then in the case of i = s, the second equation in (15) gives asj (z) = bj (z), j = 1, . . . , s − 1. By imposing i = j = 1 in the sixth equation in (6), we obtain c1 (z)γ 1 (z) = 0 which implies that cs (z)γ s (z) = g2 (z) by the third equation in (6). Hence, the method satisfies the FSAL-property. 2.4. Exponential fitting conditions Following Albrecht’s lines in [1,2], each stage of the explicit scheme (5) can be viewed as a linear multistep method on a non-equidistant grid. Exponential fitting of the modified RKN method (5) can be presented in terms of a particular type of linear operators on C2 ([t0 , T]), the space of second-order continuously differentiable functions on the interval [t0 , T]. With each stage we associate a linear operator as follows: • for the internal stages,
Li [A(z), h]y(x) = y(x + ci h) − y(x) − ci (z)γi h˙y(x) − h2
s
aij (z)¨y(x + cj (z)h),
i = 1, . . . , s;
j=1
• for the final stages,
L[b(z), h]y(x) = y(x + h) − g1 (z)y(x) − hg2 (z)˙y(x) − h2 L[b(z), h]y(x) = y˙ (x + h) − g3 (z)˙y(x) − h
s
s
bi (z)¨y(x + ci (z)h),
i=1
bi (z)¨y(x + ci (z)h).
i=1
Definition 2.3. The modified RKN method (5) is said to be exponentially fitted if the space span{ exp(± iωt)} is a subspace of the kernels kerL and kerLi , i = 1, . . . , s. By definition, requiring the modified RKN method (5) to be exponentially fitted is equivalent to requiring the operators L and Li , i = 1, . . . , s to vanish for the functions from the set { exp(± iωt)} leading to the following equations exp(±ci (z)z) = 1 ± ci (z)γi (z)z + z2 exp(±z) = g1 (z) ± g2 (z)z + z2 exp(±z) = g3 (z) ± z
s i=1
s
s
aij (z) exp(±cj (z)z),
i = 1, . . . , s,
j=1
bi (z) exp(±ci (z)z),
i=1
bi (z) exp(±ci (z)z),
z = iωh.
(17)
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
83
Using the definitions cosh(z) = (exp(z) + exp(− z))/2 and sinh(z) = (exp(z) − exp(− z))/2, the equations (17) imply that s
aij (z) cosh(cj (z)z) =
cosh(ci (z)z) − 1 , z2
aij (z) sinh(cj (z)z) =
sinh(ci (z)z) − ci (z)γi (z)z , z2
bi (z) cosh(ci (z)z) =
cosh(z) − g1 (z) , z2
j=1
s j=1 s i=1
s i=1
s
(18)
sinh(z) − zg2 (z) bi (z) sinh(ci (z)z) = , z2 bi (z) sinh(ci (z)z) =
cosh(z) − g3 (z) , z
bi (z) cosh(ci (z)z) =
sinh(z) . z
i=1
s
i = 1, . . . , s,
i=1
Eqs. (18) are referred to as the EF conditions. 2.5. Algebraic order conditions In this section, we will present algebraic order conditions for the modified RKN method (5). Like the case of a classical RKN method, for an modified RKN method, the local truncation errors in the approximations of the solution and its derivative can be expressed as ⎛ ⎞ kj p−1 e1 = y(x0 + h) − y1 = hj+1 ⎝ di (j+1) F (j) (τji )(y0 )⎠ + O(hp+1 ), i=1
j=1
e˙ 1 = y˙ (x0 + h) − y˙ 1 =
p
kj (j) j h ( d˙ i F (j) (τji )(y0 )) + O(hp+1 ),
j=1
i=1
where F(j) (τ ji )(y0 ) denotes the elementary differential associated to the Nyström tree τ ji of order j, kj is the number of (j) Nyström trees of order j, and the terms di (j+1) and d˙ i depend on the coefficients of the EFRKN method. A method (5) has order p if, for every sufficiently smooth IVP (1) and for every small step size h, the local truncation errors of the numerical solutions satisfy e1 = y(t0 + h) − y1 = O(hp+1 ), e˙ 1 = y˙ (t0 + h) − y˙ 1 = O(hp+1 ), or equivalently, di (j+1) = 0, (j) d˙ i = 0,
i = 1, . . . , kj ,
i = 1, . . . , kj ,
j = 1, . . . , p − 1,
j = 1, . . . , p.
Following the idea of Franco [14], we assume the following expansions bi (z) = bi
(0)
(2)
(4)
+ bi h2 + bi h4 + · · ·,
γi (z) = 1 + γi (2) h2 + γi (4) h4 + · · ·, g1 (z) = 1 + g1 (2) h2 + g1 (4) h4 + · · ·, g3 (z) = 1 + g3 (2) h2 + g3 (4) h4 + · · ·.
bi (z) = bi (0) + bi (2) h2 + bi (4) h4 + · · ·,
aij (z) = aij (0) + aij (2) h2 + aij (4) h4 + · · ·, g2 (z) = 1 + g2 (2) h2 + g2 (4) h4 + · · ·,
84
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
In the sequel, we only consider methods with constant nodes. Then the order conditions up to order three for the modified RKN method (5) can be listed as follows: Order 1 requires:
(1) d˙ 1 :=
bi (0) − 1 = 0.
(19)
i
Order 2 requires in addition:
d1 (2)
1 (2) = 0, d˙ 2 := g3 (2) = 0, 2 i (0) 1 := bi − = 0, d2 (2) := g1 (2) = 0. 2
(2) d˙ 1 :=
bi (0) ci −
(20)
i
Order 3 requires in addition:
(3) d˙ 1 (3) d˙ 3
1 (0) 2 1 1 (3) := bi ci − bi (0) aik (0) − = 0, = 0, d˙ 2 := 2 3 6 i i k (0) 1 := bi (2) = 0, d1 (3) := bi ci − = 0, 6 i
(21)
i
d2 (3) := g2 (2) = 0.
From Theorem 2.1 in [14], we know that the modified RKN method (5) has algebraic order at least 2.
3. Construction of symmetric and symplectic explicit EFRKN methods In this section we proceed to construct practical explicit EFRKN methods with constant nodes under the symmetry, symplecticity, exponential fitting and algebraic order conditions obtained in the previous section. We refer to these methods as SSEFRKN methods.
3.1. The case s = 2 For s = 2, the SS conditions (6) for 1 ≤ j ≤ i ≤ 2 become c1 + c2 = 1,
g1 (z) = g3 (z) = 1,
b1 (z) + b2 (z) = g2 (z)b1 (z),
c1 γ1 (z) + c2 γ2 (z) = g2 (z),
b1 (z) = c2 γ2 (z)b1 (z),
a21 (z) = c2 γ2 (z)b1 (z) − b2 (z),
b1 (z) = b2 (z),
b2 (z) = c1 γ1 (z)b2 (z),
0 = a21 (z) + c1 γ1 (z)b2 (z) − b1 (z).
(22)
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
85
On the other hand, for s = 2 the EF conditions (18) give cosh(c1 z) − 1 sinh(c1 z) − c1 γ1 (z)z = 0, = 0, 2 z z2 cosh(z) − 1 a21 (z) cosh(c1 z) = , z2 sinh(z) − c2 γ2 (z)z a21 (z) sinh(c1 z) = , z2 cosh(z) − g1 (z) b1 (z) cosh(c1 z) + b2 (z) cosh(c2 z) = , z2 sinh(z) − g2 (z) b1 (z) sinh(c1 z) + b2 (z) sinh(c2 z) = , z2 cosh(z) − g3 (z) b1 (z) cosh(c1 z) + b2 (z) cosh(c2 z) = , z sinh(z) b1 (z) cosh(c1 z) + b2 (z) cosh(c2 z) = . z
(23)
Solving (22) and (23), we obtain g1 (z) = g3 (z) = 1,
c1 = 0,
c2 = 1,
γ2 (z) =
sinh(z) , z
g2 (z) = γ2 (z),
cosh(z) − 1 , b1 (z) = a21 (z), b2 (z) = 0, z2 sinh(z) b1 (z) = , b2 (z) = b1 (z). z(cosh(z) + 1)
a21 (z) =
(24)
The two-stage method of order two is denoted by SSEFRKN2. It is exactly the symplectic exponentially fitted integrator EFRKN2 presented in [24]. However, unlike [24], we do not artificially assume c2 = 1 to find the coefficients.
3.2. The case s = 3 In a similar way we can obtain the coefficients of a three-stage SSERKN method of order two c1 = 0,
c2 =
1 , 2
g1 (z) = g3 (z) = 1,
c3 = 1, g2 (z) = γ3 (z) = γ2 (z) =
2 sinh(z/2) , z
cosh(z/2) − 1 , z2 2(cosh(z/2) − 1) a31 (z) = , a32 (z) = a31 (z), z2 b1 (z) = a31 (z), b2 (z) = a32 (z), b3 (z) = 0, a21 (z) =
b1 (z) =
cosh(z/2) − 1 , z sinh(z/2)
b2 (z) = 2b1 (z),
We denote this method as SSEFRKN3.
b3 (z) = b1 (z).
(25)
86
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
3.3. The case s = 4 Our main interest is in the case s = 4. The SS conditions (6) for 1 ≤ j ≤ i ≤ 4 become
g1 (z) = g3 (z) = 1,
c1 + c4 = 1,
c1 γ1 (z) + c4 γ4 (z) = g2 (z), b1 (z) = b4 (z),
c2 + c3 = 1,
c2 γ2 (z) + c3 γ3 (z) = g2 (z),
b2 (z) = b3 (z),
b1 (z) + b4 (z) = g2 (z)b1 (z), b1 = c4 γ4 (z)b4 (z),
b2 (z) + b3 (z) = g2 (z)b2 (z),
b2 = c3 γ3 (z)b3 (z),
b3 = c2 γ2 (z)b2 (z),
(26) b4 = c1 γ1 (z)b1 (z),
a21 (z) = c2 γ2 (z)b1 (z) − b4 (z),
a31 (z) = c3 γ3 (z)b1 (z) − b4 (z),
a32 (z) = c3 γ3 (z)b2 (z) − b3 (z),
a41 (z) = c4 γ4 (z)b1 (z) − b4 (z),
a42 (z) = c4 γ4 (z)b2 (z) − b3 (z),
a43 (z) = c3 γ4 (z)b3 (z) − b2 (z).
On the other hand, the EF conditions (18) for s = 4 can be become
cosh(c1 z) − 1 = 0, z2
sinh(c1 z) − c1 γ1 (z)z = 0, z2
a21 (z) cosh(c1 z) =
cosh(c2 z) − 1 , z2
a21 (z) sinh(c1 z) =
sinh(c2 z) − c2 γ2 (z)z , z2
a31 (z) cosh(c1 z) + a32 (z) cosh(c2 z) = a31 (z) sinh(c1 z) + a32 (z) sinh(c2 z) =
cosh(c3 (z)z) − 1 , z2 sinh(c3 z) − c3 γ3 (z)z , z2
a41 (z) cosh(c1 z) + a42 (z) cosh(c2 z) + a43 (z) cosh(c3 z) = a41 (z) sinh(c1 z) + a42 (z) sinh(c2 z) + a43 (z) sinh(c3 z) =
cosh(c4 z) − 1 , z2
sinh(c4 z) − c4 γ4 (z)z . z2
cosh(z) − 1 , z2 sinh(z) − g2 (z)z b1 (z) sinh(c1 z) + b2 (z) sinh(c2 z) + b3 (z) sinh(c3 z) = , z2 sinh(z) b1 (z) cosh(c1 z) + b2 (z) cosh(c2 z) + b3 cosh(c3 z) + b4 (z) cosh(c4 z) = , z cosh(z) − 1 b1 (z) sinh(c1 z) + b2 sinh(c2 z) + b3 (z) sinh(c3 z) + b4 (z) sinh(c4 z) = . z b1 (z) cosh(c1 z) + b2 (z) cosh(c2 z) + b3 (z) cosh(c3 z) =
(27)
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
87
We solve the system of Eqs. (26) and (27) for the coefficients with c2 as a parameter as follows c1 = 0,
c4 = 1,
c3 = 1 − c2 ,
γ2 (z) =
sinh(c2 z) , c2 z
a21 (z) =
g1 (z) = g3 (z) = 1,
cosh(c2 z) − 1 , z2
z sinh(z/2) − b1 (z)z cosh( ) a21 (z) 2 b1 (z) = , , b2 (z) = c2 γ2 (z) z cosh(z(c3 − c2 )/2) b3 (z) = c2 γ2 (z)b3 (z), N0 b1 (z)z − N1 , b1 (z)b2 (z)z sinh(c2 z) − b2 (z) cosh(c2 z) N1 − g2 (z)b2 (z) cosh(c2 z) γ3 (z) = c3 b1 (z) a31 (z) = c3 γ3 (z)b1 (z),
g2 (z) =
b3 (z) = b2 (z),
b4 (z) = b1 (z),
γ4 (z) = g2 (z),
(28)
a32 (z) = g2 (z)b2 (z) − 2c2 γ2 (z)b2 (z), b1 (z) = g2 (z)b1 (z), a41 (z) = b1 (z),
b2 (z) = c3 γ3 (z)b2 (z),
a42 (z) = b2 (z),
a43 (z) = b3 (z).
where N0 =
sinh(c3 z) + 2b3 (z) sinh(c2 z), z2
N1 =
cosh(c3 z) − 1 + 2b3 (z) cosh(c2 z). z2 √ 3
√ 3
Substituting ci , bi (z), i = 1, . . ., 4 into the first equation in (21) are satisfied, we obtain c2 = 32 + 64 + 23 ≈ 1.351207191959658. With this value of c2 , it is easy to check that (19), (20) and the other equations in (21) are satisfied. Thus the method defined by (28) is of order three. Since a symmetric method must be of an even order, then the method is actually of order 4 (see Chapter II of [15]). We denote this method as SSEFRKN4. It should be noted that the new methods SSEFRKN2 and SSEFRKN3, SSEFRKN4 derived in this section all have the FSAL-property.
4. Phase-lags and periodicity regions of the new methods Now we set out to analyze the phase and stability properties of our new methods. Stability means that the numerical solutions remain bounded as we move further away from the starting point. For classical RKN methods, the phase and stability properties are checked using the second-order linear test model y¨ (t) = −μ2 y(t),
μ > 0.
(29)
Recall that the new symmetric and symplectic exponentially fitted modified RKN methods (SEFRKN) derived in the previous section are dependent on the complex number λ = iω, where ω > 0 is an estimate of the principal frequency. Applying an s-stage SEFRKN method (5) to the test model (29) yields
y1 h˙y1
= M(η , ν ) 2
2
y0 h˙y0
,
(30)
88
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
where
M(η2 , ν2 ) =
g1 (z) − η2 b(z)T N(z)−1 e
g2 (z) − η2 b(z)T N(z)−1 c(z)γ(z)
, −η2 b(z)T N(z)−1 e g3 (z) − η2 b(z)T N(z)−1 c(z)γ(z) ν = ωh, η = μh, N(z) = I + η2 A(z), b(z)T = (b1 (z), . . . , bs (z)), bT (z) = (b1 (z), . . . , bs (z)), ⎛ ⎞ 0 0 ... 0 ⎛ ⎞ ⎛ ⎞ c (z)γ (z) 1 1 1 ⎜ a21 (z) 0 ... 0⎟ ⎜ ⎟ ⎜ ⎟ ⎜ .. ⎟ .. ⎟ , A(z) = ⎜ ⎟. e = ⎝ . ⎠ , c(z)γ(z) = ⎜ ⎜ ⎝ ⎠ . . . .. .. .. ⎟ ⎝ .. ⎠ . . 1 cs (z)γs (z) as1 (z) . . . as,s−1 (z) 0
The phase and stability behavior of the numerical solution depends on the eigenvalues or the spectrum of the stability matrix M = M(η2 , ν2 ). Eliminating y˙ 0 and y˙ 1 from (30) and the equation that is obtained from (30) by replacing the subscript 0 by 1 gives the difference equation y2 − tr(M(η2 , ν2 ))y1 + det(M(η2 , ν2 ))y0 = 0. Accordingly, the characteristic equation is given by ξ 2 − tr(M)ξ + det(M) = 0,
(31)
where tr(M) and det(M) are the trace and the determinant of M, respectively. Definition 4.1. For the characteristic Eq. (31) obtained by applying an exponentially fitted method to Eq. (29), set the ratio r = ω/μ. The two quantities tr(M) √ φ(η, r) = η − arccos , 2 det(M) √ d(η, r) = 1 − det(M) are called the dispersion error and the dissipation error of the numerical method under consideration, respectively. The method is said to be dispersive of order p and dissipative of order q, if φ(η, r) = O(ηp+1 ) and d(η, r) = O(ηq+1 ), respectively. The method is called zero-dissipative if d(η, r) = 0. For the method SSEFRKN2 method, we have r2 − 1 3 η + O(η5 ), 24 which shows that SSEFRKN2 is dispersive of order 2. For the method SSEFRKN3, we have φ(η, r) =
r2 − 1 3 η + O(η5 ), 96 which shows that SSEFRKN3 is dispersive of order 2. For the method SSEFRKN4 method, we have φ(η, r) =
12c2 4 − 18c2 3 + 6c2 − 1 2 3 r η + O(η5 ), 24 which shows that SSEFRKN4 is dispersive of order 2. It is noted that for the methods SSEFRKN2, SERRKN3 and SSEFRKN4, det(M) ≡ 1, since they are all symplectic. Therefore these methods are all zero-dissipative. The stability properties of an EFRKN method are characterized by the spectral radius ρ(M). Having in mind that the matrix M depends on the variables η and ν, the periodicity and stability intervals typical in the classical RKN methods now become two-dimensional regions in the (H, ν)-plane. So, for the new methods we use an adaptation of the terminology introduced by Coleman and Ixaru [10]: φ(η, r) =
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95 SSEFRKN2
30
SSEFRKN4
8
SSEFRKN3
30
89
25
20
20
5
15
15
4
6
ν
ν
ν
7 25
3 10
10
5
5
2
0 0
5
10
15 η
20
25
30
0 0
1 5
(a)
10
15 η
(b)
20
25
30
0 0
1
2
3
4
η
5
6
7
8
(c)
Fig. 1. Periodicity regions for the methods (a) SSEFRKN2, (b) SSEFRKN3 and (c) SSEFRKN4.
(i) (ii) (iii) (iv)
Rs = {η > 0, ν > 0|ρ(M) < 1} is the region of stability; Rp = {η > 0, ν > 0|ρ(M) = 1andtr(M)2 < 4det(M)} is the region of periodicity; If Rs = (0, ∞) × (0, ∞) except possibly for a discrete set of curves, the method is A-stable; If Rp = (0, ∞) × (0, ∞) except possibly for a discrete set of curves, the method is P-stable.
Note again that the symplecticity implies that det(M) ≡ 1, and hence symplectic methods have nonempty regions of periodicity. The regions of periodicity of the methods SSEFRKN2, SSEFRKN3 and SSEFRKN4 are depicted in Fig. 1.
5. Numerical experiments To test the numerical performance of the new modified RKN methods, we carry out experiments on five model problems to illustrate the effectiveness and efficiency of the new methods SSEFRKN3 and SSEFRKN4. The codes used in the comparisons are: • SIMOS2: The symplectic modified two-stage RKN method of order two proposed by Simos and Vigo-Aguiar [21]. • SSEFRKN2: The symplectic exponentially fitted modified two-stage RKN method of order two with FSAL-property proposed by Van de Vyver [24], that is, the method given by (24). • EFRKNS23F: The symplectic exponentially fitted modified three-stage RKN method of order two with FSALproperty proposed by Franco [14]. • EFRKNS4: The symplectic exponentially fitted modified three-stage RKN method of order four proposed by Franco [14]. • SSEFRKN3: The symmetric and symplectic exponentially fitted modified three-stage RKN method of order two with FSAL-property proposed in this paper. • SSEFRKN4: The symmetric and symplectic exponentially fitted modified four-stage RKN method of order four with FSAL-property proposed in this paper. The criterion used in the numerical comparisons is efficiency in terms of the decimal logarithm of the maximum global error (log10 (MGE)) versus the computational effort measured in the number of function evaluations required by each method. Problem 1. Consider the “almost periodic” orbit problem (studied by Franco and Palacios [12]) z¨ + z = ε exp(iψx),
z(0) = 1,
z˙ (0) = i,
90
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
Fig. 2. Efficiency curves for solution in (a) Problem 1, (b) Problem 2 and (c) Problem 3.
where ε = 0.001 and ψ = 0.01. The analytical solution z(x) = u(x) + iv(x) is given by ε 1 − ε − ψ2 cos(x) + cos(ψx), 2 1−ψ 1 − ψ2 1 − εψ − ψ2 ε v(x) = sin(x) + sin(ψx). 2 1−ψ 1 − ψ2 u(x) =
The problem is a nonhomogeneous Hamiltonian system H(p, q) =
1 1 (p1 2 + p2 2 ) + (q1 2 + q2 2 ) − ε(cos(ψt)q1 + sin(ψt)q2 ), 2 2
where q = (q1 , q2 )T = (u, v)T , p = (p1 , p2 )T = (˙q1 , q˙ 2 )T . The initial values are q(0) = (1, 0)T , q˙ (0) = (0, 1)T . We integrate the problem on the interval [0, 1000] with fitting frequency ω = 1. The stepsizes are taken as h = 1/2j , j = 1, . . ., 5 for the methods SIMOS2, SSEFRKN2, SSEFRKN3, and 3h/2 for the other three methods EFRKNS23F, EFRKNS4 and SSERKN4. We will follow this convention in the following experiments. The numerical results are displayed in Fig. 2(a).
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
91
Problem 2. Consider the nonhomogeneous Duffing equation (studied in [21]) q + q¨ + q3 = 0.002 cos(1.01t),
q(0) = 0.200426728069666,
q˙ (0) = 0.
The Hamiltonian of this system is 1 2 1 2 1 4 p + q + q − 0.002 cos(1.01t)q. 2 2 4 An “exact” solution was first given by Van Dooren [11], but a more accurate solution is given by Zhao et al. [27] as H(p, q) =
q(t) =
4
A2i+1 cos[(2i + 1)1.01t],
i=0
with (A1 , A3 , A5 , A7 , A9 ) = (0.20017947753661502, 2.46946143255559 × 10−4 , 3.0401498519692437 × 10−7 , 3.743490701609247 × 10−10 , 4.609682949622697 × 10−13 ). The equation is integrated on the interval [0, 1000] with ω = 1, h = 1/2j , j = 0, . . ., 4. The numerical results are displayed in Fig. 2(b). Problem 3. Consider the perturbed orbit problem (studied in [14]) q¨ 1 = − q¨ 1 = −
q1 (q1 2 + q2 q1
2 )3/2
(q1 + q2
2 )3/2
2
− −
(2ε + ε2 )q1 5/2
(q1 2 + q2 2 ) (2ε + ε2 )q1 (q1 + q2 2
2 )5/2
,
q1 (0) = 1,
q1 (0) = 0,
,
q2 (0) = 0,
q2 (0) = 1 + ε.
with ε = 0.001 . The Hamiltonian of this problem is given by H(p, q) =
1 1 2ε + ε2 (p1 2 + p2 2 ) − − . 1/2 3/2 2 (q1 2 + q2 2 ) 3(q1 2 + q2 2 )
The analytic solution of the problem is given by q1 (t) = cos((1 + ε)t),
q2 (t) = sin((1 + ε)t).
representing a periodic motion with two dominant frequencies and a small perturbation of low frequency. We take the fitting frequency ω = 1. The system is integrated on the interval [0, 100] with the step sizes h = 1/2j , j = 4, . . ., 7. The numerical results are displayed in Fig. 2(c). In the following two experiments, since the exact solutions are not available in a closed form, we compare the errors of the Hamiltonian produced by the six methods. Problem 4. Consider the pendulum with the Hamiltonian (studied in [19]) given by p2 − a cos(q), a > 0. 2 The Hamiltonian system is equivalent to the following second-order equation H(p, q) =
q¨ = −a sin(q). The initial values are q(0) = 0, q˙ (0) = 1. For these initial conditions the constant value of the Hamiltonian is 3 5 H0 = (1/2) − a. Since for |q| < 1√the equation can be expanded as q¨ + aq = a( q3! − q5! + · · ·), a reasonable choice for the fitted frequency is ω = a. We take the parameter a = 0.2. The problem is integrated on the interval [0, 5000] and the step size h = 1/4. Fig. 3 presents the time evolution of numerical Hamiltonian for each methods on the interval [4930, 5000]. From this figure we see that the symmetric methods SSEFRKN2, SSEFRKN3 and SSEFRKN4 are much more precise in conserving the Hamiltonian than the non-symmetric methods SIMOS2, EFRKNS23F and EFRKNS4. Then we integrate the system on the interval [0, 100] with the step h = 1/2j , j = 2, . . ., 5. Fig. 4 presents the efficiency of each method in the preservation of the Hamiltonian H(p, q).
92
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95 Problem 4: Time−evolution of Hamiltonian for SSEFRKN2 0.55
0.5
0.5
0.45
0.45
0.4
0.4
0.35
0.35
Hamiltonian
Hamiltonian
Problem 4: Time−evolution of Hamiltonian for SIMOS2 0.55
0.3 0.25
0.3 0.25
0.2
0.2
0.15
0.15
0.1 0.05 4930
0.1 4940
4950
4960
4970
4980
4990
0.05 4930
5000
4940
4950
4960
Time
(a) 0.5
0.5
0.45
0.45
0.4
0.4
0.35
0.35
0.3 0.25
0.25
0.2
0.2 0.15
0.1
0.1 4960
4970
4980
4990
0.05 4930
5000
4940
4950
4960
Time
(c)
4980
4990
5000
(d) Problem 4: Time−evolution of Hamiltonian for SSEFRKN4
0.55
0.55
0.5
0.5
0.45
0.45
0.4
0.4
0.35
0.35
Hamiltonian
Hamiltonian
4970 Time
Problem 4: Time−evolution of Hamiltonian for EFRKNS4
0.3 0.25
0.3 0.25
0.2
0.2
0.15
0.15
0.1 0.05 4930
5000
0.3
0.15
4950
4990
Problem 4: Time−evolution of Hamiltonian for EFRKNS23F 0.55
Hamiltonian
Hamiltonian
Problem 4: Time−evolution of Hamiltonian for SSEFRKN3
4940
4980
(b)
0.55
0.05 4930
4970 Time
0.1 4940
4950
4960
4970
4980
4990
5000
0.05 4930
4940
4950
Time
4960
4970
4980
4990
5000
Time
(e)
(f)
Fig. 3. Numerical Hamiltonian for the methods (a) SIMOS2, (b) SSEFRKN2, (c) SSEFRKN3, (d) EFRKNS23F, (e) EFRKNS4 and (f) SSEFRKN4 in Problem 4.
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
93
Problem 4: Efficiency curves for Hamiltonian −2 SIMOS2 SSERKN2 EFRKN23F EFRKNS4 SSEFRKN3 SSEFRKN4
log10(Error of Hamiltonian)
−3
−4
−5
−6
−7
−8
−9
0
1000
2000 3000 4000 5000 Number of function evaluations
6000
7000
Fig. 4. Efficiency curve for Hamiltonian in Problem 4.
Problem 5. The two-body problem (studied in Hairer et al. [15]) consists of finding the positions and velocities of two massive bodies that attract each other gravitationally, given their masses, positions, and velocities at some initial time. The first body is located in the origin while the second body is located in the plane with Cartesian coordinates (q1 , q2 ). The velocity of the second body is given by (˙q1 , q˙ 2 ) = (p1 , p2 ). The Hamiltonian is given by H = T + V,
T =
1 (p1 2 + p2 2 ), 2
1 V = − . q1 2 + q 2 2
The Hamiltonian system is equivalent to the following second-order problem q¨ 1 + q¨ 2 +
q1 3/2
= 0,
3/2
= 0.
(q1 + q2 2 ) q2 2
(q1 2 + q2 2 )
For 0 < e < 1, the initial values are given as follows: q1 (0) = 1 − e,
q2 (0) = 0,
q˙ 1 (0) = 0,
q˙ 2 (0) =
1+e . 1−e
(32)
The exact solution of the problem is q1 = cos(E) − e, q2 = 1 − e2 sin(E). where e is the eccentricity of the orbit and the eccentric anomaly E is expressed as an implicit function of the independent variable t by Kepler’s equation t = E − e sin(E). The system has the conservative energy H and the conservative angular √ momentum M = q1 p2 − q2 p1 . The initial energy H0 = − (1/2(1 − e)) and the initial momentum M0 = 1 − e2 . In this experiment, we choose e = 10−2 . Fig. 5 presents the efficiency curves for the Hamiltonian and for the momentum on the interval [0, 100] with the stepsizes h = 1/2i , i = 3, . . ., 6. From Fig. 2, we can see that the symmetric methods SSEFRKN2, SSSEFRKN3 and SSEFRKN4 are more precise than the non-symmetric methods of the same algebraic order when they are applied to Problems 1 to 3. In Figs. 4 and 5(a), we see the same result in energy-preservation for Problems 4 and 5. Among all the methods, the SSEFRKN4 is the most efficient. It is noted that, in Problems 1 and 4, the second-order methods SSEFRKN2 and SSEFRKN3 is even more efficient than the fourth-order method EFRKNS4. Fig. 5(b) shows that, for Problem 5, among all the methods we select, the new second-order method SSEFRKN3 is the most efficient in preserving the momentum.
94
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
Fig. 5. Efficiency curve for (a) Hamiltonian and (b) momentum in Problem 5.
6. Conclusions A new approach to constructing RKN-type methods for oscillatory Hamiltonian systems is presented in this paper. This approach is based on a combination of the symmetry, symplecticity and exponential fitting conditions. Like the existing RKN integrators adapted to oscillators (see [14], for example), the coefficients of the new SSEFRKN methods depend on the product of the dominant frequency ω and the step size h. When the parameter z approaches zero, the SSEFRKN methods reduce to classical RKN methods, and symmetric RKN methods and the symmetry and symplecticity conditions for SSERKN methods reduce to those for RKN methods. The numerical experiments carried out on several practical oscillatory Hamiltonian systems show that the new SSEFRKN methods are more efficient than non-symmetric symplectic integrators of the same order. In every experiment carried out, the method SSEFRKN4 is shown to be the most efficient among the methods we use for comparison. The FSAL property also contributes to the high efficiency of the new methods SSEFRKN3 and SSEFRKN4. Following the approach of this paper, higher order and multidimensional symmetric and symplectic exponentially fitted RKN methods can be constructed (see Section 3.2 of the monograph by Wu et al. [25]). On the other hand, the development of SSEFRKN methods with variable nodes and embedded pairs will be interesting topics for future work. Acknowledgements The authors are grateful to the anonymous referees for their careful reading of the manuscript and for their invaluable comments and suggestions, which largely help to improve this paper. This research was jointly supported by National Natural Science Foundation of China under Grant No. 11171155, the Fundamental Research Fund for the Central Universities under Grant No. Y0201100265. References [1] P. Albrecht, The extension of the theory of A-methods to RK methods., in: K. Strehmel (Ed.), Numerical Treatment of Differential Equations, Proceedings of the 4th Seminar NUMDIFF-4, Tuebner-Texte Zur Mathematik, Tuebner, Leipzig, 1987, pp. 8–18. [2] P. Albrecht, A new theoretical approach to Runge Kutta methods, SIAM Journal on Numerical Analysis 24 (1987) 391–406. [3] G. Vanden Berghe, M. Van Daele, H. Van de Vyver, Exponential fitted Runge–Kutta methods of collocation type: fixed or variable knot points? Journal of Computational and Applied Mathematics 159 (2003) 217–239. [4] G. Vanden Berghe, H. De Meyer, M. Van Daele, T. Van Hecke, Exponentially fitted Runge–Kutta methods, Journal of Computational and Applied Mathematics 125 (2000) 107–115. [5] J.C. Butcher, Numerical Methods for Ordinary Differential Equations, Second Edition, John Wiley & Sons Ltd, Chichester, 2008. [6] M.P. Calvo, J.M. Sanz-Serna, High-order symplectic Runge–Kutta–Nyström methods, SIAM Journal on Scientific Computing 14 (1993) 1237–1252.
X. You, B. Chen / Mathematics and Computers in Simulation 94 (2013) 76–95
95
[7] J. Candy, W. Rozmus, A symplectic integration algorithm for separable Hamiltonian functions, Journal of Computational Physics 92 (1991) 230–256. [8] Z. Chen, X. You, W. Shi, Z. Liu, Symmetric and symplectic ERKN methods for oscillatory Hamiltonian systems, Computer Physics Communications 183 (2012) 86–98. [9] C. Burnton, R. Scherer, Gauss–Runge–Kutta–Nyström methods, BIT 38 (1) (1998) 12–21. [10] J.P. Coleman, L.Gr. Ixaru, P-stability and exponential-fitting methods for y = f(x, y), IMA Journal of Numerical Analysis 16 (1996) 179–199. [11] R. Van Dooren, Stabilization of Cowell’s classical finite difference method for numerical integration, Journal of Computational Physics 16 (1974) 186–192. [12] J.M. Franco, M. Palacios, High-order P-stable multistep methods, Journal of Computational and Applied Mathematics 30 (1990) 1–10. [13] J.M. Franco, Exponentially fitted explicit Runge–Kutta–Nyström methods, Journal of Computational and Applied Mathematics 167 (2004) 1–19. [14] J.M. Franco, Exponentially fitted symplectic integrators of RKN type for solving oscillatory problems, Computer Physics Communications 177 (2007) 479–492. [15] E. Hairer, S.P. Nørsett, G. Wanner, Solving Ordinary Differential Equations I, Nonstiff Problems, Second revised edition, Springer-Verlag, Berlin, 1993. [16] L.Gr. Ixaru, G. Vanden Berghe, Exponential Fitting, Kluwer, Boston, Dordrecht, London, 2004. [17] B. Paternoster, Runge–Kutta(–Nyström) methods for ODEs with periodic solutions based on trigonometric polynomials, Applied Numerical Mathematics 28 (1998) 401–412. [18] J.M. Sanz-Serna, Symplectic integrators for Hamiltonian problems: an overview, Acta Numerica 1 (1992) 243–286. [19] J.M. Sanz-Serna, M.P. Calvo, Numerical Hamiltonian Problems, Chapman and Hall, London, 1994. [20] T.E. Simos, An exponentially-fitted Runge–Kutta method for the numerical integration of initial-value problems with periodic or oscillating solutions, Computer Physics Communications 115 (1998) 1–8. [21] T.E. Simos, J. Vigo-Aguiar, Exponentially fitted symplectic integrator, Physical Review E 67 (2003) 1–7. [22] A. Tocino, J. Vigo-Aguiar, Symplectic Conditions for Exponential Fitting Runge–Kutta–Nyström Methods, Mathematical and Computer Modelling 42 (2005) 873–876. [23] H. Van de Vyver, A fourth-order symplectic exponentially fitted integrator, Computer Physics Communications 174 (2006) 255–262. [24] H. Van de Vyver, A symplectic exponentially fitted modified Runge–Kutta–Nyström method for the numerical integration of orbital problems, New Astronomy 10 (2005) 261–269. [25] X. Wu, X. You, B. Wang, Structure-Preserving Algorithms for Oscillatory Differential Equations, Springer, Berlin, Heidelberg, 2013. [26] H. Yang, X. Wu, X. You, Y. Fang, Extended RKN-type methods for numerical integration of perturbed oscillators, Computer Physics Communications 180 (2009) 1777–1794. [27] D. Zhao, Z. Wang, Y. Dai, Importance of the first-order derivative formula in the Obrechkoff methods, Computer Physics Communications 167 (2005) 65–75.