- Email: [email protected]
Special multistep methods based on numerical differentiation for solving the initial value problem
Applied Mathematics and Computation 181 (2006) 500–510 www.elsevier.com/locate/amc
Special multistep methods based on numerical differentiation for solving the initial value problem P.S. Rama Chandra Rao Department of Mathematics, Kakatiya Institute of Technology and Science, Warangal 506015, Andhra Pradesh, India
Abstract In this paper, methods based on numerical differentiation and their regions of absolute stability of some classes of special multistep methods are derived. Some of the methods given in Henrici [P. Henrici, Discrete Variable Methods in Ordinary Differential Equations, Wiley, 1962] and some methods derived in this paper are applied to solve a second-order initial value problem hence establishing the superiority of the methods derived in this paper. Ó 2006 Elsevier Inc. All rights reserved. Keywords: Multistep methods; Special multistep methods; Numerical differentiation; Initial value problem; Absolute stability
1. Introduction Special multistep methods based on numerical integration such as Adams-Bash forth methods, AdamsMoulton methods and methods based on numerical differentiation for solving first-order differential equations have been derived in Henrici [3] and Gear [1]. The methods based on numerical differentiation for first-order differential equations have been shown to be stiffly stable by Gear [1], High order stiffly stable methods were considered by Jain [4]. Numerical differentiation methods with an off-step point were considered by Gragg and Statter [2]. High order numerical differentiation type formulas with one or more off -step points were derived in [5]. The motivation for the work carried out in this paper arises from the methods based on numerical differentiation for the first-order differential equations, Special multistep methods based on numerical integration for the solution of the special second-order differential equations have been derived in Henrici [3]. Method of reduction of order for solving singularly perturbed two point boundary value problems was considered in [7]. Generalized multistep methods for special second-order differential equations were considered in [6]. The methods based on numerical differentiation for the special second-order differential equation have been derived in this paper. In Henrici [3] methods based on Numerical Integration have been derived by integrating Y00 = f(x, y) twice and replacing the function f(x, y) by an interpolating polynomial. Some of the special methods that arise are Stormer’s methods, Cowell’s methods, etc. Special multistep methods have been derived by replacing Y(x) on the left hand side of Y00 (x) = f(x, y) by an interpolating polynomial and differentiating it E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.12.063
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
501
twice. We have studied a class of explicit methods and the resulting implicit methods. It is found that the explicit and implicit methods have order (k 1). The regions of absolute stability of the methods are derived. 2. General linear multistep methods for special second-order differential equations The special second-order differential equation Y 00 ¼ f ðx; yÞ;
Y 00 ð0Þ ¼ Y 00 .
Y ð0Þ ¼ Y 0 ;
ð1Þ
Often arises in a number of applications such as mechanical problems without dissipation. A general linear multistep method of step number k for the numerical solution of (1) is give by Y nþ1 ¼
k X
aj y nþ1j þ h2
j¼1
k X
ð2Þ
bj y nþ1j ;
j¼0
where aj, bj are constants and h is the step size. Introducing the polynomials qðnÞ ¼ nk
k X
aj nk1
and
rðnÞ ¼
j¼1
k X
bj nkj .
ð3Þ
j¼1
Eq. (2) can be written as qðEÞy nkþ1 h2 rðEÞy 00nkþ1 ¼ 0;
ð4Þ
where E is the shift operator defined by E(yn) = yn+1. Applying (4) to y00 = ky, we get the characteristic equation qðnÞ hrðnÞ ¼ 0;
where h ¼ kh2 .
ð5Þ
h are in general, complex and the region of absolute stability The roots ni of the characteristic equation (5) and is defined to be the region of the complex h-plane such that the roots of the characteristic equation (5) lie within the unit circle whenever h lies in the interior of the region. If we denote the region of absolute stability of R and its boundary by oR, then the locus of oR is given by hðhÞ ¼ qðeih Þ=rðeih Þ;
0 6 h 6 2p.
ð6Þ
3. Derivation of the methods Let p(x) be the backward difference interpolating polynomial of Y(x) at the (k + 1) abscissas xn+1, xn, . . . , xnk+1. Then p(x) is given by pðxÞ ¼
k X m ð1Þm ðs m Þr y nþ1 ;
s ¼ ðx xnþ1 Þ=h.
ð7Þ
m¼0
Differentiating (7) twice with respect to x, we get X k 1 d2 m 00 p ðxÞ ¼ ½ð1Þm ðs m Þr y nþ1 . h2 m¼0 ds2 Replacing y00 (x) by p00 (x) in Eq. (1) and putting x = xn+1r, i.e. s = r, we get k X
dr;m rm y nþ1 ¼ h2 fnþ1r ;
ð8Þ
d2 m s ð1Þ ðm Þ . 2 ds
ð9Þ
m¼0
where dr;m ¼
502
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
4. Generating function for the coefficients dr,m Let Dr;t ¼
1 X
dr;m tm .
ð10Þ
m¼0
Substituting (9) in (10) and simplifying, we get s
2
Dr;t ¼ ð1 tÞ ½logð1 tÞ at s ¼ r; 1 X s 2 ) dr;m tm ¼ ð1 tÞ ½logð1 tÞ at s ¼ r.
ð11Þ
m¼0
5. Explicit methods Putting r = 1 in (8), we get a class of explicit methods given by k X
d1;m rm y nþ1 ¼ h2 fn .
ð12Þ
m¼0
From (11) it follows that d1,m is the coefficient of tm in the expansion of (1 t)[Log(1 t)]2 in powers of t. The coefficients d1,m are tabulated in Table 1. Differences in (12) are expressed in terms of function values. Eq. (12) takes the form k X
aj y nþ1j ¼ h2 fn .
ð13Þ
j¼0
The coefficients aj are tabulated in Table 2. It can be seen that the local truncation error of the formula (13) is given by LTE ¼ d1;kþ1 hkþ1 y ðkþ1Þ ðgÞ.
ð14Þ
Table 1 Coefficients d1,m : m = 0 (1) 7 m
0
1
2
3
4
5
6
7
d1,m
0
0
1
0
1/12
1/12
13/180
11/180
Table 2 Coefficients aj : j = 0 (1) k, k = 2 (1) 7 k
j 0
1
2
3
4
5
6
7
2 3
1 1 11 12 10 12 137 180 126 180
2 2 20 12 15 12 147 180 70 180
1 0 6 12 4 12 255 180 486 180
4 12 14 12 470 180 855 180
1 12 6 12 285 180 670 180
1 12 93 180 324 180
13 180 90 180
11 180
4 5 6 7
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
503
It follows that the k-step method (14) has order (k 1). When k = 2, 3, the method (14) reduces to the Stormer’s method y nþ1 2y n þ y n1 ¼ h2 fn ;
ð15Þ
which is absolutely stable for h 2 [4, 0]. The regions of absolute stability for k = 4, 5, 6 and 7 are shown in Figs. 1–4, respectively. (Taking real part on x-axis and imaginary part on y-axis.) The region of absolute stability is the region lying outside the boundary.
Fig. 1. r = 1, k = 4.
Fig. 2. r = 1, k = 5.
504
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
Fig. 3. r = 1, k = 6.
Fig. 4. r = 1, k = 7.
6. Implicit methods Putting r = 0 in (8), we get the implicit method k X d0;m rm y nþ1 ¼ h2 fnþ1 ;
ð16Þ
m¼0
where from (16), d0,m is the coefficient of tm in the expansion of [log(1 t)]2 in powers of t. The coefficient d0,m are tabulated in Table 3.
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
505
Table 3 Coefficients d0,m; m = 0 (1) 8 m
0
1
2
3
4
5
6
7
8
d0,m
0
0
1
1
11/12
5/6
137/180
7/10
1089/1680
Differences in (16) are expressed in terms of the function values, (16) takes the form k X
aj y nþ1j ¼ h2 fnþ1 .
ð17Þ
j¼0
The coefficients aj are tabulated in Table 4. The coefficients aj in a slightly different form are given in [1] for k = 2 (1) 6. It can be seen that the local truncation error of the formula (17) is given by LTE ¼ d0;kþ1 hkþ1 y ðkþ1Þ ðgÞ.
ð18Þ
Table 4 Coefficients aj, j = 0 (1) k, k = 2 (1) 7 k
j 0
1
2
3
4
5
6
7
2 3
1 2 35 12 45 12 812 180 938 180
2 5 104 12 154 12 3132 180 4014 180
1 4 114 12 214 12 5265 180 7911 180
1 56 12 156 12 5080 180 9490 180
11 12 61 12 2970 180 7380 180
10 12 972 180 3618 180
137 180 1019 180
126 180
4 5 6 7
Fig. 5. r = 0, k = 2.
506
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
It follows that the k-step method (17) has order (k 1). For the method (17), we have qðnÞ ¼
k X
aj nkj
and
rðnÞ ¼ nk .
ð19Þ
j¼0
The regions of absolute stability for k = 2 (1) 7 are shown in Figs. 5–10. (Taking real part on x-axis and imaginary part on y-axis.) The region of absolute stability is the region lying outside the boundary.
Fig. 6. r = 0, k = 3.
Fig. 7. r = 0, k = 4.
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
507
Fig. 8. r = 0, k = 5.
Fig. 9. r = 0, k = 6.
7. Numerical experiment In this section, we have applied two of the methods given in Henrici [3] and two of the methods derived in this paper to solve the differential equation y 00 ¼ 2ex þ y;
yð0Þ ¼ 0;
y 0 ð0Þ ¼ 1
ð20Þ
508
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
Fig. 10. r = 0, k = 7.
in the interval [0, 2] with h = 0.01 and h = 0.02 and the results are given in Tables 5–8. The methods considered are: (1) The explicit third-order Stormer’s method [3]: y nþ1 2y n þ y n1 ¼ ðh2 =12Þð13f n 2f n1 þ fn2 Þ.
ð21Þ
Table 5 Solution by the third-order Stormer and explicit third-order ND with h = 0.01 x
Exact solution
Numerical solution by third-order Stormer
Absolute error
Numerical solution by explicit third-order ND
Absolute error
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0
0.89816689E01 0.15839171E 00 0 0.20408261E 00 0.22477448E 00 0.21782982E 00 0.18003601E 00 0.10754061E 00 0.42206645E 00 0.16060901E 00 0.36787885E 00 0.63328600E 00 0.96521336E 00 0.13733168E 01 0.18686733E 01 0.24639683E 01 0,31737099E 01 0.40144348E 01 0.50050077E 01 0.61668558E 01 0.75243769E 01
0.89754999E01 0.15825027E 00 0.20385867E 00 0.22446233E 00 0.21742368E 00 0.17952093E00 0.10692543E 00 0.49525425E02 0.16146421E 00 0.36886311E 00 0.63440281E 00 0.96643871E 00 0.13745842E 01 0.18699493E 01 0.24652100E 01 0.31748810E 01 0.40155220E 01 0.50059834E 01 0.61676874E 01 0.75250063E 01
0.61690807E04 0.14144182E03 0.22393465E03 0.31214952E03 0.40614605E03 0.50693750E03 0.61517954E03 0.73187798E03 0.85520744E03 0.98425150E03 0.11168122E02 0.12253523E02 0.12674332E02 0.12760162E02 0.12416840E02 0.11711121E02 0.10871887E02 0.97560883E03 0.83160400E03 0.62942505E03
0.89813292E01 0.15836817E 00 0.20401984E 00 0.22465330E 00 0.21763104E 00 0.17973918E 00 0.10712641E 00 0.47717355E02 0.16130012E 00 0.36869931E 00 0.63421696E 00 0.96618986E 00 0.13741579E 01 0.18690939E 01 0.24636812E 01 0.31724329E 01 0.40119104E 01 0.50008917E 01 0.61607552E 01 0.75157928E 01
0.37550926E05 0.23961067E04 0.62644482E04 0.12105703E03 0.19842386E03 0.29736757E03 0.41466951E03 0.55077299E03 0.69087744E03 0.82039833E03 0.93144178E03 0.97608566E03 0.84018707E03 0.41961670E03 0.28800964E03 0.12779236E02 0.25253296E02 0.4117012 E02 0.61016083E02 0.85849762E02
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
509
Table 6 Solution by the third-order Stormer and explicit third-order ND with h = 0.02 x
Exact solution
Numerical solution by third-order Stormer
Absolute error
Numerical solution by explicit third-order ND
Absolute error
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0
0.089816689E01 0.15839171E00 0.20408255E00 0.22477496E 00 0.21782970E 00 0.18003589E 00 0.10754061E 00 0.42206645 E02 0.16060930E00 0.36787927E00 0.63328600E00 0.96521682E 00 0.13733168E 01 0.18686733E 01 0.24639683E 01 0.31737099E 01 0.40144444E 01 0.50050077E 01 0.61668682E 01 0.75243769E 01
0.89714527E01 0.15811831E 00 0.20363516E 00 0.22414744E 00 0.21701640E 00 0.17902732E 00 0.10632610E 00 0.56536905E02 0.16227430E 00 0.36979210E 00 0.63546371E 00 0.96766591E 00 0.13760481E 01 0.18717012E 01 0.24673147E 01 0.31773930E 01 0.40184813E 01 0.50094366E 01 0.61717129E 01 0.75296831E 01
0.10216236E03 0.27340651E03 0.44739246E03 0.62751770E03 0.81330538E03 0.10085702E 02 0.12145042E02 0.14330260E02 0.16649961E02 0.19128323E02 0.21777153E02 0.24490952E02 0.27313232E02 0.30279160E02 0.33464432E02 0.36830902E02 0.40369034E02 0.44288635E02 0.48446655E02 0.53062439E02
0.89814484E01 0.15838271E00 0.20406234E00 0.22473764E00 0.21777177E00 0.17995220E00 0.10742712E 00 0.43677911E02 0.16079080E 00 0.36809212E 00 0.63352615E 00 0.96546608E 00 0.13735361E 01 0.18637868E 01 0.24639130E 01 0.31734142E 01 0.40138359E 01 0.50040121E 01 0.61653681E 01 0.75222578E 01
0.25629997E05 0.93579292E05 0.20027161E04 0.37610531E04 0.57578087E04 0.84340572E04 0.11396408E03 0.14682859E03 0.18155575E03 021320581E03 0.24062395E03 0.24849176E03 0.21839142E03 0.11253357E03 0.56266785E04 0.29659271E03 0.60939789E03 0.99658966E02 0.15010834E02 0.21200180E02
Table 7 Solution by the fourth-order Numerov and the implicit third-order (k = 4) ND with h = 0.01 x
Exact solution
Numerical solution by the Numerov method
Absolute error
Numerical solution by implicit third-order ND
Absolute error
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0
0.89816689E01 0.15839171E 00 0.20408261E 00 0.22477448E 00 0.21782982E 00 0.18003601E 00 0.10754061E 00 0.42206645E02 0.16060901E 00 0.36787885E 00 0.63328600E 00 0.96521336E 00 0.13733168E 01 0.18686733E 01 0.24639683E 01 0.317370993 01 0.40144348E 01 0.50050077E 01 0.61668558E 01 0.75243769E 01
0.89810729E01 0.15837002E 00 0.20403171E 00 0.22467715E 00 0.21766943E 00 0.17979515E 00 0.10720360E 00 0.46671517E 00 0.1611667E 00 0.36853814E 00 0.63402724E 00 0.96597791E 00 0.13739767E 01 0.18690596E 01 0.24639025E 01 0.31729794E 01 0.40127935E 01 0.50021467E 01 0.61624279E 01 0.75179682E 01
0.59604645E05 0.21696091E04 0.50902367E04 0.97334385E04 0.16039610E03 0.24086237E03 0.33700466E03 0.44648722E03 0.5566106 E03 0.65928698E03 0.74124336E03 0.76454878E03 0.65994263E03 0.38623810E03 0.65803528E04 0.73051453E03 0.16412735E02 0.28610229E02 0.44279099E02 0.64086914E02
0.89816146E01 0.15838894E+00 0.20407532E+00 0.22475812E+00 0.21780018E+00 0.17998917E+00 0.10747372E+00 0.43100547E02 0.16072173E+00 0.36801278E+00 0.63343662E+00 0.96537262E+00 0.13734691E+01 0.18687968E+01 0.24640326E+01 0.31736782E+01 0.40142522E+01 0.50045738E+01 0.61660881E+01 0.752231447E+01
0.26077032E06 0.25033951E05 0.76144934E05 0.16659498E04 0.29787421E04 0.46730042E04 0.66801906E04 0.89339912E04 0.11242926E03 0.13333559E03 0.14889240E03 0.15527010E03 0.14853477E03 0.12004375E04 0.57935715E04 0.37193298E04 0.19407272E03 0.44250488E03 0.78535080E03 0.124644523E02
(2) The explicit third-order numerical differentiation (ND) method (13) with k = 4: y nþ1 ¼ ð20=11Þy n ð6=11Þy n1 ð4=11Þy n2 þ ð1=11Þy n3 þ ð12=11Þh2 fn .
ð22Þ
derived in this paper. (3) The implicit fourth-order Numerov’s method [3]: y nþ1 2y n þ y n1 ¼ ðh2 =12Þðfnþ1 þ 10f n þ fn1 Þ.
ð23Þ
510
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
Table 8 Solution by the fourth-order Numerov and the implicit third-order (k = 4) ND with h = 0.02 X
Exact solution
Numerical solution by the Numerov method
Absolute error
Numerical solution by implicit third-order ND
Absolute error
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0
0.89816689E01 0.15839171E 00 0.20408255E 00 0.22477496E 00 0.21782970E 00 0.18003589E 00 0.10754061E 00 0.42206645E02 0.16060930E 00 0.36787927E 00 0.63328600E 00 0.96521682E 00 0.13733168E 01 0.18686733E 01 0.24639683E 01 0.31737099E 01 0.40144444E 01 0.50050077E 01 0.61668682E 01 0.75243769E 01
0.8981746E01 0.15839088E 00 0.20407718E 00 0.22475988E 00 0.21780157E 00 0.17998946E 00 0.10747236E 00 0.43136887E 00 0.16072696E 00 0.36801898E 00 0.63344175E 00 0.96536684E 00 0.13734283E 01 0.18686972E 01 0.24638605E 01 0.31734152E 01 0.40138912E 01 0.50041208E 01 0.61655321E 01 0.75224791E 01
0.77486038E06 0.83446503E06 0.53644180E05 0.15079975E04 0.28133392E04 0.46432018E04 0,68247318E04 0.93024224E04 0.11765957E03 0.13971329E03 0.15574694E03 0.15002489E03 0.11157990E03 0.23841858E04 0.10776520E03 0.29468536E03 0.55313110E03 0.88691711E03 0.13360977E02 0.18978119E02
0.89816376E01 0.15839095E+00 0.20408171E+00 0.22477235E+00 0.21782596E+00 0.18002991E+00 0.10753208E+00 0.42324257E02 0.16062501E+00 0.36790019E+00 0.63331491E+00 0.96525234E+00 0.13733662E+01 0.18687357E+01 0.24640474E+01 0.31738043E+01 0.40145464E+01 0.50051322E+01 0.61669965E+01 0.75245180E+01
0.29802322E07 0.49173832E06 0.12218952E05 0.2428893E05 0.40084124E05 0.59902668E05 0.84415078E05 0.11710916E04 0.15705824E04 0.20742416E04 0.27179718E04 0.34987926E04 0.45657158E04 0.58889389E04 0.72717667E04 0.8893013 E04 0.10013580E03 0.11587143E03 0.12302399E03 0.12683868E03
(4) The implicit third-order ND method (17) with k = 4: y nþ1 ¼ ð104=35Þy n ð114=35Þy n1 þ ð56=35Þy n2 ð11=35Þy n3 þ ð12=35Þh2 fnþ1 .
ð24Þ
derived in this paper.
8. Discussion and conclusion The methods based on numerical integration are found to have closed regions of absolute stability, the methods based on numerical differentiation are found to be absolutely stable outside some closed boundaries. The numerical solution obtained by the numerical differentiation methods derived in this paper is more accurate than those obtained by the corresponding methods given in Henrici [3]. Acknowledgements The author has been thankful to Prof. S.K. Lakshmana Rao for his encouragement for the past two and half decades. References [1] [2] [3] [4] [5]
C.W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice Hall, 1971. W.B. Gragg, H.J. Statter, Generalized multistep predictor–corrector methods, Journal of the ACM 11 (1964) 188–209. P. Henrici, Discrete Variable Methods in Ordinary Differential Equations, Wiley, 1962. M.K. Jain, Numerical Solution of Differential Equations, Wiley Eastern Limited, 1984. C. Prabhakar Rao, S.R.K. Iyengar, High order numerical differentiation type formulas with and off-set point for stiff ordinary differential equations, Journal of the Institute of Mathematics and its Applications 17 (1976) 281–293. [6] C. Prabhakar Rao, S.R.K. Iyengar, A note on generalized multistep methods for special second order differential equations, Journal of Mathematical and Physical sciences 10 (1976) 457–460. [7] Y.N. Reddy, P. Pramod Chakravarthy, Method of reduction of order for solving singularly perturbed two-point boundary value problems, Applied Mathematics and Computation 136 (2003) 27–45.
Special multistep methods based on numerical differentiation for solving the initial value problem P.S. Rama Chandra Rao Department of Mathematics, Kakatiya Institute of Technology and Science, Warangal 506015, Andhra Pradesh, India
Abstract In this paper, methods based on numerical differentiation and their regions of absolute stability of some classes of special multistep methods are derived. Some of the methods given in Henrici [P. Henrici, Discrete Variable Methods in Ordinary Differential Equations, Wiley, 1962] and some methods derived in this paper are applied to solve a second-order initial value problem hence establishing the superiority of the methods derived in this paper. Ó 2006 Elsevier Inc. All rights reserved. Keywords: Multistep methods; Special multistep methods; Numerical differentiation; Initial value problem; Absolute stability
1. Introduction Special multistep methods based on numerical integration such as Adams-Bash forth methods, AdamsMoulton methods and methods based on numerical differentiation for solving first-order differential equations have been derived in Henrici [3] and Gear [1]. The methods based on numerical differentiation for first-order differential equations have been shown to be stiffly stable by Gear [1], High order stiffly stable methods were considered by Jain [4]. Numerical differentiation methods with an off-step point were considered by Gragg and Statter [2]. High order numerical differentiation type formulas with one or more off -step points were derived in [5]. The motivation for the work carried out in this paper arises from the methods based on numerical differentiation for the first-order differential equations, Special multistep methods based on numerical integration for the solution of the special second-order differential equations have been derived in Henrici [3]. Method of reduction of order for solving singularly perturbed two point boundary value problems was considered in [7]. Generalized multistep methods for special second-order differential equations were considered in [6]. The methods based on numerical differentiation for the special second-order differential equation have been derived in this paper. In Henrici [3] methods based on Numerical Integration have been derived by integrating Y00 = f(x, y) twice and replacing the function f(x, y) by an interpolating polynomial. Some of the special methods that arise are Stormer’s methods, Cowell’s methods, etc. Special multistep methods have been derived by replacing Y(x) on the left hand side of Y00 (x) = f(x, y) by an interpolating polynomial and differentiating it E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.12.063
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
501
twice. We have studied a class of explicit methods and the resulting implicit methods. It is found that the explicit and implicit methods have order (k 1). The regions of absolute stability of the methods are derived. 2. General linear multistep methods for special second-order differential equations The special second-order differential equation Y 00 ¼ f ðx; yÞ;
Y 00 ð0Þ ¼ Y 00 .
Y ð0Þ ¼ Y 0 ;
ð1Þ
Often arises in a number of applications such as mechanical problems without dissipation. A general linear multistep method of step number k for the numerical solution of (1) is give by Y nþ1 ¼
k X
aj y nþ1j þ h2
j¼1
k X
ð2Þ
bj y nþ1j ;
j¼0
where aj, bj are constants and h is the step size. Introducing the polynomials qðnÞ ¼ nk
k X
aj nk1
and
rðnÞ ¼
j¼1
k X
bj nkj .
ð3Þ
j¼1
Eq. (2) can be written as qðEÞy nkþ1 h2 rðEÞy 00nkþ1 ¼ 0;
ð4Þ
where E is the shift operator defined by E(yn) = yn+1. Applying (4) to y00 = ky, we get the characteristic equation qðnÞ hrðnÞ ¼ 0;
where h ¼ kh2 .
ð5Þ
h are in general, complex and the region of absolute stability The roots ni of the characteristic equation (5) and is defined to be the region of the complex h-plane such that the roots of the characteristic equation (5) lie within the unit circle whenever h lies in the interior of the region. If we denote the region of absolute stability of R and its boundary by oR, then the locus of oR is given by hðhÞ ¼ qðeih Þ=rðeih Þ;
0 6 h 6 2p.
ð6Þ
3. Derivation of the methods Let p(x) be the backward difference interpolating polynomial of Y(x) at the (k + 1) abscissas xn+1, xn, . . . , xnk+1. Then p(x) is given by pðxÞ ¼
k X m ð1Þm ðs m Þr y nþ1 ;
s ¼ ðx xnþ1 Þ=h.
ð7Þ
m¼0
Differentiating (7) twice with respect to x, we get X k 1 d2 m 00 p ðxÞ ¼ ½ð1Þm ðs m Þr y nþ1 . h2 m¼0 ds2 Replacing y00 (x) by p00 (x) in Eq. (1) and putting x = xn+1r, i.e. s = r, we get k X
dr;m rm y nþ1 ¼ h2 fnþ1r ;
ð8Þ
d2 m s ð1Þ ðm Þ . 2 ds
ð9Þ
m¼0
where dr;m ¼
502
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
4. Generating function for the coefficients dr,m Let Dr;t ¼
1 X
dr;m tm .
ð10Þ
m¼0
Substituting (9) in (10) and simplifying, we get s
2
Dr;t ¼ ð1 tÞ ½logð1 tÞ at s ¼ r; 1 X s 2 ) dr;m tm ¼ ð1 tÞ ½logð1 tÞ at s ¼ r.
ð11Þ
m¼0
5. Explicit methods Putting r = 1 in (8), we get a class of explicit methods given by k X
d1;m rm y nþ1 ¼ h2 fn .
ð12Þ
m¼0
From (11) it follows that d1,m is the coefficient of tm in the expansion of (1 t)[Log(1 t)]2 in powers of t. The coefficients d1,m are tabulated in Table 1. Differences in (12) are expressed in terms of function values. Eq. (12) takes the form k X
aj y nþ1j ¼ h2 fn .
ð13Þ
j¼0
The coefficients aj are tabulated in Table 2. It can be seen that the local truncation error of the formula (13) is given by LTE ¼ d1;kþ1 hkþ1 y ðkþ1Þ ðgÞ.
ð14Þ
Table 1 Coefficients d1,m : m = 0 (1) 7 m
0
1
2
3
4
5
6
7
d1,m
0
0
1
0
1/12
1/12
13/180
11/180
Table 2 Coefficients aj : j = 0 (1) k, k = 2 (1) 7 k
j 0
1
2
3
4
5
6
7
2 3
1 1 11 12 10 12 137 180 126 180
2 2 20 12 15 12 147 180 70 180
1 0 6 12 4 12 255 180 486 180
4 12 14 12 470 180 855 180
1 12 6 12 285 180 670 180
1 12 93 180 324 180
13 180 90 180
11 180
4 5 6 7
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
503
It follows that the k-step method (14) has order (k 1). When k = 2, 3, the method (14) reduces to the Stormer’s method y nþ1 2y n þ y n1 ¼ h2 fn ;
ð15Þ
which is absolutely stable for h 2 [4, 0]. The regions of absolute stability for k = 4, 5, 6 and 7 are shown in Figs. 1–4, respectively. (Taking real part on x-axis and imaginary part on y-axis.) The region of absolute stability is the region lying outside the boundary.
Fig. 1. r = 1, k = 4.
Fig. 2. r = 1, k = 5.
504
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
Fig. 3. r = 1, k = 6.
Fig. 4. r = 1, k = 7.
6. Implicit methods Putting r = 0 in (8), we get the implicit method k X d0;m rm y nþ1 ¼ h2 fnþ1 ;
ð16Þ
m¼0
where from (16), d0,m is the coefficient of tm in the expansion of [log(1 t)]2 in powers of t. The coefficient d0,m are tabulated in Table 3.
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
505
Table 3 Coefficients d0,m; m = 0 (1) 8 m
0
1
2
3
4
5
6
7
8
d0,m
0
0
1
1
11/12
5/6
137/180
7/10
1089/1680
Differences in (16) are expressed in terms of the function values, (16) takes the form k X
aj y nþ1j ¼ h2 fnþ1 .
ð17Þ
j¼0
The coefficients aj are tabulated in Table 4. The coefficients aj in a slightly different form are given in [1] for k = 2 (1) 6. It can be seen that the local truncation error of the formula (17) is given by LTE ¼ d0;kþ1 hkþ1 y ðkþ1Þ ðgÞ.
ð18Þ
Table 4 Coefficients aj, j = 0 (1) k, k = 2 (1) 7 k
j 0
1
2
3
4
5
6
7
2 3
1 2 35 12 45 12 812 180 938 180
2 5 104 12 154 12 3132 180 4014 180
1 4 114 12 214 12 5265 180 7911 180
1 56 12 156 12 5080 180 9490 180
11 12 61 12 2970 180 7380 180
10 12 972 180 3618 180
137 180 1019 180
126 180
4 5 6 7
Fig. 5. r = 0, k = 2.
506
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
It follows that the k-step method (17) has order (k 1). For the method (17), we have qðnÞ ¼
k X
aj nkj
and
rðnÞ ¼ nk .
ð19Þ
j¼0
The regions of absolute stability for k = 2 (1) 7 are shown in Figs. 5–10. (Taking real part on x-axis and imaginary part on y-axis.) The region of absolute stability is the region lying outside the boundary.
Fig. 6. r = 0, k = 3.
Fig. 7. r = 0, k = 4.
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
507
Fig. 8. r = 0, k = 5.
Fig. 9. r = 0, k = 6.
7. Numerical experiment In this section, we have applied two of the methods given in Henrici [3] and two of the methods derived in this paper to solve the differential equation y 00 ¼ 2ex þ y;
yð0Þ ¼ 0;
y 0 ð0Þ ¼ 1
ð20Þ
508
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
Fig. 10. r = 0, k = 7.
in the interval [0, 2] with h = 0.01 and h = 0.02 and the results are given in Tables 5–8. The methods considered are: (1) The explicit third-order Stormer’s method [3]: y nþ1 2y n þ y n1 ¼ ðh2 =12Þð13f n 2f n1 þ fn2 Þ.
ð21Þ
Table 5 Solution by the third-order Stormer and explicit third-order ND with h = 0.01 x
Exact solution
Numerical solution by third-order Stormer
Absolute error
Numerical solution by explicit third-order ND
Absolute error
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0
0.89816689E01 0.15839171E 00 0 0.20408261E 00 0.22477448E 00 0.21782982E 00 0.18003601E 00 0.10754061E 00 0.42206645E 00 0.16060901E 00 0.36787885E 00 0.63328600E 00 0.96521336E 00 0.13733168E 01 0.18686733E 01 0.24639683E 01 0,31737099E 01 0.40144348E 01 0.50050077E 01 0.61668558E 01 0.75243769E 01
0.89754999E01 0.15825027E 00 0.20385867E 00 0.22446233E 00 0.21742368E 00 0.17952093E00 0.10692543E 00 0.49525425E02 0.16146421E 00 0.36886311E 00 0.63440281E 00 0.96643871E 00 0.13745842E 01 0.18699493E 01 0.24652100E 01 0.31748810E 01 0.40155220E 01 0.50059834E 01 0.61676874E 01 0.75250063E 01
0.61690807E04 0.14144182E03 0.22393465E03 0.31214952E03 0.40614605E03 0.50693750E03 0.61517954E03 0.73187798E03 0.85520744E03 0.98425150E03 0.11168122E02 0.12253523E02 0.12674332E02 0.12760162E02 0.12416840E02 0.11711121E02 0.10871887E02 0.97560883E03 0.83160400E03 0.62942505E03
0.89813292E01 0.15836817E 00 0.20401984E 00 0.22465330E 00 0.21763104E 00 0.17973918E 00 0.10712641E 00 0.47717355E02 0.16130012E 00 0.36869931E 00 0.63421696E 00 0.96618986E 00 0.13741579E 01 0.18690939E 01 0.24636812E 01 0.31724329E 01 0.40119104E 01 0.50008917E 01 0.61607552E 01 0.75157928E 01
0.37550926E05 0.23961067E04 0.62644482E04 0.12105703E03 0.19842386E03 0.29736757E03 0.41466951E03 0.55077299E03 0.69087744E03 0.82039833E03 0.93144178E03 0.97608566E03 0.84018707E03 0.41961670E03 0.28800964E03 0.12779236E02 0.25253296E02 0.4117012 E02 0.61016083E02 0.85849762E02
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
509
Table 6 Solution by the third-order Stormer and explicit third-order ND with h = 0.02 x
Exact solution
Numerical solution by third-order Stormer
Absolute error
Numerical solution by explicit third-order ND
Absolute error
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0
0.089816689E01 0.15839171E00 0.20408255E00 0.22477496E 00 0.21782970E 00 0.18003589E 00 0.10754061E 00 0.42206645 E02 0.16060930E00 0.36787927E00 0.63328600E00 0.96521682E 00 0.13733168E 01 0.18686733E 01 0.24639683E 01 0.31737099E 01 0.40144444E 01 0.50050077E 01 0.61668682E 01 0.75243769E 01
0.89714527E01 0.15811831E 00 0.20363516E 00 0.22414744E 00 0.21701640E 00 0.17902732E 00 0.10632610E 00 0.56536905E02 0.16227430E 00 0.36979210E 00 0.63546371E 00 0.96766591E 00 0.13760481E 01 0.18717012E 01 0.24673147E 01 0.31773930E 01 0.40184813E 01 0.50094366E 01 0.61717129E 01 0.75296831E 01
0.10216236E03 0.27340651E03 0.44739246E03 0.62751770E03 0.81330538E03 0.10085702E 02 0.12145042E02 0.14330260E02 0.16649961E02 0.19128323E02 0.21777153E02 0.24490952E02 0.27313232E02 0.30279160E02 0.33464432E02 0.36830902E02 0.40369034E02 0.44288635E02 0.48446655E02 0.53062439E02
0.89814484E01 0.15838271E00 0.20406234E00 0.22473764E00 0.21777177E00 0.17995220E00 0.10742712E 00 0.43677911E02 0.16079080E 00 0.36809212E 00 0.63352615E 00 0.96546608E 00 0.13735361E 01 0.18637868E 01 0.24639130E 01 0.31734142E 01 0.40138359E 01 0.50040121E 01 0.61653681E 01 0.75222578E 01
0.25629997E05 0.93579292E05 0.20027161E04 0.37610531E04 0.57578087E04 0.84340572E04 0.11396408E03 0.14682859E03 0.18155575E03 021320581E03 0.24062395E03 0.24849176E03 0.21839142E03 0.11253357E03 0.56266785E04 0.29659271E03 0.60939789E03 0.99658966E02 0.15010834E02 0.21200180E02
Table 7 Solution by the fourth-order Numerov and the implicit third-order (k = 4) ND with h = 0.01 x
Exact solution
Numerical solution by the Numerov method
Absolute error
Numerical solution by implicit third-order ND
Absolute error
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0
0.89816689E01 0.15839171E 00 0.20408261E 00 0.22477448E 00 0.21782982E 00 0.18003601E 00 0.10754061E 00 0.42206645E02 0.16060901E 00 0.36787885E 00 0.63328600E 00 0.96521336E 00 0.13733168E 01 0.18686733E 01 0.24639683E 01 0.317370993 01 0.40144348E 01 0.50050077E 01 0.61668558E 01 0.75243769E 01
0.89810729E01 0.15837002E 00 0.20403171E 00 0.22467715E 00 0.21766943E 00 0.17979515E 00 0.10720360E 00 0.46671517E 00 0.1611667E 00 0.36853814E 00 0.63402724E 00 0.96597791E 00 0.13739767E 01 0.18690596E 01 0.24639025E 01 0.31729794E 01 0.40127935E 01 0.50021467E 01 0.61624279E 01 0.75179682E 01
0.59604645E05 0.21696091E04 0.50902367E04 0.97334385E04 0.16039610E03 0.24086237E03 0.33700466E03 0.44648722E03 0.5566106 E03 0.65928698E03 0.74124336E03 0.76454878E03 0.65994263E03 0.38623810E03 0.65803528E04 0.73051453E03 0.16412735E02 0.28610229E02 0.44279099E02 0.64086914E02
0.89816146E01 0.15838894E+00 0.20407532E+00 0.22475812E+00 0.21780018E+00 0.17998917E+00 0.10747372E+00 0.43100547E02 0.16072173E+00 0.36801278E+00 0.63343662E+00 0.96537262E+00 0.13734691E+01 0.18687968E+01 0.24640326E+01 0.31736782E+01 0.40142522E+01 0.50045738E+01 0.61660881E+01 0.752231447E+01
0.26077032E06 0.25033951E05 0.76144934E05 0.16659498E04 0.29787421E04 0.46730042E04 0.66801906E04 0.89339912E04 0.11242926E03 0.13333559E03 0.14889240E03 0.15527010E03 0.14853477E03 0.12004375E04 0.57935715E04 0.37193298E04 0.19407272E03 0.44250488E03 0.78535080E03 0.124644523E02
(2) The explicit third-order numerical differentiation (ND) method (13) with k = 4: y nþ1 ¼ ð20=11Þy n ð6=11Þy n1 ð4=11Þy n2 þ ð1=11Þy n3 þ ð12=11Þh2 fn .
ð22Þ
derived in this paper. (3) The implicit fourth-order Numerov’s method [3]: y nþ1 2y n þ y n1 ¼ ðh2 =12Þðfnþ1 þ 10f n þ fn1 Þ.
ð23Þ
510
P.S. Rama Chandra Rao / Applied Mathematics and Computation 181 (2006) 500–510
Table 8 Solution by the fourth-order Numerov and the implicit third-order (k = 4) ND with h = 0.02 X
Exact solution
Numerical solution by the Numerov method
Absolute error
Numerical solution by implicit third-order ND
Absolute error
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0
0.89816689E01 0.15839171E 00 0.20408255E 00 0.22477496E 00 0.21782970E 00 0.18003589E 00 0.10754061E 00 0.42206645E02 0.16060930E 00 0.36787927E 00 0.63328600E 00 0.96521682E 00 0.13733168E 01 0.18686733E 01 0.24639683E 01 0.31737099E 01 0.40144444E 01 0.50050077E 01 0.61668682E 01 0.75243769E 01
0.8981746E01 0.15839088E 00 0.20407718E 00 0.22475988E 00 0.21780157E 00 0.17998946E 00 0.10747236E 00 0.43136887E 00 0.16072696E 00 0.36801898E 00 0.63344175E 00 0.96536684E 00 0.13734283E 01 0.18686972E 01 0.24638605E 01 0.31734152E 01 0.40138912E 01 0.50041208E 01 0.61655321E 01 0.75224791E 01
0.77486038E06 0.83446503E06 0.53644180E05 0.15079975E04 0.28133392E04 0.46432018E04 0,68247318E04 0.93024224E04 0.11765957E03 0.13971329E03 0.15574694E03 0.15002489E03 0.11157990E03 0.23841858E04 0.10776520E03 0.29468536E03 0.55313110E03 0.88691711E03 0.13360977E02 0.18978119E02
0.89816376E01 0.15839095E+00 0.20408171E+00 0.22477235E+00 0.21782596E+00 0.18002991E+00 0.10753208E+00 0.42324257E02 0.16062501E+00 0.36790019E+00 0.63331491E+00 0.96525234E+00 0.13733662E+01 0.18687357E+01 0.24640474E+01 0.31738043E+01 0.40145464E+01 0.50051322E+01 0.61669965E+01 0.75245180E+01
0.29802322E07 0.49173832E06 0.12218952E05 0.2428893E05 0.40084124E05 0.59902668E05 0.84415078E05 0.11710916E04 0.15705824E04 0.20742416E04 0.27179718E04 0.34987926E04 0.45657158E04 0.58889389E04 0.72717667E04 0.8893013 E04 0.10013580E03 0.11587143E03 0.12302399E03 0.12683868E03
(4) The implicit third-order ND method (17) with k = 4: y nþ1 ¼ ð104=35Þy n ð114=35Þy n1 þ ð56=35Þy n2 ð11=35Þy n3 þ ð12=35Þh2 fnþ1 .
ð24Þ
derived in this paper.
8. Discussion and conclusion The methods based on numerical integration are found to have closed regions of absolute stability, the methods based on numerical differentiation are found to be absolutely stable outside some closed boundaries. The numerical solution obtained by the numerical differentiation methods derived in this paper is more accurate than those obtained by the corresponding methods given in Henrici [3]. Acknowledgements The author has been thankful to Prof. S.K. Lakshmana Rao for his encouragement for the past two and half decades. References [1] [2] [3] [4] [5]
C.W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice Hall, 1971. W.B. Gragg, H.J. Statter, Generalized multistep predictor–corrector methods, Journal of the ACM 11 (1964) 188–209. P. Henrici, Discrete Variable Methods in Ordinary Differential Equations, Wiley, 1962. M.K. Jain, Numerical Solution of Differential Equations, Wiley Eastern Limited, 1984. C. Prabhakar Rao, S.R.K. Iyengar, High order numerical differentiation type formulas with and off-set point for stiff ordinary differential equations, Journal of the Institute of Mathematics and its Applications 17 (1976) 281–293. [6] C. Prabhakar Rao, S.R.K. Iyengar, A note on generalized multistep methods for special second order differential equations, Journal of Mathematical and Physical sciences 10 (1976) 457–460. [7] Y.N. Reddy, P. Pramod Chakravarthy, Method of reduction of order for solving singularly perturbed two-point boundary value problems, Applied Mathematics and Computation 136 (2003) 27–45.