Continuous parallel-iterated RKN-type PC methods for nonstiff IVPs

Applied Numerical Mathematics 57 (2007) 1097–1107 www.elsevier.com/locate/apnum

Continuous parallel-iterated RKN-type PC methods for nonstiff IVPs Nguyen Huu Cong a,∗ , Nguyen Van Minh b a Faculty of Mathematics, Mechanics and Informatics, Hanoi University of Science, Vietnam b Faculty of Natural Science, Thai Nguyen University, Vietnam

Available online 15 November 2006

Abstract This paper investigates parallel predictor–corrector (PC) iteration schemes based on direct collocation Runge–Kutta–Nyström (RKN) corrector methods with continuous output formulas for solving nonstiff initial-value problems (IVPs) for systems of special second-order differential equations y (t) = f(t, y(t)). Consequently, the resulting parallel-iterated RKN-type PC methods are provided with continuous output formulas. The continuous numerical approximations are also used for predicting the stage values in the PC iteration processes. In this way, we obtain parallel PC methods with continuous output formulas and high-order predictors. Applications of the resulting parallel PC methods to a few widely-used test problems reveal that these new parallel PC methods are much more efficient when compared with the parallel-iterated RKN (PIRKN) methods and the sequential ODEX2 and DOPRIN codes from the literature. © 2006 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Runge–Kutta–Nyström methods; Predictor–corrector methods; Stability; Parallelism

1. Introduction The arrival of parallel computers influences the development of numerical methods for the numerical solution of nonstiff initial-value problems (IVPs) for the systems of special second-order, ordinary differential equations (ODEs)   y (t) = f t, y(t) , y(t0 ) = y0 , y (t0 ) = y0 , t0  t  T , (1.1) where y, f ∈ Rd . Among various numerical methods proposed so far, the most efficient methods for solving these problems are the explicit Runge–Kutta–Nyström (RKN) methods. In the literature, sequential explicit RKN methods up to order 10 can be found in e.g., [16–20,22,23]. In order to exploit the facilities of parallel computers, several class of parallel predictor–corrector (PC) methods based on RKN corrector methods have been investigated in e.g., [3–11,14,15,28,12,13]. A common challenge in the latter-mentioned papers is to reduce, for a given order of accuracy, the required number of sequential f-evaluations per step, using parallel processors. In the present paper, we investigate a particular class of parallel-iterated RKN-type PC methods based on direct collocation RKN corrector methods with continuous output formulas. The continuous numerical approximations also are used as starting stage values in the * Corresponding author. Current address: School of Graduate Studies, Vietnam National University, 144 Xuan Thuy, Cau Giay, Hanoi, Vietnam.

E-mail address: [email protected] (N.H. Cong). 0168-9274/$30.00 © 2006 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2006.10.002

1098

N.H. Cong, N. Van Minh / Applied Numerical Mathematics 57 (2007) 1097–1107

PC iteration process. In this way we obtain parallel PC methods that will be termed continuous parallel-iterated RKNtype PC methods (CPIRKN methods). Thus, we have achieved PC methods with dense output formulas and high-order predictors. As a consequence, the resulting new CPIRKN methods require few numbers of sequential f-evaluations per step in the PC iteration process. In Section 2, we shall consider RKN corrector methods with continuous output formulas (continuous RKN methods). Section 3 formulates and investigates the CPIRKN methods, where the order of accuracy, the rate of convergence and the stability property are considered. Furthermore, in Section 4, we present numerical comparisons of CPIRKN methods with traditional parallel-iterated RKN methods (PIRKN methods) and sequential numerical codes. 2. Continuous RKN methods A numerical method is inefficient, if the number of output points becomes very large (cf. [24, p. 188]). Therefore, in the literature, efficient numerical methods are often provided with a continuous output formula. For constructing CPIRKN methods with such a continuous output formula in Section 3, in this section, we consider a continuous extension of implicit RKN methods. Our starting point is an s-stage direct collocation (discrete) RKN method (see e.g., [4,12,25]) Yn,i = un + hci un + h2

s 

aij f(tn + cj h, Yn,j ),

i = 1, . . . , s,

(2.1a)

j =1

un+1 = un + hun + h2

s 

bj f(tn + cj h, Yn,j ),

(2.1b)

j =1

un+1 = un + h

s 

dj f(tn + cj h, Yn,j ).

(2.1c)

j =1

Let us consider a continuous output formula defined by un+ξ = un + hξ un + h2

s 

bj (ξ )f(tn + cj h, Yn,j ).

(2.1d)

j =1

Here, in (2.1), 0  ξ  2, un+ξ ≈ y(tn+ξ ), with tn+ξ = tn + ξ h, un+1 ≈ y(tn+1 ), un ≈ y(tn ), un+1 ≈ y (tn+1 ), un ≈ y (tn ) and h is the stepsize. Furthermore, Yn,i , i = 1, . . . , s, are the stage vector components representing numerical approximations to the exact solution values y(tn + ci h), i = 1, . . . , s at nth step. The s × s matrix A = (aij ), s-dimensional vectors b = (bj ), b(ξ ) = (bj (ξ )) and c = (cj ) are the method parameters in matrix or vector form. The method defined by (2.1) will be called the continuous RKN method. The step point and stage order of the (discrete) RKN method defined by (2.1a)–(2.1c) will be referred to as the step point and stage order of the continuous RKN method. By the collocation principle, the continuous RKN corrector method (2.1) is of step point order p and stage order r both at least equal s (see [25]). This continuous RKN method can be conveniently presented by the Butcher tableau (see e.g., [25,2]) c yn+1  yn+1 yn+ξ

A bT dT T b (ξ )

The matrix A and the vectors b and d are defined by the order conditions (see e.g., [12,25]). They can be explicitly expressed in terms of the collocation vector c as (cf. [12]) A = P R −1 ,

bT = gT R −1 ,

d = gˆ S −1

(2.2)

N.H. Cong, N. Van Minh / Applied Numerical Mathematics 57 (2007) 1097–1107

1099

where j +1   j −1   j −1  ci , R = (rij ) = j ci P = (pij ) = , S = (sij ) = ci , j +1     1 1 , gˆ = (gˆ i ) = , i, j = 1, . . . , s. g = (gi ) = i+1 i



The vector b(ξ ) in the continuous output formula (2.1d) is a vector function of ξ . It satisfies the continuity conditions b(0) = 0 and b(1) = b and will be determined by order conditions. For the fixed stepsize h, these order conditions can be derived by replacing un+ξ , un and Yn,j in (2.1d) with the exact solution values and by requiring that the residue is of order s + 2 in h. Using Taylor expansions for sufficiently smooth function y(t) in the neighbourhood of tn , we obtain the order conditions for determining b(ξ )   j +1 ξ (j ) T j −1 = 0, j = 1, . . . , s. (2.3a) − b (ξ )j c D = j +1 The conditions (2.3a) can be seen to be of the form

bT (ξ )R − gT diag ξ 2 , ξ 3 , . . . , ξ s+1 = 0T .

(2.3b)

From (2.3b) the explicit expression of the vector b(ξ ) is

bT (ξ ) = gT diag ξ 2 , ξ 3 , . . . , ξ s+1 R −1 .

(2.3c)

In view of (2.2) and (2.3c), it follows that the continuity conditions for the vector b(ξ ) are clearly verified. We have to note that the formula in (2.1b) is a special case of the continuous formula (2.1d) with ξ = 1. It is evident that if the conditions (2.3) are satisfied, then we have the local order relation   y(tn+ξ ) − un+ξ = O hs+2 . (2.4) For the global order of continuous approximation defined by (2.1d) (continuous order), we have the following theorem: Theorem 2.1. If the function f is Lipschitz continuous and if the continuous RKN corrector method (2.1) is of step point order p, then the continuous output formula defined by (2.1d) gives rise to a continuous approximation of order (continuous order) p ∗ = min{p, s + 2}. Proof. Let us consider the global error estimate (without the local assumption: un = y(tn ), un = y (tn )) y(tn+ξ ) − un+ξ = y(tn+ξ ) − un − hξ un − h2

s 

bj (ξ )f(tn + cj h, Yn,j )

j =1

= y(tn+ξ ) − y(tn ) − hξ y (tn ) − h2

s 

  bj (ξ )f tn + cj h, y(tn + cj h)

j =1

  + y(tn ) − un + hξ y (tn ) − un + h2

s 

  bj (ξ ) f tn + cj h, y(tn + cj h) − f(tn + cj h, Yn,j ) .

(2.5)

j =1

Since the function f is Lipschitz continuous, the following global order estimates hold: y(tn+ξ ) − y(tn ) − hξ y (tn ) − h2

s 

    bj (ξ )f tn + cj h, y(tn + cj h) = O hs+2 ,

j =1

    y(tn ) − un + hξ y (tn ) − un = O hp , h2

s  j =1

    bj (ξ ) f tn + cj h, y(tn + cj h) − f(tn + cj h, Yn,j ) = O hs+2 .

(2.6)

1100

N.H. Cong, N. Van Minh / Applied Numerical Mathematics 57 (2007) 1097–1107

Table 1 Values of NCDp |NCDp∗ for problem (2.7) obtained by various continuous RKN methods Methods

p

p∗

Nstp = 200

Nstp = 400

Nstp = 800

Nstp = 1600

Nstp = 3200

Cont.Radau Cont.Gauss Cont.Radau Cont.Gauss

3 4 5 6

3 4 5 5

1.5|1.4 2.9|2.4 4.0|3.4 5.7|3.2

2.4|2.3 4.1|3.6 5.5|4.9 7.5|4.7

3.3|3.2 5.3|4.7 7.0|6.4 9.3|6.3

4.2|4.2 6.5|5.9 8.5|7.9 11.0|7.8

5.1|5.1 7.7|7.1 10.0|9.5 12.8|9.3

The proof of Theorem 2.1 follows from (2.5) and (2.6). In view of Theorem 2.1, if the step point order p of the continuous RKN method (2.1) is not less than s + 2, then the continuous order p ∗ of the approximation defined by (2.1d) is equal to s + 2. 2 Example 2.1. In order to show the order p ∗ of the continuous approximation (continuous order) as stated in Theorem 2.1, we consider continuous RKN methods based on direct collocation Radau IIA and Gauss–Legendre methods (see [4,25]). These methods will be called continuous Radau IIA (Cont.Radau) and continuous Gauss–Legendre (Cont.Gauss). We restrict our consideration to the 2-stage and 3-stage methods and apply them to the nonlinear Fehlberg problem (cf. e.g., [16,17,19,20]): ⎛ −4t 2 − 2 2 2 ⎞ y1 (t)+y2 (t) d2 y(t) ⎝ ⎠ y(t), = 2 2  dt −4t 2 2 2 y1 (t)+y2 (t)

y(0) = (0, 1) , T

 y (0) = (−2 π/2, 0)T ,



π/2  t  10,

(2.7)

with highly oscillating exact solution given by y(t) = (cos(t 2 ), sin(t 2 ))T . The absolute global error of the (discrete) approximation of order p obtained at tn+1 = 9 and of the continuous approximation of order p ∗ obtained at tn+1.5 = tn + 1.5h = tn+1 + 0.5h are defined by y(tn+1 ) − yn+1 ∞ and y(tn+1.5 ) − yn+1.5 ∞ , respectively. Table 1 list the average numbers of correct decimal digits, i.e., the values defined by NCDp = − log10 y(tn+1 ) − yn+1 ∞ and by NCDp∗ = − log10 y(tn+1.5 ) − yn+1.5 ∞ . The values NCDp |NCDp∗ listed in Table 1 nicely show the theoretical orders p and p ∗ of the continuous RKN methods. 3. CPIRKN methods In this section, we consider a parallel PC iteration scheme based on the continuous RKN (corrector) methods. This iteration scheme is given by Yn,i = yn−1 + h(1 + ci )yn−1 + h2 (0)

s 

 (m)  bj (1 + ci )f tn−1 + cj h, Yn−1,j ,

i = 1, . . . , s,

(3.1a)

j =1

Yn,i = yn + hci yn + h2 (k)

s 

  aij f tn + cj h, Yk−1 n,j ,

i = 1, . . . , s, k = 1, . . . , m,

(3.1b)

j =1

yn+1 = yn + hyn + h2

s 

 (m)  bj f tn + cj h, Yn,j ,

(3.1c)

j =1

yn+1 = yn + h

s 

  dj f tn + cj h, Y(m) n,j ,

(3.1d)

j =1

yn+ξ = yn + hξ yn + h

s  j =1

 (m)  bj (ξ )f tn + cj h, Yn,j .

(3.1e)

N.H. Cong, N. Van Minh / Applied Numerical Mathematics 57 (2007) 1097–1107

1101

Regarding (3.1a) as a predictor method and (2.1) as a corrector method, we arrive at a PC method in PE(CE)m E (m) mode. Since the evaluations of f(tn−1 + cj h, Yn−1,j ), j = 1, . . . , s, are available from the preceding step, we have in fact, a PC method in P (CE)m E mode. In the PC method (3.1), the predictions (3.1a) are obtained by using continuous output formula (3.1e) from the (0) previous step. If in (3.1a) we set Yn,i = yn + hci yn , i = 1, . . . , s, the PC method (3.1a)–(3.1d) becomes the original parallel-iterated RKN methods (PIRKN methods) considered in [4,28]. Therefore, we call the method (3.1), a con(k−1) tinuous parallel-iterated RKN-type PC method (CPIRKN method). Notice that the s components f(tn + cj h, Yn,j ), j = 1, . . . , s, can be evaluated in parallel, provided that s processors are available, so that the number of sequential f-evaluations per step of length h in each processor equals s ∗ = m + 1. Theorem 3.1. If the function f is Lipschitz continuous and if the continuous RKN corrector method (2.1) has step point order p, then the CPIRKN method (3.1) has step point order q = min{p, 2m + s + 2} and gives rise to a continuous approximation of order (continuous order) q ∗ = min{p, s + 2}. Proof. The proof of this theorem is very simple. Suppose that f is Lipschitz continuous, yn = un = y(tn ) and yn = (0) un = y (tn ). Since Yn,i − Yn,i = O(hs+2 ) (see (2.4)) and each iteration raises the order of the iteration error by 2, we obtain the following order relations  2m+s+2  , i = 1, . . . , s, Yn,i − Y(m) n,i = O h un+1 − yn+1 = h2

s 

     bj f(tn + cj h, Yn,j ) − f tn + cj h, Y(m) = O h2m+s+4 , n,j

j =1

un+1 − yn+1 = h

s 

    (m)  dj f(tn + cj h, Yn,j ) − f tn + cj h, Yn,j = O h2m+s+3 .

(3.2)

j =1

Hence, for the local truncation error of the CPIRKN method (3.1), we may write      y(tn+1 ) − yn+1 = y(tn+1 ) − un+1 + [un+1 − yn+1 ] = O hp+1 + O h2m+s+4 ,      y (tn+1 ) − yn+1 = y (tn+1 ) − un+1 + [un+1 − yn+1 ] = O hp+1 + O h2m+s+3 .

(3.3)

The order relations (3.3) gives the step point order q as stated in Theorem 3.1 for the CPIRKN method. Furthermore, for the continuous order q ∗ of the continuous approximations defined by (3.1e), we may also write  y(tn+ξ ) − yn+ξ = y(tn+ξ ) − un+ξ + [un+ξ − yn+ξ ]  = y(tn+ξ ) − un+ξ + [un − yn ] + hξ [un − yn ] + h2

s 

  (m)  bj (ξ ) f(tn + cj h, Yn,j ) − f tn + cj h, Yn,j .

(3.4)

j =1

From (3.2), (3.3) and Theorem 2.1 we have the following global order relations   y(tn+ξ ) − un+ξ = O hmin{p,s+2} ,     un − yn = un − y(tn ) + y(tn ) − yn = O hmin{p,2m+s+2}     hξ [un − yn ] = hξ un − y (tn ) + hξ y (tn ) − yn = O hmin{p+1,2m+s+3} h2

s 

    (m)  bj (ξ ) f(tn + cj h, Yn,j ) − f tn + cj h, Yn,j = O h2m+s+3 .

(3.5)

j =1

The relations (3.4) and (3.5) then complete the proof of Theorem 3.1.

2

Remark. From Theorem 3.1, we see that by setting m = [(p − s − 1)/2] ([·] denoting the integer part function), we have a CPIRKN method of maximum step point order q = p (order of the corrector method) with minimum number of sequential f-evaluations per step s ∗ = [(p − s + 1)/2].

1102

N.H. Cong, N. Van Minh / Applied Numerical Mathematics 57 (2007) 1097–1107

3.1. Rate of convergence The rate of convergence of CPIRKN methods is defined by using the model test equation y  (t) = λy(t), where λ runs through the eigenvalues of the Jacobian matrix ∂f/∂y (cf. e.g., [4,6,7]). For this equation, we obtain the iteration error equation  (j −1) (j ) Yn − Yn = zA Yn − Yn , z := h2 λ, j = 1, . . . , m. (3.6) Hence, with respect to the model test equation, the convergence rate is determined by the spectral radius ρ(zA) of the iteration matrix zA. Requiring that ρ(zA) < 1, leads us to the convergence condition |z| <

1 ρ(A)

or h <

1 . ρ(∂f/∂y)ρ(A)

(3.7)

We shall call ρ(A) the convergence factor and 1/ρ(A) the convergence boundary of the CPIRKN method. One can exploit the freedom in the choice of the collocation vector c of continuous RKN correctors for minimizing the convergence factor ρ(A), or equivalently, for maximizing the convergence region denoted by Sconv and defined as

Sconv := z: |z| < 1/ρ(A) . (3.8) The convergence factors ρ(A) for the CPIRKN methods used in the numerical experiments can be found in Section 4. 3.2. Stability intervals The linear stability of the CPIRKN methods (3.1) is investigated by again using the model test equation y  (t) = λy(t), where λ is assumed to be negative real. By defining the matrix T  B = b(1 + c1 ), . . . , b(1 + cs ) , (0)

(0)

(0)

for the model test equation, we can present the starting vector Yn = (Yn,1 , . . . , Y1,s )T defined by (3.1a) in the form  Y(0) n = eyn−1 + h(e + c)yn−1 + zBYn−1 , (m)

where z := h2 λ. Applying (3.1a)–(3.1c) to the model test equation yields  (m−1) Y(m) n = eyn + chyn + zAYn  = I + zA + · · · + (zA)m−1 (eyn + chyn ) + (zA)m Y(0) n  (m) m+1 m m−1 eyn =z A BYn−1 + I + zA + · · · + (zA)  m−1  m m  chyn + z A eyn−1 + zm Am (e + c)hyn−1 + I + zA + · · · + (zA) ,

(3.9a)

yn+1 = yn + hyn + zbT Y(m) n



(m) = zm+2 bT Am BYn−1 + 1 + zbT I + zA + · · · + (zA)m−1 e yn 

 + 1 + zbT I + zA + · · · + (zA)m−1 c hyn + zm+1 bT Am eyn−1 + zm+1 bT Am (e + c)hyn−1 ,

(3.9b)

 = hyn + zdT Y(m) hyn+1 n

 (m) = zm+2 dT Am BYn−1 + zdT I + zA + · · · + (zA)m−1 eyn 

 + 1 + zdT I + zA + · · · + (zA)m−1 c hyn + zm+1 dT Am eyn−1 + zm+1 dT Am (e + c)hyn−1 .

From (3.9) we are led to the recursion ⎛ (m) ⎞ ⎛ (m) ⎞ Yn Yn−1 ⎜ yn+1 ⎟ ⎜ yn ⎟ ⎜  ⎟ ⎜ ⎟ ⎜ hyn+1 ⎟ = Mm (z) ⎜ hyn ⎟ , ⎝ ⎝ ⎠ ⎠ yn yn−1   hyn hyn−1

(3.9c)

(3.10a)

N.H. Cong, N. Van Minh / Applied Numerical Mathematics 57 (2007) 1097–1107

where Mm (z) is an (s + 4) × (s + 4) matrix defined by ⎛ m+1 m A B Pm−1 (z)e Pm−1 (z)c z ⎜ zm+2 bT Am B 1 + zbT Pm−1 (z)e 1 + zbT Pm−1 (z)c ⎜ Mm (z) = ⎜ zm+2 dT Am B zdT Pm−1 (z)e 1 + zdT Pm−1 (z)c ⎝ T 1 0 0 0 1 0T

1103

⎞ zm Am e zm Am (e + c) zm+1 bT Am e zm+1 bT Am (e + c) ⎟ ⎟ zm+1 dT Am e zm+1 dT Am (e + c) ⎟ , ⎠ 0 0 0 0 (3.10b)

where Pm−1 (z) = I + zA + · · · + (zA)m−1 . The matrix Mm (z) defined by (3.10) which determines the stability of the CPIRKN methods, will be called the amplification matrix, its spectral radius ρ(Mm (z)) the stability function. For a given number of iterations m, the stability interval denoted by (−βstab (m), 0) of the CPIRKN methods is defined as    

−βstab (m), 0 := z: ρ Mm (z) < 1, z  0 . We also call βstab (m) the stability boundary for a given m. The stability boundaries βstab (m) for the CPIRKN methods used in the numerical experiments can be found in Section 4. 4. Numerical experiments This section will report numerical results for the CPIRKN methods. We confine our considerations to the CPIRKN methods with direct collocation continuous RKN corrector methods based on symmetric collocation vector c investigated in [6,11]. The continuous s-stage RKN corrector methods (2.1) based on these symmetric collocation vectors have the orders p = p ∗ equal s + 1 or s depending on whether s is odd or even (cf. [6] and Theorem 2.1 in this paper). The symmetric collocation vectors were chosen such that the spectral radius ρ(A) of the RKN metrix A is minimized, so that the CPIRKN methods defined by (3.1) have “optimal” rate of convergence (see [6,11]). Table 2 below lists the stability boundaries of the CPIRKN methods with continuous RKN corrector methods based on symmetric collocation vectors considered in [6,11] with s = 3, 4, 5, 6 and with corresponding orders p = 4, 4, 6, 6. The associated CPIRKN methods based on s-stage, p-order continuous RKN corrector methods will be denoted by CPIRKNsp. For s = 3, 4, 5, 6 and p = 4, 4, 6, 6 we have the methods CPIRKN34, CPIRKN44, CPIRKN56 and CPIRKN66, respectively. We observe that the stability boundaries of these CPIRKN methods show a rather irregular behaviour, however they are sufficiently large for nonstiff IVPs of the form (1.1). In the following, we shall compare the above CPIRKN methods with explicit parallel RKN methods and sequential codes from the literature. For the CPIRKN methods, in the first step, we always use the trivial predictions given by Y0,i = y0 + hci y0 , (0)

i = 1, . . . , s.

The absolute error obtained at the end point of the integration interval is presented in the form 10−NCD (NCD may be interpreted as the average number of correct decimal digits). The computational efforts are measured by the values of Nseq denoting the total number of sequential f-evaluations required over the total number of integration steps denoted by Nstp . Ignoring load balancing factors and communication times between processors in parallel methods, the comparison of various methods in this section is based on Nseq and the obtained NCDs. The numerical experiments with small widely-used test problems taken from the literature below show a potential superiority in a parallel setting of the new CPIRKN methods over extant methods. This superiority will be significant in a parallel machine if the test problems Table 2 Stability boundaries βstab (m) for various CPIRKN methods Methods

CPIRKN34

CPIRKN44

CPIRKN56

CPIRKN66

βstab (1) βstab (2) βstab (3) βstab (4) βstab (5) βstab (6)

1.472 0.087 2.367 0.535 2.039 9.765

3.114 0.184 0.424 6.701 1.236 2.051

0.075 1.579 0.617 1.582 9.869 3.417

0.155 5.996 0.790 5.926 2.309 6.031

1104

N.H. Cong, N. Van Minh / Applied Numerical Mathematics 57 (2007) 1097–1107

are large enough and/or the f-evaluations are expensive (cf. e.g., [1]). In order to see the convergence behaviour of our CPIRKN methods, we follow a dynamical strategy in all PC methods for determining the number of iterations in the successive steps. It seems natural to require that the iteration error is of the same order in h as the local error of the corrector. This leads us to the stopping criterion (cf. e.g., [3,4,6–8,10])   (m) Y − Y(m−1)   TOL = Chp , (4.1) n

∞

n

where C is a problem- and method-dependent parameter, p is the step point order of the corrector method. All the computations were carried out on a 15-digit precision computer. An actual implementation on a parallel machine is a subject of further studies. 4.1. Comparison with parallel methods We shall report numerical results obtained by PIRKN methods, ones of the best parallel explicit RKN methods available in the literature proposed in [4,28], and the CPIRKN methods considered in this paper. We consider indirect PIRKN (Ind.PIRKN) methods investigated in [28] and direct PIRKN (Dir.PIRKN) methods investigated in [4]. We select a test set of three problems taken from the literature. These three problems possess exact solutions in closed form. Initial conditions are taken from the exact solutions. 4.1.1. Linear nonautonomous problem As a first numerical test, we apply the various p-order PC methods to the linear nonautonomous problem (cf. e.g., [4,6,7])   d2 y(t) −2α(t) + 1 −α(t) + 1 = y(t), 2(α(t) − 1) α(t) − 2 dt 2

α(t) = max 2 cos2 (t), sin2 (t) , 0  t  20, y (0) = (−1, 2)T ,

y(0) = (0, 0)T ,

(4.2)

y(t) = (− sin(t), 2 sin(t))T .

with exact solution The numerical results listed in Table 3 clearly show that the CPIRKN methods are much more efficient than the indirect and direct PIRKN methods of the same order. For this linear problem, all the CPIRKN methods need only about one iteration per step. Notice that because of round-off errors, we cannot expect 15 digits accuracy for the numerical results. As a consequence, Table 3 contains four empty spots in the last two lines where the numerical results were in the neighbourhood of the accuracy-limits of the computer and therefore considered as unreliable. 4.1.2. Nonlinear Fehlberg problem For the second numerical test, we apply the various p-order PC methods to the well-known nonlinear Fehlberg problem (2.7) considered in Section 2. The numerical results are reported in Table 4. These numerical results show that the CPIRKN methods are again by far superior to the indirect and direct PIRKN methods of the same order. For this nonlinear Fehlberg problem, the number of iterations m needed at each step for all CPIRKN methods is one or two. Table 3 Values of NCD/Nseq for problem (4.2) obtained by various p-order parallel PC methods PC methods

p

Nstp = 80

Nstp = 160

Nstp = 320

Nstp = 640

Nstp = 1280

C

Ind.PIRKN Dir.PIRKN

4 4

4.0/239 5.2/239

5.3/480 6.4/479

6.5/960 7.6/960

7.7/1920 8.8/1920

8.9/3840 10.0/3840

10−1 10−1

CPIRKN34 CPIRKN44

4 4

5.6/161 5.8/161

6.9/321 7.0/321

8.1/641 8.2/641

9.3/1281 9.4/1281

10.5/2561 10.6/2561

10−1 10−1

Ind.PIRKN Dir.PIRKN

6 6

7.4/360 8.0/354

9.2/721 9.9/710

11.0/1441 11.7/1420

12.8/2881 13.5/2839

14.6/5769 15.3/5678

10−3 10−3

CPIRKN56 CPIRKN66

6 6

9.8/173 10.2/162

11.7/322 11.9/322

13.8/642 13.9/642

10−3 10−3

N.H. Cong, N. Van Minh / Applied Numerical Mathematics 57 (2007) 1097–1107

1105

Table 4 Values of NCD/Nseq for problem (2.7) obtained by various p-order parallel PC methods PC methods

p

Nstp = 200

Nstp = 400

Nstp = 800

Nstp = 1600

Nstp = 3200

C

Ind.PIRKN Dir.PIRKN

4 4

1.7/728 2.4/722

2.8/1457 3.6/1445

4.0/2915 4.8/2889

5.2/5829 6.0/5778

6.5/11658 7.2/11555

102 102

CPIRKN34 CPIRKN44

4 4

3.3/523 3.3/473

4.6/1007 4.5/866

5.8/1942 5.7/1601

7.0/3713 6.9/3201

8.2/7033 8.1/6401

102 102

Ind.PIRKN Dir.PIRKN

6 6

4.0/900 5.0/896

5.8/1812 6.8/1807

7.6/3625 8.6/3615

9.4/7247 10.4/7230

11.2/14496 12.2/14458

103 103

CPIRKN56 CPIRKN66

6 6

6.5/526 6.7/468

8.3/999 8.5/878

10.1/1941 10.3/1611

11.9/3763 12.1/3202

13.7/7254 13.9/6402

103 103

Table 5 Values of NCD/Nseq for problem (4.3) obtained by various p-order parallel PC methods PC methods

p

Nstp = 100

Nstp = 200

Nstp = 400

Nstp = 800

Nstp = 1600

C

Ind.PIRKN Dir.PIRKN

4 4

2.9/229 2.8/229

3.7/600 4.9/600

4.9/1200 6.2/1200

6.1/2400 7.4/2400

7.3/4800 8.6/4800

101 101

CPIRKN34 CPIRKN44

4 4

4.6/201 3.3/201

5.8/401 4.5/401

6.9/801 5.7/801

8.1/1601 6.9/1601

9.3/3201 8.1/3201

101 101

Ind.PIRKN Dir.PIRKN

6 6

5.0/400 5.8/400

6.8/400 7.5/800

8.6/1600 9.3/1600

10.4/3200 11.1/3200

12.2/6400 12.9/6400

10−1 10−1

CPIRKN56 CPIRKN66

6 6

7.8/227 6.4/210

9.2/440 8.2/402

10.8/831 10.0/802

12.7/1602 11.7/1602

14.5/3202 13.6/3202

10−1 10−1

4.1.3. Newton’s equation of motion problem The third numerical example is the two-body gravitational problem for Newton’s equation of motion (see [27, p. 245]). d2 y1 (t) y1 (t) =−  , 2 dt ( y 2 (t) + y 2 (t))3

y2 (t) d2 y2 (t) =−  , 2 d t ( y 2 (t) + y 2 (t))3

y1 (0) = 1 − ε,

y1 (0) = 0,

1

2

y2 (0) = 0,

1

y2 (0) =

0  t  20,

2



1+ε . 1−ε

(4.3)

This problem can also√be found in [20] or in the test set of problems in [26]. The solution components are y1 (t) = cos(u(t))−ε, y2 (t) = (1 + ε)(1 − ε) sin(u(t)), where u(t) is the solution of Kepler’s equation t = u(t)−ε sin(u(t)) and ε denotes the eccentricity of the orbit. In this example, we set ε = 0.3. The numerical results for this problem are given in Table 5 and give rise to nearly the same conclusions as formulated in the two previous examples. 4.2. Comparison with sequential codes In Section 4.1, the CPIRKN methods were compared with indirect and direct PIRKN methods. In this section, we shall compare these CPIRKN methods with some of the best sequential codes currently available. We restricted the numerical experiments to the comparison of two of our 6-order CPIRKN56 and CPIRKN66 methods with two well-known sequential codes for the nonlinear Fehlberg problem (2.7), that is the codes DOPRIN and ODEX2 taken from [24]. We reproduced the best results obtained by these sequential codes given in the literature (cf. e.g., [28,9]) and added the results obtained by CPIRKN56 and CPIRKN66 methods. In spite of the fact that the results of the sequential codes are obtained using a stepsize strategy, whereas CPIRKN56 and CPIRKN66 methods are applied with fixed stepsizes, it is the CPIRKN56 and CPIRKN66 methods that are the most efficient (see Table 6).

1106

N.H. Cong, N. Van Minh / Applied Numerical Mathematics 57 (2007) 1097–1107

Table 6 Comparison with sequential codes for problem (2.7) Codes/Methods

Nstp

NCD

49 53 43 57 81

2.7 4.8 6.5 8.8 10.9

746 1122 1493 2039 2907

DOPRIN (from [24])

79 353 1208 4466

3.8 8.3 12.3 16.3

633 2825 9665 35729

CPIRKN56 (in this paper)

200 400 800 1600 3200

6.5 8.3 10.1 11.9 13.7

526 999 1941 3763 7254

CPIRKN66 (in this paper)

200 400 800 1600 3200

6.7 8.5 10.3 12.1 13.9

468 878 1611 3202 6402

ODEX2 (from [24])

Nseq

5. Concluding remarks In this paper, we proposed a new class of parallel PC methods called continuous parallel-iterated RKN-type PC methods (CPIRKN methods) based on continuous RKN corrector methods. Three numerical experiments showed that the CPIRKN methods are much superior to the well-known PIRKN methods and ODEX2 and DOPRIN codes available in the literature. The paper limits its focus to IVPs of the form y (t) = f(t, y(t)), y(t0 ) = y0 , y (t0 ) = y0 , however there has been proposed RKN methods for the more general problem y (x) = f(t, y(t), y (t)), y(t0 ) = y0 , y (t0 ) = y0 (see e.g., [21]). In a forthcoming paper, we will extend the ideas of this paper to this more general problem. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]

K. Burrage, Parallel and Sequential Methods for Ordinary Differential Equations, Clarendon Press, Oxford, 1995. J.C. Butcher, The Numerical Analysis of Ordinary Differential Equations, Runge–Kutta and General Linear Methods, Wiley, New York, 1987. N.H. Cong, An improvement for parallel-iterated Runge–Kutta–Nyström methods, Acta Math. Viet. 18 (1993) 295–308. N.H. Cong, Note on the performance of direct and indirect Runge–Kutta–Nyström methods, J. Comput. Appl. Math. 45 (1993) 347–355. N.H. Cong, Direct collocation-based two-step Runge–Kutta–Nyström methods, SEA Bull. Math. 19 (1995) 49–58. N.H. Cong, Explicit symmetric Runge–Kutta–Nyström methods for parallel computers, Comput. Math. Appl. 31 (1996) 111–122. N.H. Cong, Explicit parallel two-step Runge–Kutta–Nyström methods, Comput. Math. Appl. 32 (1996) 119–130. N.H. Cong, RKN-type parallel block PC methods with Lagrange-type predictors, Comput. Math. Appl. 35 (1998) 45–57. N.H. Cong, Explicit pseudo two-step RKN methods with stepsize control, Appl. Numer. Math. 38 (2001) 135–144. N.H. Cong, N.T. Hong Minh, Parallel block PC methods with RKN-type correctors and Adams-type predictors, Internat. J. Comput. Math. 74 (2000) 509–527. N.H. Cong, N.T. Hong Minh, Fast convergence PIRKN-type PC methods with Adams-type predictors, Internat. J. Comput. Math. 77 (2001) 373–387. N.H. Cong, N.T. Hong Minh, Parallel-iteratet pseudo two-step Runge–Kutta–Nyström methods for nonstiff second-order IVPs, Comput. Math. Appl. 44 (2002) 143–155. N.H. Cong, N.V. Minh, Improved parallel-iteratet pseudo two-step RKN methods for nonstiff problems, SEA Bull. Math., in press. N.H. Cong, K. Strehmel, R. Weiner, Runge–Kutta–Nyström-type parallel block predictor–corrector methods, Adv. Comput. Math. 10 (1999) 115–133. N.H. Cong, K. Strehmel, R. Weiner, A general class of explicit pseudo two-step RKN methods on parallel computers, Comput. Math. Appl. 38 (1999) 17–30.

N.H. Cong, N. Van Minh / Applied Numerical Mathematics 57 (2007) 1097–1107

1107

[16] E. Fehlberg, Klassische Runge–Kutta–Nyström Formeln mit Schrittweiten-Kontrolle für Differentialgleichungen x  = f (t, x), Computing 10 (1972) 305–315. [17] E. Fehlberg, Eine Runge–Kutta–Nyström Formel 9-ter Ordnung mit Schrittweitenkontrolle für Differentialgleichungen x  = f (t, x), Z. Angew. Math. Mech. 61 (1981) 477–485. [18] E. Fehlberg, S. Filippi, J. Gräf, Eine Runge–Kutta–Nyström Formelpaar der Ordnung 10(11) für Differentialgleichungen y  = f (t, y), Z. Angew. Math. Mech. 66 (1986) 265–270. [19] S. Filippi, J. Gräf, Ein Runge–Kutta–Nyström Formelpaar der Ordnung 11(12) für Differentialgleichungen der Form y  = f (t, y), Computing 34 (1985) 271–282. [20] S. Filippi, J. Gräf, New Runge–Kutta–Nyström formula-pairs of order 8(7), 9(8), 10(9) and 11(10) for differential equations of the form y  = f (t, y), J. Comput. Appl. Math. 14 (1986) 361–370. [21] J.M. Fine, Low order Runge–Kutta–Nyström methods, Computing 38 (1987) 281–297. [22] E. Hairer, Méthodes de Nyström pour l’équation differentielle y  (t) = f (t, y), Numer. Math. 27 (1977) 283–300. [23] E. Hairer, A one-step method of order 10 for y  (t) = f (t, y), IMA J. Numer. Anal. 2 (1982) 83–94. [24] E. Hairer, S.P. Nørsett, G. Wanner, Solving Ordinary Differential Equations, I. Nonstiff Problems, second revised ed., Springer, Berlin, 1993. [25] P.J. van der Houwen, B.P. Sommeijer, N.H. Cong, Stability of collocation-based Runge–Kutta–Nyström methods, BIT 31 (1991) 469–481. [26] T.E. Hull, W.H. Enright, B.M. Fellen, A.E. Sedgwick, Comparing numerical methods for ordinary differential equations, SIAM J. Numer. Anal. 9 (1972) 603–637. [27] L.F. Shampine, M.K. Gordon, Computer Solution of Ordinary Differential Equations, the Initial Value Problems, W.H. Freeman and Company, San Francisco, 1975. [28] B.P. Sommeijer, Explicit, high-order Runge–Kutta–Nyström methods for parallel computers, Appl. Numer. Math. 13 (1993) 221–240.