- Email: [email protected]
A Runge–Kutta–Nyström pair for the numerical integration of perturbed oscillators
Computer Physics Communications 167 (2005) 129–142 www.elsevier.com/locate/cpc
A Runge–Kutta–Nyström pair for the numerical integration of perturbed oscillators Hans Van de Vyver Received 23 June 2004; accepted 30 December 2004
Abstract New Runge–Kutta–Nyström methods especially designed for the numerical integration of perturbed oscillators are presented in this paper. They are capable of exactly integrating the harmonic or unperturbed oscillator. We construct an embedded 4(3) RKN pair that is based on the FSAL property. The new method is much more efficient than previously derived RKN methods for some subclasses of problems. 2005 Elsevier B.V. All rights reserved. PACS: 02.60.Jh; 02.70.Bf Keywords: Explicit Runge–Kutta–Nyström methods; Embedded pairs; Variable step-size algorithms; Initial value problems; Oscillating solutions
1. Introduction We are concerned with initial value problems (IVP) related to non-stiff second-order ordinary differential equations (ODE) of the form y = µ2 y + f (x, y),
y(x0 ) = y0 ,
y (x0 ) = y0 ,
|µ2 |
µ ∈ R or iR.
(1.1)
In the case that |f (x, y)| the solution exhibits a pronounced exponential (µ ∈ R) or oscillatory (µ ∈ iR) character. Equations having oscillatory solutions are of particular interest, examples occur frequently in celestial mechanics, in quantum mechanical scattering problems and elsewhere. There are many procedures in existence for the construction of Runge–Kutta(–Nyström) (RK(N)) methods for an efficient integration of oscillatory problems. An exhaustive list of references would be near impossible. For example, in the mid-eighties Van der Houwen et al. [1–3] constructed RK(N) methods with a reduced or null phase error. A first example of a phase-fitted RKN method (which means that the phase-lag order is infinity) is due to Paternoster [4]. Simos et al. [5–7] applied some modified explicit RK(N) methods (phase-fitted or with minimal phase-lag) to solve Schrödinger equations. E-mail address: [email protected] (H. Van de Vyver). 0010-4655/$ – see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2004.12.011
130
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
More recently, new approaches for constructing RK(N) methods adapted to perturbed oscillators are presented by González et al. [8], García et al. [9] and Franco [10–12]. Other optimized RK methods for oscillatory problems are designed by Tsitouras and Simos [13] and Anastassi and Simos [14]. For an efficient numerical integration of oscillatory problems is exponential fitting also an interesting option. The study of exponentially fitted RK(N) (EFRK(N)) methods is a relatively new development and we refer to [15–18] as pioneering work in this area. The first implicit EFRK(N) methods are presented by Paternoster [16]. A disadvantage of these methods is that they are fully implicit which means that a high computational effort is required. Very recently, explicit EFRKN methods are presented by Simos [19] and Franco [20]. In particular, in [20] was an embedded pair of EFRKN methods and a variable step-size algorithm developed and discussed. The pair has the RKN pair of orders 4(3) derived by Dormand et al. [21] as underlying classical pair. In this paper we present a technique for the construction of RKN methods that integrate exactly Eq. (1.1) with f (x, y) = 0. We will construct a new embedded pair of RKN methods that is also based on the RKN4(3) pair from [21]. The derivation technique can be seen as an extension of Simos’ RK approach [15] but here we start from another point of view. Numerical examples reveal that the new embedded pair is much more efficient than other ones for some subclasses of (1.1). The paper is organized as follow: in Section 2 we present some basic elements of RKN methods and embedded pairs. In Section 3 we describe the derivation of the new RKN methods and the expressions of the local truncation errors will be discussed. We will explain the difference between the RKN approach developed in this paper and EFRK(N) methods. In Section 4 some numerical experiments show the properties of the developed algorithms. Finally, in Section 5 some conclusions are drawn.
2. Embedded pairs of RKN methods For the numerical solution of the system of second-order ODEs y = f (x, y),
y(x0 ) = y0 ,
y (x0 ) = y0 ,
(2.1)
we consider the explicit Runge–Kutta–Nyström (RKN) method Yi = yn + ci hyn + h2
i−1
aij f (xn + cj h, Yj ), j =1 s b¯i f (xn + ci h, Yi ), yn+1 = yn + hyn + h2 i=1 s yn+1 = yn + h bi f (xn + ci h, Yi ). i=1
i = 1, . . . , s, (2.2)
The RKN method is completely determined by means of its Butcher tableau 0 c2 a21 c3 a31 .. .. . . cs as1 b¯1 b1
a32 .. . . . . as2 . . . b¯2 . . . b2 . . .
¯ b), or equivalent by the quadruple (c, A, b,
(2.3)
as,s−1 b¯s−1 b¯s bs−1 bs
where c, b, b¯ ∈ Rs , A ∈ Rs×s and s denotes the number of stages of the RKN method. An embedded q(p) pair of ¯ b) of order q and another RKN method (c, A, b¯ ∗ , b∗ ) of order RKN methods is based on the RKN method (c, A, b,
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
131
p < q. Embedded pairs of explicit RK(N) methods are widely used in variable step-size algorithms because they provide a cheap error estimation. A local error estimation in the integration point xn+1 = xn + h is determined by the expressions ∗ δn+1 = yn+1 − yn+1
∗ and δn+1 = yn+1 − yn+1 .
(2.4)
EST n+1 (h) = max{δn+1 ∞ , δn+1 ∞ } represents the local error estimation to control the step-size h for RKN processes by the well-known standard formula (see [22]) 1/(p+1) Tol hnew = 0.9hold (2.5) , EST n+1 (hold )
where Tol is the maximum allowable error. If EST n+1 (h) < Tol, then the step is accepted and we adopt the widely used procedure of performing local extrapolation: although, we are actually controlling an estimate of the local error in the lower order solution, it is the higher order solution we accept at each point. If EST n+1 (h) Tol then the step is rejected and we calculate hnew given by (2.5). It should be noted that several variants of local error estimation and step-size control algorithms are in existence (for example, Shampine [23–25]).
3. Explicit RKN methods fitted to perturbed oscillators 3.1. A fourth-order RKN method We consider the four-stage explicit RKN method which is in short-hand notation presented by the following Butcher-tableau: 0 1 1 4 32 7 119 7 10 1000 500
1 a41 a42 b¯1 b¯2 b1 b2
(3.1) a43 b¯3 0 b3 b4
The A-values from the first three rows are taken from the well-known fourth-order RKN method derived by Dormand et al. [21]. We are focused on the unperturbed problem: y = µ2 y,
y(x0 ) = y0 ,
y (x0 ) = y0
(µ ∈ R or µ ∈ iR).
(3.2)
The exact solution is given by: y(x) = χ1 exp(µx) + χ2 exp(−µx),
(3.3)
whereby χ1 , χ2 ∈ C are determined by the initial conditions. We remark that in practice for most of the applications is µ ∈ iR. Our aim is to determine the b- and b∗ -values in the way that the RKN method (3.1) integrates exactly problem (3.2). We solve (3.2) with the RKN method and we obtain for the numerical solution at x1 = x0 + h: 3 3 (0) (1) ψ¯ 2i µ2i y0 + 1 + ψ¯ 2i µ2i hy0 , y1 = 1 + (3.4) i=1
i=1
132
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
hy1
=
4
(0) ψ2i µ2i
y0 + 1 +
i=1
4
(1) ψ2i µ2i
hy0 .
(3.5)
i=1
¯ The ψ-values are given by: (0) ψ¯ 2 =
3
b¯i ,
(0) ψ¯ 4 =
(1) ψ¯ 2
=
(0) ψ¯ 6 = b¯3 a32 a21 ,
b¯i aij ,
i=2 j =1
i=1 3
3 i−1
b¯i ci ,
(2) ψ¯ 4
(3.6)
= b¯3 a32 c2 ,
(1) ψ¯ 6
= 0.
i=2
The ψ -values are given by: (0) ψ2
=
4
(0) ψ4
bi ,
=
i=1 (1) ψ2
=
4
bi ci ,
(1) ψ4
4 i−1
bi aij ,
i=2 j =1 4 i−1
=
(0) ψ6
=
j −1 4 i−1
bi aij aj k ,
(0)
ψ8 = b4 a43 a32 a21 ,
i=3 j =2 k=1
bi aij cj ,
(1) ψ6
= b4 a43 a32 c2 ,
(3.7) (1) ψ8
= 0.
i=3 j =2
i=2
Substituting y0 = y(x), y0 = y (x), y1 = y(x + h) and y1 = y (x + h) in (3.4) and (3.5) (y(x) given in (3.3)), one easily finds the next relations (ν = µh): 3
(0) ψ¯ 2i ν 2i +
i=1
3
3
(1) ψ¯ 2i ν 2i+1 = exp(ν) − 1 − ν,
i=1
(0) ψ¯ 2i ν 2i
−
i=1
3
(3.8) (1) ψ¯ 2i ν 2i+1 = exp(−ν) − 1 + ν,
i=1
and 4
(0)
ψ2i ν 2i +
i=1
4 i=1
4
(1)
ψ2i ν 2i+1 = ν exp(ν) − ν,
i=1
(0) ψ2i ν 2i
−
4
(3.9) (1) ψ2i ν 2i+1
= −ν exp(−ν) + ν.
i=1
We consider the A- and c-values from the RKN method (3.1). We put Eqs. (3.8) with in addition a second-order condition for RKN methods (see [22]): b¯1 + b¯2 + b¯3 = 12 ,
(3.10)
¯ and we solve the obtained linear system for the unknown b-values. The following result emerges: 1 −42ν 4 − 119 sinh(ν)ν 3 + 4 −283 + 238 cosh(ν) ν 2 b¯1 = 630 − 3420 sinh(ν)ν + 7200 −1 + cosh(ν) /ν 4 , 1 51ν 4 + 17 sinh(ν)ν 3 + 8 47 − 17 cosh(ν) ν 2 b¯2 = 90 + 560 sinh(ν)ν + 1600 1 − cosh(ν) /ν 4 , 50 −3ν 2 − sinh(ν)ν + 8 −1 + cosh(ν) /ν 4 . b¯3 = 63
(3.11)
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
133
A crucial property of the classical fourth-order RKN method of Dormand et al. [21] is that the numerical solution of the y-component at each step is equal to the last stage. In that case one can use the FSAL technique (the last evaluation at any step is the same as the first evaluation at the next step). The number of stages is four but it effectively uses only three stages at each step. Our intention is to construct a RKN method that is also based on the FSAL technique. Therefore we put now a4i = b¯i (i = 1, 2, 3). In the same spirit we put Eq. (3.9) together with second-order conditions: b1 + b2 + b3 + b4 = 1, b1 c1 + b2 c2 + b3 c3 + b4 c4 = 12 ,
(3.12)
and we solve the obtained linear system for the unknown b-values. One finds: 2 b3 = 1000 −17 sinh(ν) ν 5 + 2 sinh(ν) −61 + 17 cosh(ν) ν 4 2 + 8 281 + 382 cosh(ν) + 57 cosh(ν) ν 3 − 16 sinh(ν) 347 + 13 cosh(ν) ν 2 − 256 17 + 28 cosh(ν) −1 + cosh(ν) ν + 11520 sinh(ν) −1 + cosh(ν) /(7ν 3 N ), b4 = 17 sinh(ν) −61 + 17 cosh(ν) ν 4 − 264 43 + 17 cosh(ν) ν 3 − 16 sinh(ν) 3659 + 391 cosh(ν) ν 2 + 31680 −1 + cosh(ν) ν + 129600 sinh(ν) −1 + cosh(ν) /(νN ), b2 = 12 − b3 c3 − b4 c4 /c2 ,
(3.13)
b1 = 1 − b2 − b3 − b4 , whereby 2 2 N = 289 sinh(ν) ν 4 − 62416 + 60928 cosh(ν) + 6256 cosh(ν) ν 2 2 + 129600 sinh(ν) .
(3.14)
For smaller values of |ν| the above formulae are subject to heavy cancellations and in that case the following Taylor series expansions must be used: b¯1 = b¯2 = b¯3 =
1 2 2 61 1 2719 7957 4 6 8 10 14 + 945 ν − 226800 ν − 84672 ν − 12573792000 ν − 3432645216000 ν − · · · , 8 7 61 1 1201 179 2 4 6 8 10 27 − 1620 ν + 226800 ν + 81648 ν + 5388768000 ν + 75442752000 ν + · · · , 25 5 1 1 1 2 6 8 10 189 + 2268 ν − 2286144 ν − 150885504 ν − 18307441152 ν − · · · ,
(3.15)
1 13 2 941 69649 3619447 8743305199 4 6 8 10 14 + 2520 ν − 476280 ν − 1028764800 ν + 855520807680 ν + 30028780349568000 ν − · · · , 32 79 3823 843067 54068059 1481673821 2 4 6 8 10 81 − 21870 ν + 1180980 ν + 8928208800 ν − 7424698438080 ν − 3384505391904000 ν + · · · , 250 125 2 34885 12995 42650603 226036171 4 6 8 10 567 − 15309 ν − 23147208 ν − 1249949232 ν + 10394577813312 ν + 1824248406236256 ν − · · · , 5 193 2 677 194389 946459 607036303 4 6 8 10 54 + 29160 ν + 2755620 ν − 11904278400 ν − 899963447040 ν + 26728914377088000 ν + · · · .
(3.16)
and b1 = b2 = b3 = b4 =
The series contain terms up to ν 10 in order to obtain the coefficients with arithmetic precision of 16 digits. We ¯ observe that for ν → 0 the b- and b-values from the classical fourth-order RKN method [21] are recovered.
134
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
It is also interesting to check the algebraic fourth-order conditions for RKN methods (see [22]): b¯ T .c =
b¯ T .e = 12 , b¯ T .c2 =
1 12
+
7 2 8640 ν
+
61 4 63628800 ν
+ ···,
bT .e = 1, bT .c2 = bT .c3 =
1 3 1 4
bT .A.c =
1 6
+
b¯ T .A.e =
1 61 2 4 2160 ν + 907200 ν + · · · , 1 7 61 2 4 24 + 17280 ν + 7257600 ν + · · · ,
bT .c = 12 , +
31 1423 2 4 12960 ν − 4898880 ν + · · · , 13 2 7207 + 3456 ν − 32659200 ν4 + · · · , 1 77 17387 2 4 24 + 116640 ν − 440899200 ν + · · · .
It is clear that the new RKN method is error: y(x0 + h) − y1 = O(h5 ) and
bT .A.e =
1 6
bT (c.Ae) =
+
31 1423 2 4 25920 ν − 9797760 ν + · · · , 1 13 2 7207 4 8 + 6912 ν − 65318400 ν + · · · ,
(3.17)
algebraically of fourth order, in the sense that we can write for the local y (x0 + h) − y = O(h5 ). 1
When we apply the RKN method to problem (2.1) the principal local truncation errors (PLTE) for the y-component and its derivative are given by: h5 (fy − µ2 )y (3) , 2160
h5 31(fy − µ2 )y (4) − 12(fyy y + fxy )y (3) . PLTE(y ) = 25920 It is visible that our RKN method integrates exactly problem (3.2). In particular, for problems of the form PLTE(y) =
y = µ2 y + g(x)
(µ ∈ R or µ ∈ iR),
(3.18)
(3.19)
it can be seen that our method obtains a raise of one order of accuracy with respect to the underlying classical RKN method. We will indeed verify in Section 4 that especially for this kind of problems or similar to those, the numerical performance of the new RKN method is superior to some previously derived RKN methods for oscillatory problems. Remark. Another possible way to construct a RKN method that integrates exactly problem (3.2) is replacing the 1 and then solving the system (3.9) for the unknowns b¯2 and b¯3 . One can follow a order condition (3.10) by b¯1 = 14 similar strategy for the b-values. By this action we obtain for the second-order conditions: b¯1 + b¯2 + b¯3 =
1 2
+ O(h4 ),
b1 + b2 + b3 + b4 = 1 + O(h4 ), b1 c1 + b2 c2 + b3 c3 + b4 c4 =
1 2
+ O(h3 ).
This has the consequence that the PLTE on y and y will contain µ4 , moreover it can be verified that there is no increase of order for problems of the form (3.19). 3.2. A third-order RKN method Now our interest lies in constructing an embedded pair of explicit RKN methods based on the above presented ¯ b). To simplify the construction of the method and guided by the classical fourth-order RKN method (c, A, b, embedded RKN4(3) pair from Dormand et al. [21] we consider the RKN method (c, A, b¯ ∗ , b∗ ) whereby: b¯3∗ =
3 20 ,
1 b¯4∗ = − 20
and b4∗ = − 13 .
(3.20)
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
135
We follow a similar approach as described above and we solve the linear system (3.8) for the unknown b¯ ∗ -values: 1 483ν 5 − 250 sinh(ν)ν 4 + 200 79 + 10 cosh(ν) ν 3 − 13000 sinh(ν)ν 2 b¯1∗ = 40000 + 40000 3 + cosh(ν) ν − 160000 sinh(ν) /ν 3 , (3.21) 1 −357ν 5 − 4200ν 3 + 2000 sinh(ν)ν 2 − 40000ν + 40000 sinh(ν) /ν 3 . b¯2∗ = 10000 In the same lines we solve Eqs. (3.9) with in addition the first equation from (3.12) for the b∗ -values. The following values are found: 1 −119 sinh(ν)ν 4 + −931 + 595 cosh(ν) ν 3 − 564 sinh(ν)ν 2 b1∗ = 1890 − 180 103 + 17 cosh(ν) ν + 21600 sinh(ν) ν 3 , 1 17 sinh(ν)ν 4 + 493 − 85 cosh(ν) ν 3 + 152 sinh(ν)ν 2 b2∗ = (3.22) 270 3 + 80 59 + cosh(ν) ν − 4800 sinh(ν) /ν , 50 − sinh(ν)ν 2 + −29 + 5 cosh(ν) ν + 24 sinh(ν) /ν 3 . b3∗ = 189 For smaller values of ν, say |ν| < 0.3, the following Taylor expansion series should be preferred: 7 b¯1∗ = − 150 −
b¯2∗ =
67 150
b1∗ =
13 21
+
1801 2 3 121 23 421 4 6 8 10 120000 ν − 2800 ν − 3628800 ν − 31933440 ν − 41513472000 ν − · · · , 929 2 31 23 13 1 4 6 8 10 30000 ν + 12600 ν + 453600 ν + 19958400 ν + 176904000 ν + · · · ,
(3.23)
and b2∗ b3∗
+
137 2 247 2959 14251 23483 4 6 8 10 1890 ν + 1587600 ν − 19051200 ν − 2514758400 ν − 245188944000 ν − · · · , 221 2 713 643 6329 41453 4 6 8 10 = − 20 27 − 1620 ν − 680400 ν + 4082400 ν + 1077753600 ν + 420323904000 ν + · · · , 145 2 85 5 31 67 4 6 8 10 = 275 189 + 2268 ν + 95256 ν − 2286144 ν − 150885504 ν − 23538138624 ν − · · · .
(3.24)
For ν → 0 the RKN pair reduces to the classical RKN4(3) pair from Dormand et al. [21]. Introducing the A-, b¯ ∗ -, and b∗ -coefficients into the order condition equations for a third-order RKN method (see [22]), we have: b¯ ∗T .e =
1 2
+
383 2 24000 ν
+
1 4 720 ν
+ ···,
b∗ T .e = 1, b∗ .A.e
=
1 6
b¯ ∗T .c = b∗ T .c =
+
197 2 17280 ν
+
1349 4 7257600 ν
+ ···,
b∗ T .c2
1 6 1 2
=
+
929 31 2 4 120000 ν + 50400 ν + · · · , 23 2 47 + 2160 ν + 129600 ν4 + · · · , 1 197 2 1349 4 3 + 8640 ν + 3628800 ν + · · · .
(3.25)
It is clear that the RKN method is algebraically of third order, in the sense that we can write for the local error: y(x0 + h) − y1 = O(h4 ) and y (x0 + h) − y = O(h4 ). 1 Based on the order condition equations we obtain the PLTE for the y-component and its derivative:
383 4 PLTE(y) = 24000 h (fy − µ2 )y (2) + (fxx + 2fxy y + fyy y 2 ) ,
23 4 PLTE(y ) = 4320 h 3(y (5) − µ2 y (3) ) − (fy − µ2 )y (3) .
(3.26)
Again, we observe that the RKN method integrates exactly problem (3.2). 3.3. Exponentially fitted RK(N) methods As an example we take Simos’ and Franco’s fourth-order EFRKN methods [19,20]. Just like our fourth-order RKN method these methods are also based on the classical fourth-order RKN method of Dormand et al. [21].
136
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
An error analysis shows that both EFRKN methods have the same PLTE, so they achieve the same accuracy. The main difference between both methods is that Franco’s method shares the FSAL property. Now, we will explain the difference between the RKN methods developed above and EFRKN methods. The aim of exponential fitting is to construct numerical methods that integrate exactly a conventionally chosen mixture of powers and exponential functions. It was explained in [17] that EFRK(N) methods integrate exactly every linear combination of functions from the set {exp(±µx) | µ ∈ R or µ ∈ iR}. The fact that the considered EFRKN methods are of exponential fitting type is reflected in the next expressions of the PLTE: h5 fy (y (3) − µ2 y (1) ), PLTE [19,20], y = 2160 (3.27)
h5 (4) 2 (2) (3) 2 (1) 31fy (y − µ y ) − 12(fyy y + fxy )(y − µ y ) . PLTE [19,20], y = 25920 Observing the expressions of the PLTEs (3.18)–(3.26) it is clear that our RKN methods are not exponentially fitted so they are only suitable for the numerical integration of perturbed oscillators. In general, the region of applicability of EFRK(N) methods is larger than those of RK(N) methods adapted to perturbed oscillators. It is quite misleading to suppose that our RKN methods are of exponential fitting type. Indeed, multistep methods that integrate exactly the unperturbed problem (3.2) are always exponentially fitted. With the error analysis we have illustrated that this is not the case for RK(N) methods. It is explained in [16,17] that for the construction of EFRK(N) methods it is necessary that beside the final stage, every internal stage has also to be fitted to problem (3.2).
4. Numerical experiments 4.1. The test problems We consider some popular test problems which appear in many papers on numerical methods for oscillatory problems. Problem 1. Inhomogeneous equation y + 100y = 99 sin(x),
y(0) = 1,
y (0) = 11.
(4.1)
The equation has been solved in the interval [0, 100] with µ = 10i. The exact solution is given by: y(x) = cos(10x) + sin(10x) + sin(x). Problem 2. An “almost periodic” orbit problem studied by Stiefel and Bettis [26] z + z = 0.001eix ,
z(0) = 1,
z (0) = 0.9995i.
(4.2)
The equation has been solved in the interval [0, 1000] with µ = i. The exact solution is given by: z(x) = (1 − 0.0005ix)eix . The solution represents a motion of a perturbation of a circular orbit in the complex plane. The problem may be solved either as a single equation in complex arithmetic or as a pair of uncoupled equations. Problem 3. An “almost periodic” orbit problem studied by Franco and Palacios [27] z + z = eiψx ,
z(0) = 1,
z (0) = i,
(4.3)
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
137
where = 0.001 and ψ = 0.01. The equation has been solved in the interval [0, 1000] with µ = i. The analytical solution z(x) = u(x) + iv(x) is given by: 1 − − ψ2 cos(x) + cos(ψx), 1 − ψ2 1 − ψ2 2 1 − ψ − ψ sin(x) + sin(ψx). v(x) = 2 1−ψ 1 − ψ2
u(x) =
Problem 4. The undamped Duffing’s equation y + y + y 3 = 0.002 cos[1.01x].
(4.4)
The equation has been solved in the interval [0, 300] with µ = i. The solution is not known in a closed form. Accuracy is judged with the Galerkin approximation obtained by Van Dooren [28], expressed as y(x) =
4
a2i+1 cos 1.01(2i + 1)x ,
i=0
with a1 = 0.200179477536, a3 = 0.246946143 · 10−3 , a5 = 0.304014 · 10−6 , a7 = 0.374 · 10−9 and a9 < 10−12 . The appropriate initial conditions are y(0) = 0.200426728067 and y (0) = 0. Problem 5. An equation related to Bessel’s equation 1 y + 100 + 2 y = 0. 4x
(4.5)
The √ equation has been solved in the interval [1, 10] with µ = 10i. The exact solution is given by y(x) = xJ0 (10x), where J0 is the Bessel function of the first kind of order zero. The initial values are y(1) = J0 (10) and y (1) = 12 J0 (10) − 10J1 (10), where J1 is the Bessel function of order 1. Problem 6. A perturbed system [12] y1 + 25y1 + y12 + y22 = f1 (x), y2 + 25y2 + y12 + y22 = f2 (x),
y1 (0) = 1, y2 (0) = ,
y1 (0) = 0,
y2 (0) = 5,
(4.6)
with = 10−3 and f1 (x) = 1 + 2 + 2 sin(5x + x 2 ) + 2 cos(x 2 ) + (25 − 4x 2 ) sin(x 2 ), f2 (x) = 1 + 2 + 2 sin(5x + x 2 ) − 2 sin(x 2 ) + (25 − 4x 2 ) cos(x 2 ). The equation has been solved in the interval [0, 5] with µ = 5i. The analytical solution is given by: y1 (x) = cos(5x) + sin(x 2 ),
y2 (x) = sin(5x) + cos(x 2 ).
4.2. Comparisons with fixed step-size We will compare the numerical performance of the new fourth-order RKN method with other RKN methods specifically designed for oscillatory problems. The RKN methods used in the comparisons have been denoted by: • RKNA16: The classical fourth-order RKN method with phase-lag order 10 and interval of periodicity [0, 3.59] derived by Van der Houwen and Sommeijer [1] (p. 616, indicated as A1.6).
138
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
Table 1 Slopes of the best fitting straight lines from Fig. 2
Inhomogeneous Stiefel–Bettis Franco–Palacios
RKN4
FRKN4
RKN4:6
3.886 3.982 3.977
4.888 4.962 –
6.010 – 6.008
• EFRKN4: The trigonometric version of the fourth-order EFRKN method presented by Franco [20] (Section 3.3). The FSAL property is satisfied. • FRKN4: The new fourth-order RKN method with FSAL property derived in Section 3.1. • RKN4:6: The fourth-order RKN method that reaches order six for the unperturbed problem (3.2) proposed by García et al. [9] (Section 2.2). • PHRKN4: The fourth-order modified RKN method with phase-lag of order infinity presented by Simos and Vigo-Aguiar [7] (Section 3). • ARKN4: The fourth-order ARKN method derived by Franco [11] (Section 3). • ARKN5: The fifth-order ARKN method given by Franco [12] (Section 3). 4.2.1. Testing the efficiency In Fig. 1 we display the efficiency curves: accuracy versus the computational cost measured by the number of function evaluations required by each method. It is not surprising that the numerical performance of FRKN4 is superior to all the other RKN methods for Problems 1–3. For Duffing’s equation we observe that EFRKN4 and ARKN5 are the most efficient methods. For Bessel’s equation is ARKN5 the most efficient one and we remark that the efficiency curves for FRKN4 and EFRKN4 coincide. On the other hand, for the perturbed system is FRKN4 the most efficient method. It is difficult to choose the most suitable method for a given problem. However, we have showed that FRKN4 is suitable for the numerical integration of perturbed oscillators. 4.2.2. Verification of the error analysis Another goal with the numerical experiments in fixed step-size mode is to illustrate that an application of the new fourth-order RKN method to linear problems (3.19) will raise the order of accuracy by one unit with respect to the underlying classical RKN method. For this part of the numerical experiments we only consider the classical fourthorder RKN method of Dormand et al. [21] (denoted as RKN4), FRKN4 and RKN4:6. These methods use three function evaluations at each step. In Fig. 2 we display − log10 (h) versus − log10 (GE) where GE is the maximum global error over the whole integration interval when we apply the considered RKN methods to Problems 1–3. The slopes of the best fitting straight lines are showed in Table 1. Clearly, these numbers confirm that a raise of order is obtained for FRKN4 applied to Problems 1–3. This is a direct reflection of the expressions of the PLTEs (3.18) and the considerations made in Section 3.1. Beside this we have detected a raise of order by two units when we apply RKN4:6 to problems of the form (3.19). Although RKN4:6 reaches order six while FRKN4 reaches only order five for linear problems (3.19), the most accurate method is FRKN4. This can be easily understood by the fact that the construction of RKN4:6 is based on requirements that are softened: it does not exactly integrate the unperturbed problem (3.2). 4.3. Comparisons with variable step-size The codes used in the comparisons have been denoted by: • EFRKN4(3): The trigonometric version of the embedded 4(3) pair with FSAL property presented by Franco [20] (Sections 3.3 and 3.4). • FRKN4(3): The new pair derived in this paper. The FSAL technique is used. • RKN4:6(3:4): The embedded 4(3) pair proposed by García et al. [9] (Section 3). The lower order method reaches order four for the unperturbed problem (3.2).
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
139
Fig. 1. The number of function evaluations versus log10 (GE) where GE is the maximum global error over the whole integration interval when we apply the considered RKN methods to Problems 1–6.
140
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
Fig. 2. − log10 (h) versus − log10 (GE) where GE is the maximum global error over the whole integration interval when we apply RKN4, FRKN4 and RKN4:6 to Problems 1–3.
• ARKN4(3): The embedded 4(3) pair of ARKN methods derived by Franco [11] (Section 3). • ARKN5(3): The embedded 5(3) pair of ARKN methods given by Franco [12] (Section 3). The efficiency curves are drawn in Fig. 3. 5. Conclusions We have introduced an new technique for the construction of RKN methods fitted to perturbed oscillators. It can be seen as an extension of Simos’ RK approach [15]. The new RKN methods integrate exactly the unperturbed problem (3.2). Based on this technique we present a new embedded 4(3) pair of RKN methods that shares the FSAL property. We have showed that the new fourth-order RKN method achieves order five for linear problems of the form (3.19). An error analysis, a discussion and a comparison with some previously derived RKN methods (such as EFRKN methods) are made. Although, it is difficult to determine the most efficient RKN method for a given oscillatory problem, numerical experiments reveal an acceptable performance of the new RKN methods compared with some other RKN methods especially designed for oscillatory problems.
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
141
Fig. 3. The number of function evaluations versus log10 (GE) where GE is the maximum global error over the whole integration interval when we apply the considered embedded RKN methods to Problems 1–6.
142
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
Acknowledgements The author wish to thank the anonymous referees for their careful reading of the manuscript and their valuable comments and fruitful suggestions that greatly improved this paper.
References [1] P.J. Van der Houwen, B.P. Sommeijer, Explicit Runge–Kutta(–Nyström) methods with reduced phase errors for computing oscillating solution, SIAM J. Numer. Anal. 24 (1987) 595–617. [2] P.J. Van der Houwen, B.P. Sommeijer, Phase-lag analysis of implicit Runge–Kutta methods, SIAM J. Numer. Anal. 26 (1989) 214–229. [3] P.J. Van der Houwen, B.P. Sommeijer, Diagonally implicit Runge–Kutta–Nyström methods for oscillatory problems, SIAM J. Numer. Anal. 26 (1989) 414–429. [4] B. Paternoster, A phase-fitted collocation-based Runge–Kutta(–Nyström) method, Appl. Numer. Math. 35 (2000) 339–355. [5] T.E. Simos, E. Dimas, A.B. Sideridis, A Runge–Kutta–Nyström method for the numerical integration of special second-order periodic initial-value problems, J. Comput. Appl. Math. 51 (1994) 317–326. [6] T.E. Simos, A modified Runge–Kutta method for the numerical solution of ODE’s with oscillation solutions, Appl. Math. Lett. 6 (1996) 61–66. [7] T.E. Simos, J.V. Aguiar, A modified Runge–Kutta–Nyström method with phase-lag of order infinity for the numerical solution of the Schrödinger equation and related problems, Int. J. Mod. Phys. C 11 (2000) 1195–1208. [8] A.B. González, P. Martín, J.M. Farto, A new family of Runge–Kutta type methods for the numerical integration of perturbed oscillators, Numer. Math. 82 (1999) 635–646. [9] A. García, P. Martín, A.B. González, New methods for oscillatory problems based on classical codes, Appl. Numer. Math. 42 (2002) 141–157. [10] J.M. Franco, Runge–Kutta–Nyström methods adapted to the numerical integration of perturbed oscillators, Comput. Phys. Comm. 147 (2002) 770–787. [11] J.M. Franco, Embedded pairs of explicit ARKN methods for the numerical integration of perturbed oscillators, J. Comput. Methods Sci. Engrg. 3 (2003) 415–424. [12] J.M. Franco, A 5(3) pair of explicit ARKN methods for the numerical integration of perturbed oscillators, J. Comput. Appl. Math. 161 (2003) 283–293. [13] Ch. Tsitouras, T.E. Simos, Optimized Runge–Kutta pairs for problems with oscillating solutions, J. Comput. Appl. Math. 147 (2002) 397–409. [14] Z.A. Anastassi, T.E. Simos, Special optimized Runge–Kutta methods for IVP’s with oscillating solutions, Int. J. Mod. Phys. C 15 (2004) 1–15. [15] T.E. Simos, An exponentially-fitted Runge–Kutta method for the numerical integration of initial-value problems with periodic or oscillating solutions, Comput. Phys. Comm. 115 (1998) 1–8. [16] B. Paternoster, Runge–Kutta(–Nyström) methods for ODEs with periodic solutions based in trigonometric polynomials, Appl. Numer. Math. 28 (1998) 401–412. [17] G. Vanden Berghe, H. De Meyer, M. Van Daele, T. Van Hecke, Exponentially-fitted explicit Runge–Kutta methods, Comput. Phys. Comm. 123 (1999) 7–15. [18] J.P. Coleman, S.C. Duxbury, Mixed collocation methods for y = f (x, y), J. Comput. Appl. Math. 126 (2000) 47–75. [19] T.E. Simos, Exponentially-fitted Runge–Kutta–Nyström method for the numerical solution of initial-value problems with oscillating solutions, Appl. Math. Lett. 15 (2002) 217–225. [20] J.M. Franco, Exponentially fitted explicit Runge–Kutta–Nyström methods, J. Comput. Appl. Math. 167 (2004) 1–19. [21] J.R. Dormand, M.E.A. El-Mikkawy, P.J. Prince, Families of Runge–Kutta–Nyström formulae, IMA J. Numer. Anal. 7 (1987) 235–250. [22] E. Hairer, S.P. Nørsett, G. Wanner, Solving Ordinary Differential Equations I, Nonstiff Problems, second ed., Springer-Verlag, Berlin, 1993. [23] L.F. Shampine, Local error estimation in the solution of ordinary differential equations, Math. Comp. 27 (1973) 91–97. [24] L.F. Shampine, Local error estimation by doubling, Computing 34 (1985) 179–190. [25] L.F. Shampine, Some practical Runge–Kutta formulas, Math. Comp. 46 (1986) 135–150. [26] E. Stiefel, D.G. Bettis, Stabilization of Cowell’s method, Numer. Math 13 (1969) 154–175. [27] J.M. Franco, M. Palacios, High-order P-stable multistep methods, J. Comput. Appl. Math. 30 (1990) 1–10. [28] R. Van Dooren, Stabilization of Cowell’s classical finite difference methods for numerical integration, J. Comput. Phys. 16 (1974) 186– 192.
A Runge–Kutta–Nyström pair for the numerical integration of perturbed oscillators Hans Van de Vyver Received 23 June 2004; accepted 30 December 2004
Abstract New Runge–Kutta–Nyström methods especially designed for the numerical integration of perturbed oscillators are presented in this paper. They are capable of exactly integrating the harmonic or unperturbed oscillator. We construct an embedded 4(3) RKN pair that is based on the FSAL property. The new method is much more efficient than previously derived RKN methods for some subclasses of problems. 2005 Elsevier B.V. All rights reserved. PACS: 02.60.Jh; 02.70.Bf Keywords: Explicit Runge–Kutta–Nyström methods; Embedded pairs; Variable step-size algorithms; Initial value problems; Oscillating solutions
1. Introduction We are concerned with initial value problems (IVP) related to non-stiff second-order ordinary differential equations (ODE) of the form y = µ2 y + f (x, y),
y(x0 ) = y0 ,
y (x0 ) = y0 ,
|µ2 |
µ ∈ R or iR.
(1.1)
In the case that |f (x, y)| the solution exhibits a pronounced exponential (µ ∈ R) or oscillatory (µ ∈ iR) character. Equations having oscillatory solutions are of particular interest, examples occur frequently in celestial mechanics, in quantum mechanical scattering problems and elsewhere. There are many procedures in existence for the construction of Runge–Kutta(–Nyström) (RK(N)) methods for an efficient integration of oscillatory problems. An exhaustive list of references would be near impossible. For example, in the mid-eighties Van der Houwen et al. [1–3] constructed RK(N) methods with a reduced or null phase error. A first example of a phase-fitted RKN method (which means that the phase-lag order is infinity) is due to Paternoster [4]. Simos et al. [5–7] applied some modified explicit RK(N) methods (phase-fitted or with minimal phase-lag) to solve Schrödinger equations. E-mail address: [email protected] (H. Van de Vyver). 0010-4655/$ – see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2004.12.011
130
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
More recently, new approaches for constructing RK(N) methods adapted to perturbed oscillators are presented by González et al. [8], García et al. [9] and Franco [10–12]. Other optimized RK methods for oscillatory problems are designed by Tsitouras and Simos [13] and Anastassi and Simos [14]. For an efficient numerical integration of oscillatory problems is exponential fitting also an interesting option. The study of exponentially fitted RK(N) (EFRK(N)) methods is a relatively new development and we refer to [15–18] as pioneering work in this area. The first implicit EFRK(N) methods are presented by Paternoster [16]. A disadvantage of these methods is that they are fully implicit which means that a high computational effort is required. Very recently, explicit EFRKN methods are presented by Simos [19] and Franco [20]. In particular, in [20] was an embedded pair of EFRKN methods and a variable step-size algorithm developed and discussed. The pair has the RKN pair of orders 4(3) derived by Dormand et al. [21] as underlying classical pair. In this paper we present a technique for the construction of RKN methods that integrate exactly Eq. (1.1) with f (x, y) = 0. We will construct a new embedded pair of RKN methods that is also based on the RKN4(3) pair from [21]. The derivation technique can be seen as an extension of Simos’ RK approach [15] but here we start from another point of view. Numerical examples reveal that the new embedded pair is much more efficient than other ones for some subclasses of (1.1). The paper is organized as follow: in Section 2 we present some basic elements of RKN methods and embedded pairs. In Section 3 we describe the derivation of the new RKN methods and the expressions of the local truncation errors will be discussed. We will explain the difference between the RKN approach developed in this paper and EFRK(N) methods. In Section 4 some numerical experiments show the properties of the developed algorithms. Finally, in Section 5 some conclusions are drawn.
2. Embedded pairs of RKN methods For the numerical solution of the system of second-order ODEs y = f (x, y),
y(x0 ) = y0 ,
y (x0 ) = y0 ,
(2.1)
we consider the explicit Runge–Kutta–Nyström (RKN) method Yi = yn + ci hyn + h2
i−1
aij f (xn + cj h, Yj ), j =1 s b¯i f (xn + ci h, Yi ), yn+1 = yn + hyn + h2 i=1 s yn+1 = yn + h bi f (xn + ci h, Yi ). i=1
i = 1, . . . , s, (2.2)
The RKN method is completely determined by means of its Butcher tableau 0 c2 a21 c3 a31 .. .. . . cs as1 b¯1 b1
a32 .. . . . . as2 . . . b¯2 . . . b2 . . .
¯ b), or equivalent by the quadruple (c, A, b,
(2.3)
as,s−1 b¯s−1 b¯s bs−1 bs
where c, b, b¯ ∈ Rs , A ∈ Rs×s and s denotes the number of stages of the RKN method. An embedded q(p) pair of ¯ b) of order q and another RKN method (c, A, b¯ ∗ , b∗ ) of order RKN methods is based on the RKN method (c, A, b,
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
131
p < q. Embedded pairs of explicit RK(N) methods are widely used in variable step-size algorithms because they provide a cheap error estimation. A local error estimation in the integration point xn+1 = xn + h is determined by the expressions ∗ δn+1 = yn+1 − yn+1
∗ and δn+1 = yn+1 − yn+1 .
(2.4)
EST n+1 (h) = max{δn+1 ∞ , δn+1 ∞ } represents the local error estimation to control the step-size h for RKN processes by the well-known standard formula (see [22]) 1/(p+1) Tol hnew = 0.9hold (2.5) , EST n+1 (hold )
where Tol is the maximum allowable error. If EST n+1 (h) < Tol, then the step is accepted and we adopt the widely used procedure of performing local extrapolation: although, we are actually controlling an estimate of the local error in the lower order solution, it is the higher order solution we accept at each point. If EST n+1 (h) Tol then the step is rejected and we calculate hnew given by (2.5). It should be noted that several variants of local error estimation and step-size control algorithms are in existence (for example, Shampine [23–25]).
3. Explicit RKN methods fitted to perturbed oscillators 3.1. A fourth-order RKN method We consider the four-stage explicit RKN method which is in short-hand notation presented by the following Butcher-tableau: 0 1 1 4 32 7 119 7 10 1000 500
1 a41 a42 b¯1 b¯2 b1 b2
(3.1) a43 b¯3 0 b3 b4
The A-values from the first three rows are taken from the well-known fourth-order RKN method derived by Dormand et al. [21]. We are focused on the unperturbed problem: y = µ2 y,
y(x0 ) = y0 ,
y (x0 ) = y0
(µ ∈ R or µ ∈ iR).
(3.2)
The exact solution is given by: y(x) = χ1 exp(µx) + χ2 exp(−µx),
(3.3)
whereby χ1 , χ2 ∈ C are determined by the initial conditions. We remark that in practice for most of the applications is µ ∈ iR. Our aim is to determine the b- and b∗ -values in the way that the RKN method (3.1) integrates exactly problem (3.2). We solve (3.2) with the RKN method and we obtain for the numerical solution at x1 = x0 + h: 3 3 (0) (1) ψ¯ 2i µ2i y0 + 1 + ψ¯ 2i µ2i hy0 , y1 = 1 + (3.4) i=1
i=1
132
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
hy1
=
4
(0) ψ2i µ2i
y0 + 1 +
i=1
4
(1) ψ2i µ2i
hy0 .
(3.5)
i=1
¯ The ψ-values are given by: (0) ψ¯ 2 =
3
b¯i ,
(0) ψ¯ 4 =
(1) ψ¯ 2
=
(0) ψ¯ 6 = b¯3 a32 a21 ,
b¯i aij ,
i=2 j =1
i=1 3
3 i−1
b¯i ci ,
(2) ψ¯ 4
(3.6)
= b¯3 a32 c2 ,
(1) ψ¯ 6
= 0.
i=2
The ψ -values are given by: (0) ψ2
=
4
(0) ψ4
bi ,
=
i=1 (1) ψ2
=
4
bi ci ,
(1) ψ4
4 i−1
bi aij ,
i=2 j =1 4 i−1
=
(0) ψ6
=
j −1 4 i−1
bi aij aj k ,
(0)
ψ8 = b4 a43 a32 a21 ,
i=3 j =2 k=1
bi aij cj ,
(1) ψ6
= b4 a43 a32 c2 ,
(3.7) (1) ψ8
= 0.
i=3 j =2
i=2
Substituting y0 = y(x), y0 = y (x), y1 = y(x + h) and y1 = y (x + h) in (3.4) and (3.5) (y(x) given in (3.3)), one easily finds the next relations (ν = µh): 3
(0) ψ¯ 2i ν 2i +
i=1
3
3
(1) ψ¯ 2i ν 2i+1 = exp(ν) − 1 − ν,
i=1
(0) ψ¯ 2i ν 2i
−
i=1
3
(3.8) (1) ψ¯ 2i ν 2i+1 = exp(−ν) − 1 + ν,
i=1
and 4
(0)
ψ2i ν 2i +
i=1
4 i=1
4
(1)
ψ2i ν 2i+1 = ν exp(ν) − ν,
i=1
(0) ψ2i ν 2i
−
4
(3.9) (1) ψ2i ν 2i+1
= −ν exp(−ν) + ν.
i=1
We consider the A- and c-values from the RKN method (3.1). We put Eqs. (3.8) with in addition a second-order condition for RKN methods (see [22]): b¯1 + b¯2 + b¯3 = 12 ,
(3.10)
¯ and we solve the obtained linear system for the unknown b-values. The following result emerges: 1 −42ν 4 − 119 sinh(ν)ν 3 + 4 −283 + 238 cosh(ν) ν 2 b¯1 = 630 − 3420 sinh(ν)ν + 7200 −1 + cosh(ν) /ν 4 , 1 51ν 4 + 17 sinh(ν)ν 3 + 8 47 − 17 cosh(ν) ν 2 b¯2 = 90 + 560 sinh(ν)ν + 1600 1 − cosh(ν) /ν 4 , 50 −3ν 2 − sinh(ν)ν + 8 −1 + cosh(ν) /ν 4 . b¯3 = 63
(3.11)
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
133
A crucial property of the classical fourth-order RKN method of Dormand et al. [21] is that the numerical solution of the y-component at each step is equal to the last stage. In that case one can use the FSAL technique (the last evaluation at any step is the same as the first evaluation at the next step). The number of stages is four but it effectively uses only three stages at each step. Our intention is to construct a RKN method that is also based on the FSAL technique. Therefore we put now a4i = b¯i (i = 1, 2, 3). In the same spirit we put Eq. (3.9) together with second-order conditions: b1 + b2 + b3 + b4 = 1, b1 c1 + b2 c2 + b3 c3 + b4 c4 = 12 ,
(3.12)
and we solve the obtained linear system for the unknown b-values. One finds: 2 b3 = 1000 −17 sinh(ν) ν 5 + 2 sinh(ν) −61 + 17 cosh(ν) ν 4 2 + 8 281 + 382 cosh(ν) + 57 cosh(ν) ν 3 − 16 sinh(ν) 347 + 13 cosh(ν) ν 2 − 256 17 + 28 cosh(ν) −1 + cosh(ν) ν + 11520 sinh(ν) −1 + cosh(ν) /(7ν 3 N ), b4 = 17 sinh(ν) −61 + 17 cosh(ν) ν 4 − 264 43 + 17 cosh(ν) ν 3 − 16 sinh(ν) 3659 + 391 cosh(ν) ν 2 + 31680 −1 + cosh(ν) ν + 129600 sinh(ν) −1 + cosh(ν) /(νN ), b2 = 12 − b3 c3 − b4 c4 /c2 ,
(3.13)
b1 = 1 − b2 − b3 − b4 , whereby 2 2 N = 289 sinh(ν) ν 4 − 62416 + 60928 cosh(ν) + 6256 cosh(ν) ν 2 2 + 129600 sinh(ν) .
(3.14)
For smaller values of |ν| the above formulae are subject to heavy cancellations and in that case the following Taylor series expansions must be used: b¯1 = b¯2 = b¯3 =
1 2 2 61 1 2719 7957 4 6 8 10 14 + 945 ν − 226800 ν − 84672 ν − 12573792000 ν − 3432645216000 ν − · · · , 8 7 61 1 1201 179 2 4 6 8 10 27 − 1620 ν + 226800 ν + 81648 ν + 5388768000 ν + 75442752000 ν + · · · , 25 5 1 1 1 2 6 8 10 189 + 2268 ν − 2286144 ν − 150885504 ν − 18307441152 ν − · · · ,
(3.15)
1 13 2 941 69649 3619447 8743305199 4 6 8 10 14 + 2520 ν − 476280 ν − 1028764800 ν + 855520807680 ν + 30028780349568000 ν − · · · , 32 79 3823 843067 54068059 1481673821 2 4 6 8 10 81 − 21870 ν + 1180980 ν + 8928208800 ν − 7424698438080 ν − 3384505391904000 ν + · · · , 250 125 2 34885 12995 42650603 226036171 4 6 8 10 567 − 15309 ν − 23147208 ν − 1249949232 ν + 10394577813312 ν + 1824248406236256 ν − · · · , 5 193 2 677 194389 946459 607036303 4 6 8 10 54 + 29160 ν + 2755620 ν − 11904278400 ν − 899963447040 ν + 26728914377088000 ν + · · · .
(3.16)
and b1 = b2 = b3 = b4 =
The series contain terms up to ν 10 in order to obtain the coefficients with arithmetic precision of 16 digits. We ¯ observe that for ν → 0 the b- and b-values from the classical fourth-order RKN method [21] are recovered.
134
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
It is also interesting to check the algebraic fourth-order conditions for RKN methods (see [22]): b¯ T .c =
b¯ T .e = 12 , b¯ T .c2 =
1 12
+
7 2 8640 ν
+
61 4 63628800 ν
+ ···,
bT .e = 1, bT .c2 = bT .c3 =
1 3 1 4
bT .A.c =
1 6
+
b¯ T .A.e =
1 61 2 4 2160 ν + 907200 ν + · · · , 1 7 61 2 4 24 + 17280 ν + 7257600 ν + · · · ,
bT .c = 12 , +
31 1423 2 4 12960 ν − 4898880 ν + · · · , 13 2 7207 + 3456 ν − 32659200 ν4 + · · · , 1 77 17387 2 4 24 + 116640 ν − 440899200 ν + · · · .
It is clear that the new RKN method is error: y(x0 + h) − y1 = O(h5 ) and
bT .A.e =
1 6
bT (c.Ae) =
+
31 1423 2 4 25920 ν − 9797760 ν + · · · , 1 13 2 7207 4 8 + 6912 ν − 65318400 ν + · · · ,
(3.17)
algebraically of fourth order, in the sense that we can write for the local y (x0 + h) − y = O(h5 ). 1
When we apply the RKN method to problem (2.1) the principal local truncation errors (PLTE) for the y-component and its derivative are given by: h5 (fy − µ2 )y (3) , 2160
h5 31(fy − µ2 )y (4) − 12(fyy y + fxy )y (3) . PLTE(y ) = 25920 It is visible that our RKN method integrates exactly problem (3.2). In particular, for problems of the form PLTE(y) =
y = µ2 y + g(x)
(µ ∈ R or µ ∈ iR),
(3.18)
(3.19)
it can be seen that our method obtains a raise of one order of accuracy with respect to the underlying classical RKN method. We will indeed verify in Section 4 that especially for this kind of problems or similar to those, the numerical performance of the new RKN method is superior to some previously derived RKN methods for oscillatory problems. Remark. Another possible way to construct a RKN method that integrates exactly problem (3.2) is replacing the 1 and then solving the system (3.9) for the unknowns b¯2 and b¯3 . One can follow a order condition (3.10) by b¯1 = 14 similar strategy for the b-values. By this action we obtain for the second-order conditions: b¯1 + b¯2 + b¯3 =
1 2
+ O(h4 ),
b1 + b2 + b3 + b4 = 1 + O(h4 ), b1 c1 + b2 c2 + b3 c3 + b4 c4 =
1 2
+ O(h3 ).
This has the consequence that the PLTE on y and y will contain µ4 , moreover it can be verified that there is no increase of order for problems of the form (3.19). 3.2. A third-order RKN method Now our interest lies in constructing an embedded pair of explicit RKN methods based on the above presented ¯ b). To simplify the construction of the method and guided by the classical fourth-order RKN method (c, A, b, embedded RKN4(3) pair from Dormand et al. [21] we consider the RKN method (c, A, b¯ ∗ , b∗ ) whereby: b¯3∗ =
3 20 ,
1 b¯4∗ = − 20
and b4∗ = − 13 .
(3.20)
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
135
We follow a similar approach as described above and we solve the linear system (3.8) for the unknown b¯ ∗ -values: 1 483ν 5 − 250 sinh(ν)ν 4 + 200 79 + 10 cosh(ν) ν 3 − 13000 sinh(ν)ν 2 b¯1∗ = 40000 + 40000 3 + cosh(ν) ν − 160000 sinh(ν) /ν 3 , (3.21) 1 −357ν 5 − 4200ν 3 + 2000 sinh(ν)ν 2 − 40000ν + 40000 sinh(ν) /ν 3 . b¯2∗ = 10000 In the same lines we solve Eqs. (3.9) with in addition the first equation from (3.12) for the b∗ -values. The following values are found: 1 −119 sinh(ν)ν 4 + −931 + 595 cosh(ν) ν 3 − 564 sinh(ν)ν 2 b1∗ = 1890 − 180 103 + 17 cosh(ν) ν + 21600 sinh(ν) ν 3 , 1 17 sinh(ν)ν 4 + 493 − 85 cosh(ν) ν 3 + 152 sinh(ν)ν 2 b2∗ = (3.22) 270 3 + 80 59 + cosh(ν) ν − 4800 sinh(ν) /ν , 50 − sinh(ν)ν 2 + −29 + 5 cosh(ν) ν + 24 sinh(ν) /ν 3 . b3∗ = 189 For smaller values of ν, say |ν| < 0.3, the following Taylor expansion series should be preferred: 7 b¯1∗ = − 150 −
b¯2∗ =
67 150
b1∗ =
13 21
+
1801 2 3 121 23 421 4 6 8 10 120000 ν − 2800 ν − 3628800 ν − 31933440 ν − 41513472000 ν − · · · , 929 2 31 23 13 1 4 6 8 10 30000 ν + 12600 ν + 453600 ν + 19958400 ν + 176904000 ν + · · · ,
(3.23)
and b2∗ b3∗
+
137 2 247 2959 14251 23483 4 6 8 10 1890 ν + 1587600 ν − 19051200 ν − 2514758400 ν − 245188944000 ν − · · · , 221 2 713 643 6329 41453 4 6 8 10 = − 20 27 − 1620 ν − 680400 ν + 4082400 ν + 1077753600 ν + 420323904000 ν + · · · , 145 2 85 5 31 67 4 6 8 10 = 275 189 + 2268 ν + 95256 ν − 2286144 ν − 150885504 ν − 23538138624 ν − · · · .
(3.24)
For ν → 0 the RKN pair reduces to the classical RKN4(3) pair from Dormand et al. [21]. Introducing the A-, b¯ ∗ -, and b∗ -coefficients into the order condition equations for a third-order RKN method (see [22]), we have: b¯ ∗T .e =
1 2
+
383 2 24000 ν
+
1 4 720 ν
+ ···,
b∗ T .e = 1, b∗ .A.e
=
1 6
b¯ ∗T .c = b∗ T .c =
+
197 2 17280 ν
+
1349 4 7257600 ν
+ ···,
b∗ T .c2
1 6 1 2
=
+
929 31 2 4 120000 ν + 50400 ν + · · · , 23 2 47 + 2160 ν + 129600 ν4 + · · · , 1 197 2 1349 4 3 + 8640 ν + 3628800 ν + · · · .
(3.25)
It is clear that the RKN method is algebraically of third order, in the sense that we can write for the local error: y(x0 + h) − y1 = O(h4 ) and y (x0 + h) − y = O(h4 ). 1 Based on the order condition equations we obtain the PLTE for the y-component and its derivative:
383 4 PLTE(y) = 24000 h (fy − µ2 )y (2) + (fxx + 2fxy y + fyy y 2 ) ,
23 4 PLTE(y ) = 4320 h 3(y (5) − µ2 y (3) ) − (fy − µ2 )y (3) .
(3.26)
Again, we observe that the RKN method integrates exactly problem (3.2). 3.3. Exponentially fitted RK(N) methods As an example we take Simos’ and Franco’s fourth-order EFRKN methods [19,20]. Just like our fourth-order RKN method these methods are also based on the classical fourth-order RKN method of Dormand et al. [21].
136
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
An error analysis shows that both EFRKN methods have the same PLTE, so they achieve the same accuracy. The main difference between both methods is that Franco’s method shares the FSAL property. Now, we will explain the difference between the RKN methods developed above and EFRKN methods. The aim of exponential fitting is to construct numerical methods that integrate exactly a conventionally chosen mixture of powers and exponential functions. It was explained in [17] that EFRK(N) methods integrate exactly every linear combination of functions from the set {exp(±µx) | µ ∈ R or µ ∈ iR}. The fact that the considered EFRKN methods are of exponential fitting type is reflected in the next expressions of the PLTE: h5 fy (y (3) − µ2 y (1) ), PLTE [19,20], y = 2160 (3.27)
h5 (4) 2 (2) (3) 2 (1) 31fy (y − µ y ) − 12(fyy y + fxy )(y − µ y ) . PLTE [19,20], y = 25920 Observing the expressions of the PLTEs (3.18)–(3.26) it is clear that our RKN methods are not exponentially fitted so they are only suitable for the numerical integration of perturbed oscillators. In general, the region of applicability of EFRK(N) methods is larger than those of RK(N) methods adapted to perturbed oscillators. It is quite misleading to suppose that our RKN methods are of exponential fitting type. Indeed, multistep methods that integrate exactly the unperturbed problem (3.2) are always exponentially fitted. With the error analysis we have illustrated that this is not the case for RK(N) methods. It is explained in [16,17] that for the construction of EFRK(N) methods it is necessary that beside the final stage, every internal stage has also to be fitted to problem (3.2).
4. Numerical experiments 4.1. The test problems We consider some popular test problems which appear in many papers on numerical methods for oscillatory problems. Problem 1. Inhomogeneous equation y + 100y = 99 sin(x),
y(0) = 1,
y (0) = 11.
(4.1)
The equation has been solved in the interval [0, 100] with µ = 10i. The exact solution is given by: y(x) = cos(10x) + sin(10x) + sin(x). Problem 2. An “almost periodic” orbit problem studied by Stiefel and Bettis [26] z + z = 0.001eix ,
z(0) = 1,
z (0) = 0.9995i.
(4.2)
The equation has been solved in the interval [0, 1000] with µ = i. The exact solution is given by: z(x) = (1 − 0.0005ix)eix . The solution represents a motion of a perturbation of a circular orbit in the complex plane. The problem may be solved either as a single equation in complex arithmetic or as a pair of uncoupled equations. Problem 3. An “almost periodic” orbit problem studied by Franco and Palacios [27] z + z = eiψx ,
z(0) = 1,
z (0) = i,
(4.3)
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
137
where = 0.001 and ψ = 0.01. The equation has been solved in the interval [0, 1000] with µ = i. The analytical solution z(x) = u(x) + iv(x) is given by: 1 − − ψ2 cos(x) + cos(ψx), 1 − ψ2 1 − ψ2 2 1 − ψ − ψ sin(x) + sin(ψx). v(x) = 2 1−ψ 1 − ψ2
u(x) =
Problem 4. The undamped Duffing’s equation y + y + y 3 = 0.002 cos[1.01x].
(4.4)
The equation has been solved in the interval [0, 300] with µ = i. The solution is not known in a closed form. Accuracy is judged with the Galerkin approximation obtained by Van Dooren [28], expressed as y(x) =
4
a2i+1 cos 1.01(2i + 1)x ,
i=0
with a1 = 0.200179477536, a3 = 0.246946143 · 10−3 , a5 = 0.304014 · 10−6 , a7 = 0.374 · 10−9 and a9 < 10−12 . The appropriate initial conditions are y(0) = 0.200426728067 and y (0) = 0. Problem 5. An equation related to Bessel’s equation 1 y + 100 + 2 y = 0. 4x
(4.5)
The √ equation has been solved in the interval [1, 10] with µ = 10i. The exact solution is given by y(x) = xJ0 (10x), where J0 is the Bessel function of the first kind of order zero. The initial values are y(1) = J0 (10) and y (1) = 12 J0 (10) − 10J1 (10), where J1 is the Bessel function of order 1. Problem 6. A perturbed system [12] y1 + 25y1 + y12 + y22 = f1 (x), y2 + 25y2 + y12 + y22 = f2 (x),
y1 (0) = 1, y2 (0) = ,
y1 (0) = 0,
y2 (0) = 5,
(4.6)
with = 10−3 and f1 (x) = 1 + 2 + 2 sin(5x + x 2 ) + 2 cos(x 2 ) + (25 − 4x 2 ) sin(x 2 ), f2 (x) = 1 + 2 + 2 sin(5x + x 2 ) − 2 sin(x 2 ) + (25 − 4x 2 ) cos(x 2 ). The equation has been solved in the interval [0, 5] with µ = 5i. The analytical solution is given by: y1 (x) = cos(5x) + sin(x 2 ),
y2 (x) = sin(5x) + cos(x 2 ).
4.2. Comparisons with fixed step-size We will compare the numerical performance of the new fourth-order RKN method with other RKN methods specifically designed for oscillatory problems. The RKN methods used in the comparisons have been denoted by: • RKNA16: The classical fourth-order RKN method with phase-lag order 10 and interval of periodicity [0, 3.59] derived by Van der Houwen and Sommeijer [1] (p. 616, indicated as A1.6).
138
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
Table 1 Slopes of the best fitting straight lines from Fig. 2
Inhomogeneous Stiefel–Bettis Franco–Palacios
RKN4
FRKN4
RKN4:6
3.886 3.982 3.977
4.888 4.962 –
6.010 – 6.008
• EFRKN4: The trigonometric version of the fourth-order EFRKN method presented by Franco [20] (Section 3.3). The FSAL property is satisfied. • FRKN4: The new fourth-order RKN method with FSAL property derived in Section 3.1. • RKN4:6: The fourth-order RKN method that reaches order six for the unperturbed problem (3.2) proposed by García et al. [9] (Section 2.2). • PHRKN4: The fourth-order modified RKN method with phase-lag of order infinity presented by Simos and Vigo-Aguiar [7] (Section 3). • ARKN4: The fourth-order ARKN method derived by Franco [11] (Section 3). • ARKN5: The fifth-order ARKN method given by Franco [12] (Section 3). 4.2.1. Testing the efficiency In Fig. 1 we display the efficiency curves: accuracy versus the computational cost measured by the number of function evaluations required by each method. It is not surprising that the numerical performance of FRKN4 is superior to all the other RKN methods for Problems 1–3. For Duffing’s equation we observe that EFRKN4 and ARKN5 are the most efficient methods. For Bessel’s equation is ARKN5 the most efficient one and we remark that the efficiency curves for FRKN4 and EFRKN4 coincide. On the other hand, for the perturbed system is FRKN4 the most efficient method. It is difficult to choose the most suitable method for a given problem. However, we have showed that FRKN4 is suitable for the numerical integration of perturbed oscillators. 4.2.2. Verification of the error analysis Another goal with the numerical experiments in fixed step-size mode is to illustrate that an application of the new fourth-order RKN method to linear problems (3.19) will raise the order of accuracy by one unit with respect to the underlying classical RKN method. For this part of the numerical experiments we only consider the classical fourthorder RKN method of Dormand et al. [21] (denoted as RKN4), FRKN4 and RKN4:6. These methods use three function evaluations at each step. In Fig. 2 we display − log10 (h) versus − log10 (GE) where GE is the maximum global error over the whole integration interval when we apply the considered RKN methods to Problems 1–3. The slopes of the best fitting straight lines are showed in Table 1. Clearly, these numbers confirm that a raise of order is obtained for FRKN4 applied to Problems 1–3. This is a direct reflection of the expressions of the PLTEs (3.18) and the considerations made in Section 3.1. Beside this we have detected a raise of order by two units when we apply RKN4:6 to problems of the form (3.19). Although RKN4:6 reaches order six while FRKN4 reaches only order five for linear problems (3.19), the most accurate method is FRKN4. This can be easily understood by the fact that the construction of RKN4:6 is based on requirements that are softened: it does not exactly integrate the unperturbed problem (3.2). 4.3. Comparisons with variable step-size The codes used in the comparisons have been denoted by: • EFRKN4(3): The trigonometric version of the embedded 4(3) pair with FSAL property presented by Franco [20] (Sections 3.3 and 3.4). • FRKN4(3): The new pair derived in this paper. The FSAL technique is used. • RKN4:6(3:4): The embedded 4(3) pair proposed by García et al. [9] (Section 3). The lower order method reaches order four for the unperturbed problem (3.2).
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
139
Fig. 1. The number of function evaluations versus log10 (GE) where GE is the maximum global error over the whole integration interval when we apply the considered RKN methods to Problems 1–6.
140
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
Fig. 2. − log10 (h) versus − log10 (GE) where GE is the maximum global error over the whole integration interval when we apply RKN4, FRKN4 and RKN4:6 to Problems 1–3.
• ARKN4(3): The embedded 4(3) pair of ARKN methods derived by Franco [11] (Section 3). • ARKN5(3): The embedded 5(3) pair of ARKN methods given by Franco [12] (Section 3). The efficiency curves are drawn in Fig. 3. 5. Conclusions We have introduced an new technique for the construction of RKN methods fitted to perturbed oscillators. It can be seen as an extension of Simos’ RK approach [15]. The new RKN methods integrate exactly the unperturbed problem (3.2). Based on this technique we present a new embedded 4(3) pair of RKN methods that shares the FSAL property. We have showed that the new fourth-order RKN method achieves order five for linear problems of the form (3.19). An error analysis, a discussion and a comparison with some previously derived RKN methods (such as EFRKN methods) are made. Although, it is difficult to determine the most efficient RKN method for a given oscillatory problem, numerical experiments reveal an acceptable performance of the new RKN methods compared with some other RKN methods especially designed for oscillatory problems.
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
141
Fig. 3. The number of function evaluations versus log10 (GE) where GE is the maximum global error over the whole integration interval when we apply the considered embedded RKN methods to Problems 1–6.
142
H. Van de Vyver / Computer Physics Communications 167 (2005) 129–142
Acknowledgements The author wish to thank the anonymous referees for their careful reading of the manuscript and their valuable comments and fruitful suggestions that greatly improved this paper.
References [1] P.J. Van der Houwen, B.P. Sommeijer, Explicit Runge–Kutta(–Nyström) methods with reduced phase errors for computing oscillating solution, SIAM J. Numer. Anal. 24 (1987) 595–617. [2] P.J. Van der Houwen, B.P. Sommeijer, Phase-lag analysis of implicit Runge–Kutta methods, SIAM J. Numer. Anal. 26 (1989) 214–229. [3] P.J. Van der Houwen, B.P. Sommeijer, Diagonally implicit Runge–Kutta–Nyström methods for oscillatory problems, SIAM J. Numer. Anal. 26 (1989) 414–429. [4] B. Paternoster, A phase-fitted collocation-based Runge–Kutta(–Nyström) method, Appl. Numer. Math. 35 (2000) 339–355. [5] T.E. Simos, E. Dimas, A.B. Sideridis, A Runge–Kutta–Nyström method for the numerical integration of special second-order periodic initial-value problems, J. Comput. Appl. Math. 51 (1994) 317–326. [6] T.E. Simos, A modified Runge–Kutta method for the numerical solution of ODE’s with oscillation solutions, Appl. Math. Lett. 6 (1996) 61–66. [7] T.E. Simos, J.V. Aguiar, A modified Runge–Kutta–Nyström method with phase-lag of order infinity for the numerical solution of the Schrödinger equation and related problems, Int. J. Mod. Phys. C 11 (2000) 1195–1208. [8] A.B. González, P. Martín, J.M. Farto, A new family of Runge–Kutta type methods for the numerical integration of perturbed oscillators, Numer. Math. 82 (1999) 635–646. [9] A. García, P. Martín, A.B. González, New methods for oscillatory problems based on classical codes, Appl. Numer. Math. 42 (2002) 141–157. [10] J.M. Franco, Runge–Kutta–Nyström methods adapted to the numerical integration of perturbed oscillators, Comput. Phys. Comm. 147 (2002) 770–787. [11] J.M. Franco, Embedded pairs of explicit ARKN methods for the numerical integration of perturbed oscillators, J. Comput. Methods Sci. Engrg. 3 (2003) 415–424. [12] J.M. Franco, A 5(3) pair of explicit ARKN methods for the numerical integration of perturbed oscillators, J. Comput. Appl. Math. 161 (2003) 283–293. [13] Ch. Tsitouras, T.E. Simos, Optimized Runge–Kutta pairs for problems with oscillating solutions, J. Comput. Appl. Math. 147 (2002) 397–409. [14] Z.A. Anastassi, T.E. Simos, Special optimized Runge–Kutta methods for IVP’s with oscillating solutions, Int. J. Mod. Phys. C 15 (2004) 1–15. [15] T.E. Simos, An exponentially-fitted Runge–Kutta method for the numerical integration of initial-value problems with periodic or oscillating solutions, Comput. Phys. Comm. 115 (1998) 1–8. [16] B. Paternoster, Runge–Kutta(–Nyström) methods for ODEs with periodic solutions based in trigonometric polynomials, Appl. Numer. Math. 28 (1998) 401–412. [17] G. Vanden Berghe, H. De Meyer, M. Van Daele, T. Van Hecke, Exponentially-fitted explicit Runge–Kutta methods, Comput. Phys. Comm. 123 (1999) 7–15. [18] J.P. Coleman, S.C. Duxbury, Mixed collocation methods for y = f (x, y), J. Comput. Appl. Math. 126 (2000) 47–75. [19] T.E. Simos, Exponentially-fitted Runge–Kutta–Nyström method for the numerical solution of initial-value problems with oscillating solutions, Appl. Math. Lett. 15 (2002) 217–225. [20] J.M. Franco, Exponentially fitted explicit Runge–Kutta–Nyström methods, J. Comput. Appl. Math. 167 (2004) 1–19. [21] J.R. Dormand, M.E.A. El-Mikkawy, P.J. Prince, Families of Runge–Kutta–Nyström formulae, IMA J. Numer. Anal. 7 (1987) 235–250. [22] E. Hairer, S.P. Nørsett, G. Wanner, Solving Ordinary Differential Equations I, Nonstiff Problems, second ed., Springer-Verlag, Berlin, 1993. [23] L.F. Shampine, Local error estimation in the solution of ordinary differential equations, Math. Comp. 27 (1973) 91–97. [24] L.F. Shampine, Local error estimation by doubling, Computing 34 (1985) 179–190. [25] L.F. Shampine, Some practical Runge–Kutta formulas, Math. Comp. 46 (1986) 135–150. [26] E. Stiefel, D.G. Bettis, Stabilization of Cowell’s method, Numer. Math 13 (1969) 154–175. [27] J.M. Franco, M. Palacios, High-order P-stable multistep methods, J. Comput. Appl. Math. 30 (1990) 1–10. [28] R. Van Dooren, Stabilization of Cowell’s classical finite difference methods for numerical integration, J. Comput. Phys. 16 (1974) 186– 192.