- Email: [email protected]
Construction of high-order quadratically stable second-derivative general linear methods for the numerical integration of stiff ODEs
Journal of Computational and Applied Mathematics 303 (2016) 218–228
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Construction of high-order quadratically stable second-derivative general linear methods for the numerical integration of stiff ODEs A. Abdi Faculty of Mathematical Sciences, University of Tabriz, Tabriz, Iran
article
info
Article history: Received 29 September 2015 Received in revised form 28 February 2016 Keywords: Stiff differential equations General linear methods Second derivative methods A- and L- stability Quadratic stability
abstract Theory of general linear methods (GLMs) for the numerical solution of autonomous system of ordinary differential equations of the form y′ = f (y) is extended to include the second derivative y′′ = g (y) := f ′ (y)f (y). This extension of GLMs is called second derivative general linear methods (SGLMs). In this paper we will construct two-stage A- and L-stable SGLMs of order p and stage order q = p up to six with low error constants. We will show the efficiency of the proposed methods by implementing on some well-known stiff problems. © 2016 Elsevier B.V. All rights reserved.
1. Introduction The characterization and discovery of methods suitable for the approximate numerical integration of stiff autonomous differential equations y′ (x) = f (y(x)), y(x0 ) = y0 ,
x ∈ [x0 , x],
(1)
where f : Rm → Rm and m is the dimensionality of the system, has attracted the interest of many researchers. Ideally such methods should be A-stable. But the order of an A-stable linear multistep method cannot exceed two (Dahlquist’s second barrier [1]). On the other hand, although there exist A-stable Runge–Kutta methods of arbitrarily high order, however the order of a diagonally implicit Runge–Kutta (DIRK) method with s stages and even an s-parallel SDIRK (singly diagonally implicit Runge–Kutta) method cannot exceed s + 1 [2] (see also [3]). Consequently research interest has been focused on the exploration of other paths. The idea of combining traditional one-step methods with traditional multistep methods is very natural and has been followed by many people over the last 50 years. Examples are given by the modified multistep methods [4–6]. All these methods belong to the class of general linear methods (GLMs). GLMs were introduced by Butcher [7] (see also [8]) and further investigated, for instance, in [9–19]. On the other hand, since an increase in the order of accuracy can be achieved by using the higher derivatives of the solution, especially second derivative, several authors have introduced methods for the solution of stiff systems which incorporate the higher derivatives in the integration formulas (e.g., [20–27]). For system (1), let g (y) := y′′ = f ′ (y)f (y). For problems in which g can be calculated along with f , at a moderate additional cost, second derivative methods become feasible. Although the mentioned GLMs include linear multistep methods, Runge Kutta and many other standard methods, but they were extended to second derivative general linear methods (SGLMs) to cover second derivative methods, too. These methods were introduced by Butcher and Hojjati in [28] and E-mail address: [email protected]. http://dx.doi.org/10.1016/j.cam.2016.02.054 0377-0427/© 2016 Elsevier B.V. All rights reserved.
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
219
investigated more in [29–33]. Abdi and Hojjati have divided SGLMs into four types [29,30], together with their intended applications (for nonstiff or stiff ODEs) and architectures (sequential or parallel). In [30,32], they have also obtained some upper bounds on the order of SGLMs of different types. These barriers on the order were found under the assumption of the Runge–Kutta stability (RKS) property (see [8,19]). Also, order arrows are used to obtain simpler and more illuminating proofs for these barriers in [34]. References are also made to [35,10] which include a discussion of order arrows. Construction and implementation of Nordsieck SGLMs for the numerical integration of stiff systems of ODEs are discussed by Abdi and Hojjati in [31]. In [36], to find more efficient methods with extensive absolute stability region for stiff systems, considering modified extended second derivative multistep methods [27] (MESDMMs) in SGLMs structure, perturbations of these methods were introduced which preserve their orders and improve their stability properties. The organization of the paper is as follows. After a short review of SGLMs, especially order conditions and types of these methods in Section 2, construction of A- and L-stable SGLMs with quadratic stability (QS) property of order p and stage order q = p up to six with low error constants is given in Section 3. In Section 4 some numerical experiments are presented to confirm the efficiency of the constructed methods. Finally, in Section 5 some concluding remarks are given and plans for future research are outlined. 2. A short introduction to SGLMs SGLMs for the numerical solution (1) are multivalue and multistage methods of the form
r s s [n] [n−1] [n] [n] aij g (Yj ) + uij yj , aij f (Yj ) + h2 Yi = h
i = 1, 2, . . . , s,
s s r [n] [n] [n] 2 ) + h b f ( Y = h b g ( Y ) + vij yj[n−1] , y ij ij j j i
i = 1, 2, . . . , r ,
j=1
j =1
j=1
j =1
j =1
(2)
j =1
where n = 1, 2, . . . , N , Nh = x − x0 . Here, introducing xn := x0 + nh, the internal stages Yi , i = 1, 2, . . . , s, are approximations of stage order q to y(xn−1 + ci h), i.e., [n]
= y(xn−1 + ci h) + O (hq+1 ),
[n]
Yi
and the external stages yi , i = 1, 2, . . . , r, are approximations of order p to the linear combinations of scaled derivatives of y(xn ), i.e., [n]
[n]
yi
= αi0 y(xn ) + αi1 hy′ (xn ) + αi2 h2 y′′ (xn ) + · · · + αip hp y(p) (xn ) + O (hp+1 ),
for some real numbers αi1 , αi2 , . . . , αip (cf. [29,32]). The vector y[n] denotes the computed quantities at step n, which also represent the incoming values for step n + 1. These methods are characterized by the abscissa vector c = [c1 c2 · · · cs ]T , the coefficient matrices A = [aij ] ∈ Rs×s , B = [bij ] ∈ Rr ×s ,
A = [aij ] ∈ Rs×s , B = [bij ] ∈ Rr ×s ,
U = [uij ] ∈ Rs×r , V = [vij ] ∈ Rr ×r ,
and the four integers: p, the order of the method; q, the stage order; r, the number of external stages; and s, the number of internal stages. They can be represented by a partitioned (s + r ) × (2s + r ) matrix of the form
A
A
U
B
B
V
.
It follows from [32] that the SGLM (2) has order p and stage order q = p if and only if exp(cz ) = zA exp(cz ) + z 2 A exp(cz ) + UWZ + O(z p+1 ), exp(z )WZ = zB exp(cz ) + z B exp(cz ) + VWZ + O(z 2
p+1
).
Here, exp(cz ) = [exp(c1 z ) exp(c2 z ) · · · exp(cs z )]T , Z = [1 z · · · z p ]T and
α10 α20 W = .. . αr0
α11 α21 .. . αr1
··· ··· .. . ···
α1p α2p .. . . αrp
(3) (4)
220
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
Expanding exp(cz ) and exp(z ) in (3) and (4) into power series around z = 0 and comparing the coefficients of z k , k = 0, 1, . . . , p, in the resulting expressions, the stage order (3) and order (4) can be reformulated in the form
e − U α0 = 0, α0 − V α0 = 0, c − Ae − U α1 = 0, α1 + α0 − Be − V α1 = 0, k−2 k−1 ck Ac Ac − − − U αk = 0, k ! ( k − 1 )! ( k − 2)! k Bc k−1 Bc k−2 αk−l − − − V αk = 0, l! (k − 1)! (k − 2)! l=0
(5)
for k = 2, 3, . . . , p and e = [1 1 · · · 1] ∈ Rs . The first and second equations of (5) are called pre-consistency and consistency conditions, respectively [30]. The most convenient choice for U, in the case of r = s, is U = I, where I is the identity matrix of appropriate dimension. This choice is valid since a transformation can be found to force U = I while keeping the structure of the other matrices unchanged. The matrix W is uniquely determined by the coefficients of the method if r = s and U = I. In this case it follows from order conditions (3) that [32] W = C − ACK − ACK 2 , s×(p+1)
where C = (Cij ) ∈ R Cij =
j −1 ci
(j − 1)!
,
(6)
is the Vandermonde matrix with components
1 ≤ i ≤ s, 1 ≤ j ≤ p + 1 ,
and K ∈ R(p+1)×(p+1) is the shifting matrix defined by K = [0 e1 · · · ep ] with ej as the jth unit vector. It was proved in [29] that the method (2) with p = q = r = s, U = I and Ve = e has order p and stage order p if and only if B = B0 − AB1 − AB2 − VB3 − (B − V A)B4 + VA, where
1+c
φj (x)dx B0 = φj (cj ) c s i φj (x)dx 0 B3 = φj (cj ) i
s
,
0
B1 =
i,j=1
,
B4 =
i ,j = 1
φj (1 + ci ) φj (cj )
φj′ (ci )
s
φj (cj )
s
,
B2 =
i,j=1
φj′ (1 + ci ) φj (cj )
s
,
i,j=1
,
i,j=1
and s
φi (x) =
(x − cj ),
i = 1, 2, . . . , s.
j=1,j̸=i
In the recent papers on SGLMs, these methods were investigated under the assumption that the coefficient matrices A and A have the form
λ
a21 A= .. .
λ .. .
as1
as2
,
..
. ···
µ
a21 A= .. .
µ .. .
as1
as2
λ
,
..
. ···
µ
with λ ≥ 0 and µ ≤ 0, where the explicit methods corresponding to λ = µ = 0 are called type 1 methods and the implicit formulas corresponding to λ > 0, µ < 0 are called type 2 methods. Type 3 and type 4 methods for which A and A have the form A = diag {λ, λ, . . . , λ} , A = diag {µ, µ, . . . , µ} , with λ = µ = 0 and λ > 0, µ < 0, respectively. Order barriers for different types of these methods under the assumption of the Runge–Kutta stability (RKS) property are obtained in [30,32]. Examples of many these methods of various types, which are appropriate for nonstiff or stiff differential systems in a sequential or parallel computing environment, are given in [29,30,32,28]. Also, as explained in these papers the stability properties of the method (2) are determined by the stability matrix M (z ) = V + (zB + z 2 B)(I − zA − z 2 A)−1 U , and the corresponding stability function p(w, z ) = det(w I − M (z )), where w and z are complex numbers.
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
221
3. Construction of two stage methods Since we are interested in developing a code suitable for solving stiff problems, A-stability should be necessary. In this section, we construct type 2 methods with r = s = 2 of order and stage order 3 up to 6 which are A-stable. It will be assumed that V is a rank-one matrix of the form
V =
v1 v1
v2 , v2
with v1 + v2 = 1. This choice implies that V is power bounded. Consequently, these methods are zero-stable. Moreover, it will also be assumed that U = I. Coefficients partitioned matrix for these methods is given by
A
A
U
B
B
V
λ
0
µ
0
1
0
λ
µ
0
b12
1−v
1
v
b22
1−v
v
a 21 = b 11
b12
a21 b11
b21
b22
b21
with λ > 0 and µ < 0. Stability properties of these methods are determined by quadratic stability functions, i.e., the stability function p(w, z ) of such methods has the form p(w, z ) = (1 − λz − µz 2 )2 w 2 − p1 (z )w + p2 (z ),
(7)
where p1 (z ) and p2 (z ) are polynomials of degree less than or equal to 2s = 4. For these methods to be A-stable, the Schur criterion [37] (see also [38]) can be used. This states that, a quadratic polynomial d0 w 2 + d1 w + d2 has both its roots in the open unit disc if and only if
|d0 |2 − |d2 |2 > 0,
(8)
(|d0 |2 − |d2 |2 )2 − |d0 d1 − d2 d1 |2 > 0.
(9)
and
In the stability function (7), p1 (z ) = (1 − λz − µz 2 )2 tr(M (z )), and p2 (z ) = (1 − λz − µz 2 )2 det(M (z )), are polynomials of degree at most four. So, assuming A-stability of such methods, the L-stability requirement is equivalent to lim
z →∞
p1 ( z )
(1 − λz − µ )
z2 2
= 0,
lim
z →∞
p2 (z )
(1 − λz − µz 2 )2
= 0,
which imply that the coefficients of z 4 in the numerator of tr(M (z )) and det(M (z )) must be zero. By simple algebraic calculations, these lead to
µ2 − tr(B)µ + a21 b12 det(B − µV ) + a21 v(b12 − b22 )
= =
0, 0.
(10)
Denoting the error constant of the method of order p by Ep (̸= 0), the stability function p(w, z ) in (7) satisfies p(exp(z ), z ) = Ep z p+1 + O (z p+2 ).
(11)
In the construction of the methods, we will use the expressions ‘‘L-conditions’’ and ‘‘E-condition’’ to mean satisfying (10) and (11), respectively. 3.1. Methods of order 3 In this subsection, we construct L-stable methods of order p = q = 2s − 1 = 3, r = s and abscissa vector c = [0 1]T with the error constant E3 = 10−5 . Solving stage order and order conditions (3) and (4), L-conditions and E-condition, we obtain √ a four-parameter family of methods of order and stage order p = q = 3 depending on a21 , a21 , λ, and µ with v = (α ± β)/γ where α, β and γ are complicated expressions in λ, µ and a21 . Then we apply the Schur criterion, i.e. conditions (8) and (9), to find A-stable (and so L-stable) methods in the parameter space (λ, µ) for selected values of the parameters a21 and a21 and β ≥ 0. Restricting
222
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
Fig. 1. L-stability regions in the (λ, µ)-plane for quadratically stable SGLMs with s = 2, p = q = 3 and the error constant E3 = 10−5 for specific values of the parameters.
λ and µ to 0 < λ ≤ 2 and −2 ≤ µ < 0, pairs of (λ, µ) values giving L-stability are presented in Fig. 1. The four cases in Fig. 1 correspond to: √ α+ β . γ √ α− β v= γ . √ v = α+γ β . √ v = α−γ β .
• Case 1: a21 = 12 , a21 = − 12 , v = • Case 2: a21 = 12 , a21 = − 12 , • Case 3: a21 = 2, a21 = −2, • Case 4: a21 = 2, a21 = −2,
We present here just a single example, from family of Case 1 methods, characterized by (λ, µ) = ( 34 , − 51 ). The coefficients for the method are
3 4
0
− 51
0
1
0
1 2
3 4
− 21
− 51
0
−801+940 v−475 v 2
21−250 v+325 v 2 60 d
89+127 v−120 v 2 30 d
−1+29 v−40 v 2
60 d
1−v
−96+103 v−95 v 2
−138+65 v+65 v 2
−8−6 v−40 v 2
12 d
12 d
83+11 v−60 v 2 15 d
, v
where d = −13 + 5v and v =
15 d 15 d
1−v
1
(12)
v
√ 104509+ 42225283081 . 360000
3.2. Methods of order 4 In this subsection, we construct L-stable methods of order p = q = 2s = 4, r = s and abscissa vector c = [0 1]T with the error constant E4 = 10−5 . Solving stage order and order conditions (3) and (4), L-conditions and E-condition, we√obtain a two-parameter family of methods of order and stage order p = q = 4 depending on a21 and µ with λ = (α ± β)/(1 + 12µ) where α and β are expressions in µ. Then we apply the Schur criterion to find A-stable (and so L-stable) methods in the parameter space (a21 , µ) for β ≥ 0. Restricting a21 and µ to −1 ≤ a21 ≤ 3 and −1 ≤ µ < 0, pairs of (a21 , µ) values giving L-stability are presented in Fig. 2.
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
223
Fig. 2. L-stability region in the (a21 , µ)-plane for quadratically stable SGLMs with s = 2, p = q = 4 and the error constant E4 = 10−5 .
Fig. 3. A-stability region in the (λ, c1 )-plane for quadratically stable SGLMs with s = 2, p = q = 5 and the error constant E5 = 10−5 for specific values of the parameter µ.
√
We present here just a single example, characterized by λ = (α+ β)/(1+12µ) and (a21 , µ) = ( 14 , − 14 ). The coefficients for the method are
λ
0
− 14
0
1
0
1 4
λ
−205+486 λ−288 λ2
− 14
0
1
−83−156 λ+288 λ2
−9(7−22 λ+16 λ2 )
−18953+64500 λ−73152 λ2 +27648 λ3
16 d2
8 d2
72 d1 d2
3 d1 8 d2
9 d1 4 d2
4 d2 −9 d1 4 d2
−443+676 λ−192 λ2
87−182 λ+96 λ2 8 d2
−14333+48636 λ−55008 λ2 +20736 λ3
1103−2520 λ+1440 λ2 72 d1 d2
9 d1 4 d2
4 d2 −9 d1 4 d2
16 d2
where d1 = 8λ − 7, d2 = 6λ − 5 and λ =
9 d1
72 d1 d2
(13)
√ 567509− 807465081 . 600000
3.3. Methods of order 5 In this subsection, we construct A-stable methods of order p = q = 2s + 1 = 5, r = s and abscissa vector c = [c1 1]T with c1 as a free parameter and the error constant E5 = 10−5 . Solving stage order and order conditions (3) and (4), we obtain a four-parameter family of methods of order and stage order p = q = 5 depending on λ, µ, and c1 ̸= 0, 1. Then we apply the Schur criterion to find A-stable methods in the parameter space (λ, c1 ) for selected values of the parameter µ. Restricting λ and c1 to 0 < λ ≤ 2 and c1 ∈ [−1, 1) \ {0}, pairs of (λ, c1 ) values giving A-stability are presented in Fig. 3 for selected values of the parameter µ. 3 187 We present here just a single example characterized by µ = − 20 and (λ, c1 ) = ( 250 , − 12 ). The coefficients for the method are
187 250
0
3 − 20
0
1
0
98481 62000
187 250
10713 12400
3 − 20
0
1
55804 50625
31376 − 253125
5381 − 16875
31 675
2521 2025
496 − 2025
.
41448379 25110000
− 11938687 15693750
− 1265651 8370000
38813 − 334800
2521 2025
496 − 2025
(14)
224
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
Fig. 4. A-stability region in the (λ, c1 )-plane for quadratically stable SGLMs with s = 2, p = q = 6.
3.4. Methods of order 6 In this subsection, we construct A-stable methods of order p = q = 2s + 2 = 6, r = s and abscissa vector c = [c1 1]T with c1 as a free parameter. Solving stage order and order conditions (3) and (4), we obtain a two-parameter family of methods of order and stage order p = q = 6 depending on λ and c1 ̸= 0, 1. Then we apply the Schur criterion to find A-stable methods in the parameter space (λ, c1 ). Restricting λ and c1 to 0 < λ ≤ 2 and c1 ∈ [−1, 1) \ {0}, pairs of (λ, c1 ) values giving A-stability are presented in Fig. 4. By searching for the smallest absolute value of error constant within these A-acceptable values for (λ, c1 ), we obtain 179 , −1) and the error constant E6 ≈ −0.899 × 10−3 . (λ, c1 ) = ( 250 179 Setting (λ, c1 ) = ( 250 , −1), the coefficients matrices of the method take the following forms
179 250
0
− 16
0
1
0
502 375
179 250
139 125
− 61
0
1
666071 600000
627 200000
64011 − 200000
37 8000
937 800
− 137 800
.
383357 200000
85319 − 600000
13611 − 200000
257 1600
937 800
− 137 800
(15)
4. Numerical experiments In this section we present some numerical results which aim to indicate accuracy and efficiency of the proposed methods in Section 3 in the integration of stiff differential systems. We use the following nonlinear stiff problems: P1. The non-linear stiff system of ODEs
y′1 (x) = −1002y1 (x) + 1000y22 (x), y′2 (x) = y1 (x) − y2 (x)(1 + y2 (x)),
y1 (0) = 1, y2 (0) = 1,
where the exact solution is [y1 (x), y2 (x)]T = [exp(−2x), exp(−x)]T and x ∈ [0, X ]. P2. The non-linear stiff system of ODEs
y′1 (x) = −10004y1 (x) + 10000y42 (x), y′2 (x) = y1 (x) − y2 (x)(1 + y32 (x)),
y1 (0) = 1, y2 (0) = 1,
where x ∈ [0, X ] and the exact solution is y1 (x) = exp(−4 x) and y2 (x) = exp(−x). P3. The non-linear stiff system of ODEs
′ 3 6 3 6 y1 (x) = −1000y1 y2 − cos (x)sin (x) − sin(x), y′2 (x) = −1000 y52 y43 − sin9 (x) + cos(x) ′ y3 (x) = −1000 y21 y33 − cos2 (x) sin3 (x) + cos(x)
y1 (0) = 1, y2 (0) = 0, y3 (0) = 0,
where the exact solution is [y1 (x), y2 (x), y3 (x)]T = [cos(x), sin(x), sin(x)]T and x ∈ [0, 1].
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
225
Table 1 The global error for SGLMs (12) and (13) at the end of the interval of integration [0, 1] applied to the problem P1. k
Method of order 3
Method of order 4
∥eh (X )∥ 1.02 × 10−9 6.91 × 10−11 6.89 × 10−12 7.72 × 10−13
5 6 7 8
p
∥eh (X )∥
p
3.88 3.27 3.16
2.94 × 10−10 9.49 × 10−12 3.09 × 10−13 1.97 × 10−14
4.95 4.94 3.97
P4. A system of differential equations, called CUSP, resulting from discretization of the diffusion terms by the method of lines of the periodic boundary-value problem [26]
1 ∂ 2y ∂y = − (y3 + ay + b) + σ 2 , ∂t ε ∂x ∂a ∂ 2a = b + 0.07v + σ 2 , ∂t ∂x ∂ b ∂ 2b 2 = (1 − a )b − a − 0.4y + 0.035v + σ 2 , ∂t ∂x where
v=
u u + 0.1
,
u = (y − 0.7)(y − 1.3).
This problem takes the form
y˙ i = −ε −1 (y3i + ai yi + bi ) + D(yi−1 − 2yi + yi+1 ), a˙ i = bi + 0.07vi + D(ai−1 − 2ai + ai+1 ), i = 1, 2, . . . , Nx ˙ bi = (1 − a2i )bi − ai − 0.4yi + 0.035vi + D(bi−1 − 2bi + bi+1 ), where
vi =
ui ui + 0.1
,
ui = (yi − 0.7)(yi − 1.3),
D = σ Nx2 ,
with periodic boundary conditions y0 := yNx , yNx +1 := y1 ,
a0 := aNx , aNx +1 := a1 ,
b0 := bNx , bNx +1 := b1 .
We take σ = 1/144, ε = 10−4 , Nx = 32 and the initial values as yi (0) = 0,
ai (0) = −2 cos
2iπ Nx
,
bi (0) = 2 sin
2iπ , Nx
i = 1, 2, . . . , Nx ,
with tout = 1.1 [26]. The implementation of the proposed methods in Section 3, needs a suitable starting procedure that gives the initial input vector y[0] . Here, similar to the introduced starting procedure in [31], we use the Runge–Kutta methods of suitable order, which give sufficient output information to provide an approximation to y[0] . Also, since for the all proposed methods we [n] have cs = c2 = 1, finishing procedure does not need and we use Ys[n] = Y2 as an approximation of order q = p to y(xn ). We have applied the methods (12) and (13) of orders 3 and 4, respectively, to the problem P1 with X = 1 and the methods (14) and (15) of orders 5 and 6, respectively, to the problem P2 with X = 2. We have implemented these methods with fixed stepsize h = X /2k , in correspondence to several integer values of k. The results of numerical experiments are presented in Tables 1 and 2. In these tables, we have listed norm of error ∥eh (X )∥ at the endpoint of integration X and numerical estimate to the order of convergence, p, computed by the formula p = log2 (∥eh (X )∥/∥eh/2 (X )∥), where eh (X ) and eh/2 (X ) are errors corresponding to stepsizes h and h/2. We can observe that the proposed methods are convergent with expected order p = 3, 4, 5 and 6. To illustrate the efficiency of the proposed methods in this paper, we have applied the Runge–Kutta–Gauss (RK–Gauss) methods of orders p = 4 and p = 6 [8] and the constructed SGLMs in Section 3 of the same orders to the problem P3. In the implementation of implicit methods, such as RK–Gauss methods, to approximate stage values, the Jacobian matrix ∂ f /∂ y must be found. So, the function g for the second derivative methods can be calculated by (∂ f /∂ y)f without any additional computational cost. Also, in order to save computational effort, we use (∂ f /∂ y)2 (cf. [20,24]) as an approximation
226
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
Table 2 The global error for SGLMs (14) and (15) at the end of the interval of integration [0, 2] applied to the problem P2. k
Method of order 5
∥eh (X )∥ 3 4 5 6
9.13 × 10−11 6.30 × 10−12 1.68 × 10−13 5.72 × 10−15
Method of order 6 p
∥eh (X )∥
p
3.86 5.23 4.88
1.29 × 10−9 1.98 × 10−11 3.02 × 10−13 6.33 × 10−15
6.03 6.03 5.58
Fig. 5. Comparison of the number of function evaluations for the methods of order 4 (top) and order 6 (bottom) applied to the problem P3.
to the Jacobian matrix ∂ g /∂ y. This approximation is a piecewise constant approximation to the Jacobian matrix. In Fig. 5, we compare the number of function evaluations, nfe, versus accuracy of the methods. This figure confirms the efficiency of the proposed SGLMs. Also, to compare the performance of the new methods and SGLMs previously constructed, we have applied the constructed methods in [31] and in the present paper of orders p = 3 and p = 4 to the problem P4. The number of function evaluations versus the error of the methods are plotted in Fig. 6 which illustrates superiority of the new methods from computational cost point of view. In this figure, error is the difference between the solution reported in [39] as a reference solution and the numerical results of the applied SGLMs at the endpoint of integration tout = 1.1.
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
227
Fig. 6. Performance comparison of the methods (12) and (13) with the SGLMs constructed in [31] of orders 3 and 4.
5. Conclusion The construction of the class of second derivative diagonally implicit multistage integration methods (SDIMSIMs), as a subclass of SGLMs, with p = q = r = s and RKS property has been described in [29], and some examples of these methods of all four types up to the order p ≤ 4 have been given. But, it is no longer possible to solve the nonlinear systems of equations for satisfying RKS property by symbolic manipulation packages to derive methods of higher orders (p ≥ 5). So, in this paper, we described the construction of A- and L-stable SGLMs with low error constants and QS properties, i.e., methods for which stability properties are determined by quadratic stability polynomials. The methods which were derived have the number of internal stages s = 2 and the stage order q equal to the order p which are presented up to the order six. The numerical experiments confirm that the new constructed methods are competitive with existing similar methods of the same order. It is our intention to investigate implementation of these methods in a variable stepsize environment. The implementation issues including the choice of appropriate starting procedures and stage predictors, strategies for efficient solution of system of nonlinear equations, estimation of local truncation error, design of stepsize changing strategies and comparisons with other methods are the subjects of future work. References [1] G. Dahlquist, A special stability problem for linear multistep methods, BIT 3 (1963) 27–43. [2] A. Iserles, S.P. Nørsett, On the theory of parallel Runge–Kutta methods, IMA J. Numer. Anal. 10 (1990) 463–488. [3] K.R. Jackson, S.P. Nørsett, The potential for parallelism in Runge–Kutta methods. Part 1: RK formulas in standard form, SIAM J. Numer. Anal. 32 (1995) 49–82. [4] J.C. Butcher, A modified multistep method for the numerical integration of ordinary differential equations, J. Assoc. Comput. Mach. 12 (1965) 124–135. [5] C.W. Gear, Hybrid methods for initial value problems in ordinary differential equations, SIAM J. Numer. Anal. 2 (1965) 69–86. [6] W.B. Gragg, H.J. Stetter, Generalized multistep predictor corrector methods, J. Assoc. Comput. Mach. 11 (1964) 188–209. [7] J.C. Butcher, On the convergence of numerical solutions to ordinary differential equations, Math. Comp. 20 (1966) 1–10. [8] J.C. Butcher, Numerical Methods for Ordinary Differential Equations, Wiley, New York, 2008. [9] J.C. Butcher, Diagonally-implicit multistage integration methods, Appl. Numer. Math. 11 (1993) 347–363. [10] J.C. Butcher, General linear methods, Acta Numer. 15 (2006) 157–256. [11] J.C. Butcher, Z. Jackiewicz, Implementation of diagonally implicit multistage integration methods for ordinary differential equations, SIAM J. Numer. Anal. 34 (1997) 2119–2141. [12] J.C. Butcher, Z. Jackiewicz, Construction of high order diagonally implicit multistage integration methods for ordinary differential equations, Appl. Numer. Math. 27 (1998) 1–12. [13] J.C. Butcher, Z. Jackiewicz, A reliable error estimation for diagonally implicit multistage integration methods, BIT 41 (2001) 656–665. [14] J.C. Butcher, Z. Jackiewicz, H.D. Mittelmann, A nonlinear optimization approach to the construction of general linear methods of high order, J. Comput. Appl. Math. 81 (1997) 181–196. [15] J.C. Butcher, H. Podhaisky, On error estimation in general linear methods for stiff ODEs, Appl. Numer. Math. 56 (2006) 345–357. [16] J.C. Butcher, W.M. Wright, The construction of practical general linear methods, BIT 43 (2003) 695–721. [17] A. Cardone, Z. Jackiewicz, Explicit Nordsieck methods with quadratic stability, Numer. Algorithms 60 (2012) 1–25. [18] L.L. Hewitt, A.T. Hill, Algebraically stable general linear methods and the G-matrix, BIT 49 (2009) 93–111. [19] Z. Jackiewicz, General Linear Methods for Ordinary Differential Equations, Wiley, New Jersey, 2009. [20] J.R. Cash, Second derivative extended backward differentiation formulas for the numerical integration of stiff systems, SIAM J. Numer. Anal. 18 (1981) 21–36. [21] P.C. Chakravarti, M.S. Kamel, Stiffly stable second derivative multistep methods with higher order and improved stability regions, BIT 23 (1983) 75–83. [22] T.M.H. Chan, R.P.K. Chan, A simplified approach to the order conditions of integration methods, Computing 77 (2006) 237–252. [23] R.P.K. Chan, A.Y.J. Tsai, On explicit two-derivative Runge–Kutta methods, Numer. Algorithms 53 (2010) 171–194. [24] W.H. Enright, Second derivative multistep methods for stiff ordinary differential equations, SIAM J. Numer. Anal. 11 (1974) 321–331. [25] G.K. Gupta, Implementing second-derivative multistep methods using the Nordsieck polynomial representation, Math. Comp. 32 (1978) 13–18. [26] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential–Algebraic Problems, Springer, Berlin, 1996. [27] G. Hojjati, M.Y. Rahimi Ardabili, S.M. Hosseini, New second derivative multistep methods for stiff systems, Appl. Math. Model. 30 (2006) 466–476.
228
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
[28] J.C. Butcher, G. Hojjati, Second derivative methods with RK stability, Numer. Algorithms 40 (2005) 415–429. [29] A. Abdi, M. Braś, G. Hojjati, On the construction of second derivative diagonally implicit multistage integration methods, Appl. Numer. Math. 76 (2014) 1–18. [30] A. Abdi, G. Hojjati, An extension of general linear methods, Numer. Algorithms 57 (2011) 149–167. [31] A. Abdi, G. Hojjati, Implementation of Nordsieck second derivative methods for stiff ODEs, Appl. Numer. Math. 94 (2015) 241–253. [32] A. Abdi, G. Hojjati, Maximal order for second derivative general linear methods with Runge–Kutta stability, Appl. Numer. Math. 61 (2011) 1046–1058. [33] A.K. Ezzeddine, G. Hojjati, A. Abdi, Sequential second derivative general linear methods for stiff systems, Bull. Iranian Math. Soc. 40 (2014) 83–100. [34] A. Abdi, J.C. Butcher, Order bounds for second derivative approximations, BIT 52 (2012) 273–281. [35] A. Abdi, J.C. Butcher, Applications of order arrows, Appl. Numer. Math. 62 (2012) 556–566. [36] A.K. Ezzeddine, G. Hojjati, A. Abdi, Perturbed second derivative multistep methods, J. Numer. Math. 23 (2015) 235–245. [37] J. Schur, Über Potenzreihen die im Innern des Einheitskreises beschränkt sind, J. Reine Angew. Math. 147 (1916) 205–232. [38] J.D. Lambert, Computational Methods in Ordinary Differential Equations, Wiley, New York, 1973. [39] E. Hairer, http://www.unige.ch/~hairer/testset/testset.html.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Construction of high-order quadratically stable second-derivative general linear methods for the numerical integration of stiff ODEs A. Abdi Faculty of Mathematical Sciences, University of Tabriz, Tabriz, Iran
article
info
Article history: Received 29 September 2015 Received in revised form 28 February 2016 Keywords: Stiff differential equations General linear methods Second derivative methods A- and L- stability Quadratic stability
abstract Theory of general linear methods (GLMs) for the numerical solution of autonomous system of ordinary differential equations of the form y′ = f (y) is extended to include the second derivative y′′ = g (y) := f ′ (y)f (y). This extension of GLMs is called second derivative general linear methods (SGLMs). In this paper we will construct two-stage A- and L-stable SGLMs of order p and stage order q = p up to six with low error constants. We will show the efficiency of the proposed methods by implementing on some well-known stiff problems. © 2016 Elsevier B.V. All rights reserved.
1. Introduction The characterization and discovery of methods suitable for the approximate numerical integration of stiff autonomous differential equations y′ (x) = f (y(x)), y(x0 ) = y0 ,
x ∈ [x0 , x],
(1)
where f : Rm → Rm and m is the dimensionality of the system, has attracted the interest of many researchers. Ideally such methods should be A-stable. But the order of an A-stable linear multistep method cannot exceed two (Dahlquist’s second barrier [1]). On the other hand, although there exist A-stable Runge–Kutta methods of arbitrarily high order, however the order of a diagonally implicit Runge–Kutta (DIRK) method with s stages and even an s-parallel SDIRK (singly diagonally implicit Runge–Kutta) method cannot exceed s + 1 [2] (see also [3]). Consequently research interest has been focused on the exploration of other paths. The idea of combining traditional one-step methods with traditional multistep methods is very natural and has been followed by many people over the last 50 years. Examples are given by the modified multistep methods [4–6]. All these methods belong to the class of general linear methods (GLMs). GLMs were introduced by Butcher [7] (see also [8]) and further investigated, for instance, in [9–19]. On the other hand, since an increase in the order of accuracy can be achieved by using the higher derivatives of the solution, especially second derivative, several authors have introduced methods for the solution of stiff systems which incorporate the higher derivatives in the integration formulas (e.g., [20–27]). For system (1), let g (y) := y′′ = f ′ (y)f (y). For problems in which g can be calculated along with f , at a moderate additional cost, second derivative methods become feasible. Although the mentioned GLMs include linear multistep methods, Runge Kutta and many other standard methods, but they were extended to second derivative general linear methods (SGLMs) to cover second derivative methods, too. These methods were introduced by Butcher and Hojjati in [28] and E-mail address: [email protected]. http://dx.doi.org/10.1016/j.cam.2016.02.054 0377-0427/© 2016 Elsevier B.V. All rights reserved.
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
219
investigated more in [29–33]. Abdi and Hojjati have divided SGLMs into four types [29,30], together with their intended applications (for nonstiff or stiff ODEs) and architectures (sequential or parallel). In [30,32], they have also obtained some upper bounds on the order of SGLMs of different types. These barriers on the order were found under the assumption of the Runge–Kutta stability (RKS) property (see [8,19]). Also, order arrows are used to obtain simpler and more illuminating proofs for these barriers in [34]. References are also made to [35,10] which include a discussion of order arrows. Construction and implementation of Nordsieck SGLMs for the numerical integration of stiff systems of ODEs are discussed by Abdi and Hojjati in [31]. In [36], to find more efficient methods with extensive absolute stability region for stiff systems, considering modified extended second derivative multistep methods [27] (MESDMMs) in SGLMs structure, perturbations of these methods were introduced which preserve their orders and improve their stability properties. The organization of the paper is as follows. After a short review of SGLMs, especially order conditions and types of these methods in Section 2, construction of A- and L-stable SGLMs with quadratic stability (QS) property of order p and stage order q = p up to six with low error constants is given in Section 3. In Section 4 some numerical experiments are presented to confirm the efficiency of the constructed methods. Finally, in Section 5 some concluding remarks are given and plans for future research are outlined. 2. A short introduction to SGLMs SGLMs for the numerical solution (1) are multivalue and multistage methods of the form
r s s [n] [n−1] [n] [n] aij g (Yj ) + uij yj , aij f (Yj ) + h2 Yi = h
i = 1, 2, . . . , s,
s s r [n] [n] [n] 2 ) + h b f ( Y = h b g ( Y ) + vij yj[n−1] , y ij ij j j i
i = 1, 2, . . . , r ,
j=1
j =1
j=1
j =1
j =1
(2)
j =1
where n = 1, 2, . . . , N , Nh = x − x0 . Here, introducing xn := x0 + nh, the internal stages Yi , i = 1, 2, . . . , s, are approximations of stage order q to y(xn−1 + ci h), i.e., [n]
= y(xn−1 + ci h) + O (hq+1 ),
[n]
Yi
and the external stages yi , i = 1, 2, . . . , r, are approximations of order p to the linear combinations of scaled derivatives of y(xn ), i.e., [n]
[n]
yi
= αi0 y(xn ) + αi1 hy′ (xn ) + αi2 h2 y′′ (xn ) + · · · + αip hp y(p) (xn ) + O (hp+1 ),
for some real numbers αi1 , αi2 , . . . , αip (cf. [29,32]). The vector y[n] denotes the computed quantities at step n, which also represent the incoming values for step n + 1. These methods are characterized by the abscissa vector c = [c1 c2 · · · cs ]T , the coefficient matrices A = [aij ] ∈ Rs×s , B = [bij ] ∈ Rr ×s ,
A = [aij ] ∈ Rs×s , B = [bij ] ∈ Rr ×s ,
U = [uij ] ∈ Rs×r , V = [vij ] ∈ Rr ×r ,
and the four integers: p, the order of the method; q, the stage order; r, the number of external stages; and s, the number of internal stages. They can be represented by a partitioned (s + r ) × (2s + r ) matrix of the form
A
A
U
B
B
V
.
It follows from [32] that the SGLM (2) has order p and stage order q = p if and only if exp(cz ) = zA exp(cz ) + z 2 A exp(cz ) + UWZ + O(z p+1 ), exp(z )WZ = zB exp(cz ) + z B exp(cz ) + VWZ + O(z 2
p+1
).
Here, exp(cz ) = [exp(c1 z ) exp(c2 z ) · · · exp(cs z )]T , Z = [1 z · · · z p ]T and
α10 α20 W = .. . αr0
α11 α21 .. . αr1
··· ··· .. . ···
α1p α2p .. . . αrp
(3) (4)
220
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
Expanding exp(cz ) and exp(z ) in (3) and (4) into power series around z = 0 and comparing the coefficients of z k , k = 0, 1, . . . , p, in the resulting expressions, the stage order (3) and order (4) can be reformulated in the form
e − U α0 = 0, α0 − V α0 = 0, c − Ae − U α1 = 0, α1 + α0 − Be − V α1 = 0, k−2 k−1 ck Ac Ac − − − U αk = 0, k ! ( k − 1 )! ( k − 2)! k Bc k−1 Bc k−2 αk−l − − − V αk = 0, l! (k − 1)! (k − 2)! l=0
(5)
for k = 2, 3, . . . , p and e = [1 1 · · · 1] ∈ Rs . The first and second equations of (5) are called pre-consistency and consistency conditions, respectively [30]. The most convenient choice for U, in the case of r = s, is U = I, where I is the identity matrix of appropriate dimension. This choice is valid since a transformation can be found to force U = I while keeping the structure of the other matrices unchanged. The matrix W is uniquely determined by the coefficients of the method if r = s and U = I. In this case it follows from order conditions (3) that [32] W = C − ACK − ACK 2 , s×(p+1)
where C = (Cij ) ∈ R Cij =
j −1 ci
(j − 1)!
,
(6)
is the Vandermonde matrix with components
1 ≤ i ≤ s, 1 ≤ j ≤ p + 1 ,
and K ∈ R(p+1)×(p+1) is the shifting matrix defined by K = [0 e1 · · · ep ] with ej as the jth unit vector. It was proved in [29] that the method (2) with p = q = r = s, U = I and Ve = e has order p and stage order p if and only if B = B0 − AB1 − AB2 − VB3 − (B − V A)B4 + VA, where
1+c
φj (x)dx B0 = φj (cj ) c s i φj (x)dx 0 B3 = φj (cj ) i
s
,
0
B1 =
i,j=1
,
B4 =
i ,j = 1
φj (1 + ci ) φj (cj )
φj′ (ci )
s
φj (cj )
s
,
B2 =
i,j=1
φj′ (1 + ci ) φj (cj )
s
,
i,j=1
,
i,j=1
and s
φi (x) =
(x − cj ),
i = 1, 2, . . . , s.
j=1,j̸=i
In the recent papers on SGLMs, these methods were investigated under the assumption that the coefficient matrices A and A have the form
λ
a21 A= .. .
λ .. .
as1
as2
,
..
. ···
µ
a21 A= .. .
µ .. .
as1
as2
λ
,
..
. ···
µ
with λ ≥ 0 and µ ≤ 0, where the explicit methods corresponding to λ = µ = 0 are called type 1 methods and the implicit formulas corresponding to λ > 0, µ < 0 are called type 2 methods. Type 3 and type 4 methods for which A and A have the form A = diag {λ, λ, . . . , λ} , A = diag {µ, µ, . . . , µ} , with λ = µ = 0 and λ > 0, µ < 0, respectively. Order barriers for different types of these methods under the assumption of the Runge–Kutta stability (RKS) property are obtained in [30,32]. Examples of many these methods of various types, which are appropriate for nonstiff or stiff differential systems in a sequential or parallel computing environment, are given in [29,30,32,28]. Also, as explained in these papers the stability properties of the method (2) are determined by the stability matrix M (z ) = V + (zB + z 2 B)(I − zA − z 2 A)−1 U , and the corresponding stability function p(w, z ) = det(w I − M (z )), where w and z are complex numbers.
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
221
3. Construction of two stage methods Since we are interested in developing a code suitable for solving stiff problems, A-stability should be necessary. In this section, we construct type 2 methods with r = s = 2 of order and stage order 3 up to 6 which are A-stable. It will be assumed that V is a rank-one matrix of the form
V =
v1 v1
v2 , v2
with v1 + v2 = 1. This choice implies that V is power bounded. Consequently, these methods are zero-stable. Moreover, it will also be assumed that U = I. Coefficients partitioned matrix for these methods is given by
A
A
U
B
B
V
λ
0
µ
0
1
0
λ
µ
0
b12
1−v
1
v
b22
1−v
v
a 21 = b 11
b12
a21 b11
b21
b22
b21
with λ > 0 and µ < 0. Stability properties of these methods are determined by quadratic stability functions, i.e., the stability function p(w, z ) of such methods has the form p(w, z ) = (1 − λz − µz 2 )2 w 2 − p1 (z )w + p2 (z ),
(7)
where p1 (z ) and p2 (z ) are polynomials of degree less than or equal to 2s = 4. For these methods to be A-stable, the Schur criterion [37] (see also [38]) can be used. This states that, a quadratic polynomial d0 w 2 + d1 w + d2 has both its roots in the open unit disc if and only if
|d0 |2 − |d2 |2 > 0,
(8)
(|d0 |2 − |d2 |2 )2 − |d0 d1 − d2 d1 |2 > 0.
(9)
and
In the stability function (7), p1 (z ) = (1 − λz − µz 2 )2 tr(M (z )), and p2 (z ) = (1 − λz − µz 2 )2 det(M (z )), are polynomials of degree at most four. So, assuming A-stability of such methods, the L-stability requirement is equivalent to lim
z →∞
p1 ( z )
(1 − λz − µ )
z2 2
= 0,
lim
z →∞
p2 (z )
(1 − λz − µz 2 )2
= 0,
which imply that the coefficients of z 4 in the numerator of tr(M (z )) and det(M (z )) must be zero. By simple algebraic calculations, these lead to
µ2 − tr(B)µ + a21 b12 det(B − µV ) + a21 v(b12 − b22 )
= =
0, 0.
(10)
Denoting the error constant of the method of order p by Ep (̸= 0), the stability function p(w, z ) in (7) satisfies p(exp(z ), z ) = Ep z p+1 + O (z p+2 ).
(11)
In the construction of the methods, we will use the expressions ‘‘L-conditions’’ and ‘‘E-condition’’ to mean satisfying (10) and (11), respectively. 3.1. Methods of order 3 In this subsection, we construct L-stable methods of order p = q = 2s − 1 = 3, r = s and abscissa vector c = [0 1]T with the error constant E3 = 10−5 . Solving stage order and order conditions (3) and (4), L-conditions and E-condition, we obtain √ a four-parameter family of methods of order and stage order p = q = 3 depending on a21 , a21 , λ, and µ with v = (α ± β)/γ where α, β and γ are complicated expressions in λ, µ and a21 . Then we apply the Schur criterion, i.e. conditions (8) and (9), to find A-stable (and so L-stable) methods in the parameter space (λ, µ) for selected values of the parameters a21 and a21 and β ≥ 0. Restricting
222
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
Fig. 1. L-stability regions in the (λ, µ)-plane for quadratically stable SGLMs with s = 2, p = q = 3 and the error constant E3 = 10−5 for specific values of the parameters.
λ and µ to 0 < λ ≤ 2 and −2 ≤ µ < 0, pairs of (λ, µ) values giving L-stability are presented in Fig. 1. The four cases in Fig. 1 correspond to: √ α+ β . γ √ α− β v= γ . √ v = α+γ β . √ v = α−γ β .
• Case 1: a21 = 12 , a21 = − 12 , v = • Case 2: a21 = 12 , a21 = − 12 , • Case 3: a21 = 2, a21 = −2, • Case 4: a21 = 2, a21 = −2,
We present here just a single example, from family of Case 1 methods, characterized by (λ, µ) = ( 34 , − 51 ). The coefficients for the method are
3 4
0
− 51
0
1
0
1 2
3 4
− 21
− 51
0
−801+940 v−475 v 2
21−250 v+325 v 2 60 d
89+127 v−120 v 2 30 d
−1+29 v−40 v 2
60 d
1−v
−96+103 v−95 v 2
−138+65 v+65 v 2
−8−6 v−40 v 2
12 d
12 d
83+11 v−60 v 2 15 d
, v
where d = −13 + 5v and v =
15 d 15 d
1−v
1
(12)
v
√ 104509+ 42225283081 . 360000
3.2. Methods of order 4 In this subsection, we construct L-stable methods of order p = q = 2s = 4, r = s and abscissa vector c = [0 1]T with the error constant E4 = 10−5 . Solving stage order and order conditions (3) and (4), L-conditions and E-condition, we√obtain a two-parameter family of methods of order and stage order p = q = 4 depending on a21 and µ with λ = (α ± β)/(1 + 12µ) where α and β are expressions in µ. Then we apply the Schur criterion to find A-stable (and so L-stable) methods in the parameter space (a21 , µ) for β ≥ 0. Restricting a21 and µ to −1 ≤ a21 ≤ 3 and −1 ≤ µ < 0, pairs of (a21 , µ) values giving L-stability are presented in Fig. 2.
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
223
Fig. 2. L-stability region in the (a21 , µ)-plane for quadratically stable SGLMs with s = 2, p = q = 4 and the error constant E4 = 10−5 .
Fig. 3. A-stability region in the (λ, c1 )-plane for quadratically stable SGLMs with s = 2, p = q = 5 and the error constant E5 = 10−5 for specific values of the parameter µ.
√
We present here just a single example, characterized by λ = (α+ β)/(1+12µ) and (a21 , µ) = ( 14 , − 14 ). The coefficients for the method are
λ
0
− 14
0
1
0
1 4
λ
−205+486 λ−288 λ2
− 14
0
1
−83−156 λ+288 λ2
−9(7−22 λ+16 λ2 )
−18953+64500 λ−73152 λ2 +27648 λ3
16 d2
8 d2
72 d1 d2
3 d1 8 d2
9 d1 4 d2
4 d2 −9 d1 4 d2
−443+676 λ−192 λ2
87−182 λ+96 λ2 8 d2
−14333+48636 λ−55008 λ2 +20736 λ3
1103−2520 λ+1440 λ2 72 d1 d2
9 d1 4 d2
4 d2 −9 d1 4 d2
16 d2
where d1 = 8λ − 7, d2 = 6λ − 5 and λ =
9 d1
72 d1 d2
(13)
√ 567509− 807465081 . 600000
3.3. Methods of order 5 In this subsection, we construct A-stable methods of order p = q = 2s + 1 = 5, r = s and abscissa vector c = [c1 1]T with c1 as a free parameter and the error constant E5 = 10−5 . Solving stage order and order conditions (3) and (4), we obtain a four-parameter family of methods of order and stage order p = q = 5 depending on λ, µ, and c1 ̸= 0, 1. Then we apply the Schur criterion to find A-stable methods in the parameter space (λ, c1 ) for selected values of the parameter µ. Restricting λ and c1 to 0 < λ ≤ 2 and c1 ∈ [−1, 1) \ {0}, pairs of (λ, c1 ) values giving A-stability are presented in Fig. 3 for selected values of the parameter µ. 3 187 We present here just a single example characterized by µ = − 20 and (λ, c1 ) = ( 250 , − 12 ). The coefficients for the method are
187 250
0
3 − 20
0
1
0
98481 62000
187 250
10713 12400
3 − 20
0
1
55804 50625
31376 − 253125
5381 − 16875
31 675
2521 2025
496 − 2025
.
41448379 25110000
− 11938687 15693750
− 1265651 8370000
38813 − 334800
2521 2025
496 − 2025
(14)
224
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
Fig. 4. A-stability region in the (λ, c1 )-plane for quadratically stable SGLMs with s = 2, p = q = 6.
3.4. Methods of order 6 In this subsection, we construct A-stable methods of order p = q = 2s + 2 = 6, r = s and abscissa vector c = [c1 1]T with c1 as a free parameter. Solving stage order and order conditions (3) and (4), we obtain a two-parameter family of methods of order and stage order p = q = 6 depending on λ and c1 ̸= 0, 1. Then we apply the Schur criterion to find A-stable methods in the parameter space (λ, c1 ). Restricting λ and c1 to 0 < λ ≤ 2 and c1 ∈ [−1, 1) \ {0}, pairs of (λ, c1 ) values giving A-stability are presented in Fig. 4. By searching for the smallest absolute value of error constant within these A-acceptable values for (λ, c1 ), we obtain 179 , −1) and the error constant E6 ≈ −0.899 × 10−3 . (λ, c1 ) = ( 250 179 Setting (λ, c1 ) = ( 250 , −1), the coefficients matrices of the method take the following forms
179 250
0
− 16
0
1
0
502 375
179 250
139 125
− 61
0
1
666071 600000
627 200000
64011 − 200000
37 8000
937 800
− 137 800
.
383357 200000
85319 − 600000
13611 − 200000
257 1600
937 800
− 137 800
(15)
4. Numerical experiments In this section we present some numerical results which aim to indicate accuracy and efficiency of the proposed methods in Section 3 in the integration of stiff differential systems. We use the following nonlinear stiff problems: P1. The non-linear stiff system of ODEs
y′1 (x) = −1002y1 (x) + 1000y22 (x), y′2 (x) = y1 (x) − y2 (x)(1 + y2 (x)),
y1 (0) = 1, y2 (0) = 1,
where the exact solution is [y1 (x), y2 (x)]T = [exp(−2x), exp(−x)]T and x ∈ [0, X ]. P2. The non-linear stiff system of ODEs
y′1 (x) = −10004y1 (x) + 10000y42 (x), y′2 (x) = y1 (x) − y2 (x)(1 + y32 (x)),
y1 (0) = 1, y2 (0) = 1,
where x ∈ [0, X ] and the exact solution is y1 (x) = exp(−4 x) and y2 (x) = exp(−x). P3. The non-linear stiff system of ODEs
′ 3 6 3 6 y1 (x) = −1000y1 y2 − cos (x)sin (x) − sin(x), y′2 (x) = −1000 y52 y43 − sin9 (x) + cos(x) ′ y3 (x) = −1000 y21 y33 − cos2 (x) sin3 (x) + cos(x)
y1 (0) = 1, y2 (0) = 0, y3 (0) = 0,
where the exact solution is [y1 (x), y2 (x), y3 (x)]T = [cos(x), sin(x), sin(x)]T and x ∈ [0, 1].
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
225
Table 1 The global error for SGLMs (12) and (13) at the end of the interval of integration [0, 1] applied to the problem P1. k
Method of order 3
Method of order 4
∥eh (X )∥ 1.02 × 10−9 6.91 × 10−11 6.89 × 10−12 7.72 × 10−13
5 6 7 8
p
∥eh (X )∥
p
3.88 3.27 3.16
2.94 × 10−10 9.49 × 10−12 3.09 × 10−13 1.97 × 10−14
4.95 4.94 3.97
P4. A system of differential equations, called CUSP, resulting from discretization of the diffusion terms by the method of lines of the periodic boundary-value problem [26]
1 ∂ 2y ∂y = − (y3 + ay + b) + σ 2 , ∂t ε ∂x ∂a ∂ 2a = b + 0.07v + σ 2 , ∂t ∂x ∂ b ∂ 2b 2 = (1 − a )b − a − 0.4y + 0.035v + σ 2 , ∂t ∂x where
v=
u u + 0.1
,
u = (y − 0.7)(y − 1.3).
This problem takes the form
y˙ i = −ε −1 (y3i + ai yi + bi ) + D(yi−1 − 2yi + yi+1 ), a˙ i = bi + 0.07vi + D(ai−1 − 2ai + ai+1 ), i = 1, 2, . . . , Nx ˙ bi = (1 − a2i )bi − ai − 0.4yi + 0.035vi + D(bi−1 − 2bi + bi+1 ), where
vi =
ui ui + 0.1
,
ui = (yi − 0.7)(yi − 1.3),
D = σ Nx2 ,
with periodic boundary conditions y0 := yNx , yNx +1 := y1 ,
a0 := aNx , aNx +1 := a1 ,
b0 := bNx , bNx +1 := b1 .
We take σ = 1/144, ε = 10−4 , Nx = 32 and the initial values as yi (0) = 0,
ai (0) = −2 cos
2iπ Nx
,
bi (0) = 2 sin
2iπ , Nx
i = 1, 2, . . . , Nx ,
with tout = 1.1 [26]. The implementation of the proposed methods in Section 3, needs a suitable starting procedure that gives the initial input vector y[0] . Here, similar to the introduced starting procedure in [31], we use the Runge–Kutta methods of suitable order, which give sufficient output information to provide an approximation to y[0] . Also, since for the all proposed methods we [n] have cs = c2 = 1, finishing procedure does not need and we use Ys[n] = Y2 as an approximation of order q = p to y(xn ). We have applied the methods (12) and (13) of orders 3 and 4, respectively, to the problem P1 with X = 1 and the methods (14) and (15) of orders 5 and 6, respectively, to the problem P2 with X = 2. We have implemented these methods with fixed stepsize h = X /2k , in correspondence to several integer values of k. The results of numerical experiments are presented in Tables 1 and 2. In these tables, we have listed norm of error ∥eh (X )∥ at the endpoint of integration X and numerical estimate to the order of convergence, p, computed by the formula p = log2 (∥eh (X )∥/∥eh/2 (X )∥), where eh (X ) and eh/2 (X ) are errors corresponding to stepsizes h and h/2. We can observe that the proposed methods are convergent with expected order p = 3, 4, 5 and 6. To illustrate the efficiency of the proposed methods in this paper, we have applied the Runge–Kutta–Gauss (RK–Gauss) methods of orders p = 4 and p = 6 [8] and the constructed SGLMs in Section 3 of the same orders to the problem P3. In the implementation of implicit methods, such as RK–Gauss methods, to approximate stage values, the Jacobian matrix ∂ f /∂ y must be found. So, the function g for the second derivative methods can be calculated by (∂ f /∂ y)f without any additional computational cost. Also, in order to save computational effort, we use (∂ f /∂ y)2 (cf. [20,24]) as an approximation
226
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
Table 2 The global error for SGLMs (14) and (15) at the end of the interval of integration [0, 2] applied to the problem P2. k
Method of order 5
∥eh (X )∥ 3 4 5 6
9.13 × 10−11 6.30 × 10−12 1.68 × 10−13 5.72 × 10−15
Method of order 6 p
∥eh (X )∥
p
3.86 5.23 4.88
1.29 × 10−9 1.98 × 10−11 3.02 × 10−13 6.33 × 10−15
6.03 6.03 5.58
Fig. 5. Comparison of the number of function evaluations for the methods of order 4 (top) and order 6 (bottom) applied to the problem P3.
to the Jacobian matrix ∂ g /∂ y. This approximation is a piecewise constant approximation to the Jacobian matrix. In Fig. 5, we compare the number of function evaluations, nfe, versus accuracy of the methods. This figure confirms the efficiency of the proposed SGLMs. Also, to compare the performance of the new methods and SGLMs previously constructed, we have applied the constructed methods in [31] and in the present paper of orders p = 3 and p = 4 to the problem P4. The number of function evaluations versus the error of the methods are plotted in Fig. 6 which illustrates superiority of the new methods from computational cost point of view. In this figure, error is the difference between the solution reported in [39] as a reference solution and the numerical results of the applied SGLMs at the endpoint of integration tout = 1.1.
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
227
Fig. 6. Performance comparison of the methods (12) and (13) with the SGLMs constructed in [31] of orders 3 and 4.
5. Conclusion The construction of the class of second derivative diagonally implicit multistage integration methods (SDIMSIMs), as a subclass of SGLMs, with p = q = r = s and RKS property has been described in [29], and some examples of these methods of all four types up to the order p ≤ 4 have been given. But, it is no longer possible to solve the nonlinear systems of equations for satisfying RKS property by symbolic manipulation packages to derive methods of higher orders (p ≥ 5). So, in this paper, we described the construction of A- and L-stable SGLMs with low error constants and QS properties, i.e., methods for which stability properties are determined by quadratic stability polynomials. The methods which were derived have the number of internal stages s = 2 and the stage order q equal to the order p which are presented up to the order six. The numerical experiments confirm that the new constructed methods are competitive with existing similar methods of the same order. It is our intention to investigate implementation of these methods in a variable stepsize environment. The implementation issues including the choice of appropriate starting procedures and stage predictors, strategies for efficient solution of system of nonlinear equations, estimation of local truncation error, design of stepsize changing strategies and comparisons with other methods are the subjects of future work. References [1] G. Dahlquist, A special stability problem for linear multistep methods, BIT 3 (1963) 27–43. [2] A. Iserles, S.P. Nørsett, On the theory of parallel Runge–Kutta methods, IMA J. Numer. Anal. 10 (1990) 463–488. [3] K.R. Jackson, S.P. Nørsett, The potential for parallelism in Runge–Kutta methods. Part 1: RK formulas in standard form, SIAM J. Numer. Anal. 32 (1995) 49–82. [4] J.C. Butcher, A modified multistep method for the numerical integration of ordinary differential equations, J. Assoc. Comput. Mach. 12 (1965) 124–135. [5] C.W. Gear, Hybrid methods for initial value problems in ordinary differential equations, SIAM J. Numer. Anal. 2 (1965) 69–86. [6] W.B. Gragg, H.J. Stetter, Generalized multistep predictor corrector methods, J. Assoc. Comput. Mach. 11 (1964) 188–209. [7] J.C. Butcher, On the convergence of numerical solutions to ordinary differential equations, Math. Comp. 20 (1966) 1–10. [8] J.C. Butcher, Numerical Methods for Ordinary Differential Equations, Wiley, New York, 2008. [9] J.C. Butcher, Diagonally-implicit multistage integration methods, Appl. Numer. Math. 11 (1993) 347–363. [10] J.C. Butcher, General linear methods, Acta Numer. 15 (2006) 157–256. [11] J.C. Butcher, Z. Jackiewicz, Implementation of diagonally implicit multistage integration methods for ordinary differential equations, SIAM J. Numer. Anal. 34 (1997) 2119–2141. [12] J.C. Butcher, Z. Jackiewicz, Construction of high order diagonally implicit multistage integration methods for ordinary differential equations, Appl. Numer. Math. 27 (1998) 1–12. [13] J.C. Butcher, Z. Jackiewicz, A reliable error estimation for diagonally implicit multistage integration methods, BIT 41 (2001) 656–665. [14] J.C. Butcher, Z. Jackiewicz, H.D. Mittelmann, A nonlinear optimization approach to the construction of general linear methods of high order, J. Comput. Appl. Math. 81 (1997) 181–196. [15] J.C. Butcher, H. Podhaisky, On error estimation in general linear methods for stiff ODEs, Appl. Numer. Math. 56 (2006) 345–357. [16] J.C. Butcher, W.M. Wright, The construction of practical general linear methods, BIT 43 (2003) 695–721. [17] A. Cardone, Z. Jackiewicz, Explicit Nordsieck methods with quadratic stability, Numer. Algorithms 60 (2012) 1–25. [18] L.L. Hewitt, A.T. Hill, Algebraically stable general linear methods and the G-matrix, BIT 49 (2009) 93–111. [19] Z. Jackiewicz, General Linear Methods for Ordinary Differential Equations, Wiley, New Jersey, 2009. [20] J.R. Cash, Second derivative extended backward differentiation formulas for the numerical integration of stiff systems, SIAM J. Numer. Anal. 18 (1981) 21–36. [21] P.C. Chakravarti, M.S. Kamel, Stiffly stable second derivative multistep methods with higher order and improved stability regions, BIT 23 (1983) 75–83. [22] T.M.H. Chan, R.P.K. Chan, A simplified approach to the order conditions of integration methods, Computing 77 (2006) 237–252. [23] R.P.K. Chan, A.Y.J. Tsai, On explicit two-derivative Runge–Kutta methods, Numer. Algorithms 53 (2010) 171–194. [24] W.H. Enright, Second derivative multistep methods for stiff ordinary differential equations, SIAM J. Numer. Anal. 11 (1974) 321–331. [25] G.K. Gupta, Implementing second-derivative multistep methods using the Nordsieck polynomial representation, Math. Comp. 32 (1978) 13–18. [26] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential–Algebraic Problems, Springer, Berlin, 1996. [27] G. Hojjati, M.Y. Rahimi Ardabili, S.M. Hosseini, New second derivative multistep methods for stiff systems, Appl. Math. Model. 30 (2006) 466–476.
228
A. Abdi / Journal of Computational and Applied Mathematics 303 (2016) 218–228
[28] J.C. Butcher, G. Hojjati, Second derivative methods with RK stability, Numer. Algorithms 40 (2005) 415–429. [29] A. Abdi, M. Braś, G. Hojjati, On the construction of second derivative diagonally implicit multistage integration methods, Appl. Numer. Math. 76 (2014) 1–18. [30] A. Abdi, G. Hojjati, An extension of general linear methods, Numer. Algorithms 57 (2011) 149–167. [31] A. Abdi, G. Hojjati, Implementation of Nordsieck second derivative methods for stiff ODEs, Appl. Numer. Math. 94 (2015) 241–253. [32] A. Abdi, G. Hojjati, Maximal order for second derivative general linear methods with Runge–Kutta stability, Appl. Numer. Math. 61 (2011) 1046–1058. [33] A.K. Ezzeddine, G. Hojjati, A. Abdi, Sequential second derivative general linear methods for stiff systems, Bull. Iranian Math. Soc. 40 (2014) 83–100. [34] A. Abdi, J.C. Butcher, Order bounds for second derivative approximations, BIT 52 (2012) 273–281. [35] A. Abdi, J.C. Butcher, Applications of order arrows, Appl. Numer. Math. 62 (2012) 556–566. [36] A.K. Ezzeddine, G. Hojjati, A. Abdi, Perturbed second derivative multistep methods, J. Numer. Math. 23 (2015) 235–245. [37] J. Schur, Über Potenzreihen die im Innern des Einheitskreises beschränkt sind, J. Reine Angew. Math. 147 (1916) 205–232. [38] J.D. Lambert, Computational Methods in Ordinary Differential Equations, Wiley, New York, 1973. [39] E. Hairer, http://www.unige.ch/~hairer/testset/testset.html.