- Email: [email protected]
Symplectic Partitioned Runge–Kutta methods with minimal phase-lag
Computer Physics Communications 181 (2010) 1251–1254
Contents lists available at ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
Symplectic Partitioned Runge–Kutta methods with minimal phase-lag Th. Monovasilis a , Z. Kalogiratou b , T.E. Simos c,∗,1 a b c
Department of International Trade, Technological Educational Institution of Western Macedonia at Kastoria, Kastoria, P.O. Box 30, 52100 Greece Department of Informatics and Computer Technology, Technological Educational Institution of Western Macedonia at Kastoria, Kastoria, P.O. Box. 30, 52100 Greece Department of Computer Science and Technology, Faculty of Science and Technology, University of Peloponnessos, Greece
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 21 May 2009 Received in revised form 19 March 2010 Accepted 24 March 2010 Available online 30 March 2010
Symplectic Partitioned Runge–Kutta (SPRK) methods with minimal phase-lag are derived. Specifically two new symplectic methods are constructed of second and third order with fifth phase-lag order. The methods are tested on the numerical integration of Hamiltonian problems and the Schrödinger equation. © 2010 Elsevier B.V. All rights reserved.
Keywords: Symplectic Partitioned Runge–Kutta methods Hamiltonian problems Schrödinger equation Phase-lag
1. Introduction In the last decades a lot of research has been performed in the area of numerical integration of Hamiltonian systems. Hamiltonian systems appear in many areas of mechanics, physics, chemistry, and elsewhere. They are first-order ODEs that can be expressed as
p = −
∂H ( p , q), ∂q
q =
∂H ( p , q) ∂p
where p , q ∈ Rd and H ( p , q, x) is a sufficiently smooth real function called the Hamiltonian of the system [1]. Symplecticity is a characteristic property of Hamiltonian systems and many authors developed and applied symplectic schemes for the numerical integration of such systems. Additionally the solution of Hamiltonian systems often has an oscillatory behavior and have been solved in the literature with methods which take into account the nature of the problem. There are two categories of such methods, with constant coefficients and with coefficients depending on the frequency. In this work we shall consider methods of the first category. These methods can be applied to any oscillatory problem even if the frequency is not initially known. Many authors developed such methods with specific properties as P-stability or minimal phase-lag. Van de Vyrer [2] constructed a symplectic Runge–Kutta– Nystrom method with minimal phase-lag. In this work Symplectic
* 1
Corresponding author. E-mail address: [email protected] (T.E. Simos). Active member of the European Academy of Sciences and Arts.
0010-4655/$ – see front matter doi:10.1016/j.cpc.2010.03.013
©
2010 Elsevier B.V. All rights reserved.
Partitioned Runge–Kutta (SPRK) methods with minimal phase-lag are derived. Specifically two new symplectic methods are constructed of second and third order with fifth phase-lag order. In Sections 2 and 3 we give some basic theory concerning SPRK methods and the phase-lag property. The construction of the new methods is presented in Section 4. Numerical illustrations are given in Section 5. 2. Symplectic Partitioned Runge–Kutta methods In this work we shall consider systems with separable Hamiltonian of the form
H ( p , q, x) = T ( p ) + V (q, x),
T ( p) =
1 2
p T p.
In this case the Hamiltonian system is
p = f (q, x),
q = p ,
where
f (q, x) = −
∂V (q, x). ∂q
A numerical method is symplectic if
dpn+1 ∧ dqn+1 = dpn ∧ dqn . Separable Hamiltonian systems can be solved with SPRK methods. The advantage of using SPRK is that there exist explicit SPRK methods, while SRK methods cannot be explicit.
1252
Th. Monovasilis et al. / Computer Physics Communications 181 (2010) 1251–1254
A PRK scheme is specified by two tableaux
C 1 a11
.. .
... .. . ... ...
.. .
C s as1 c1
a1s
D 1 A 11
.. .
.. .
ass
D s A s1
.. .
cs
d1
... .. . ... ...
Ci =
M=
,
hp 0
As (v 2)
B s (v 2)
C s (v 2)
D s (v 2)
,
v = wh.
.. .
The eigenvalues of the M are called amplification factors of the method and are the roots of the characteristic equation
A ss
ξ 2 − tr M v 2 ξ + det M v 2 = 0.
ds
ai j ,
Di =
and
s
φ( v ) = v − arccos
Ai j .
j =1
The first tableau is used for the integration of p components and the second tableau is used for the integration of the q components as follows:
P i = pn + h
s
Q i = qn + h
s
Ai j P j ,
p n +1 = p n + h
s
c i f ( Q i , x + D i h),
c i A i j + d j a ji − c i d j = 0,
Ci =
i
c j,
j =1
Di =
i −1
d j,
i = 1, 2, . . . , s .
j =1
The fact that the tableaux of the SPRK method are constant along columns imply a favorable implementation of the method using only two d-dimensional vectors since P i +1 and Q i +1 can be overwritten on P i and Q i . The computation proceeds in the following form:
P 0 = pn , Q 1 = qn , for i = 1, . . . , s, P i = P i −1 + hc i f ( Q i , xn + D i h), Q i +1 = Q i + hdi P i ,
= P s,
n +1
= Q s +1 .
c 2 d1 + c 3 (d1 + d2 ) = 1/2 (second order), c 2 d21 + c 3 (d1 + d2 )2 = 1/3 (third order), d3 + d2 (c 1 + c 2 )2 + d1 c 12 = 1/3.
i , j = 1, 2, . . . , s .
As it was mentioned previously there exist explicit PRK methods. Assume the following explicit form ai j = 0 for i < j and A i j = 0 for i j. Then due to the symplecticness requirement the coefficients ai j and A i j (and consequently C i and D i ) are fully determined in terms of the coefficients c i and di .
n +1
α ( v ) = 1 − det M v 2 .
d 1 + d 2 + d 3 = 1,
The above PRK method is symplectic when the following equalities are satisfied
Ai j = d j ,
and the dissipation (amplification error) is
c 1 + c 2 + c 3 = 1 (first order),
di P i .
i =1
ai j = c j ,
,
In this work we consider methods with three stages and orders two and three. The symplecticness conditions are
i =1
qn+1 = qn + h
2 det( M ( v 2 ))
4. Construction of the specific methods
i = 1, 2, . . . , s ,
j =1 s
tr( M ( v 2 ))
For a symplectic PRK method the determinant of the amplification matrix is zero, so the methods we construct here are zero dissipative.
ai j f ( Q j , x + D j h),
j =1
q
q0
The phase-lag (dispersion) of the method is
j =1
p
= Mn
hpn
A 1s
where s
qn
For the second order method we set c 1 = c 3 and solve the order conditions leaving c 1 and d3 as free parameters. We substitute into φ( v ) and take the Taylor expansion. The coefficient of v 3 is
pl3 =
24(d3 − 1)c 13 − 12(d23 + d3 − 2)c 12 + 4(3d23 − 2)c 1 + 1 48c 1 − 24 5
and the coefficient of v is
pl5 = t1/t2, where
t1 = 320(d3 − 1)2 c 16 + 320(d3 − 2)(d3 − 1)2 c 15
+ 80 d43 − 6d33 + 17d23 + 17d23 − 16d3 + 4 c 14 − 80 2d43 − 2d33 + 8d23 − 5d3 − 1 c 13 + 8 10d43 + 25d23 − 5d3 − 14 c 12 − 8 5d23 − 4 c 1 − 3,
t2 = 640(1 − 2c 1 )2 . Since We solve the equations pl3 = 0, pl5 = 0 and obtain the following coefficients
c 1 = 0.5974665433347971,
c 2 = −0.1949330866695942,
c 3 = 0.5974665433347971,
d1 = −0.18554773759667065,
3. Dispersion and dissipation in PRK methods
d2 = 0.961876740574417,
d3 = 0.2236709955392291.
Phase-lag analysis of numerical methods for second order equations is based on the scalar test equation q = − w 2 q, where w is a real constant [3]. For the numerical solution of this equation we can write
We refer to this method as New23. For the third order method the coefficient of v 3 in the Taylor expansion of φ( v ) is zero. In order to achieve fifth phase lag order the additional condition that the coefficient of v 5 is zero is added
Th. Monovasilis et al. / Computer Physics Communications 181 (2010) 1251–1254
1253
Table 1 Two-body with integration interval [0, 1000].
h = 1/2 h = 1/4 h = 1/8
Symp23
Ruth
Symp33
New23
New33
3.4 × 10−4 2.3 × 10−5 1.5 × 10−6
6.2 × 10−5 4.4 × 10−7 3.0 × 10−9
4.1 × 10−5 5.1 × 10−7 7.3 × 10−9
9.6 × 10−5 2.1 × 10−6 7.0 × 10−8
7.5 × 10−6 5.5 × 10−8 5.7 × 10−10
Table 2 Stiefel–Bettis with integration interval is [0, 1000].
h = 1/4 h = 1/8 h = 1/16
Symp23
Ruth
Symp33
New23
New33
0.65 0.14 3.0 × 10−2
1.4 × 10−3 6.6 × 10−5 1.6 × 10−6
9.4 × 10−4 6.5 × 10−5 4.8 × 10−6
1.0 × 10−3 2.4 × 10−4 5.7 × 10−5
2.0 × 10−4 2.6 × 10−5 3.2 × 10−6
to the system of the order conditions and the following coefficients are obtained.
c 1 = 0.2603116924199056,
c 2 = 1.0941427983167422,
c 3 = −0.35445449073664803,
d1 = 0.6308476929866689,
d2 = −0.0941427983167424,
d3 = 0.4632951053300734.
We refer to this method as New33. 5. Numerical results The new methods (New23, New33) are tested on the numerical integration of Hamiltonian problems and the Schrödinger equation. For comparison we use the optimal second and third order SPRK methods Symp23, Symp33 [4], the well-known third and fourth order SPRK methods of Ruth Ruth [5]. 5.1. The two-body problem The following system of equations is known as the two-body problem and is a standard symplectic test case:
p 1 = −
q1
(q21
p 2 = −
q2
(q21
,
+ q22 )3
E0 E 10 E l5 E 20 E 30
Symp23
Y4
Ruth
Symp33
New23
New33
90 – – – –
4 – – – –
1 214 661 1508 5113
0 83 285 692 2301
0 8 16 18 121
0 7 15 18 120
Table 4 Absolute error (×10−6 ) of the eigenvalues of the harmonic oscillator with step size h = 0.05.
E 40 E 50 E 60 E 70 E 80 E 100
Symp23
Y4
Ruth
Symp33
New23
New33
– – – – – –
– – – – – –
691 1348 2335 3724 5592 –
328 639 1105 1756 2627 5155
6 13 26 47 78 190
6 13 26 47 78 189
5.3. The Schrödinger equation We shall use our new methods for the computation of the eigenvalues of the one-dimensional time-independent Schrödinger equation. The Schrödinger equation may be written in the form
q1 = p 1 ,
,
+ q22 )3
Table 3 Absolute error (×10−6 ) of the eigenvalues of the harmonic oscillator with step size h = 0.1.
q2 = p 2
1
− y + V (x) y = E y , 2
x ∈ [a, b], y (a) = y (b) = 0
with initial conditions p 1 (0) = 0, q1 (0) = 1, p 2 (0) = 1, q2 (0) = 0. The exact solution is q1 (x) = cos x, q2 (x) = sin x. In Table 1 we give the maximum Euclidean norm of the global error. The integration interval is [0, 1000].
where E is the energy eigenvalue, V (x) the potential, and y (x) the wave function. The problems used are the harmonic oscillator and the doubly anharmonic oscillator.
5.2. Stiefel–Bettis
5.3.1. The harmonic oscillator The potential of the one-dimensional harmonic oscillator is
A classical test case is the following orbital problem studied by Stiefel and Bettis.
V (x) =
y 1 = − y 1 + 0.001 cos(x),
y 1 (0) = 1,
y 1 (0) = 0,
y 2 = − y 2 + 0.001 sin(x),
y 2 (0) = 0,
y 2 (0) = 0.9995,
1 2
we consider k = 1. The integration interval is [− R , R ]. The exact eigenvalues are given by
with exact solution
y 1 (x) = cos(x) + 0.0005x sin(x), y 2 (x) = sin(x) − 0.0005x cos(x). In Table 2 we give the maximum error with integration interval [0, 1000] and several step size.
kx2
En = n +
1 2
,
n = 0, 1 , 2 , . . . .
In Tables 3 and 4 we give the absolute error of several eigenvalues up to E 100 computed with step sizes h = 0.1 and h = 0.05. The integration interval ranges from R = 5 to R = 20. Both new methods give very accurate eigenvalues.
1254
Th. Monovasilis et al. / Computer Physics Communications 181 (2010) 1251–1254
Table 5 Absolute error (×10−6 ) of the eigenvalues of the doubly anharmonic oscillator with step size h = 0.05. Exact
Symp23
Y4
Ruth
Symp33
New23
New33
0.807447 5.553677 12.534335 21.118364 31.030942 42.104446 54.222484 67.29805 81.262879 96.061534 111.647831
93 2427 – – – – – – – – –
2 317 3787 – – – – – – – –
1 17 85 293 796 1837 3757 7025 – – –
0 4 1 41 193 559 1300 2621 – – –
0 4 13 32 63 112 186 297 470 730 1126
0 4 12 30 59 107 178 286 455 711 1102
V (x) =
2
In this paper two new symplectic partitioned Runge–Kutta methods with minimum phase lag are presented. Both methods have three stages and algebraic orders two and three with phase lag order five. The new methods have been tested against other three and four stages well-known methods. The evidence of the numerical results is that the new methods are far more efficient. Since all methods tested (apart from the fourth order method of Yoshida) have three stages the computational cost is identical. Acknowledgements The authors wish to thank the anonymous referee for his/her careful reading of the manuscript and his/her fruitful comments and suggestions.
5.3.2. Doubly anharmonic oscillator The potential is
1
6. Conclusion
2
4
6
x + λ1 x + λ2 x
we take λ1 = λ2 = 1/2. The integration interval is [− R , R ]. In Table 5 we give the computed eigenvalues up to E 20 with step h = 0.05 and integration interval [−3, 3]. The performance of all methods considered is similar with this the harmonic oscillator.
References [1] J.M. Sanz-Serna, M.P. Calvo, Numerical Hamiltonian Problem, Chapman and Hall, London, 1994. [2] H. Van de Vyrer, A symplectic Runge–Kutta–Nystrom method with minimal phase lag, Physics Letters A 367 (2007) 16–24. [3] P.J. Van Der Houwen, B.P. Sommeijer, Explicit Runge–Kutta(–Nystrom) methods with reduced phase errors for computing oscillating solutions, SIAM Journal of Numerical Analysis 24 (1987) 595–617. [4] R. McLachlan, P. Atela, The accuracy of symplectic integrators, Nonlinearity 5 (1992) 541–562. [5] R.D. Ruth, A canonical integration technique, IEEE Transactions on Nuclear Science, NS 30 (1983) 2669–2671.
Contents lists available at ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
Symplectic Partitioned Runge–Kutta methods with minimal phase-lag Th. Monovasilis a , Z. Kalogiratou b , T.E. Simos c,∗,1 a b c
Department of International Trade, Technological Educational Institution of Western Macedonia at Kastoria, Kastoria, P.O. Box 30, 52100 Greece Department of Informatics and Computer Technology, Technological Educational Institution of Western Macedonia at Kastoria, Kastoria, P.O. Box. 30, 52100 Greece Department of Computer Science and Technology, Faculty of Science and Technology, University of Peloponnessos, Greece
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 21 May 2009 Received in revised form 19 March 2010 Accepted 24 March 2010 Available online 30 March 2010
Symplectic Partitioned Runge–Kutta (SPRK) methods with minimal phase-lag are derived. Specifically two new symplectic methods are constructed of second and third order with fifth phase-lag order. The methods are tested on the numerical integration of Hamiltonian problems and the Schrödinger equation. © 2010 Elsevier B.V. All rights reserved.
Keywords: Symplectic Partitioned Runge–Kutta methods Hamiltonian problems Schrödinger equation Phase-lag
1. Introduction In the last decades a lot of research has been performed in the area of numerical integration of Hamiltonian systems. Hamiltonian systems appear in many areas of mechanics, physics, chemistry, and elsewhere. They are first-order ODEs that can be expressed as
p = −
∂H ( p , q), ∂q
q =
∂H ( p , q) ∂p
where p , q ∈ Rd and H ( p , q, x) is a sufficiently smooth real function called the Hamiltonian of the system [1]. Symplecticity is a characteristic property of Hamiltonian systems and many authors developed and applied symplectic schemes for the numerical integration of such systems. Additionally the solution of Hamiltonian systems often has an oscillatory behavior and have been solved in the literature with methods which take into account the nature of the problem. There are two categories of such methods, with constant coefficients and with coefficients depending on the frequency. In this work we shall consider methods of the first category. These methods can be applied to any oscillatory problem even if the frequency is not initially known. Many authors developed such methods with specific properties as P-stability or minimal phase-lag. Van de Vyrer [2] constructed a symplectic Runge–Kutta– Nystrom method with minimal phase-lag. In this work Symplectic
* 1
Corresponding author. E-mail address: [email protected] (T.E. Simos). Active member of the European Academy of Sciences and Arts.
0010-4655/$ – see front matter doi:10.1016/j.cpc.2010.03.013
©
2010 Elsevier B.V. All rights reserved.
Partitioned Runge–Kutta (SPRK) methods with minimal phase-lag are derived. Specifically two new symplectic methods are constructed of second and third order with fifth phase-lag order. In Sections 2 and 3 we give some basic theory concerning SPRK methods and the phase-lag property. The construction of the new methods is presented in Section 4. Numerical illustrations are given in Section 5. 2. Symplectic Partitioned Runge–Kutta methods In this work we shall consider systems with separable Hamiltonian of the form
H ( p , q, x) = T ( p ) + V (q, x),
T ( p) =
1 2
p T p.
In this case the Hamiltonian system is
p = f (q, x),
q = p ,
where
f (q, x) = −
∂V (q, x). ∂q
A numerical method is symplectic if
dpn+1 ∧ dqn+1 = dpn ∧ dqn . Separable Hamiltonian systems can be solved with SPRK methods. The advantage of using SPRK is that there exist explicit SPRK methods, while SRK methods cannot be explicit.
1252
Th. Monovasilis et al. / Computer Physics Communications 181 (2010) 1251–1254
A PRK scheme is specified by two tableaux
C 1 a11
.. .
... .. . ... ...
.. .
C s as1 c1
a1s
D 1 A 11
.. .
.. .
ass
D s A s1
.. .
cs
d1
... .. . ... ...
Ci =
M=
,
hp 0
As (v 2)
B s (v 2)
C s (v 2)
D s (v 2)
,
v = wh.
.. .
The eigenvalues of the M are called amplification factors of the method and are the roots of the characteristic equation
A ss
ξ 2 − tr M v 2 ξ + det M v 2 = 0.
ds
ai j ,
Di =
and
s
φ( v ) = v − arccos
Ai j .
j =1
The first tableau is used for the integration of p components and the second tableau is used for the integration of the q components as follows:
P i = pn + h
s
Q i = qn + h
s
Ai j P j ,
p n +1 = p n + h
s
c i f ( Q i , x + D i h),
c i A i j + d j a ji − c i d j = 0,
Ci =
i
c j,
j =1
Di =
i −1
d j,
i = 1, 2, . . . , s .
j =1
The fact that the tableaux of the SPRK method are constant along columns imply a favorable implementation of the method using only two d-dimensional vectors since P i +1 and Q i +1 can be overwritten on P i and Q i . The computation proceeds in the following form:
P 0 = pn , Q 1 = qn , for i = 1, . . . , s, P i = P i −1 + hc i f ( Q i , xn + D i h), Q i +1 = Q i + hdi P i ,
= P s,
n +1
= Q s +1 .
c 2 d1 + c 3 (d1 + d2 ) = 1/2 (second order), c 2 d21 + c 3 (d1 + d2 )2 = 1/3 (third order), d3 + d2 (c 1 + c 2 )2 + d1 c 12 = 1/3.
i , j = 1, 2, . . . , s .
As it was mentioned previously there exist explicit PRK methods. Assume the following explicit form ai j = 0 for i < j and A i j = 0 for i j. Then due to the symplecticness requirement the coefficients ai j and A i j (and consequently C i and D i ) are fully determined in terms of the coefficients c i and di .
n +1
α ( v ) = 1 − det M v 2 .
d 1 + d 2 + d 3 = 1,
The above PRK method is symplectic when the following equalities are satisfied
Ai j = d j ,
and the dissipation (amplification error) is
c 1 + c 2 + c 3 = 1 (first order),
di P i .
i =1
ai j = c j ,
,
In this work we consider methods with three stages and orders two and three. The symplecticness conditions are
i =1
qn+1 = qn + h
2 det( M ( v 2 ))
4. Construction of the specific methods
i = 1, 2, . . . , s ,
j =1 s
tr( M ( v 2 ))
For a symplectic PRK method the determinant of the amplification matrix is zero, so the methods we construct here are zero dissipative.
ai j f ( Q j , x + D j h),
j =1
q
q0
The phase-lag (dispersion) of the method is
j =1
p
= Mn
hpn
A 1s
where s
qn
For the second order method we set c 1 = c 3 and solve the order conditions leaving c 1 and d3 as free parameters. We substitute into φ( v ) and take the Taylor expansion. The coefficient of v 3 is
pl3 =
24(d3 − 1)c 13 − 12(d23 + d3 − 2)c 12 + 4(3d23 − 2)c 1 + 1 48c 1 − 24 5
and the coefficient of v is
pl5 = t1/t2, where
t1 = 320(d3 − 1)2 c 16 + 320(d3 − 2)(d3 − 1)2 c 15
+ 80 d43 − 6d33 + 17d23 + 17d23 − 16d3 + 4 c 14 − 80 2d43 − 2d33 + 8d23 − 5d3 − 1 c 13 + 8 10d43 + 25d23 − 5d3 − 14 c 12 − 8 5d23 − 4 c 1 − 3,
t2 = 640(1 − 2c 1 )2 . Since We solve the equations pl3 = 0, pl5 = 0 and obtain the following coefficients
c 1 = 0.5974665433347971,
c 2 = −0.1949330866695942,
c 3 = 0.5974665433347971,
d1 = −0.18554773759667065,
3. Dispersion and dissipation in PRK methods
d2 = 0.961876740574417,
d3 = 0.2236709955392291.
Phase-lag analysis of numerical methods for second order equations is based on the scalar test equation q = − w 2 q, where w is a real constant [3]. For the numerical solution of this equation we can write
We refer to this method as New23. For the third order method the coefficient of v 3 in the Taylor expansion of φ( v ) is zero. In order to achieve fifth phase lag order the additional condition that the coefficient of v 5 is zero is added
Th. Monovasilis et al. / Computer Physics Communications 181 (2010) 1251–1254
1253
Table 1 Two-body with integration interval [0, 1000].
h = 1/2 h = 1/4 h = 1/8
Symp23
Ruth
Symp33
New23
New33
3.4 × 10−4 2.3 × 10−5 1.5 × 10−6
6.2 × 10−5 4.4 × 10−7 3.0 × 10−9
4.1 × 10−5 5.1 × 10−7 7.3 × 10−9
9.6 × 10−5 2.1 × 10−6 7.0 × 10−8
7.5 × 10−6 5.5 × 10−8 5.7 × 10−10
Table 2 Stiefel–Bettis with integration interval is [0, 1000].
h = 1/4 h = 1/8 h = 1/16
Symp23
Ruth
Symp33
New23
New33
0.65 0.14 3.0 × 10−2
1.4 × 10−3 6.6 × 10−5 1.6 × 10−6
9.4 × 10−4 6.5 × 10−5 4.8 × 10−6
1.0 × 10−3 2.4 × 10−4 5.7 × 10−5
2.0 × 10−4 2.6 × 10−5 3.2 × 10−6
to the system of the order conditions and the following coefficients are obtained.
c 1 = 0.2603116924199056,
c 2 = 1.0941427983167422,
c 3 = −0.35445449073664803,
d1 = 0.6308476929866689,
d2 = −0.0941427983167424,
d3 = 0.4632951053300734.
We refer to this method as New33. 5. Numerical results The new methods (New23, New33) are tested on the numerical integration of Hamiltonian problems and the Schrödinger equation. For comparison we use the optimal second and third order SPRK methods Symp23, Symp33 [4], the well-known third and fourth order SPRK methods of Ruth Ruth [5]. 5.1. The two-body problem The following system of equations is known as the two-body problem and is a standard symplectic test case:
p 1 = −
q1
(q21
p 2 = −
q2
(q21
,
+ q22 )3
E0 E 10 E l5 E 20 E 30
Symp23
Y4
Ruth
Symp33
New23
New33
90 – – – –
4 – – – –
1 214 661 1508 5113
0 83 285 692 2301
0 8 16 18 121
0 7 15 18 120
Table 4 Absolute error (×10−6 ) of the eigenvalues of the harmonic oscillator with step size h = 0.05.
E 40 E 50 E 60 E 70 E 80 E 100
Symp23
Y4
Ruth
Symp33
New23
New33
– – – – – –
– – – – – –
691 1348 2335 3724 5592 –
328 639 1105 1756 2627 5155
6 13 26 47 78 190
6 13 26 47 78 189
5.3. The Schrödinger equation We shall use our new methods for the computation of the eigenvalues of the one-dimensional time-independent Schrödinger equation. The Schrödinger equation may be written in the form
q1 = p 1 ,
,
+ q22 )3
Table 3 Absolute error (×10−6 ) of the eigenvalues of the harmonic oscillator with step size h = 0.1.
q2 = p 2
1
− y + V (x) y = E y , 2
x ∈ [a, b], y (a) = y (b) = 0
with initial conditions p 1 (0) = 0, q1 (0) = 1, p 2 (0) = 1, q2 (0) = 0. The exact solution is q1 (x) = cos x, q2 (x) = sin x. In Table 1 we give the maximum Euclidean norm of the global error. The integration interval is [0, 1000].
where E is the energy eigenvalue, V (x) the potential, and y (x) the wave function. The problems used are the harmonic oscillator and the doubly anharmonic oscillator.
5.2. Stiefel–Bettis
5.3.1. The harmonic oscillator The potential of the one-dimensional harmonic oscillator is
A classical test case is the following orbital problem studied by Stiefel and Bettis.
V (x) =
y 1 = − y 1 + 0.001 cos(x),
y 1 (0) = 1,
y 1 (0) = 0,
y 2 = − y 2 + 0.001 sin(x),
y 2 (0) = 0,
y 2 (0) = 0.9995,
1 2
we consider k = 1. The integration interval is [− R , R ]. The exact eigenvalues are given by
with exact solution
y 1 (x) = cos(x) + 0.0005x sin(x), y 2 (x) = sin(x) − 0.0005x cos(x). In Table 2 we give the maximum error with integration interval [0, 1000] and several step size.
kx2
En = n +
1 2
,
n = 0, 1 , 2 , . . . .
In Tables 3 and 4 we give the absolute error of several eigenvalues up to E 100 computed with step sizes h = 0.1 and h = 0.05. The integration interval ranges from R = 5 to R = 20. Both new methods give very accurate eigenvalues.
1254
Th. Monovasilis et al. / Computer Physics Communications 181 (2010) 1251–1254
Table 5 Absolute error (×10−6 ) of the eigenvalues of the doubly anharmonic oscillator with step size h = 0.05. Exact
Symp23
Y4
Ruth
Symp33
New23
New33
0.807447 5.553677 12.534335 21.118364 31.030942 42.104446 54.222484 67.29805 81.262879 96.061534 111.647831
93 2427 – – – – – – – – –
2 317 3787 – – – – – – – –
1 17 85 293 796 1837 3757 7025 – – –
0 4 1 41 193 559 1300 2621 – – –
0 4 13 32 63 112 186 297 470 730 1126
0 4 12 30 59 107 178 286 455 711 1102
V (x) =
2
In this paper two new symplectic partitioned Runge–Kutta methods with minimum phase lag are presented. Both methods have three stages and algebraic orders two and three with phase lag order five. The new methods have been tested against other three and four stages well-known methods. The evidence of the numerical results is that the new methods are far more efficient. Since all methods tested (apart from the fourth order method of Yoshida) have three stages the computational cost is identical. Acknowledgements The authors wish to thank the anonymous referee for his/her careful reading of the manuscript and his/her fruitful comments and suggestions.
5.3.2. Doubly anharmonic oscillator The potential is
1
6. Conclusion
2
4
6
x + λ1 x + λ2 x
we take λ1 = λ2 = 1/2. The integration interval is [− R , R ]. In Table 5 we give the computed eigenvalues up to E 20 with step h = 0.05 and integration interval [−3, 3]. The performance of all methods considered is similar with this the harmonic oscillator.
References [1] J.M. Sanz-Serna, M.P. Calvo, Numerical Hamiltonian Problem, Chapman and Hall, London, 1994. [2] H. Van de Vyrer, A symplectic Runge–Kutta–Nystrom method with minimal phase lag, Physics Letters A 367 (2007) 16–24. [3] P.J. Van Der Houwen, B.P. Sommeijer, Explicit Runge–Kutta(–Nystrom) methods with reduced phase errors for computing oscillating solutions, SIAM Journal of Numerical Analysis 24 (1987) 595–617. [4] R. McLachlan, P. Atela, The accuracy of symplectic integrators, Nonlinearity 5 (1992) 541–562. [5] R.D. Ruth, A canonical integration technique, IEEE Transactions on Nuclear Science, NS 30 (1983) 2669–2671.