c 1999 Society for Industrial and Applied Mathematics °
SIAM J. SCI. COMPUT. Vol. 21, No. 2, pp. 747–763
HIGH PHASE-LAG-ORDER RUNGE–KUTTA AND ¨ NYSTROM PAIRS∗ S. N. PAPAKOSTAS† AND CH. TSITOURAS† Abstract. We exploit the freedom in the selection of the free parameters of one family of eighth-algebraic-order Runge–Kutta (RK) pairs and of three families of fourth-, sixth-, and eighthorder RK Nystr¨ om (RKN) pairs with the purpose of obtaining specific pairs of the highest possible phase-lag order, which are also characterized by minimized principal truncation error coefficients. We present a method for the analytic derivation of the dissipation-order conditions for RK methods and the phase-lag- and dissipation-order conditions for Nystr¨ om methods. The RK pairs we study here are based on a one-parameter generalization of some older families of pairs. An algorithm and specific optimized 8(6) Nystr¨ om pairs are also provided. For a class of initial value problems, whose solution is known to be described by free oscillations or free oscillations of low frequency with forced oscillations of high frequency superimposed, over long integration intervals, these new pairs seem to offer some advantages with respect to some older pairs. The latter are of the same algebraic orders as the new ones but are characterized by the minimal phase-lag order according to their algebraic order and number of stages. Key words. Runge–Kutta, Nystr¨ om methods, periodic initial value problems, pairs of embedded methods, hyperbolic equations, phase-lag, dispersion order AMS subject classification. 65L05 PII. S1064827597315509
1. Introduction. In a previous work [10] (and its enhanced version [11]) the authors considered the problem of the construction of specially devoted methods for periodic initial value problems of the form (1.1)
y 0 = f (x, y) , y (x0 ) = y0 , x ∈ [x0 , xe ] , f : R × Rm → Rm ,
whenever f is supposed to be sufficiently smooth and the problem (1.1) is known to be characterized by a solution described by free oscillations or free oscillations of low frequency with forced oscillations of high frequency superimposed, over long integration intervals. Runge–Kutta (RK) methods for problems of this type were introduced by van der Houwen and Sommeijer in [8]. One method for solving (1.1) numerically, particularly for problems with a low cost in function evaluations, is by means of explicit high-order RK pairs. An s-stage RK pair is characterized by a set of coefficients in A, b, ˆb, c (A ∈ Rs×s , bT , ˆbT , c ∈ Rs ), and it is applied on problem (1.1) in the following way: Given the approximate solution yn ' y (xn ) at the point xn , an approximation yn+1 ' y (xn+1 ) is computed at xn+1 from the formulas yn+1 = yn + hn
s X
bi ki ,
i=1
where hn = xn+1 − xn , ∗ Received by the editors January 24, 1997; accepted for publication (in revised form) May 21, 1998; published electronically October 14, 1999. http://www.siam.org/journals/sisc/21-2/31550.html † National Technical University of Athens, Department of Mathematics, Zografou Campus, Athens 157 80, Greece (
[email protected],
[email protected]).
747
748
S. N. PAPAKOSTAS AND CH. TSITOURAS
k1 = f (xn , yn ) , Ã kj = f
xn + hn cj , yn + hn
s−1 X
! aji ki
,
j = 2, 3, . . . , s,
i=1
and the error quantity en+1 = hn
s X
bi − ˆbi ki
i=1
is used for the purposes of step-size control (see Shampine [16]). An RK pair is said to be of an algebraic order p(q) (p > q), whenever both its associated RK methods (characterized by A, b, c and A, ˆb, c, respectively) satisfy a set of polynomial order conditions, being in one-to-one correspondence with the set of (rooted) trees T (see, for example, Hairer, N¨ orsett, and Wanner [7] or Butcher [1]). The solution of this system of order conditions is usually carried out through the use of various types of simplifying assumptions (the most commonly used are those of the type discussed by Curtis [2] and Butcher [1] and theoretically justified in [12]). All known RK methods and pairs of orders exceeding four satisfy the simplifying assumption T
Ae = c, e ={1, . . . , 1} . | {z } s
The application of an RK method to the test problem (1.2)
y 0 = iωy,
ω ∈ R,
leads to the numerical scheme yn+1 = P (ωhn ) yn , hn = xn+1 −xn , where the function P∞ −1 P (v) = P (ωh) satisfies the relation P (v) = 1 + vb (I − vA) e = i=0 ti v i , and for i ≥ 1, ti = bAi−1 e, t0 = 1. The numbers ti depend only on the coefficients of the method. It must be observed that for explicit methods (that is, for A lower triangular), the summation in the determination of P (v) above is finite (specifically, i runs from 0 through s). The phase-lag (or dispersion) order of an RK method is defined as the order of approximation of the argument of the exponential function by the argument of P along the imaginary axis. Symbolically, the phase-lag order of a method is q, whenever ¡ δ (v) = O v q+1 , for δ (v) = v − arg (P (v)). For RK methods this notion has been introduced in [8]. The imaginary stability interval of an RK method II = (0, v0 ) is defined by the relations |P (v)| < 1 and |P (v0 + θ)| > 1 for every v ∈ II and every suitably small positive θ. A method characterized by a nonvanishing imaginary stability interval is called dissipative. Although for an RK method the phase-lag property is defined for the special problem (1.2), as it was shown by the numerical tests presented in [10], RK pairs of high phase-lag order exhibit a remarkable numerical performance on a much wider class of test problems. It seems that for a certain class of initial value problems (as those whose solutions are described by free oscillations or free oscillations of low frequency with forced oscillations of high frequency superimposed, over long integration intervals), it might be advantageous to use pairs of methods of high phase-lag order which minimize their leading truncation error coefficients, instead of pairs of the same
¨ PAIRS HIGH PHASE-LAG-ORDER RUNGE–KUTTA AND NYSTROM
749
algebraic order as the latter, but with a phase-lag order equal to the minimal allowed by the number of stages and by their algebraic order. In this article we intend to further explore these assertions by obtaining high algebraic-order, high phase-lag-order RK pairs. We also intend to study the numerical performance of RK Nystr¨ om (RKN) pairs applied on a special second-order initial value problem of the form y 00 = f (x, y) , y (x0 ) = y0 , y 0 (x0 ) = y00 , x ∈ [x0 , xe ] , f : R × Rm → Rm ,
(1.3)
which is known to be characterized by a periodic solution. Nystr¨ om pairs are characterized by the coefficients in A, b, b0 , ˆb, bˆ0 , and c, where A ∈ Rs×s and bT , b0T , ˆbT , bˆ0 T , c ∈ Rs . If yn ' y (xn ), y 0 ' y 0 (xn ) is the solution, and its derivative, comn puted at the point xn , then these pairs are applied on (1.3) in order to approximate 0 ' y 0 (xn+1 ) at xn+1 according to the formulas yn+1 ' y (xn+1 ), yn+1 yn+1 = yn +
hn yn0
0 yn+1 = yn0 + hn
+
s X
h2n
s X
bi ki ,
i=1
b0i ki ,
i=1
where k1 = f (xn , yn ) , Ã kj = f
xn + hn cj , yn +
hn yn0
+
h2n
s−1 X
! aji ki
,
j = 2, 3, . . . , s,
i=1
and some norm of the vector (eTn+1 , (e0n+1 )T )T , with en+1 = h2n
s X
bi − ˆbi ki ,
i=1
e0n+1 = hn
s X
b0i − ˆb0i ki ,
i=1
is used for feeding the step-size change control algorithm. For higher-order pairs it is advantageous to consider the following set of simplifying assumptions: c2 , 2 0 b = b (I − C) , ˆb = bˆ0 (I − C) ,
Ae =
where C = diag (c). None of these conditions is absolutely necessary (see, for example, the low algebraic-order methods of van der Houwen and Sommeijer [8]), but their application, especially for higher-order methods and pairs, greatly simplifies the system of nonlinear order conditions.
750
S. N. PAPAKOSTAS AND CH. TSITOURAS
In the next section we shall seek explicit formulas for the determination of the phase-lag and dissipation order of RK and Nystr¨ om methods. In the third section we shall describe the derivation of special eighth-algebraic-order RK pairs, which belong to a three-parameter (b13 , ˆb12 , and ˆb13 ) extension of the family proposed by Verner in [18] and studied theoretically in its full generality, leading to the only known efficiently implementable algorithms, in [12]. The type of simplifying assumptions used for the derivation in [18] (as well as in other newer ones) are based on those proposed by Fehlberg [6]. In particular, the implied derivation by Prince and Dormand in [15] (which does not contain any algorithms) seems to miss the free parameter ˆb13 . High phase-lag-order Nystr¨ om pairs of algebraic orders 4(3), 6(4), and 8(6) are derived in the fourth section. As we have found in practice (see also [10]), for the types of periodic initial value problems studied here, it seems to be more advantageous to consider pairs of the highest achievable phase-lag-order, instead of pairs of high dissipative order as well. In the last section we evaluate the numerical performance of all new RK and Nystr¨om pairs proposed by our study, with respect to other existing pairs of the same algebraic orders, which, however, were not specifically designed for periodic initial value problems. 2. Phase-lag- and dissipation-order of Runge–Kutta and Nystr¨ om methods. The dissipation order of an RK method (see van der Houwen and Sommeijer [8]) is defined as the order of approximation of the modulus of the exponential function by the modulus of the characteristic function P (v) of the along the imaginary ¡ method axis. That is, the dissipation order is q, iff α (v) = O v q+1 , for α (v) = 1 − |P (v)|. In practical situations one is interested in estimating both the phase-lag and dissipation order of a pth algebraic-order RK method. Let ti = bAi−1 e and t−1 = 0. Explicit formulas for both these quantities are offered by the following theorem. Theorem 2.1 (phase-lag- and dissipation-order conditions for RK methods). An RK method is of phase-lag order 2rp iff for every k = 1, . . . , rp , Xp (k) = 0 and Xp (rp + 1) 6= 0, where Xp (k) =
¡ k X 22(k−n+1) 22(k−n+1) − 1 B2(k−n+1) t2n−2 − t2k−1 , (2 (k − n + 1))! n=1
and B2n = B2n (0) are the Bernoulli numbers. Moreover, a method is of dissipation order rd ≥ 2 iff for every k = 2, . . . , rd , ˜ d (rd + 1) 6= 0, where ˜ d (k) = 0 and X X ˜ d (k) = X
k X ¡
t2n−2 t2(k−n) + t2n−1 t2(k−n)−1 .
n=1
Proof. The proof concerning the phase-lag-order conditions has been given in [10]. For the dissipation-order conditions one should consider the formal expansion Ã∞ !2 à ∞ !2 X X 2 i+1 i+1 2i−2 2i−1 (−1) t2i−2 v + (−1) t2i−1 v −1 |P (v)| − 1 = i=1
=
∞ X i=2
i=1 i X j=1
¡
t2j−2 t2(i−j) + t2j−1 t2(i−j)−1 v 2i−2 .
¨ PAIRS HIGH PHASE-LAG-ORDER RUNGE–KUTTA AND NYSTROM
751
If a specific method is of algebraic-order p, then one should take into account in the previous theorem that from the algebraic-order conditions it is ti = i!1 for 0 ≤ i ≤ p. Next consider the test problem (2.1)
y 00 = −λ2 y,
y (0) = 1,
y 0 (0) = iλ,
λ ∈ R,
whose exact solution is y¯ = eiλx .
(2.2)
When (2.1) is solved numerically by a Nystr¨ om method, the following recursive relation is obtained: Yn+1 = R (zn ) Yn ,
zn = −Hn2 ,
Hn = λ hn ,
T
where Yn = [yn , hn yn0 ] and R (z) =
−1
−1
1 + zb (I − zA) e 1 + zb (I − zA) c . −1 −1 zb0 (I − zA) e 1 + zb0 (I − zA) c
T y (xn ) , hn y¯0 (xn )] of (2.2) satisfies the The theoretical solution Y¯n = Y¯ (xn ) = [¯ recursive relation
¯ (Hn ) Y¯n , Y¯n+1 = R where ¯ (Hn ) = R
cos Hn −Hn sin Hn
sin Hn Hn
cos Hn
.
The phase-lag (or dispersion) order of a Nystr¨ om method (see [8]) is defined as the ¯ by the argument of the order of approximation of the arguments of the eigenvalues of R ¡ eigenvalues of R. Clearly this is expressed by the number q, where δ˜ (v) = O v q+1 p and δ˜ (v) = v − arccos(trace(R(z)/(2 det R(z)))). 0 Let, for i ≥ 1, σ2 i−1 = bAi−1 e, σ2 i = bAi−1 c, σ20 i−1 = b0 Ai−1 e, σ2i = b0 Ai−1 c, 0 0 and σ−1 = 0, σ0 = σ−1 = σ0 = 1. Applying the formal Neumann expansion of −1 (I − zA) (see Ortega [9, p. 201, relation 5.3.18]) we may formally write ∞ P σ2 i−1 z i 1 + σ 1 z + σ 3 z 2 + · · · 1 + σ2 z + σ 4 z 2 + · · · i=0 = P R (z) = ∞ σ10 z + σ30 z 2 + · · · 1 + σ20 z + σ40 z 2 + · · · σ20 i−1 z i i=0
∞ P i=0 ∞ P i=0
σ2 i z i 0 i σ2i z
.
The interval of periodicity of a Nystr¨ om method IP = (0, v0 ), applied on problem (2.1), is defined to satisfy det (R (z)) ≡ 1 and |trace (R (z))| ≤ 2, for every z ∈ IP and |trace (R (z + θ))| > 2, for some suitable small θ. Methods with an extended interval of periodicity might be especially suited for mildly stiff, oscillatory problems. Theorem 2.2 (phase-lag-order conditions for Nystr¨ om methods). A Nystr¨ om method is of phase-lag order 2rN p iff for any k = 1, . . . , rN p , X (k) = 0 and X (rN p + 1) 6= 0,
752
S. N. PAPAKOSTAS AND CH. TSITOURAS
where XN p (k) = 2πk − ρk +
k X 22n+1 πk−n , (2n)! n=0
and πk = ρk =
k P j=0 k P
0 0 σ2k−2j−1 σ2j − σ2k−2j−1 σ2 j ,
j=0
0 σ2k−2j−1 + σ2k−2j
¡ 0 σ2j−1 + σ2j .
Proof. Let Tp = tr (R) and Dp = det (R). We may write (2.3)
Tp =
∞ X
(σ2i−1 + σ2i ) z i ,
i=0
and after performing the multiplications ∞ i ∞ i−1 X X X X 0 σ2i−2j−1 σ2j 0 z i − σ2i−2j−1 σ2j z i . Dp = i=0
j=0
i=1
j=0
Since σ 0 −1 = 0, we find Dp =
(2.4)
∞ X
πi z i .
i=0
After carrying out the necessary multiplication, we may estimate from (2.3) Tp2 =
(2.5)
∞ X
ρi z i .
i=0
Substituting (2.4), (2.5) in the relation 4Dp cos2 (z) − Tp2 and expanding cos2 (z) = P∞ 22i i 1 i=0 (2i)! z ) we find 2 (1 + 4Dp cos2 (z) − Tp2 = 2
∞ X i=0
πi z i +
!Ã ∞ ! ∞ X 22i X zi − 2πi z i ρi z i , (2i)! i=0 i=0 i=0
̰ X
from which we may estimate the required equations as the multipliers in the series ∞ i 2j+1 X X 2 2πi − ρi + πi−j z i . 4Dp cos2 (z) − Tp2 = (2j)! i=0 j=0 3. Some high phase-lag-order RK pairs. The solution of the order conditions for the construction of high algebraic-order RK pairs is achieved by means of employing a set of simplifying assumptions. According to the number and type of simplifying assumptions they satisfy, RK pairs belong to certain families of solutions of the algebraic-order conditions, characterized by a set of free parameters.
1550002331 7531441519 0 0 0 0 0 0 0 0 0 0 0 0
1 39
− 819919454 7954880963
166589347 4324203682
1434700181 4526297859
216865477 4503783464
764510756 13040330925
199659293 12702785561
793314604 13976557807
1266074123 6281065917
− 3403758010 22022585129
− 657508583 7361604142
5170243423 15755825726
169842013 3932943796
522276081 5878413005
1 39
1551985603 15106987586
4655956809 30213975172
514748023 1400737386
753846329 1732179602
81892839 604421644
2 7
6192277787 12208988746
10 13
6 7
1
1
b
ˆ b
0
0
0
0
0
0
0
0
0
0
0
−17244495369 14827907888
2244683731 19421949501
0
0
62486325091 20061313835
22432995201 7608186638
− 2357986081 765762355
958376589 5714262451
− 1119897545 20078489442
873774614 42183754199
− 6187148133 7800444566
9127994833 10060356632
5335414531 10119211873
− 1974918656 10673172099
551682764 839617083
0
0
3845693529 21142204130
− 10426657685 − 1159962323 4368315947 11041458463
− 538046068 6364136977
− 1751809551 − 125019663783 − 12406223420 4639249355 25359617453 6830353731
1115105471 6543055469
49525568891 22957107026 5017941652 6447774631
5262730239 11622135734
− 2937544370 5006731581
14028594 1101021281
657452440 5173579471
− 895931653 22522139750
3736705947 22745211070
4666993327 14783147739
− 4581355363 3046479129
2985235807 7278171282
− 7555597260 16483583449
452445933 4716444208
2363387954 10609390473
14672042545 12090788322
19094602291 8082377928
4105332252 3814764721
1806141113 11369404446
77341300113 24385589366
− 9226813649 3827536113
7536718113 25936600354
19903609720 8774029877
1227998211 3278129960
37789086389 9899353391
0
− 11 18
11 18
1126082023 − 51833701 1 6379136305 5104603713 18 − 7363979475 5002757021 14423836544 8792250993
701266701 6738259450
− 5956267453 2693279968 11979576404 5123250599
967235657 − 5946035525 − 387575575 3942473888 11448823629 2463131197
− 33791676863 − 22523121910 14898822572 13774108167
25789073016 9056069051
405612119 9044847963
Table 3.1 Coefficients of NEW8(7)P, selection (10) (rational approximations accurate to 21 significant digits).
754
S. N. PAPAKOSTAS AND CH. TSITOURAS Table 3.2 The main characteristics of the RK and Nystr¨ om pairs compared here.
Type RK RKN
Pairs of alg. orders p(q) NEW8(7)P PD8(7) NEW4(3)P DEP4(3) NEW6(4)P DEP6(4) NEW8(6) NEW8(6)P DEP8(6)
Effective Phase-lag number order of stages 16 13 8 13 8 3 4 3 10 5 6 5 8 8 16 8 8 8
° ° ° ° °T (p+1) ° °T 0(p+1) ° 2
5.83 · 10−6 4.51 · 10−6 1.1 · 10−2 4.6 · 10−4 8.4 · 10−5 8.7 · 10−5 1.7 · 10−7 2.7 · 10−6 8.3 · 10−7
IR
2
7.4 · 10−3 1.8 · 10−3 6.7 · 10−5 7.7 · 10−5 1.6 · 10−7 2.6 · 10−6 8.2 · 10−7
IIM
D∞
(−6.55, 0) (0, 0) 4.9 (−5.16, 0) (1.51, 3.70) 16.7 — ∗ 1.5 — ∗ 6.1 5.61 ∗ 0.69 — ∗ 1.1 3.39 ∗ 22.3 — ∗ 12.9 3.29 ∗ 9.7
IIM : Imaginary stability interval (not applicable to Nystr¨ om methods). IR : Real°stability interval (for Nystr¨ om methods this is simply the stability interval). ° °T (p+1) ° : Euclidean norm of the vector consisting of all the principal truncation error coefficients 2 of the higher-order method of a pair (regarding the solution y). ° ° °T 0(p+1) ° : As before but regarding y0. 2 D∞ : The maximum modulus of all coefficients of a pair.
In the following we shall restrict our attention to the family of type Ia (as characterized in [12]) of algebraic order 8(7) with 13 stages. A theoretical study of this family, as well as details of specific algorithms, which may be used for the construction of specific pairs, is given in [12]. These pairs are characterized by the condition a13,12 = 0. Consequently t13 = 0. According to Theorem 2.1, this restricts the maximal phase-lag order to 16. The phase-lag-order conditions in this case result in the following set of additional conditions: (3.1)
t9 =
1 , 9!
t10 =
71 , 259459200
t11 =
1 , 43243200
t12 =
1 . 778377600
An alternative approach is to try to increase the phase-lag order of the seventh algebraic-order formula as well. We do not, however, support this strategy, because many numerical tests we have performed indicate that, in general, it is more advantageous to use the otherwise restricted parameters of the family in order to further minimize the principal truncation error coefficients of the pth-order method of the pair. The free parameters of the 8(7) family of RK pairs are c2 , c5 , c6 , c7 , c8 , c10 , c11 , a87 , b13 , ˆb12 , and ˆb13 (one more parameter than those identified in [15] and three 1 more than those initially proposed in [18]). Since t9 = 9! , we see that by default these pairs are of a tenth dissipation order. For the solution of (3.1), along with the algorithms mentioned before, we used both a symbolic algebra manipulation package and its FORTRAN version in conjunction with a numerical minimization to find those values of the free parameters ° ° routine that lead to a minimized value of °T (9) °2 . This may be achieved by making use of the only widely available (fast) algorithm for this case, based on the theoretical justification proposed in [12]. For incorporating the phase-lag-order conditions we proceeded as follows. From (3.1), using a symbolic algebra manipulation package, we express the values of c5 , c6 , c7 , and a87 in terms of the remaining free parameters, leaving as free parameters only c2 , c8 , c10 , c11 , b13 , ˆb12 , and ˆb13 . We found an optimal pair of
¨ PAIRS HIGH PHASE-LAG-ORDER RUNGE–KUTTA AND NYSTROM
755
phase-lag order 16 (NEW8(7)P) to be described by the values of the free parameters (3.2)
c2 =
1 39 ,
c5 =
c7 =
81892839 604421644 ,
c8 = 27 ,
a87 =
551682764 839617083 ,
b13 =
514748023 1400737386 ,
c6 =
753846329 1732179602 ,
c10 =
10 13 ,
ˆb12 = − 11 , 18
1 18 ,
c11 = 67 , ˆb13 =
11 18 .
The parameters of this family are presented in Table 3.1. All its major characteristics may be found in Table 3.2. The extended real stability interval was a coincidental side effect of the minimization. 4. High phase-lag-order Nystr¨ om pairs. We shall first investigate the possibility of constructing high phase-lag-order Nystr¨ om pairs for the families of algebraicorders 4(3) and 6(4) suggested and studied by Dormand, El-Mikkawy, and Prince [3], because these seem to be fairly general families of pairs based on a reasonably small number of simplifying assumptions. According to our numerical tests we decided to sacrifice the possibility of obtaining pairs with a nonvanishing interval of periodicity in favor of finding pairs with maximal phase-lag orders. All the pairs below satisfy the first step as last (FSAL) device. Table 4.1 Coefficients of NEW4(3)P (rational approximations accurate to 20 digits).
0 5280320246 8739982487 2488284716 4549257065
1 b 0
b ˆb
ˆb0
379540831 2079644486 4090178078 30784316359 1382866212 8058611935 1382866212 8058611935 1382866212 8058611935 2788682442 15739802773 1153046258 2574473725
724764885 43347833381 958895179 − 4232828403 958895179 − 4232828403 − 5103277582 8917268805 1966645306 8825920383 27711564209 4540932828
3246882200 5850906049 3246882200 5850906049 12283001905 10027502947 3 20 − 31509757256 6039640179
0 2640160123 15021458620 1 − 20 − 13
Table 4.2 Coefficients of NEW6(4)P (rational approximations accurate to 20 digits).
0 76064096 555208869 61651457 172436989 473 677 1521284172 2494038851
1 b b0 ˆb bˆ0
148104835 15781657211 83570507 15621004272 313507335 5002407628 497059253 7416116119 1104491309 17344385380 1104491309 17344385380 1104491309 17344385380 390850314 4665518297 390850314 4665518297
1008730685 17224387836 232561219 6632504445 31814195 − 4301239521 2297852298 21296988487 2297852298 21296988487 928753894 7428602053 879866760 14012015573 831255784 11424277409
2487592367 16999280447 666859859 4681708182 1270882233 4862169760 1270882233 4862169760 1088487657 2675475233 1237986347 4111942715 5386054494 11493559817
164695106 − 10269973361 1745513301 8827769149 1745513301 8827769149 2281030107 3476164510 1838896521 8824986790 5380034471 7780066871
854905921 − 6541617807 854905921 − 6541617807 5717085047 − 17062458528
1 12
1945509358 − 12470194255
− 25
1 12
576314810 3215962383
2236434251 13285504895 10010095879 36964622736
− 171339392336 7672895439 0 0 0
− 32801447959 18176875798 62469663917 6900212338 257873323 6918876884 257873323 6918876884 257873323 6918876884 479568319 9904172323 416578791 5890491199
43 46
1 1 b b0 ˆb ˆb0 0
2055300167 20115712915
1252106627 8413973203
1885846298 9184313637
1503948753 8413843957
7953080697 15591954814
1332906447 6116615518
1069201912 15512587877
2236434251 13285504895
262962495363 17824923050
31592171746 6958893399
10259024870 9108477419
17 22
0
1069201912 15512587877
− 22108842829 16963055973
− 111550006196 40089394711
− 25149249362 9340973033
− 1732991908 3477246155
8 13
1503948753 8413843957
661764535 1698709821
3451154231 7987305225
4686267579 2513053636
20699018807 16215676961
286557584 4330809711
14 37
− 3349780053 9500680168
4718243 6317185962
68790340 8728368029
326162972 − 7839732939
8943798416 − 12438207277
912003620 − 7090326959
81843218 29290841163
12794 98811
556579829 13434269006
711229321 5458138039
2215175292 16525689869
80991615 14493030629
19183400 9154053911
6397 98811
0
3783517924 4663518415
2506941335 28742592539
378512797 2226489968
980034039 25364950097
980034039 25364950097
238225641 − 2934789434
123716081 2797556961
13956454 1038655275
0
2818275669 − 3297774331 − 14061235138 8428609817
39547221 − 12689185081
1523915682 12294818705
92941557 11497613663
92941557 11497613663
260644226 13286668711
Table 4.3 Coefficients of NEW8(6). (Rational approximations accurate to 20 significant digits.)
9 20
756 S. N. PAPAKOSTAS AND CH. TSITOURAS
− 38441458363 13770476059
− 219790214441 16972291569 95607773388 19456972039
− 7932773132 4596414763 71049658140 6236153987
− 46100359290 11392487707 0 0 0
231129203 16522527624 5224366214 9902435809
− 29750076627 9320686309 9692589089 7956379280 439359615 8227388038 439359615 8227388038 439359615 8227388038 3034171983 6757641101 1526973169 13603364219
3808925972 8563443359
31153630152 38816419001
3518845173 3981361687
1 1 b b0 ˆb ˆb0 0
0
73428054989 9858275849
87422991625 38663097972
35551163 784327566
145826941 9208283543
4696262543 13668949958
7219910059 10269302009
16933304293 7508686580
10288661125 27344034528
2104854608 6975879173
− 2181985157 4065431411
− 10492266293 7445206961
7840618031 − 23104956506
4854425756 − 21792413805
4854425756 − 21792413805
− 6139426159 4480524683
333163109 14131828198
104556807 5769254975
50438778 7712600855
1098568487 5545869701
2104854608 6975879173
109492285 6831475017
195051773 7783045689
138671903 10602170394
31506868 6423626337
0
1098568487 11091739402
768256318 391757638 15932639971 64586127263 1779344314 1385261623 7284722091 26530729094
5105325445 16293370058 2932309469 5195834421
5474459843 5430985344
1422715268 235010273 2452894439 12882148881
7697231240 59303295 16805986777 6995595491
768256318 391757638 15932639971 64586127263
5105325445 16293370058
19233236849 8591554205
575023724 51984500 15170477089 7207732003
712628830 12795996721
5332359172 4592881589
− 23775370867 10053479602
6036627635 9590311458
Table 4.4 Coefficients of NEW8(6)P . (Rational approximations accurate to 20 significant digits.)
− 1511865776 8178732255
407851393 8348978584
0
3 10
¨ PAIRS HIGH PHASE-LAG-ORDER RUNGE–KUTTA AND NYSTROM
757
758
S. N. PAPAKOSTAS AND CH. TSITOURAS
These 4(3) pairs are characterized by the free parameters c2 , c3 , ˆb3 , ˆb4 , and bˆ04 , of conditions which only ˆb3 , ˆb4 , and bˆ04 remain free after imposing the two phase-lag-order ° ° for order 8 resulting from Theorem 2.2. However, the value of °T (5) °2 does not depend on the latter parameters (which for this reason were selected as in DEP4(3)). The coefficients of one of °the two ° uniquely defined pairs (namely, that corresponding to the minimal value of °T (5) °2 ) may be found in Table 4.1. Similarly, in Table 4.2 may be found a 6(4) pair with minimized truncation error coefficients of phase-lag order 10. In this case from the free parameters c2 , c3 , c4 , bˆ05 , and b06 = bˆ06 , only c4 , bˆ05 , and b06 remain free after imposing the phase-lag-order conditions. By inspecting the principal characteristics of this new pair from Table 3.2, we see that NEW6(4)P is comparable with DEP6(4) as a general purpose pair as well. For the 8(6) pairs we applied the simplifying assumptions used in [4] and we developed the explicit algorithm described in the Appendix. The following theorem may be proved along the lines of Proposition 1 and Lemma 1 of [13] and Propositions 1 and 2 of [14] concerning RK methods. Theorem 4.1. The coefficients of any Nystr¨ om pair derived by applying the algorithm in the Appendix characterize a pair of algebraic orders 8(6). We should note that although no algorithms or hints of a possible solution were presented in [4], the pair DEP8(6) of [4] is now seen to belong in the exact family studied here. We present here two pairs from this family after properly implementing the algorithm in the Appendix. First we derived the optimized general purpose pair NEW8(6) (see Table 4.3). Pairs with smaller truncation error coefficients can be produced at a c5 = 43/59, c6 = 82/95, and cost of increased values of °D∞ . Choosing c4° = 146/345, ° ° c7 = 73/74 we get °T (p+1) °2 = 2.4 × 10−8 , °T 0(p+1) °2 = 2.5 × 10−8 , and D∞ = 49.6. The latter pair is expected to perform better at stringent tolerances, say, smaller than 10−11 . For finding a sixteenth phase-lag order, 8(6) algebraic-order pair, we used a symbolic algebra manipulation package and solved the four explicitly evaluated phaselag-order conditions with respect to the four free nodes, using 35 digits of accuracy. The coefficients of the pair NEW8(6)P, which was obtained by this procedure, may be found in Table 4.4. 5. Numerical results. We compared NEW8(7) with respect to PD8(7) [15] on the same set of test problems as in [10], namely, G1, G2, G3, G4, and G5, under om pairs were error-per-step control and for tolerances 10−3 , 10−4 ,. . .,10−9 . The Nystr¨ applied on problems G2, G4, G5, but this time we used in place of G1 the problem G0 1: y 00 (x) = −25 y (x) ,
y (0) = 1,
y 0 (0) = 0,
x ∈ [0, 1000] ,
and instead of G3 we used the wave equation described in [8] with x ∈ [0, 1000] (referª enced here as G0 3). The class of problems G0 is thus defined as G0 = G0 1, G2, G0 3, G4, G5 . We have found these problems to be representative of the type of periodic problems we are concerned with here. The compilation of the results of these tests is presented here (Tables 5.1, 5.2, 5.3, 5.4, and 5.5) in the tabular format used also in [17] and [14] (we refer to the second of these references for a more detailed explanation of how exactly these numbers are derived). This format is based on the work of Enright and Pryce [5]. Studying the respective tables we see that all new pairs exhibit better performance than older pairs. Moreover, in all cases high phase-lag-order pairs seem to perform,
¨ PAIRS HIGH PHASE-LAG-ORDER RUNGE–KUTTA AND NYSTROM
759
on this set of test problems, better than the general-purpose pairs (either old or new ones). Efficiency gains tables (Tables 5.1–5.5). Unity represents 1%. Numbers have been rounded to the nearest digit. Positive numbers mean that the first method is superior. The final row gives the mean value of efficiency gain for all tolerances in a problem. The rightmost lower number is the average efficiency gain for all problems. Empty places in the tables are due to the unavailability of data for the respective tolerances.
log global error
Table 5.1 NEW8(7)P against PD8(7). −1 −2 −3 −4 −5 −6 −7 −8 −9
G1
G2
G3
G4
28 32 35 37 38
42 38 39 37
−2 −4 10 25
36 44 30 30 40
34
39
7
36
G5
7 11 10 20 32 16
26.4%
log global error
Table 5.2 NEW6(4)P against DEP6(4). −1 −2 −3 −4 −5 −6 −7 −8 −9
G1 34 37 38 35 29 23
G2 35 30 24 18 13 11
G3
33
22
19
28 19 20 19 16 13
G4 29 42 35 25 15 7
26
G5 32 33 0 −2 −5 −5 −4 7
21.4%
log global error
Table 5.3 NEW4(3)P against DEP4(3). −1 −2 −3 −4 −5 −6 −7 −8 −9
G1 24 19 18 18 18 18
G2
G3
G4
140 137 141 146 151
29 33 31 28
131 126 131 142 153
19
143
30
136
G5 −45 −48 −39 −33 −32 −30 −40
57.6%
760
S. N. PAPAKOSTAS AND CH. TSITOURAS
log global error
Table 5.4 NEW8(6) against DEP8(6). −1 −2 −3 −4 −5 −6 −7 −8 −9 −10
D1
D2
D3
−7 −2 7 16 16
−12 4 13 12 14 18
6 8 10 17 18
6
8
12
D4 2 4 7 14 21
10
D5 −2 1 5 0 19
5
8.2%
log global error
Table 5.5 NEW8(6)P against NEW8(6). −1 −2 −3 −4 −5 −6 −7 −8 −9 −10
G0 1
G2
11 32 10 8 9 10
43 52 39 28 17 18
13
33
G0 3
G4
32 30 33 36 40 44
37 23 25 27 22
36
27
G5
−54 −50 −44 −38 −35 −44
13%
6. Appendix: Algorithm for the construction of FSAL Nystr¨ om pairs of orders 8(6). Choose c4 , c5 , c6 , c7 , and bˆ09 as free parameters. Exclude those cases that could lead to indeterminate values on some of the following rational expressions. Set c8 = 1, b2 = ˆb2 = b02 = bˆ02 = 0. Estimate the remaining coefficients in the following order. Note that in some cases parallel paths may also be followed. First set d1 := (7c5 (c6 (5c7 − 3) − 3c7 + 2) − 7c6 (3c7 − 2) + 2 (7c7 − 5)) , d2 := (c5 (5c6 (2c7 − 1) − 5c7 + 3) + c6 (3 − 5c7 ) + 3c7 − 2) , and estimate 2c4 d1 − 2c5 (7c6 (3c7 − 2) − 2 (7c7 − 5)) + 4c6 (7c7 − 5) − 5 (4c7 − 3) , 2 (7c4 d2 − 7c5 (c6 (5c7 − 3) − 3c7 + 2) + 7c6 (3c7 − 2) − 2 (7c7 − 5)) c3 c2 = . 2
c3 =
Next define f1 (x) := (x (5c6 (2c7 − 1) − 5c7 + 3) + c6 (3 − 5c7 ) + 3c7 − 2) , and compute the following three weights: b03 = −
7c4 f1 (c5 ) − 7c5 (c6 (5c7 − 3) − 3c7 + 2) + 7c6 (3c7 − 2) − 2 (7c7 − 5) , 420c3 (c3 − 1) (c3 − c4 ) (c3 − c5 ) (c3 − c6 ) (c3 − c7 )
¨ PAIRS HIGH PHASE-LAG-ORDER RUNGE–KUTTA AND NYSTROM
b04 =
761
7c3 f1 (c5 ) − 7c5 (c6 (5c7 − 3) − 3c7 + 2) + 7c6 (3c7 − 2) − 2 (7c7 − 5) , 420c4 (c3 − c4 ) (c4 − 1) (c4 − c5 ) (c4 − c6 ) (c4 − c7 )
b05 = −
7c3 f1 (c4 ) − 7c4 (c6 (5c7 − 3) − 3c7 + 2) + 7c6 (3c7 − 2) − 2 (7c7 − 5) . 420c5 (c3 − c5 ) (c4 − c5 ) (c5 − 1) (c5 − c6 ) (c5 − c7 )
Define f2 (x) := 12x (x − c7 ) (x − 1) , and find b06 = −
b03 f2 (c3 ) + b04 f2 (c4 ) + b05 f2 (c5 ) − 2c7 + 1 , 12c6 (c6 − 1) (c6 − c7 )
b07 = −
6b03 c3 (c3 − 1) + 6b04 c4 (c4 − 1) + 6b05 c5 (c5 − 1) + 6b06 c6 (c6 − 1) + 1 , 6c7 (c7 − 1)
b08 = −
2b03 c3 + 2b04 c4 + 2b05 c5 + 2b06 c6 + 2b07 c7 − 1 , 2c8
b01 = 1 − b03 − b04 − b05 − b06 − b07 − b08 , b1 = b01 ,
b3 = b03 (1 − c3 ) ,
b7 = b07 (1 − c7 ) , a32 = a65 =
a75
c33
6c2
,
b4 = b04 (1 − c4 ) ,
b5 = b05 (1 − c5 ) ,
b6 = b06 (1 − c6 ) ,
b8 = b9 = 0,
a42 = −
c34 (2c3 − c4 ) , 12c2 (c2 − c3 )
a43 =
c34 (2c2 − c4 ) , 12c3 (c2 − c3 )
9 − 20c3 − 20c4 + 56c3 c4 − 12c7 + 28c3 c7 + 28c4 c7 − 84c3 c4 c7 , 10080(c3 − c5 )c5 (−c4 + c5 )(−1 + c6 )(c6 − c7 )b06
−9c5 + 20c3 c5 + 20c4 c5 − 56c3 c4 c5 + 15c6 − 32c3 c6 − 32c4 c6 + 84c3 c4 c6 −12c26 + 28c3 c26 + 28c4 c26 − 84c3 c4 c26 − 6c7 + 12c3 c7 + 12c4 c7 − 28c3 c4 c7 +12c5 c7 − 28c3 c5 c7 − 28c4 c5 c7 + 84c3 c4 c5 c7 = , 10080(c3 − c5 )c5 (−c4 + c5 )(c5 − c6 )(c6 − c7 )(−1 + c7 )b07 3 − 6c3 − 6c4 + 14c3 c4 − 6c5 + 14c3 c5 + 14c4 c5 − 42c3 c4 c5 , 5040(c3 − c6 )c6 (−c4 + c6 )(−c5 + c6 )(−1 + c7 )b07 ¡ 2a65 b06 + 2a75 b07 − b05 c25 − 2c5 + 1 =− , 2b08 ¡ 2 2a76 b07 − b06 c26 − 2c6 + 1 b07 (c7 − 1) =− , a = . 87 2b08 2b08
a76 = a85 a86
Setting m1 m2 m3 m4
= a65 (c3 − c5 )(c4 − c5 )c5 , = a75 (c3 − c5 )(c4 − c5 )c5 + a76 (c3 − c6 )(c4 − c6 )c6 , = a85 (c3 − c5 )(c4 − c5 )c5 + a86 (c3 − c6 )(c4 − c6 )c6 + a87 (c3 − c7 )(c4 − c7 )c7 , = (1 − c5 )c5 (−c3 + c5 )(−c4 + c5 )b05 + (1 − c6 )c6 (−c3 + c6 )(−c4 + c6 )b06 + (1 − c7 )c7 (−c3 + c7 )(−c4 + c7 )b07 ,
762
S. N. PAPAKOSTAS AND CH. TSITOURAS
solve the linear system
c3 c23 c33 c43 c53 0
c4 c24 c34 c44 c54 0
c5 c25 c35 c45 c55 0
c6 c26 c36 c46 c56 m1
c7 c27 c37 c47 c57 m2
1 1 1 1 1 m3
bb0 3 bb0 4 bb0 5 bb0 6 bb0 7 bb0
=
8
1 120
− bb09 − bb09 − bb09 − bb09 − bb09 1 1 + 24 c3 c4 − 60 (c3 + c4 ) − m4 bˆ09 1 2 1 3 1 4 1 5 1 6
Pi=9 for bb0i , i = 3, 4, . . . , 8; then bb01 = 1 − i=3 bb0i and ˆb1 = bˆ01 , ˆb3 = bˆ03 (1 − c3 ), ˆb4 = bˆ04 (1 − c4 ), ˆb5 = bˆ05 (1 − c5 ), ˆb6 = bˆ06 (1 − c6 ), ˆb7 = bˆ07 (1 − c7 ), ˆb8 = ˆb9 = 0. Evaluate a52 , a62 , a72 , a82 , solving the linear system 0 b06 b07 b08 b5 −b03 a32 − b04 a42 a52 b05 c5 b06 c6 b07 c7 b08 a62 −b03 c3 a32 − b04 c4 a42 0 2 b5 c5 b06 c26 b07 c27 b08 a72 = −b03 c23 a32 − b04 c24 a42 . ˆb5 ˆb6 ˆb7 ˆb8 a82 −ˆb3 a32 − ˆb4 a42 Then estimate the following coefficients: a53 = − a54 =
12a52 c2 (c2 − c3 ) + c35 (2c3 − c5 ) , 12c4 (c3 − c4 )
a63 = − a64 =
12a52 c2 (c2 − c4 ) + c35 (2c4 − c5 ) , 12c3 (c3 − c4 )
12a62 c2 (c2 − c4 ) − 12a65 c5 (c4 − c5 ) + c36 (2c4 − c6 ) , 12c3 (c3 − c4 )
12a62 c2 (c2 − c3 ) − 12a65 c5 (c3 − c5 ) + c36 (2c3 − c6 ) , 12c4 (c3 − c4 )
a73 = −
12a72 c2 (c2 − c4 ) − 12a75 c5 (c4 − c5 ) − 12a76 c6 (c4 − c6 ) + c37 (2c4 − c7 ) , 12c3 (c3 − c4 )
12a72 c2 (c2 − c3 ) − 12a75 c5 (c3 − c5 ) − 12a76 c6 (c3 − c6 ) + c37 (2c3 − c7 ) , 12c4 (c3 − c4 ) ¡ 2a43 b04 + 2a53 b05 + 2a63 b06 + 2a73 b07 − b03 c23 − 2c3 + 1 =− , 2b08 ¡ 2a54 b05 + 2a64 b06 + 2a74 b07 − b04 c24 − 2c4 + 1 =− . 2b08
a74 = a83 a84
Finally, modify bb, bb0 , using the device bb = λ1 b+(1−λ1 )bb, and bb0 = λ2 b0 +(1−λ2 )bb0 , λ1 6= 1, λ2 6= 1, in order to satisfy the error estimation criteria set by Dormand, ElMikkawy, and Prince [3]. REFERENCES [1] J. C. Butcher, The Numerical Analysis of Ordinary Differential Equations, John Wiley and Sons, Chichester, UK, 1987. [2] A. R. Curtis, High order explicit Runge-Kutta formulae, their uses and their limitations, J. Inst. Math. Appl., 16 (1975), pp. 35–55.
¨ PAIRS HIGH PHASE-LAG-ORDER RUNGE–KUTTA AND NYSTROM
763
[3] J. R. Dormand, M. E. El-Mikkawy, and P. J. Prince, Families of Runge-Kutta-Nystrom formulae, IMA J. Numer. Anal., 7 (1987), pp. 235–250. [4] J. R. Dormand, M. E. El-Mikkawy, and P. J. Prince, High order embedded Runge-KuttaNystrom formulae, IMA J. Numer. Anal., 7 (1987), pp. 423–430. [5] W. H. Enright and J. D. Pryce, Two FORTRAN packages for assessing initial value methods, ACM Trans. Math. Software, 13 (1987), pp. 1–27. [6] E. Fehlberg, Classical Fifth, Sixth, Seventh, and Eighth Order Runge-Kutta Formulas with Stepsize Control, Tech. report TR R–287, NASA, George C. Marshall Space Flight Center, Huntsville, AL, 1968. ¨ rsett, and G. Wanner, Solving Ordinary Differential Equations I, 2nd [7] E. Hairer, S. P. No ed., Springer-Verlag, Berlin, 1993. [8] P. J. van der Houwen and B. P. Sommeijer, Explicit Runge–Kutta (–Nystr¨ om) methods with reduced phase errors for computing oscillating solutions, SIAM J. Numer. Anal., 24 (1987), pp. 595–617. [9] J. Ortega, Matrix Theory, a Second Course, The University Series in Mathematics, Plenum Press, New York, 1987. [10] G. Papageorgiou, Ch. Tsitouras, and S. N. Papakostas, Runge–Kutta pairs for periodic initial value problems, Computing, 51 (1993), pp. 151–163. [11] G. Papageorgiou, Ch. Tsitouras, and S. N. Papakostas, Runge–Kutta Pairs for Periodic Initial Value Problems, Tech. report NA 93–1, Dept. Math., Nat. Tech. Univ. Athens, 1993. [12] S. N. Papakostas, On a class of families of high order Runge–Kutta methods and pairs, unpublished manuscript. [13] S. N. Papakostas and G. Papageorgiou, A family of fifth order Runge–Kutta pairs, Math. Comp., 65 (1996), pp. 1165–1181. [14] S. N. Papakostas, Ch. Tsitouras, and G. Papageorgiou, A general family of explicit Runge– Kutta pairs of orders 6(5), SIAM J. Numer. Anal., 33 (1996), pp. 917–936. [15] P. J. Prince and J. R. Dormand, High order embedded Runge–Kutta formulae, J. Comput. Appl. Math., 7 (1981), pp. 67–75. [16] L. F. Shampine, Some practical Runge–Kutta formulas, Math. Comp., 46 (1986), pp. 135–150. [17] P. Sharp, Numerical comparisons of some explicit Runge–Kutta pairs of orders 4 through 8, ACM Trans. Math. Software, 17 (1991), pp. 387–409. [18] J. H. Verner, Explicit Runge–Kutta methods with estimates of the local truncation error, SIAM J. Numer. Anal., 15 (1978), pp. 772–790.
