- Email: [email protected]
New second derivative multistep methods for stiff systems
Applied Mathematical Modelling 30 (2006) 466–476 www.elsevier.com/locate/apm
New second derivative multistep methods for stiff systems G. Hojjati
a,*
, M.Y. Rahimi Ardabili a, S.M. Hosseini
b
a
b
Faculty of Mathematics, Tabriz University, 29 Bahman Street, 5166614766 Tabriz, Iran Department of Mathematics, Tarbiat Modarres University, P.O. Box 14115-175, Tehran, Iran Received 1 May 2004; received in revised form 1 May 2005; accepted 27 June 2005 Available online 24 August 2005
Abstract In this paper we present details of a new class of second derivative multistep methods. Stability analysis is discussed and an improvement in stability region is obtained. With a simple modification we take advantage of calling for the solution of algebraic equations with the same coefficient matrix in each step. Moreover, using IOM to solve the resulting linear systems, the coefficient matrix is not needed explicitly. Some numerical experiments and comparison with several popular codes are given, showing strong superiority of this new class of methods. 2005 Elsevier Inc. All rights reserved. Keywords: Stiff ODEs; Multistep and multiderivative methods; Stability aspects; IOM
1. Introduction In recent years, numerous works have focused on the development of more advanced and efficient methods for stiff problems. A potentially good numerical method for the solution of stiff systems of ODEs must have good accuracy and some reasonably wide region of absolute stability [1]. A-stability requirement puts a severe limitation on the choice of suitable methods for stiff problems. The search for higher order A-stable multistep methods is carried out in the two main directions: *
Corresponding author. Tel.: +98 411 3312649; fax: +98 411 3342102. E-mail address: [email protected] (G. Hojjati).
0307-904X/$ - see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2005.06.007
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
467
• use higher derivatives of the solutions; • throw in additional stages, off-step points, super-future points and like. This leads into the large field general linear methods [2]. In this paper we introduce a new class of second derivative multistep methods (SDMMnew) that have good stability properties. This new class of methods is equipped with an efficient computational algorithm in which for all iterations the final matrix is not changed. Using IOM [3] to solving linear systems in each iteration one dose not require the storage of coefficient matrix in any form. In the same lines of MF-MEBDF [4], the MF-SDMMnew is also introduced which in turn is formed to be very effective and efficient for stiff ODEs.
2. Second derivative multistep methods (SDMM) Let us consider the stiff initial value problem y 0 ðxÞ ¼ f ðx; yðxÞÞ;
yðx0 Þ ¼ y 0 ;
ð2:1Þ m
m
m
on the finite interval I = [x0, xN] where y : ½x0 ; xN ! R and f : ½x0 ; xN R ! R is continuous and differentiable. The general SDMM can be written in the form: k X
aj y nþj ¼ h
j¼0
k X j¼0
bj fnþj þ h2
k X
cj gnþj ;
ð2:2Þ
j¼0
where aj, bj, cj are parameters to be determined and y00 = fx + fy Æ f = : g(x, y). If either bk or ck is nonzero, the formula is implicit. Taylor expansion shows that the method (2.2) is of order p if and only if k X j¼0
aj jq ¼ q
k X
bj jq1 þ qðq þ 1Þ
j¼0
k X
cj jq2 ;
ð2:3Þ
j¼0
with 0 6 q 6 p. Some known important SDMM schemes that will be used for comparison are as follows: • The Enright [5] k-step formulas of order k + 2 which takes the following form: y nþ1 y n ¼ h
k X
bj fnþjkþ1 þ h2 ck gnþ1 .
j¼0
• Second derivative BDF (SDBDF) [2] with the formula ! ! k k k X X X h2 1 1 ri y nþ1 hf nþ1 g ¼ ; j j 2 nþ1 i j¼1 i¼1 j¼1 where $yn+1 = yn+1 yn. This is a class of k-step formulas of order k + 1. • Second derivative extended backward differentiation formulas (E2BD), that was introduced by Cash [6], with the following forms:
468
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
Class 1 predictor : y nþk y nþk1 ¼ h
k X
bj fnþj þ h2 ck gnþk ;
j¼0
corrector : y nþk y nþk1 ¼ h
kþ1 X
fnþj þ h2 ðc g þ c g b j k nþk kþ1 nþkþ1 Þ;
j¼0
and Class 2 predictor : the same as above; corrector : y nþk y nþk1 ¼ h
kþ1 X
fnþj þ h2c g . b j k nþk
j¼0
• Special class of SDMM, introduced by Ismail and Ibrahim [7] of the form k X
aj y nþj ¼ hbk ðfnþk b fnþk1 Þ þ h2 ck ðgnþk c gnþk1 Þ.
j¼0
For b* = 0, c* = 0 this is the same as the SDBDF method.
3. The new SDMM We are going to introduce a new SDMM that takes the following general form: k X
^ fnþk þ h2 ð^c g ^c g ^ aj y nþj ¼ hb k k nþk kþ1 nþkþ1 Þ;
ð3:1Þ
j¼0
ak ¼ 1 and the other coefficients are chosen so that (3.1) has where g(x, y) = y00 = fx + fyf and ^ order k + 2. The coefficients of k-step methods of class (3.1) are given in Table 1, for k = 1, 2, . . ., 8. In this form also we use the one super-future point technique. We design the above form such that in addition to the method have good stability properties and higher order, it is suitable to apply matrix-free technique as we will show in Section 5. Assuming that the solution values yn, yn+1, . . ., yn+k1 are available, the way in which (3.1) is used in practice is by carry out the following computations: Stage 1. Compute y nþk as the solution of k X
aj y nþj ¼ hbk fnþk þ h2 ck gnþk ;
ð3:2Þ
j¼0
where ak = 1 and the other coefficients are chosen so that (3.2) has order k + 1. The coefficients of k-step methods of class (3.2) are given in Table 2, for k = 1, 2, . . ., 8.
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
469
Table 1 Coefficients in (3.1) k
1
2
3
4
5
6
7
8
d ^ b k ^ck ^ckþ1 ^ a0 ^ a1 ^ a2 ^ a3 ^ a4 ^ a5 ^ a6 ^ a7
6 1 4/d 1/d 1
87 78/d 34/d 4/d 9/d 96/d
5239 4350/d 1530/d 108/d 136/d 1161/d 6264/d
149,809 117300/d 35928/d 1728/d 1431/d 12,325/d 52,920/d 191,808/d
7,364,773 5,500,740/d 1,524,600/d 54,000/d 32,148/d 295,875/d 1,279,000/d 3,663,000/d 10,012,500/d
6,399,049 4,595,220/d 1,179,000/d 32,400/d 14,680/d 146124/d 672,975/d 1,941,800/d 4,279,600/d 9,178,380/d
29,993,974,247 20,828,002,860/d 5,020,961,400/d 111,132,000/d 39,877,200/d 430,293,500/d 2,142,550,872/d 6,582,641,625/d 14,259,196,000/d 24,594,437,700/d 45,159,723,000/d
139,632,996,029 94,170,512,520/d 21,552,904,800/d 395,136,000/d 115,643,675/d 1,350,144,000/d 7,286,417,600/d 24,196,635,904/d 55,785,177,000/d 96,751,244,800/d 139,255,256,000/d 219,777,465,600/d
Table 2 Coefficients in (3.2) k
1
2
3
4
5
6
7
8
d bk ck a0 a1 a2 a3 a4 a5 a6 a7
2 1 1/d 1
7 6/d 2/d 1/d 8/d
85 66/d 18/d 4/d 27/d 108/d
415 300/d 72/d 9/d 64/d 216/d 576/d
12,019 8220/d 1800/d 144/d 1125/d 4000/d 9000/d 18,000/d
13,489 8820/d 1800/d 100/d 864/d 3375/d 8000/d 13,500/d 21,600/d
726,301 457,380/d 88,200/d 3600/d 34,300/d 148,276/d 385,875/d 686,000/d 926,100/d 1,234,800/d
3,144,919 1,917,720/d 352,800/d 11,025/d 115,200/d 548,800/d 1,580,544/d 3,087,000/d 4,390,400/d 4,939,200/d 5,644,800/d
Stage 2. Compute y nþkþ1 as the solution of k X
aj y nþjþ1 ¼ hbk fnþkþ1 þ h2 ck gnþkþ1 .
ð3:3Þ
j¼0
Stage 3. Evaluate gnþkþ1 ¼ gðxnþkþ1 ; y nþkþ1 Þ. Stage 4. Compute yn+k as the solution of k X j¼0
^ fnþk þ h2 ð^c g ^c g ^ aj y nþj ¼ hb k k nþk kþ1 nþkþ1 Þ.
ð3:4Þ
470
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
Lemma 3.1. Given that (i) formula (3.2) is of order k + 1, (ii) formula (3.1) is of order k + 2, (iii) the implicit algebraic equations defining y nþk and y nþkþ1 are solved exactly, then scheme (3.4) has order k + 2.
Proof. Suppose the values yn,yn+1, . . ., yn+k1 be exact. From (3.2) we have yðxnþk Þ y nþk ¼ C 1 hkþ2 y ðkþ2Þ ðxnþk Þ þ Oðhkþ3 Þ.
ð3:5Þ
From (3.5) for one super-future point and if we suppose that y(xn+k) = yn+k we have yðxnþkþ1 Þ y nþkþ1 ¼ C 1 hkþ2 y ðkþ2Þ ðxnþkþ1 Þ þ Oðhkþ3 Þ ¼ C 1 hkþ2 y ðkþ2Þ ðxnþk Þ þ C 1 hkþ3 y ðkþ3Þ ðxnþk Þ þ Oðhkþ3 Þ. But since in (3.3) we apply y nþk , we must add the error of ak1 ðyðxnþk Þ y nþk Þ to the above expression. Hence yðxnþkþ1 Þ y nþkþ1 ¼ C 1 ð1 ak1 Þhkþ2 y ðkþ2Þ ðxnþk Þ þ Oðhkþ3 Þ.
ð3:6Þ
If now C2hk+3y(k+3)(xn+k) + O(hk+4) is the defect of formula (3.1), replacing g(xn+k+1, y(xn+k+1)) by gðxnþkþ1 ; y nþkþ1 Þ adds the expression obtained in (3.6) to this error and we obtain yðxnþk Þ y nþk ¼ C 2 hkþ3 y ðkþ3Þ ðxnþk Þ h2ckþ1 gðxnþkþ1 ; yðxnþkþ1 ÞÞ gðxnþkþ1 ; y nþkþ1 Þ . From (3.6) og ðgÞðyðxnþkþ1 Þ y nþkþ1 Þ oy og ¼ ðgÞC 1 ð1 ak1 Þhkþ2 y ðkþ2Þ ðxnþk Þ þ Oðhkþ3 Þ. oy
gðxnþkþ1 ; yðxnþkþ1 ÞÞ gðxnþkþ1 ; y nþkþ1 Þ ¼
This yields
og ðkþ2Þ yðxnþk Þ y nþk ¼ h C2y ðxnþk Þ ðgÞC 1 ð1 ak1 Þh^ckþ1 y ðxnþk Þ þ Oðhkþ4 Þ; oy which shows that the order of scheme (4) is k + 2. h kþ3
ðkþ3Þ
4. Stability analysis We now examine the stability behavior of our approach. If we apply (3.2) and (3.4) to the test problem y 0 = ky for which y00 = k2y, we get k X j¼0
C j ð hÞy nþj ¼ 0;
ð4:1Þ
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
471
where 2
^ Ck ¼ 1 hb h ^ck ; k 2
Cj ¼ aj þ h ^ckþ1 d j ; a0 ak1 d0 ¼ ; A2 aj ak1 aj1 ; dj ¼ A A2
j ¼ 1; 2; . . . ; k 1;
and h ¼ kh;
2 h ck . A¼1 hbk
Therefore, the corresponding characteristic equation of the kth order difference equation of the method is pðn; hÞ ¼
k X
C j nj ¼ 0.
ð4:2Þ
j¼0
To see the zero-stability of this new method, one can easily show that by substituting h ¼ kh ¼ 0 in (4.6) the resulting characteristic polynomial satisfies the root condition and so the method is zero-stable, for example see [9]. To obtain the region of absolute stability we use the boundary locus method. By collecting coefficients of different powers of h in (4.6), we obtain 6 5 4 3 2 h þ A5 h þ A4 h þ A3 h þ A2 h þ A1 h þ A0 ¼ 0; A6
ð4:3Þ ih
where A0, A1, A2, A3, A4, A5, A6 are functions of n. Inserting n = e , Eq. (4.7) gives us six roots hi ðhÞ; i ¼ 1; 2; . . . ; 6 which describe the stability domain. The corresponding (approximate) regions of A(a)-stability found using a numerical approach are given in Table 3. Comparison of Tables 3 and 4 show that regions of A(a)-stability for our new method is larger than those of the other mentioned methods except of E2BDF-Class 1. But the modification of the new method that will be explained in Section 5, called MF-SDMMnew, is overall better than E2BDF.
Table 3 The A(a)-stability of some mentioned methods k
1 2 3 4 5 6
Enright methods
SDBDF methods
E2BDFClass 1
E2BDFClass 2
Ismail methods
p
a ()
p
a ()
p
a ()
p
a ()
p
amax ()
3 4 5 6 7 8
90 90 87.88 82.03 73.10 59.95
2 3 4 5 6 7
90 90 90 89.36 86.35 80.32
4 5 6 7 8 9
90 90 90 90 90 89
4 5 6 7 8 9
90 90 90 89 87 83
2 3 4 5 6 7
90 90 90 89.9 87.3 84.2
472
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
Table 4 A(a)-stability of new SDMM k
p
a ()
k
p
a ()
1 2 3 4 5 6
3 4 5 6 7 8
90 90 90 90 89.8 88.3
7 8 9 10 11 12
9 10 11 12 13 14
85.3 80.5 73.5 61.9 50.3 29.9
L-Stability: From (4.6) the characteristic equation takes the form: c0 y n þ c1 y nþ1 þ þ ck y nþk ¼ 0. That is ck1 ck2 c1 c0 y nþk1 y nþk2 y nþ1 y n . ck ck ck c1 ! 1 that means the method is L-stable or L(a)-stable according to We observe that yn+k ! 0 as h whether the method is A-stable or A(a)-stable, respectively. y nþk ¼
5. Matrix free SDMM To obtain this class of methods, which we call MF-SDMMnew, a hybrid application of the IOM [3] algorithm and the new SDMM will be followed in an efficient way. Now, we examine some of the more important computational aspects associated with our algorithm. We note that in implementing stages (1), (2) and (4) to integrate a nonlinear initial value problem, it is necessary to solve a system of nonlinear algebraic equation for each of the required solutions y nþk , y nþkþ1 and yn+k. In each case, these algebraic equations are solved using a modified form of Newton iteration iterated to convergence. h2 ck og and for Stage (4) the Jacobian In Stages (1) and (2) the Jacobian matrix is I hbk of oy oy of 2 og ^ matrix is I hbk oy h ^ck oy . we consider the Stage (4) as follows: ^ vk Þfnþk h2 ð^c uk Þg ¼ y nþk hðb k k nþk
k1 X
a^j y nþj þ hvk f nþk þ h2 uk gnþk h2^ckþ1 gnþkþ1 .
j¼0
The choose of vk and uk is voluntary and the order of formula is independent of vk and uk. If now ^ b and uk ¼ ^c c , then for Stage (4) we have: we choose vk ¼ b k k k k Stage 4*. Correct yn+k by solution of y nþk hbk fnþk h2 ck gnþk ¼
k1 X
2 ^ b Þf þ h2 ð^c c Þ ^ aj y nþj þ hðb ckþ1 gnþkþ1 . k k k k g nþk h ^ nþk
j¼0
Therefore, the Jacobian matrix in each of the 3 stages (1), (2) and (4*) is the same as h2 ck og . I hbk of oy oy
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
473
In each iterate a linear system Ax = b must be solved in which the matrix A is the same as h2 ck og . If we use IOM to solve this systems, the matrix of A is not needed explicitly, I hbk of oy oy only the action of A times a vector v is necessary. For details see [4]. Let F = bkf + hckg, we have , this means that we require value of product oF v. For this value, approximation is A ¼ I h oF oy oy made using the difference quotient oF F ðxn ; y n þ rvÞ F ðxn ; y n Þ v’ oy r for suitably chosen scalar r. Note that if f(xn, yn) and g(xn, yn) have been saved, then this only requires one additional evaluation of f and g. Comparing with E2BD formulas one can see that the E2BD formulas call for many Jacobian evaluations [8], whereas MF-SDMMnew formulas do not need any Jacobian evaluations. Because the modification used to form MF-SDMMnew from SDMM has an effect on stability, it is necessary to investigate whether or not the stability region for each value of k is made smaller. In fact, for order 1, the original method is A(a)-stable with a is almost 90 whereas the modified method MF-SDMMnew is A-stable. Even though full A-stability is not always attained, there is actually an improvement for all values of k (see Tables 5 and 6).
Table 5 The results of Example 1 x 0.4
40
4000
yi
MF-SDMMnew
SDMMnew 1
9.85172113860776 · 10 3.38639537898003 · 105 1.47940221854341 · 102 7.15827068719389 · 101 9.18553476455750 · 106 2.84163745745846 · 101 1.83202257776709 · 101 8.94237125277590 · 107 8.16796847986166 · 101
y1 y2 y3 y1 y2 y3 y1 y2 y3
9.85172113860778 · 101 3.38639737898005 · 105 1.47940221854326 · 102 7.15827068719401 · 101 9.18553476455793 · 106 2.84163745745835 · 101 1.83202041870106 · 101 8.94235840009368 · 107 8.16797063894054 · 101
Table 6 The results of Example 2 x
yi
MF-SDMMnew
SDMMnew
1
y1 y2 y1 y2 y1 y2 y1 y2
1.865095096003 0.7524845333162 1.898512786032 0.7289766061956 1.786196470324 0.8154281628737 1.504881815732 1.189933309100
1.865095095942 0.7524845333605 1.898529645404 0.7289651615856 1.786316509872 0.8153233443068 1.504845868983 1.190006686423
5 10 20
474
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
6. Numerical results In this section we present some numerical results to compare our new methods SDMMnew and MF-SDMMnew with that of other second derivative multistep methods. The aim of this investigation is to demonstrate numerically the generally superior performance of our new methods and in particular, MF-SDMMnew over SDMMnew with the same stepnumber. Example 1. Consider the Robertson problem: y 01 ¼ 0.04y 1 þ 104 y 2 y 3 ; y 02 ¼ 0.04y 1 104 y 2 y 3 3 107 y 22 ; y 03 ¼ 3 107 y 22 ; with initial value Y(0) = (1, 0, 0)T. We compared the results of SDMMnew and MF-SDMMnew on this problem. Example 2. Consider the van der Pols equation: y 01 ¼ y 2 ; y 02 ¼ l2 ðð1 y 21 Þy 2 y 1 Þ; with initial value y(0) = (2, 0)T. We choose l = 500. We present the numerical solution of this problem using MF-SDMMnew and SDMMnew at some selected points. Example 3. The following stiff initial value problem arose from a chemistry problem: y 01 ¼ 0.013y 2 1000y 1 y 2 2500y 1 y 3 ; y 02 ¼ 0.013y 2 1000y 1 y 2 ; y 03 ¼ 2500y 1 y 3 ; with initial value y(0) = (0, 1, 1)T. We solve this problem at x = 2 and compare the results with those of Ismail methods [7] and SDBDF [2]. The results tabulate in Table 7. The superior performance of the new method is clear. Example 4. Consider the following stiff initial value problem: y 01 ¼ ay 1 by 2 þ ða þ b 1Þex ; y 02 ¼ by 1 ay 2 ða b 1Þex ; with initial value y(0) = (1, 1)T. In order to make this system homogeneous, we introduce an additional variable y3 such that y 03 ¼ 1;
y 3 ð0Þ ¼ 0.
The eigenvalues of the Jacobian associated with the resulting system are a ± ib, 0 and the required solution is y1(x) = y2(x) = ex, y3(x) = x.
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
475
Table 7 The results of Example 3 x
yi
Exact solution
Error in MF-SDMMnew
Error in SDMMnew
Error in Ismail methods
Error in SDBDF
2.0
y1 y2 y3
0.3616933169289 · 105 0.9815029948230 1.018493388244
0.14 · 1018 0.23 · 1013 0.19 · 1012
0.11 · 1018 0.31 · 1013 0.20 · 1012
0.82 · 1010 0.61 · 105 0.57 · 105
0.31 · 108 0.18 · 105 0.57 · 105
Table 8 The results of Example 4 for the case a = 1, b = 15 x
yi
Error in MF-SDMMnew
Error in E2BD
5
y1 y2 y1 y2 y1 y2 y1 y2
0.10 · 1016 0.94 · 1017 0.12 · 1018 0.15 · 1018 0.87 · 1021 0.17 · 1020 0.11 · 1023 0.13 · 1022
0.12 · 109 0.17 · 109 0.30 · 1011 0.10 · 1011 0.84 · 1013 0.12 · 1013 0.15 · 1014 0.32 · 1015
10 15 20
Table 9 The results of Example 4, for the case a = 1, b = 30 x
yi
Error in MF-SDMMnew
Error in E2BD-Class 1
Error in E2BD-Class 2
4.5
y1 y2 y1 y2 y1 y2 y1 y2
0.3 · 1015 0.3 · 1015 0.3 · 1018 0.3 · 1018 0.4 · 1019 0.3 · 1019 0.6 · 1022 0.7 · 1022
<0.1 · 1010 <0.1 · 1010 <0.1 · 1012 <0.1 · 1012 <0.1 · 1015 <0.1 · 1015 <0.1 · 1017 <0.1 · 1017
<0.1 · 1010 <0.1 · 1010 <0.1 · 1012 <0.1 · 1012 0.8 · 1011 0.6 · 1011 0.1 · 1011 0.1 · 1011
9 13.5 18
In Table 8, we give the results obtained for the solution of this problem for the case a = 1, b = 15 at some x with MF-SDMMnew of order 4 and compare the results with second derivative extended backward differentiation formulas (E2BD) that introduced by Cash. In Table 9, the results of this problem are shown with a = 1, b = 30 using the MF-SDMMnew of order 5 and compared with the results of E2BD-Class 1 and E2BD-Class 2 reported by Cash.
Acknowledgement This paper was completed while the corresponding author was visiting The University of Auckland.
476
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
References [1] G. Dahlquist, A special stability problem for linear multistep methods, BIT 3 (1963) 27–43. [2] E. Hairer, G. Wanner, Solving Ordinary Differential Equation II: Stiff and Differential–Algebraic Problems, Springer, Berlin, 1996. [3] Y. Saad, Krylov subspace methods for solving large unsymmetric linear systems, Math. Comp. 37 (1981) 105–126. [4] S.M. Hosseini, G. Hojjati, Matrix free MEBDF method for numerical solution of systems of ODEs, Math. Comp. Model. 29 (1999) 67–77. [5] W.H. Enright, Second derivative multistep methods for stiff ordinary differential equations, SIAM J. Numer. Anal. 11 (1974) 321–331. [6] J.R. Cash, On the integration of stiff systems of ODEs using extended backward differentiation formula, Numer. Math. 34 (2) (1980) 235–246. [7] G. Ismail, I. Ibrahim, New efficient second derivative multistep methods for stiff systems, Appl. Math. Model. 23 (1999) 279–288. [8] J.R. Cash, Second derivative extended backward differentiation formulas for the numerical integration of stiff systems, SIAM J. Numer. Anal. 18 (2) (1981) 21–36. [9] J.D. Lambert, Computational Methods in Ordinary Differential Equations, John Wiley and Sons, 1972.
New second derivative multistep methods for stiff systems G. Hojjati
a,*
, M.Y. Rahimi Ardabili a, S.M. Hosseini
b
a
b
Faculty of Mathematics, Tabriz University, 29 Bahman Street, 5166614766 Tabriz, Iran Department of Mathematics, Tarbiat Modarres University, P.O. Box 14115-175, Tehran, Iran Received 1 May 2004; received in revised form 1 May 2005; accepted 27 June 2005 Available online 24 August 2005
Abstract In this paper we present details of a new class of second derivative multistep methods. Stability analysis is discussed and an improvement in stability region is obtained. With a simple modification we take advantage of calling for the solution of algebraic equations with the same coefficient matrix in each step. Moreover, using IOM to solve the resulting linear systems, the coefficient matrix is not needed explicitly. Some numerical experiments and comparison with several popular codes are given, showing strong superiority of this new class of methods. 2005 Elsevier Inc. All rights reserved. Keywords: Stiff ODEs; Multistep and multiderivative methods; Stability aspects; IOM
1. Introduction In recent years, numerous works have focused on the development of more advanced and efficient methods for stiff problems. A potentially good numerical method for the solution of stiff systems of ODEs must have good accuracy and some reasonably wide region of absolute stability [1]. A-stability requirement puts a severe limitation on the choice of suitable methods for stiff problems. The search for higher order A-stable multistep methods is carried out in the two main directions: *
Corresponding author. Tel.: +98 411 3312649; fax: +98 411 3342102. E-mail address: [email protected] (G. Hojjati).
0307-904X/$ - see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2005.06.007
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
467
• use higher derivatives of the solutions; • throw in additional stages, off-step points, super-future points and like. This leads into the large field general linear methods [2]. In this paper we introduce a new class of second derivative multistep methods (SDMMnew) that have good stability properties. This new class of methods is equipped with an efficient computational algorithm in which for all iterations the final matrix is not changed. Using IOM [3] to solving linear systems in each iteration one dose not require the storage of coefficient matrix in any form. In the same lines of MF-MEBDF [4], the MF-SDMMnew is also introduced which in turn is formed to be very effective and efficient for stiff ODEs.
2. Second derivative multistep methods (SDMM) Let us consider the stiff initial value problem y 0 ðxÞ ¼ f ðx; yðxÞÞ;
yðx0 Þ ¼ y 0 ;
ð2:1Þ m
m
m
on the finite interval I = [x0, xN] where y : ½x0 ; xN ! R and f : ½x0 ; xN R ! R is continuous and differentiable. The general SDMM can be written in the form: k X
aj y nþj ¼ h
j¼0
k X j¼0
bj fnþj þ h2
k X
cj gnþj ;
ð2:2Þ
j¼0
where aj, bj, cj are parameters to be determined and y00 = fx + fy Æ f = : g(x, y). If either bk or ck is nonzero, the formula is implicit. Taylor expansion shows that the method (2.2) is of order p if and only if k X j¼0
aj jq ¼ q
k X
bj jq1 þ qðq þ 1Þ
j¼0
k X
cj jq2 ;
ð2:3Þ
j¼0
with 0 6 q 6 p. Some known important SDMM schemes that will be used for comparison are as follows: • The Enright [5] k-step formulas of order k + 2 which takes the following form: y nþ1 y n ¼ h
k X
bj fnþjkþ1 þ h2 ck gnþ1 .
j¼0
• Second derivative BDF (SDBDF) [2] with the formula ! ! k k k X X X h2 1 1 ri y nþ1 hf nþ1 g ¼ ; j j 2 nþ1 i j¼1 i¼1 j¼1 where $yn+1 = yn+1 yn. This is a class of k-step formulas of order k + 1. • Second derivative extended backward differentiation formulas (E2BD), that was introduced by Cash [6], with the following forms:
468
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
Class 1 predictor : y nþk y nþk1 ¼ h
k X
bj fnþj þ h2 ck gnþk ;
j¼0
corrector : y nþk y nþk1 ¼ h
kþ1 X
fnþj þ h2 ðc g þ c g b j k nþk kþ1 nþkþ1 Þ;
j¼0
and Class 2 predictor : the same as above; corrector : y nþk y nþk1 ¼ h
kþ1 X
fnþj þ h2c g . b j k nþk
j¼0
• Special class of SDMM, introduced by Ismail and Ibrahim [7] of the form k X
aj y nþj ¼ hbk ðfnþk b fnþk1 Þ þ h2 ck ðgnþk c gnþk1 Þ.
j¼0
For b* = 0, c* = 0 this is the same as the SDBDF method.
3. The new SDMM We are going to introduce a new SDMM that takes the following general form: k X
^ fnþk þ h2 ð^c g ^c g ^ aj y nþj ¼ hb k k nþk kþ1 nþkþ1 Þ;
ð3:1Þ
j¼0
ak ¼ 1 and the other coefficients are chosen so that (3.1) has where g(x, y) = y00 = fx + fyf and ^ order k + 2. The coefficients of k-step methods of class (3.1) are given in Table 1, for k = 1, 2, . . ., 8. In this form also we use the one super-future point technique. We design the above form such that in addition to the method have good stability properties and higher order, it is suitable to apply matrix-free technique as we will show in Section 5. Assuming that the solution values yn, yn+1, . . ., yn+k1 are available, the way in which (3.1) is used in practice is by carry out the following computations: Stage 1. Compute y nþk as the solution of k X
aj y nþj ¼ hbk fnþk þ h2 ck gnþk ;
ð3:2Þ
j¼0
where ak = 1 and the other coefficients are chosen so that (3.2) has order k + 1. The coefficients of k-step methods of class (3.2) are given in Table 2, for k = 1, 2, . . ., 8.
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
469
Table 1 Coefficients in (3.1) k
1
2
3
4
5
6
7
8
d ^ b k ^ck ^ckþ1 ^ a0 ^ a1 ^ a2 ^ a3 ^ a4 ^ a5 ^ a6 ^ a7
6 1 4/d 1/d 1
87 78/d 34/d 4/d 9/d 96/d
5239 4350/d 1530/d 108/d 136/d 1161/d 6264/d
149,809 117300/d 35928/d 1728/d 1431/d 12,325/d 52,920/d 191,808/d
7,364,773 5,500,740/d 1,524,600/d 54,000/d 32,148/d 295,875/d 1,279,000/d 3,663,000/d 10,012,500/d
6,399,049 4,595,220/d 1,179,000/d 32,400/d 14,680/d 146124/d 672,975/d 1,941,800/d 4,279,600/d 9,178,380/d
29,993,974,247 20,828,002,860/d 5,020,961,400/d 111,132,000/d 39,877,200/d 430,293,500/d 2,142,550,872/d 6,582,641,625/d 14,259,196,000/d 24,594,437,700/d 45,159,723,000/d
139,632,996,029 94,170,512,520/d 21,552,904,800/d 395,136,000/d 115,643,675/d 1,350,144,000/d 7,286,417,600/d 24,196,635,904/d 55,785,177,000/d 96,751,244,800/d 139,255,256,000/d 219,777,465,600/d
Table 2 Coefficients in (3.2) k
1
2
3
4
5
6
7
8
d bk ck a0 a1 a2 a3 a4 a5 a6 a7
2 1 1/d 1
7 6/d 2/d 1/d 8/d
85 66/d 18/d 4/d 27/d 108/d
415 300/d 72/d 9/d 64/d 216/d 576/d
12,019 8220/d 1800/d 144/d 1125/d 4000/d 9000/d 18,000/d
13,489 8820/d 1800/d 100/d 864/d 3375/d 8000/d 13,500/d 21,600/d
726,301 457,380/d 88,200/d 3600/d 34,300/d 148,276/d 385,875/d 686,000/d 926,100/d 1,234,800/d
3,144,919 1,917,720/d 352,800/d 11,025/d 115,200/d 548,800/d 1,580,544/d 3,087,000/d 4,390,400/d 4,939,200/d 5,644,800/d
Stage 2. Compute y nþkþ1 as the solution of k X
aj y nþjþ1 ¼ hbk fnþkþ1 þ h2 ck gnþkþ1 .
ð3:3Þ
j¼0
Stage 3. Evaluate gnþkþ1 ¼ gðxnþkþ1 ; y nþkþ1 Þ. Stage 4. Compute yn+k as the solution of k X j¼0
^ fnþk þ h2 ð^c g ^c g ^ aj y nþj ¼ hb k k nþk kþ1 nþkþ1 Þ.
ð3:4Þ
470
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
Lemma 3.1. Given that (i) formula (3.2) is of order k + 1, (ii) formula (3.1) is of order k + 2, (iii) the implicit algebraic equations defining y nþk and y nþkþ1 are solved exactly, then scheme (3.4) has order k + 2.
Proof. Suppose the values yn,yn+1, . . ., yn+k1 be exact. From (3.2) we have yðxnþk Þ y nþk ¼ C 1 hkþ2 y ðkþ2Þ ðxnþk Þ þ Oðhkþ3 Þ.
ð3:5Þ
From (3.5) for one super-future point and if we suppose that y(xn+k) = yn+k we have yðxnþkþ1 Þ y nþkþ1 ¼ C 1 hkþ2 y ðkþ2Þ ðxnþkþ1 Þ þ Oðhkþ3 Þ ¼ C 1 hkþ2 y ðkþ2Þ ðxnþk Þ þ C 1 hkþ3 y ðkþ3Þ ðxnþk Þ þ Oðhkþ3 Þ. But since in (3.3) we apply y nþk , we must add the error of ak1 ðyðxnþk Þ y nþk Þ to the above expression. Hence yðxnþkþ1 Þ y nþkþ1 ¼ C 1 ð1 ak1 Þhkþ2 y ðkþ2Þ ðxnþk Þ þ Oðhkþ3 Þ.
ð3:6Þ
If now C2hk+3y(k+3)(xn+k) + O(hk+4) is the defect of formula (3.1), replacing g(xn+k+1, y(xn+k+1)) by gðxnþkþ1 ; y nþkþ1 Þ adds the expression obtained in (3.6) to this error and we obtain yðxnþk Þ y nþk ¼ C 2 hkþ3 y ðkþ3Þ ðxnþk Þ h2ckþ1 gðxnþkþ1 ; yðxnþkþ1 ÞÞ gðxnþkþ1 ; y nþkþ1 Þ . From (3.6) og ðgÞðyðxnþkþ1 Þ y nþkþ1 Þ oy og ¼ ðgÞC 1 ð1 ak1 Þhkþ2 y ðkþ2Þ ðxnþk Þ þ Oðhkþ3 Þ. oy
gðxnþkþ1 ; yðxnþkþ1 ÞÞ gðxnþkþ1 ; y nþkþ1 Þ ¼
This yields
og ðkþ2Þ yðxnþk Þ y nþk ¼ h C2y ðxnþk Þ ðgÞC 1 ð1 ak1 Þh^ckþ1 y ðxnþk Þ þ Oðhkþ4 Þ; oy which shows that the order of scheme (4) is k + 2. h kþ3
ðkþ3Þ
4. Stability analysis We now examine the stability behavior of our approach. If we apply (3.2) and (3.4) to the test problem y 0 = ky for which y00 = k2y, we get k X j¼0
C j ð hÞy nþj ¼ 0;
ð4:1Þ
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
471
where 2
^ Ck ¼ 1 hb h ^ck ; k 2
Cj ¼ aj þ h ^ckþ1 d j ; a0 ak1 d0 ¼ ; A2 aj ak1 aj1 ; dj ¼ A A2
j ¼ 1; 2; . . . ; k 1;
and h ¼ kh;
2 h ck . A¼1 hbk
Therefore, the corresponding characteristic equation of the kth order difference equation of the method is pðn; hÞ ¼
k X
C j nj ¼ 0.
ð4:2Þ
j¼0
To see the zero-stability of this new method, one can easily show that by substituting h ¼ kh ¼ 0 in (4.6) the resulting characteristic polynomial satisfies the root condition and so the method is zero-stable, for example see [9]. To obtain the region of absolute stability we use the boundary locus method. By collecting coefficients of different powers of h in (4.6), we obtain 6 5 4 3 2 h þ A5 h þ A4 h þ A3 h þ A2 h þ A1 h þ A0 ¼ 0; A6
ð4:3Þ ih
where A0, A1, A2, A3, A4, A5, A6 are functions of n. Inserting n = e , Eq. (4.7) gives us six roots hi ðhÞ; i ¼ 1; 2; . . . ; 6 which describe the stability domain. The corresponding (approximate) regions of A(a)-stability found using a numerical approach are given in Table 3. Comparison of Tables 3 and 4 show that regions of A(a)-stability for our new method is larger than those of the other mentioned methods except of E2BDF-Class 1. But the modification of the new method that will be explained in Section 5, called MF-SDMMnew, is overall better than E2BDF.
Table 3 The A(a)-stability of some mentioned methods k
1 2 3 4 5 6
Enright methods
SDBDF methods
E2BDFClass 1
E2BDFClass 2
Ismail methods
p
a ()
p
a ()
p
a ()
p
a ()
p
amax ()
3 4 5 6 7 8
90 90 87.88 82.03 73.10 59.95
2 3 4 5 6 7
90 90 90 89.36 86.35 80.32
4 5 6 7 8 9
90 90 90 90 90 89
4 5 6 7 8 9
90 90 90 89 87 83
2 3 4 5 6 7
90 90 90 89.9 87.3 84.2
472
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
Table 4 A(a)-stability of new SDMM k
p
a ()
k
p
a ()
1 2 3 4 5 6
3 4 5 6 7 8
90 90 90 90 89.8 88.3
7 8 9 10 11 12
9 10 11 12 13 14
85.3 80.5 73.5 61.9 50.3 29.9
L-Stability: From (4.6) the characteristic equation takes the form: c0 y n þ c1 y nþ1 þ þ ck y nþk ¼ 0. That is ck1 ck2 c1 c0 y nþk1 y nþk2 y nþ1 y n . ck ck ck c1 ! 1 that means the method is L-stable or L(a)-stable according to We observe that yn+k ! 0 as h whether the method is A-stable or A(a)-stable, respectively. y nþk ¼
5. Matrix free SDMM To obtain this class of methods, which we call MF-SDMMnew, a hybrid application of the IOM [3] algorithm and the new SDMM will be followed in an efficient way. Now, we examine some of the more important computational aspects associated with our algorithm. We note that in implementing stages (1), (2) and (4) to integrate a nonlinear initial value problem, it is necessary to solve a system of nonlinear algebraic equation for each of the required solutions y nþk , y nþkþ1 and yn+k. In each case, these algebraic equations are solved using a modified form of Newton iteration iterated to convergence. h2 ck og and for Stage (4) the Jacobian In Stages (1) and (2) the Jacobian matrix is I hbk of oy oy of 2 og ^ matrix is I hbk oy h ^ck oy . we consider the Stage (4) as follows: ^ vk Þfnþk h2 ð^c uk Þg ¼ y nþk hðb k k nþk
k1 X
a^j y nþj þ hvk f nþk þ h2 uk gnþk h2^ckþ1 gnþkþ1 .
j¼0
The choose of vk and uk is voluntary and the order of formula is independent of vk and uk. If now ^ b and uk ¼ ^c c , then for Stage (4) we have: we choose vk ¼ b k k k k Stage 4*. Correct yn+k by solution of y nþk hbk fnþk h2 ck gnþk ¼
k1 X
2 ^ b Þf þ h2 ð^c c Þ ^ aj y nþj þ hðb ckþ1 gnþkþ1 . k k k k g nþk h ^ nþk
j¼0
Therefore, the Jacobian matrix in each of the 3 stages (1), (2) and (4*) is the same as h2 ck og . I hbk of oy oy
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
473
In each iterate a linear system Ax = b must be solved in which the matrix A is the same as h2 ck og . If we use IOM to solve this systems, the matrix of A is not needed explicitly, I hbk of oy oy only the action of A times a vector v is necessary. For details see [4]. Let F = bkf + hckg, we have , this means that we require value of product oF v. For this value, approximation is A ¼ I h oF oy oy made using the difference quotient oF F ðxn ; y n þ rvÞ F ðxn ; y n Þ v’ oy r for suitably chosen scalar r. Note that if f(xn, yn) and g(xn, yn) have been saved, then this only requires one additional evaluation of f and g. Comparing with E2BD formulas one can see that the E2BD formulas call for many Jacobian evaluations [8], whereas MF-SDMMnew formulas do not need any Jacobian evaluations. Because the modification used to form MF-SDMMnew from SDMM has an effect on stability, it is necessary to investigate whether or not the stability region for each value of k is made smaller. In fact, for order 1, the original method is A(a)-stable with a is almost 90 whereas the modified method MF-SDMMnew is A-stable. Even though full A-stability is not always attained, there is actually an improvement for all values of k (see Tables 5 and 6).
Table 5 The results of Example 1 x 0.4
40
4000
yi
MF-SDMMnew
SDMMnew 1
9.85172113860776 · 10 3.38639537898003 · 105 1.47940221854341 · 102 7.15827068719389 · 101 9.18553476455750 · 106 2.84163745745846 · 101 1.83202257776709 · 101 8.94237125277590 · 107 8.16796847986166 · 101
y1 y2 y3 y1 y2 y3 y1 y2 y3
9.85172113860778 · 101 3.38639737898005 · 105 1.47940221854326 · 102 7.15827068719401 · 101 9.18553476455793 · 106 2.84163745745835 · 101 1.83202041870106 · 101 8.94235840009368 · 107 8.16797063894054 · 101
Table 6 The results of Example 2 x
yi
MF-SDMMnew
SDMMnew
1
y1 y2 y1 y2 y1 y2 y1 y2
1.865095096003 0.7524845333162 1.898512786032 0.7289766061956 1.786196470324 0.8154281628737 1.504881815732 1.189933309100
1.865095095942 0.7524845333605 1.898529645404 0.7289651615856 1.786316509872 0.8153233443068 1.504845868983 1.190006686423
5 10 20
474
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
6. Numerical results In this section we present some numerical results to compare our new methods SDMMnew and MF-SDMMnew with that of other second derivative multistep methods. The aim of this investigation is to demonstrate numerically the generally superior performance of our new methods and in particular, MF-SDMMnew over SDMMnew with the same stepnumber. Example 1. Consider the Robertson problem: y 01 ¼ 0.04y 1 þ 104 y 2 y 3 ; y 02 ¼ 0.04y 1 104 y 2 y 3 3 107 y 22 ; y 03 ¼ 3 107 y 22 ; with initial value Y(0) = (1, 0, 0)T. We compared the results of SDMMnew and MF-SDMMnew on this problem. Example 2. Consider the van der Pols equation: y 01 ¼ y 2 ; y 02 ¼ l2 ðð1 y 21 Þy 2 y 1 Þ; with initial value y(0) = (2, 0)T. We choose l = 500. We present the numerical solution of this problem using MF-SDMMnew and SDMMnew at some selected points. Example 3. The following stiff initial value problem arose from a chemistry problem: y 01 ¼ 0.013y 2 1000y 1 y 2 2500y 1 y 3 ; y 02 ¼ 0.013y 2 1000y 1 y 2 ; y 03 ¼ 2500y 1 y 3 ; with initial value y(0) = (0, 1, 1)T. We solve this problem at x = 2 and compare the results with those of Ismail methods [7] and SDBDF [2]. The results tabulate in Table 7. The superior performance of the new method is clear. Example 4. Consider the following stiff initial value problem: y 01 ¼ ay 1 by 2 þ ða þ b 1Þex ; y 02 ¼ by 1 ay 2 ða b 1Þex ; with initial value y(0) = (1, 1)T. In order to make this system homogeneous, we introduce an additional variable y3 such that y 03 ¼ 1;
y 3 ð0Þ ¼ 0.
The eigenvalues of the Jacobian associated with the resulting system are a ± ib, 0 and the required solution is y1(x) = y2(x) = ex, y3(x) = x.
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
475
Table 7 The results of Example 3 x
yi
Exact solution
Error in MF-SDMMnew
Error in SDMMnew
Error in Ismail methods
Error in SDBDF
2.0
y1 y2 y3
0.3616933169289 · 105 0.9815029948230 1.018493388244
0.14 · 1018 0.23 · 1013 0.19 · 1012
0.11 · 1018 0.31 · 1013 0.20 · 1012
0.82 · 1010 0.61 · 105 0.57 · 105
0.31 · 108 0.18 · 105 0.57 · 105
Table 8 The results of Example 4 for the case a = 1, b = 15 x
yi
Error in MF-SDMMnew
Error in E2BD
5
y1 y2 y1 y2 y1 y2 y1 y2
0.10 · 1016 0.94 · 1017 0.12 · 1018 0.15 · 1018 0.87 · 1021 0.17 · 1020 0.11 · 1023 0.13 · 1022
0.12 · 109 0.17 · 109 0.30 · 1011 0.10 · 1011 0.84 · 1013 0.12 · 1013 0.15 · 1014 0.32 · 1015
10 15 20
Table 9 The results of Example 4, for the case a = 1, b = 30 x
yi
Error in MF-SDMMnew
Error in E2BD-Class 1
Error in E2BD-Class 2
4.5
y1 y2 y1 y2 y1 y2 y1 y2
0.3 · 1015 0.3 · 1015 0.3 · 1018 0.3 · 1018 0.4 · 1019 0.3 · 1019 0.6 · 1022 0.7 · 1022
<0.1 · 1010 <0.1 · 1010 <0.1 · 1012 <0.1 · 1012 <0.1 · 1015 <0.1 · 1015 <0.1 · 1017 <0.1 · 1017
<0.1 · 1010 <0.1 · 1010 <0.1 · 1012 <0.1 · 1012 0.8 · 1011 0.6 · 1011 0.1 · 1011 0.1 · 1011
9 13.5 18
In Table 8, we give the results obtained for the solution of this problem for the case a = 1, b = 15 at some x with MF-SDMMnew of order 4 and compare the results with second derivative extended backward differentiation formulas (E2BD) that introduced by Cash. In Table 9, the results of this problem are shown with a = 1, b = 30 using the MF-SDMMnew of order 5 and compared with the results of E2BD-Class 1 and E2BD-Class 2 reported by Cash.
Acknowledgement This paper was completed while the corresponding author was visiting The University of Auckland.
476
G. Hojjati et al. / Applied Mathematical Modelling 30 (2006) 466–476
References [1] G. Dahlquist, A special stability problem for linear multistep methods, BIT 3 (1963) 27–43. [2] E. Hairer, G. Wanner, Solving Ordinary Differential Equation II: Stiff and Differential–Algebraic Problems, Springer, Berlin, 1996. [3] Y. Saad, Krylov subspace methods for solving large unsymmetric linear systems, Math. Comp. 37 (1981) 105–126. [4] S.M. Hosseini, G. Hojjati, Matrix free MEBDF method for numerical solution of systems of ODEs, Math. Comp. Model. 29 (1999) 67–77. [5] W.H. Enright, Second derivative multistep methods for stiff ordinary differential equations, SIAM J. Numer. Anal. 11 (1974) 321–331. [6] J.R. Cash, On the integration of stiff systems of ODEs using extended backward differentiation formula, Numer. Math. 34 (2) (1980) 235–246. [7] G. Ismail, I. Ibrahim, New efficient second derivative multistep methods for stiff systems, Appl. Math. Model. 23 (1999) 279–288. [8] J.R. Cash, Second derivative extended backward differentiation formulas for the numerical integration of stiff systems, SIAM J. Numer. Anal. 18 (2) (1981) 21–36. [9] J.D. Lambert, Computational Methods in Ordinary Differential Equations, John Wiley and Sons, 1972.