- Email: [email protected]
Explicit high order methods for the numerical integration of periodic initial-value problems
APPLIED MAThEmATiCS AND 'CO~1~3,"I'ATtON ELSEVIER
Applied Mathematics and Computation 95 (1998) 15 26
Explicit high order methods for the numerical integration of periodic initial-value problems Ch. Tsitouras a T.E. S i m o s b,. Section of Mathematics, General Sciences Department, National Technical University o/ Athens, Zografou Campus, GR-15780 Athens, Greece b Laboratory of Applied Mathematics and Computers, Department qfSciences, Technical University of Crete, Kounoupidiana, GR-73100 Chania, Crete, Greece
Abstract Two explicit two-step hybrid methods of order 7 and 8 for the numerical integration of second order periodic initial-value problems are developed in this paper. The first of them has a large interval of periodicity and the other a minimal phase-lag. Each of them has seven stages per iteration. Numerical and theoretical results obtained for several well-known problems show the efficiency of the new methods. © 1998 Elsevier Science Inc. All rights reserved. Keywords: Second order differential equations; Oscillating solutions; Interval of periodicity: Explicit methods; Phase-lag; Orbital problems; High-order method
I. Introduction In the last ten years there has been m u c h activity in the a r e a o f n u m e r i c a l s o l u t i o n o f the second o r d e r initial-value p r o b l e m y"=f(x,y),
y(xo)=Yo,
y'(xo)=Yo
(1)
involving second o r d e r o r d i n a r y differential e q u a t i o n s in which the first derivative does n o t a p p e a r explicitly. In this p a p e r we investigate the class o f the a b o v e p r o b l e m s with p e r i o d i c a l solutions. The results o f this activity are
*Corresponding author. Corresponding address: Amfithea - Paleon Faliron, 26 Menelaou Street, GR 175 64 Athens, Greece. E-mail: [email protected]. 0096-3003198l$19.00 © 1998 Elsevier Science Inc. All rights reserved. PII: S0096-3003(97) 10086-8
16
Ch. Tsitouras, T.E. Simos / Appl. Math. Comput. 95 (1998) 15~6
methods which can be applied to many problems in celestial mechanics, quantum mechanical scattering theory, theoretical physics and chemistry, and electronics (see [1,2]). The most important characteristics of an efficient method for the solution of the problem (1) are accuracy and computational efficiency. The development of methods with the above mentioned characteristics is an open problem. One of the most important properties for the numerical solution of general second order differential equations with periodical solution is the algebraic order of the method. Another important new insight is the phase-lag first introduced by Brusa and Nigro [3]. The most widely used technique for the numerical integration of (1) is the Numerov's method, with interval of periodicity (0, 6) and phase-lag of order 4. Many authors [zD10] have developed methods with minimal phase-lag for the solution of general second order differential equations with periodical solutions. All these methods have algebraic orders 4 and 6. The purpose of this paper is to introduce two explicit eighth algebraic order methods. One of them has a large interval of periodicity and the other a minimal phase-lag. Each of them has ten stages per iteration. In Section 2 we will describe the basic theory of the stability and phase-lag of symmetric multistep methods. In Section 3 we will develop the explicit methods with ten stages of order 8 or 7. Finally, in Section 4 an application of the new methods to some well known problems is presented. Theoretical and numerical results show that these new methods are more efficient than the other well known explicit methods.
2. Basic theory To investigate the stability properties of methods for solving the initial-value problem (1) Lambert and Watson [11] introduce the scalar test equation
y" = -w2y
(2)
and the interval of periodicity. When we apply a symmetric two-step method to the scalar test (2) we obtain a difference equation of the form Yn+l - 2C(H)y, + y,, t = O,
(3)
where H = wh, h is the step length, C(H) = B(H)/A(H) where A(H) and B(H) are polynomials in H and y, is the computed approximation to y(nh), n - 0, 1,2,... For the explicit methods, A(H) = 1.
Ch. Tsitouras, TE. Simos / AppL Math. Comput. 95 (1998) 15-26
17
The characteristic equation associated with Eq. (3) is s 2 - 2C(H)s + 1 = O.
(4)
Brusa and Nigro [3] introduced the frequency distortion as an important property of a method for solving special second order initial-value problems. For frequency distortion other authors [4-7,10] use the terms of phase-lag, phase error or dispersion. From now on we use the term phase-lag. Following Coleman [8] when we apply a symmetric two-step method to the scalar test equation y" = - w 2 y then we have the difference equation (3). The characteristic equation associated with (3) is given by (4). The roots of the characteristic equation (4) are denoted as sj and s2, We have the following definitions.
Definition 1 (See [12] and [13]). The method (3) is unconditionally stable if ]Sll ~< 1 and Is2] ~< 1 for all values of wh. Definition 2. Following Lambert and Watson [1 1] we say that the numerical method (3) has an interval of periodicity (0, H02), if, for all H 2 C (0,H~), Sl and s2 satisfy: Sl = e i°(H) and s2 = e i0(H)
(5)
where O(H) is a real function of H. For any method corresponding to the characteristic equation (4) the phase-lag is defined (see [14]) as the leading term in the expansion of t = H - O(H) = H - cos ' [C(H)].
(6)
If the quantity t = O ( H q+l ) as H --+ 0, the order of phase-lag is q. Definition 3 ([15]). The method (3) is P-stable if its interval o f periodicity is
(0, Theorem 1 (For the proof see [13]). A method which has the characteristic equation (4), has an interval o f periodicity (O,H~2), if f o r all H 2 E (O, H2)IC(H)] < 1. Theorem 2 (For the proof see [13]). About the method which has an interval o f periodicity (0, H 2) we can write: cos [0(H)] = C(H), where H 2 E (0, H2). Based on this Coleman [8] arrived at the following remark.
(7)
18
Ch. Tsitouras, T.E. Simos / Appl. Math. Comput. 95 (1998) 15-26
Remark 1. If the phase-lag order is q = 2r then we have: t = c H 2"+1 + O(H 2r+3) ::~ cos(H) - C ( H ) = cos(H) - cos(H - t)
: cH 2r+2 q_ O(H2r+4).
(8)
3. Construction of the new scheme
The construction of the new schemes is based on the following hypotheses. Firstly we put f . = f ( x . , y . ) and we construct an O(h 4) order approximation of Y.+I. Secondly we give an O(h 6) order approximation of y.+l. Then, we construct approximations to Y.±4 which depend on the O(h 6) order approximation to y,,.l and are of order O(h3). Then, we obtain approximations to y.±], which depend on the approximations to y.±l and to y.+l. These approximations are of order O(hT). We produce also approximations to y.±~, which depend on the approximations to y.±,_, y.±3_ and to y~+~. These approximations are also of order O(h7). For the metl~od ~ith extended interval of periodicity and based on the approximations to y.±_3 and y.±: we obtain the O(h 7) order approximation to 4 . . 5 . • 3;, ~. For the method with minimal phase-lag and based on the approximations to v,,+3 and y,,±~ and . on .the O(h 6) order approximation to y.+l we obtain an ~ . O(h 8) order approximation to y.+l. Finally, the methods depend on the approximations to V,,±a and y,,±_~ and to appropriate order approximation to v,,+l. If we have an ~(h 7) order approximation to y.+l then the final method is of seventh order and it was chosen so in order to achieve an enlarged interval of periodicity (Case 1). If we have an O(h s) order approximation to y.+l then the final method is of eighth algebraic order with 14th phase-lag order (Case 2). Consider now the new family of ten stage methods: f~ = f ( x . , y . ) , Y.+I = 2y,, - y,,-i + h2f~, Y,,+l = 2y. - Y.-I + - ( ~
Y,,+l/'> =
a l y . + a2y., , +
+1 + 10J~ + f . - 1
h2(a3 ?n+l ~- a4fn
= bly. + b2y._, +
-t- a5fn--I
f . + , + b4S. +
=
Yn--3/4 : dlYn -~- d2Yn , -t-
, , .
C3
h2(d3
?n+l -it- d4f " + dsfn_ ' -~- d6L+ll2-t-d7~_ll2),
Yn-215 = kiyn + k2yn-I + hZ(k3fn+i + k4f,~ + ksf.-i
Ch. Tsitouras, T.E. Simos / Appl. Math. Comput. 95 (1998) 15-26 Yn+2/5 ~-
19
qlY. + q2Yn 1 + h2(q3 9~n+l + q4f. + qsf.-I
+ q6L+½+ q7L fsn = Slyn + s2yn-I
+ h2(s3 )~n+l+
½+qsj~+3 +
q9L_3),
S4fn + s5fn I
+ s6J~+] + svj~_]+ssJ~+} + s9f,__~), y,+t - 2y, + y, 1 ----h2[r0L+l + rlf, + rof, I (9) Notice that Yn+l is of order O(h 6) because of the so chosen points. Based on the hypotheses described in the first paragraph of this section we have to solve two constraint optimization problems (maximization of the order of the local truncation error (L.T.E.) with the constraints described in the first paragraph of this section). The solution of each problem is the required method. For each case the evaluation of L.T.E. needs tedious algebra. The hybrid Numerov-type methods that have been used in the last 15 years for the construction of high order methods are special forms of two-step R u n g e - K u t t a Nystrom methods. Unfortunately only two-step Runge-Kutta methods have been studied seriously until now in this direction, see Jackiewitz and Tracogna [16]. Since a comprehensive presentation of the related theory is a serious undertaking, beyond the scope of this article, we give some restricted results that apply in our case. Suppose that we have to solve the autonomous system y " - - f ( y ) , where y E R". There is no need for the independent variable because we include it as a parameter in the system with the additional equation x" = 0. Now according to the Runge-Kutta-Nystrom theory, see Hairer et al. [17], we expand all the expressions of the method with respect to h around the central point xn. If the order of the method is p, taking into account that y" = f , y / 3 ) = ( O f / O x ) = ( O f / O y ) ( O y / O x ) = f ' y ' , etc., we arrive at L.T.E. of the form hP+2(ClFl + c2F2 + . . . + csF~) + O(h p+3) where c are numbers and Felementary differentials involving only y' and partial derivatives of f with respect to y. We must note here that assuming a scalar case, compression of different F,'s may occur. For example, F " f ' y '2 for a system of equations since f ' , f " , y' are matrices, but this is not true if we have a simple scalar differential equation. Finally for reasons of simplicity the L.T.E. of each method given in this paper follows for the scalar case. Case 1: The resulting method of the constraint optimization problems (maximization of the order of the L.T.E. with the constraint of maximal non-empty interval of periodicity) is given by
20
Ch. Tsitouras. T.E. Simos / Appl. Math. Comput. 95 (1998) 15-26
Y.+I = 2y,,-
Y,,-I + h2.)¢5,,
y,,+, = 2y,, - y,, . + ~
Y,,-Cl/'2) =
1
,
~(Y'+ Y"
.1],,+1+ 10./],,+/~,
,
h2 l)+3-84 ( 5 "?''+l - 34f,- 19f, l),
1 3
h2
=
Y,,+<~/2/=~( Y,,+Y,, ~)+~-~ ( - f,,+, +42f,, +'U;,_~), 1 h2 37,, /3/4) = ~(Y,, + 3y,, ~) ~ 61440 × (79.f,,+1 + 6541;,- 901f, 1-5176£, ,1'21- 416f,+<1/2)), )3n+(3/4) =
1 a(7yn -
3y. l) + ~
h2
( - 133 f,.~l + 2062> (10)
+847f, t+11872f-,_/~/2 ) + 7112f,,+¢~/2)) , 1 h2 ( 2037055677f.+ 1 Y. (2/51 = ~(3~, +2y. ~) + - 965178655367
1861136501255f, 55037424366211
34855979517/;t 530119492499~_/1/2 ) 139390705562J~+ij/2 ) 15751970659675 8841236437941 28970075634617 36875009511f, i3/4)+ 28242531293f.+(3/4)) 1078401743656 5443566587179 ' 1 7 ~,+(2/5) = 7 ( ~, --
2yn_~) +
h2 (30212893675f,,+ 1 5764995265451.~ 6663260377838 ~ 29306186848430
84986591148j~ 1 118413779356f._(i/2 ) 2425582963f,+(~/2) 5443566587179 5723915870177 10971141089971 399048287965Ji, 13/41 99940981675f.+(3/4)) 4946740775993 1 ~ J' ~.+~ =2~,-N,_~ + h e (5039283792736535~ 224415064868754
70170004296827J~_~ 12993021541469
q 105626254763288~_/3/4 / + 6039158492667f.+i3/4 ) 5994444036373 2329060521985 566592868426776~_c2/5~_ 34143682776539f~+iz/s) ) 22243851111835 3162516129838 '
21
Ch. Tsitouras, T.E. Simos / AppL Math. Comput. 95 (1998) 1576
Y,,+I - 2y. +Y.-i
[59064(fn+i+ f. l)+6381074f.+1286144
h2
-
17040240
(f.+(3/4) +J~-(3/,))+ 3984375(f,+(2/5)+J~
t. (h)
(2/5))}.
= L.T.E = h 9 ( - 0 . 0 0 l f H 2 j f (2) - 0 . 0 0 0 3 4 f ' j 3 f (2)2 - 0.00 l f 2 f ' j f _ O.O0067ft2yt3f(3) _ O.O0067fftj3f(4) --
(3)
O.O00067f'y'Sf(51).
If we apply the resulting method (10) to the scalar test equation (2) we obtain (3) with
A(H)
H2
= 1,B(H)=
l
H4
- ~-~! + 4!
H6 H8 f 6! 8!
821H 10 969H 12 ÷ 844 10! 1193 12!
447H 14 1084 14!"
(11)
Based on Definitions 2 and 3 on Theorem 1 it is easy to see that IC(H)I = IB(H)I < 1
(12)
for all H E (0, 10.4). Case 2: The resulting method of the constraint optimization problems (maximization of the order of the L.T.E. with the constraint of maximal phase-lag order) is given by
r,n+ I = 2 y . - y .
j+h2f,,.
Y,,+-I = 2y. - y . _ l +]-~
+1 + 10J~ + f . - I
1 h2 Yn (1/2)= ~(Yn ~-Yn 1 ) - ~ 3 ~ ( 1
1
)~n-(3/4)
=
5 ?n÷l-
,
34fn - 19fn-1),
h2 h2
7(Y. + 3y. ,) + 6144-----6 ( 7 9 ) . + , + 654f. - 901f. I - 51761~_(1/2/1
416f.+(,/21),
h2
Y.+(3/4) =~. ( % - 3 y . - i ) + 6144~ ( - 133 f.+l + 20622f. + 847f. 1 + 11872f._(1/2) + 7112f.+(1/2/),
Ch. Tsitouras, T.E. Simos I Appl. Math. Comput. 95 (1998) 15-26
22
1 7
Yn+(2/5) = 5 (
(
Yn -- 5 y n - , ) ÷ h2
-
24375737 f.+~ 45937500
81395963f. 29531250
11367553f,, I 5066326J~ (~/2/ 2986558f.+(1/2) 19687500 1640625 984375 +
Yn-(2/5)
=
254618368j~_/3/4 / 33818368Jn+(3/4)) 103359375 + 14765625
1
,2
{ 50099241170683f~+j
~ (3y. + 2y. ,) + n ~ 280726123412437f. 108632178137518
95~-~6]-64--6-ff~
22073902520131fn_1 38729857234995
+ 201956076335441fn_/1/21 66000541848045
263030085667820fn+t~/21 87506559488363
/
52836780893543f,_(3/4 ) 390533175453863J~+/3/41 1 21708860881282 172741337100507 Y.+l =
2y.-y._j- h2 ( -
373483676123807f. 80981605341810
q 44972104451237~ ~+ y.+~) 73238607345031 208224080396663(f._/3/4 ) + 95084935602207
fn+(3/4))
732254540265175(j~-t2/s / ÷ f.+<2/5)) ) q 16711257797879 '
h2 I59064 (f.+~ + f . i ) + 6381074f,, + 1286144 y,,+l - 2y. +yn 1 - 17040240 ~ (~+,~,4,+~ ,~,4,)+ ~984~(j~+,2., +~ ,2.,)] (13)
Ch. Tsitouras, T.E. Simos / Appl. Math. Comput. 95 (1998) 15~6
23
tn (h) = L.T.E = h I° (-0.0000 l f Z f ' 2 f (2) -- 0.000016f'3y'2f (2) - 2.7 x 10-7f3f (2)2 - O.O0036ff'y'2f (2)2 - 0.0001 ly'4f (2)3
-- O.O000043 f 3 f ' f(3) -- O.O00034 ff'2y'2 f(3) _ O.O003 5f 2y'2 f(2) f(3) -- O.O0024f'J4f(E)f (3) -- 0.000001 l f J 4 f (3)2 -- 3.4 x 10 7f4f(4) -- O.O00015 f 2f'y'2 f(4) _ O.O000081f'2y'4 f(4) _ O.O0023 fy'4 f(2) f(4) -- 2.7 x 10 7yt6f(3) f(4) _ O.O000013 f 3 y,2 f(5) _ O.O000057 ff,y,4 f(5)
- O.O00023y'6f(2)f (5) - 6.7 x 10 7f2yt4f(6)
_
_
4.2 X lO-7ftyt6f (6)
-- 9.0 X l o - S f y ' 6 f (71 -- 3.2 X 10-9j8f(8)).
(14)
Notice that special care was taken so the m a x i m u m coefficient is only 3.6 x 10 -4. If we apply the resulting method (10) to the scalar test equation (2) we obtain (3) with H2
A(H)=
1,
B(H)=I--~(-t
H4
4!
H6
H8
6! -~ 8!
H 10
HI2
10! ÷ 12!
4142090H 14 390259648662642' (15) Based on Definitions 2 and 3 and on Theorem 1 it is easy to see that IC(H)I--IB(H)[ < 1
(16)
for all H E (0, 2.89). F r o m Definition 2 and Theorem 3 we have that t = (2.12 x 1 0 - 8 ) H 14
(17)
i.e., the method has a phase-lag of order 14. This is the higher order phase-lag that it can be obtained for this family of methods. Notice that both methods are also the highest algebraic order explicit methods in the literature. 4. Numerical
results - conclusion
In this section we illustrate the new eighth order method by applying it to the numerical solution of three problems. The first is an inhomogeneous problem, the second is the well known Bessel equation and the third is a nonlinear problem.
4.1. I n h o m o g e n e o u s
equation
We consider the following problem y"= 100y+99sinx,
y ( 0 ) = 1,
y ' ( 0 ) = 11,
whose theoretical solution is y(x) = cos (10x) + sin (10x) + sin (x).
(18)
Ch. Tsitouras, ZE. Simos / Appl. Math. Comput. 95 (1998) 15~6
24
Table 1 C o m p a r i s o n o f the e n d - p o i n t g l o b a l e r r o r s in the a p p r o x i m a t i o n s o b t a i n e d b y u s i n g M e t h o d s 1~, (F.Ev. = F u n c t i o n e v a l u a t i o n s ) F.Ev.
Method 1
Method 2
Method 3
Method 4
4000 6000 8000 10 000
2.3 8.3 8.0 1.3
3.7 1.4 1.4 2.3
2.2 8.3 8.3 1.4
2.0 1.5 4.5 3.8
x x × x
10 -7 10 -9 10 -I° 10 10
x x x x
10 10 10 10
7 ~ ~ lit
x x × x
10 v 10 9 10 -I° 10 -Ill
x x x ×
10 -7 10 9 10 -ll 10 -]2
Eq. (18) has been solved numerically for 0 ~
4.2. The Bessel equation We consider the following problem
y" = -
100 +
,
y(1) = -0.2459357644513483,
y'(1) = -0.5576953439142885
(19)
whose theoretical solution is: y(x)---- V / ~ J 0 ( 1 0 x ) .
(20)
Eq. (19) has been solved numerically in order to find the 10th root of the solution which is equal to 32.59406213134967. We use the methods mentioned in the previous example. In Table 2 we present the global error at the 100th root. The steps are chosen so that the function evaluations are of equal size.
4.3. A nonlinear example [4] We consider the following nonlinear problem studied by Chawla and Rao [4] y " + lOOy = siny,
y(O) =Y0 = O,
y'(O) =Y0 -- 1.
(21)
In this example, since there is not an analytical solution, we measure the global errors in computation of y(2Orc) = 0.000392823991.
Ch. Tsitouras, TE. Simos / Appl. Math. Comput. 95 (1998) 15-26
25
Table 2 Comparison of the end-point global errors in the approximations obtained by using Methods 1 4 (F.Ev. = Function evaluations) F.Ev.
Method 1
Method 2
Method 3
Method 4
6000 7000 8000 9000
1.6 4.8 1.7 6.5
2.6 7.1 2.2 8.0
2.2 6.4 2.2 8.5
4.4 7.1 1.5 3.8
x x × ×
10 9 10 Io 10 -t° 10 tl
x x × ×
10 10 10 10
9 Io Ill tl
x x × x
10 9 10 -I° 10 -t° 10-j~
x × × ×
10 l0 10 u 10 -I1 10 I:
Table 3 Comparison of the end-point global errors in the approximations obtained by using Methods 1 4 (F.Ev. = Function evaluations) F.Ev.
Method 1
Method 2
Method 3
Method 4
10 12 14 16
4.9 1.2 3.4 1.2
1.4 3.2 9.2 3.2
1.1 2.9 9.7 3.8
2.4 2.6 3.6 5.8
000 000 000 000
z x x x
10 o 10 9 10-1° 10 Io
x x x x
10 10 10 10
s ~ ~o Io
x x x x
10 10 10 10
s 9 io to
x x x ×
10 -9 10 -II 10 Ji 10 -12
Eq. (21) has been solved numerically for 0 ~
References [I] L.D. Landau, F.M. Lifshitz, Q u a n t u m Mechanics, Pergamon Press, New York, 1965. [2] R.L. Liboff, Introductory Q u a n t u m Mechanics, Addison Wesley, Reading, MA, 1980. [3] L. Brusa, L. Nigro, A one-step method for direct integration of structural dynamic equations, Int. J. Numer. Methods Eng. 15 (1980) 685-699. [4] M.M. Chawla, P.S. Rao, Numerov-type method with minimal phase-lag for the integration of second order periodic initial-value problems II. Explicit method, J. Comput. Appl. Math. 15 (1986) 329 337.
26
Ch. Tsitouras, T.E. Simos / Appl. Math. Comput. 95 (1998) 15-26
[5] M.M. Chawla, P.S. Rao, An explicit sixth-order method with phase-lag of order eight for y" -- f(t,y), J. Comput. Appl. Math. 17 (1987) 365. [6] T.E. Simos, An explicit almost P-stable two-step method with phase-lag or order infinity for the numerical integration of second order initial value problems, Appl. Math. Comput. 49 (1992) 261-268. [7] M.M. Chawla, P.S. Rao, Numerov-type method with minimal phase-lag for the integration of second order periodic initial-value problems, J. Comput. Appl. Math. 11 (1984) 277-281. [8] J.P. Coleman, Numerical methods for y " = f(x,y), via rational approximation for the cosine, IMA J. Numer. Anal. 9 (1989) 145-165. [9] P.J. van der Houwen, B.P. Sommeijer, Explicit Runge-Kutta (-Nystrom) methods with reduced phase errors for computing oscillating solutions, SIAM J. Numer. Anal. 24 (1987) 595 617. [10[ T.E. Simos, Explicit two-step methods with minimal phase-lag for the numerical integration of special second-order initial-value problems and their application to the one-dimensional Schrodinger equation, J. Comput. Appl. Math. 39 (1992) 89-94. [11] J.D. Lambert, i.A. Watson, Symmetric multistep methods for periodic initial-value problems, J. Inst. Math. Applic. 18 (1976) 189 202. [12] R.M. Thomas, Phase properties of high order almost P-stable formulae, BIT 24 (1984) 225238. [13] T.E. Simos, G. Mousadis, Some new Numerov-type methods with minimal phase-lag for the numerical integration of the radial Schrodinger equation, Comput. Math. Appl. 83 (1994) 1145 1153. [14] J.R. Dormand, M. El-Mikkawy, P.J. Prince, High order embedded Runge-Kutta-Nystrom formulae, IMA J. Numer. Anal. 7 (1987) 423~430. [15] E. Hairer, Unconditionally stable methods for the second order differential equations, Numer. Math. 37 (1984) 355. [16] Z. Jackiewicz, S. Tracogna, A general class of two step Runge-Kutta methods for ordinary differential equations, SIAM J. Numer. Anal. 32 (1995) 1390 1427. [17] E. Hairer, S.P. Norsett, G. Wanner, Solving Ordinary Differential Equations I (Nonstiff Problems), second edition, Springer, Berlin, 1993.
Applied Mathematics and Computation 95 (1998) 15 26
Explicit high order methods for the numerical integration of periodic initial-value problems Ch. Tsitouras a T.E. S i m o s b,. Section of Mathematics, General Sciences Department, National Technical University o/ Athens, Zografou Campus, GR-15780 Athens, Greece b Laboratory of Applied Mathematics and Computers, Department qfSciences, Technical University of Crete, Kounoupidiana, GR-73100 Chania, Crete, Greece
Abstract Two explicit two-step hybrid methods of order 7 and 8 for the numerical integration of second order periodic initial-value problems are developed in this paper. The first of them has a large interval of periodicity and the other a minimal phase-lag. Each of them has seven stages per iteration. Numerical and theoretical results obtained for several well-known problems show the efficiency of the new methods. © 1998 Elsevier Science Inc. All rights reserved. Keywords: Second order differential equations; Oscillating solutions; Interval of periodicity: Explicit methods; Phase-lag; Orbital problems; High-order method
I. Introduction In the last ten years there has been m u c h activity in the a r e a o f n u m e r i c a l s o l u t i o n o f the second o r d e r initial-value p r o b l e m y"=f(x,y),
y(xo)=Yo,
y'(xo)=Yo
(1)
involving second o r d e r o r d i n a r y differential e q u a t i o n s in which the first derivative does n o t a p p e a r explicitly. In this p a p e r we investigate the class o f the a b o v e p r o b l e m s with p e r i o d i c a l solutions. The results o f this activity are
*Corresponding author. Corresponding address: Amfithea - Paleon Faliron, 26 Menelaou Street, GR 175 64 Athens, Greece. E-mail: [email protected]. 0096-3003198l$19.00 © 1998 Elsevier Science Inc. All rights reserved. PII: S0096-3003(97) 10086-8
16
Ch. Tsitouras, T.E. Simos / Appl. Math. Comput. 95 (1998) 15~6
methods which can be applied to many problems in celestial mechanics, quantum mechanical scattering theory, theoretical physics and chemistry, and electronics (see [1,2]). The most important characteristics of an efficient method for the solution of the problem (1) are accuracy and computational efficiency. The development of methods with the above mentioned characteristics is an open problem. One of the most important properties for the numerical solution of general second order differential equations with periodical solution is the algebraic order of the method. Another important new insight is the phase-lag first introduced by Brusa and Nigro [3]. The most widely used technique for the numerical integration of (1) is the Numerov's method, with interval of periodicity (0, 6) and phase-lag of order 4. Many authors [zD10] have developed methods with minimal phase-lag for the solution of general second order differential equations with periodical solutions. All these methods have algebraic orders 4 and 6. The purpose of this paper is to introduce two explicit eighth algebraic order methods. One of them has a large interval of periodicity and the other a minimal phase-lag. Each of them has ten stages per iteration. In Section 2 we will describe the basic theory of the stability and phase-lag of symmetric multistep methods. In Section 3 we will develop the explicit methods with ten stages of order 8 or 7. Finally, in Section 4 an application of the new methods to some well known problems is presented. Theoretical and numerical results show that these new methods are more efficient than the other well known explicit methods.
2. Basic theory To investigate the stability properties of methods for solving the initial-value problem (1) Lambert and Watson [11] introduce the scalar test equation
y" = -w2y
(2)
and the interval of periodicity. When we apply a symmetric two-step method to the scalar test (2) we obtain a difference equation of the form Yn+l - 2C(H)y, + y,, t = O,
(3)
where H = wh, h is the step length, C(H) = B(H)/A(H) where A(H) and B(H) are polynomials in H and y, is the computed approximation to y(nh), n - 0, 1,2,... For the explicit methods, A(H) = 1.
Ch. Tsitouras, TE. Simos / AppL Math. Comput. 95 (1998) 15-26
17
The characteristic equation associated with Eq. (3) is s 2 - 2C(H)s + 1 = O.
(4)
Brusa and Nigro [3] introduced the frequency distortion as an important property of a method for solving special second order initial-value problems. For frequency distortion other authors [4-7,10] use the terms of phase-lag, phase error or dispersion. From now on we use the term phase-lag. Following Coleman [8] when we apply a symmetric two-step method to the scalar test equation y" = - w 2 y then we have the difference equation (3). The characteristic equation associated with (3) is given by (4). The roots of the characteristic equation (4) are denoted as sj and s2, We have the following definitions.
Definition 1 (See [12] and [13]). The method (3) is unconditionally stable if ]Sll ~< 1 and Is2] ~< 1 for all values of wh. Definition 2. Following Lambert and Watson [1 1] we say that the numerical method (3) has an interval of periodicity (0, H02), if, for all H 2 C (0,H~), Sl and s2 satisfy: Sl = e i°(H) and s2 = e i0(H)
(5)
where O(H) is a real function of H. For any method corresponding to the characteristic equation (4) the phase-lag is defined (see [14]) as the leading term in the expansion of t = H - O(H) = H - cos ' [C(H)].
(6)
If the quantity t = O ( H q+l ) as H --+ 0, the order of phase-lag is q. Definition 3 ([15]). The method (3) is P-stable if its interval o f periodicity is
(0, Theorem 1 (For the proof see [13]). A method which has the characteristic equation (4), has an interval o f periodicity (O,H~2), if f o r all H 2 E (O, H2)IC(H)] < 1. Theorem 2 (For the proof see [13]). About the method which has an interval o f periodicity (0, H 2) we can write: cos [0(H)] = C(H), where H 2 E (0, H2). Based on this Coleman [8] arrived at the following remark.
(7)
18
Ch. Tsitouras, T.E. Simos / Appl. Math. Comput. 95 (1998) 15-26
Remark 1. If the phase-lag order is q = 2r then we have: t = c H 2"+1 + O(H 2r+3) ::~ cos(H) - C ( H ) = cos(H) - cos(H - t)
: cH 2r+2 q_ O(H2r+4).
(8)
3. Construction of the new scheme
The construction of the new schemes is based on the following hypotheses. Firstly we put f . = f ( x . , y . ) and we construct an O(h 4) order approximation of Y.+I. Secondly we give an O(h 6) order approximation of y.+l. Then, we construct approximations to Y.±4 which depend on the O(h 6) order approximation to y,,.l and are of order O(h3). Then, we obtain approximations to y.±], which depend on the approximations to y.±l and to y.+l. These approximations are of order O(hT). We produce also approximations to y.±~, which depend on the approximations to y.±,_, y.±3_ and to y~+~. These approximations are also of order O(h7). For the metl~od ~ith extended interval of periodicity and based on the approximations to y.±_3 and y.±: we obtain the O(h 7) order approximation to 4 . . 5 . • 3;, ~. For the method with minimal phase-lag and based on the approximations to v,,+3 and y,,±~ and . on .the O(h 6) order approximation to y.+l we obtain an ~ . O(h 8) order approximation to y.+l. Finally, the methods depend on the approximations to V,,±a and y,,±_~ and to appropriate order approximation to v,,+l. If we have an ~(h 7) order approximation to y.+l then the final method is of seventh order and it was chosen so in order to achieve an enlarged interval of periodicity (Case 1). If we have an O(h s) order approximation to y.+l then the final method is of eighth algebraic order with 14th phase-lag order (Case 2). Consider now the new family of ten stage methods: f~ = f ( x . , y . ) , Y.+I = 2y,, - y,,-i + h2f~, Y,,+l = 2y. - Y.-I + - ( ~
Y,,+l/'> =
a l y . + a2y., , +
+1 + 10J~ + f . - 1
h2(a3 ?n+l ~- a4fn
= bly. + b2y._, +
-t- a5fn--I
f . + , + b4S. +
=
Yn--3/4 : dlYn -~- d2Yn , -t-
, , .
C3
h2(d3
?n+l -it- d4f " + dsfn_ ' -~- d6L+ll2-t-d7~_ll2),
Yn-215 = kiyn + k2yn-I + hZ(k3fn+i + k4f,~ + ksf.-i
Ch. Tsitouras, T.E. Simos / Appl. Math. Comput. 95 (1998) 15-26 Yn+2/5 ~-
19
qlY. + q2Yn 1 + h2(q3 9~n+l + q4f. + qsf.-I
+ q6L+½+ q7L fsn = Slyn + s2yn-I
+ h2(s3 )~n+l+
½+qsj~+3 +
q9L_3),
S4fn + s5fn I
+ s6J~+] + svj~_]+ssJ~+} + s9f,__~), y,+t - 2y, + y, 1 ----h2[r0L+l + rlf, + rof, I (9) Notice that Yn+l is of order O(h 6) because of the so chosen points. Based on the hypotheses described in the first paragraph of this section we have to solve two constraint optimization problems (maximization of the order of the local truncation error (L.T.E.) with the constraints described in the first paragraph of this section). The solution of each problem is the required method. For each case the evaluation of L.T.E. needs tedious algebra. The hybrid Numerov-type methods that have been used in the last 15 years for the construction of high order methods are special forms of two-step R u n g e - K u t t a Nystrom methods. Unfortunately only two-step Runge-Kutta methods have been studied seriously until now in this direction, see Jackiewitz and Tracogna [16]. Since a comprehensive presentation of the related theory is a serious undertaking, beyond the scope of this article, we give some restricted results that apply in our case. Suppose that we have to solve the autonomous system y " - - f ( y ) , where y E R". There is no need for the independent variable because we include it as a parameter in the system with the additional equation x" = 0. Now according to the Runge-Kutta-Nystrom theory, see Hairer et al. [17], we expand all the expressions of the method with respect to h around the central point xn. If the order of the method is p, taking into account that y" = f , y / 3 ) = ( O f / O x ) = ( O f / O y ) ( O y / O x ) = f ' y ' , etc., we arrive at L.T.E. of the form hP+2(ClFl + c2F2 + . . . + csF~) + O(h p+3) where c are numbers and Felementary differentials involving only y' and partial derivatives of f with respect to y. We must note here that assuming a scalar case, compression of different F,'s may occur. For example, F " f ' y '2 for a system of equations since f ' , f " , y' are matrices, but this is not true if we have a simple scalar differential equation. Finally for reasons of simplicity the L.T.E. of each method given in this paper follows for the scalar case. Case 1: The resulting method of the constraint optimization problems (maximization of the order of the L.T.E. with the constraint of maximal non-empty interval of periodicity) is given by
20
Ch. Tsitouras. T.E. Simos / Appl. Math. Comput. 95 (1998) 15-26
Y.+I = 2y,,-
Y,,-I + h2.)¢5,,
y,,+, = 2y,, - y,, . + ~
Y,,-Cl/'2) =
1
,
~(Y'+ Y"
.1],,+1+ 10./],,+/~,
,
h2 l)+3-84 ( 5 "?''+l - 34f,- 19f, l),
1 3
h2
=
Y,,+<~/2/=~( Y,,+Y,, ~)+~-~ ( - f,,+, +42f,, +'U;,_~), 1 h2 37,, /3/4) = ~(Y,, + 3y,, ~) ~ 61440 × (79.f,,+1 + 6541;,- 901f, 1-5176£, ,1'21- 416f,+<1/2)), )3n+(3/4) =
1 a(7yn -
3y. l) + ~
h2
( - 133 f,.~l + 2062> (10)
+847f, t+11872f-,_/~/2 ) + 7112f,,+¢~/2)) , 1 h2 ( 2037055677f.+ 1 Y. (2/51 = ~(3~, +2y. ~) + - 965178655367
1861136501255f, 55037424366211
34855979517/;t 530119492499~_/1/2 ) 139390705562J~+ij/2 ) 15751970659675 8841236437941 28970075634617 36875009511f, i3/4)+ 28242531293f.+(3/4)) 1078401743656 5443566587179 ' 1 7 ~,+(2/5) = 7 ( ~, --
2yn_~) +
h2 (30212893675f,,+ 1 5764995265451.~ 6663260377838 ~ 29306186848430
84986591148j~ 1 118413779356f._(i/2 ) 2425582963f,+(~/2) 5443566587179 5723915870177 10971141089971 399048287965Ji, 13/41 99940981675f.+(3/4)) 4946740775993 1 ~ J' ~.+~ =2~,-N,_~ + h e (5039283792736535~ 224415064868754
70170004296827J~_~ 12993021541469
q 105626254763288~_/3/4 / + 6039158492667f.+i3/4 ) 5994444036373 2329060521985 566592868426776~_c2/5~_ 34143682776539f~+iz/s) ) 22243851111835 3162516129838 '
21
Ch. Tsitouras, T.E. Simos / AppL Math. Comput. 95 (1998) 1576
Y,,+I - 2y. +Y.-i
[59064(fn+i+ f. l)+6381074f.+1286144
h2
-
17040240
(f.+(3/4) +J~-(3/,))+ 3984375(f,+(2/5)+J~
t. (h)
(2/5))}.
= L.T.E = h 9 ( - 0 . 0 0 l f H 2 j f (2) - 0 . 0 0 0 3 4 f ' j 3 f (2)2 - 0.00 l f 2 f ' j f _ O.O0067ft2yt3f(3) _ O.O0067fftj3f(4) --
(3)
O.O00067f'y'Sf(51).
If we apply the resulting method (10) to the scalar test equation (2) we obtain (3) with
A(H)
H2
= 1,B(H)=
l
H4
- ~-~! + 4!
H6 H8 f 6! 8!
821H 10 969H 12 ÷ 844 10! 1193 12!
447H 14 1084 14!"
(11)
Based on Definitions 2 and 3 on Theorem 1 it is easy to see that IC(H)I = IB(H)I < 1
(12)
for all H E (0, 10.4). Case 2: The resulting method of the constraint optimization problems (maximization of the order of the L.T.E. with the constraint of maximal phase-lag order) is given by
r,n+ I = 2 y . - y .
j+h2f,,.
Y,,+-I = 2y. - y . _ l +]-~
+1 + 10J~ + f . - I
1 h2 Yn (1/2)= ~(Yn ~-Yn 1 ) - ~ 3 ~ ( 1
1
)~n-(3/4)
=
5 ?n÷l-
,
34fn - 19fn-1),
h2 h2
7(Y. + 3y. ,) + 6144-----6 ( 7 9 ) . + , + 654f. - 901f. I - 51761~_(1/2/1
416f.+(,/21),
h2
Y.+(3/4) =~. ( % - 3 y . - i ) + 6144~ ( - 133 f.+l + 20622f. + 847f. 1 + 11872f._(1/2) + 7112f.+(1/2/),
Ch. Tsitouras, T.E. Simos I Appl. Math. Comput. 95 (1998) 15-26
22
1 7
Yn+(2/5) = 5 (
(
Yn -- 5 y n - , ) ÷ h2
-
24375737 f.+~ 45937500
81395963f. 29531250
11367553f,, I 5066326J~ (~/2/ 2986558f.+(1/2) 19687500 1640625 984375 +
Yn-(2/5)
=
254618368j~_/3/4 / 33818368Jn+(3/4)) 103359375 + 14765625
1
,2
{ 50099241170683f~+j
~ (3y. + 2y. ,) + n ~ 280726123412437f. 108632178137518
95~-~6]-64--6-ff~
22073902520131fn_1 38729857234995
+ 201956076335441fn_/1/21 66000541848045
263030085667820fn+t~/21 87506559488363
/
52836780893543f,_(3/4 ) 390533175453863J~+/3/41 1 21708860881282 172741337100507 Y.+l =
2y.-y._j- h2 ( -
373483676123807f. 80981605341810
q 44972104451237~ ~+ y.+~) 73238607345031 208224080396663(f._/3/4 ) + 95084935602207
fn+(3/4))
732254540265175(j~-t2/s / ÷ f.+<2/5)) ) q 16711257797879 '
h2 I59064 (f.+~ + f . i ) + 6381074f,, + 1286144 y,,+l - 2y. +yn 1 - 17040240 ~ (~+,~,4,+~ ,~,4,)+ ~984~(j~+,2., +~ ,2.,)] (13)
Ch. Tsitouras, T.E. Simos / Appl. Math. Comput. 95 (1998) 15~6
23
tn (h) = L.T.E = h I° (-0.0000 l f Z f ' 2 f (2) -- 0.000016f'3y'2f (2) - 2.7 x 10-7f3f (2)2 - O.O0036ff'y'2f (2)2 - 0.0001 ly'4f (2)3
-- O.O000043 f 3 f ' f(3) -- O.O00034 ff'2y'2 f(3) _ O.O003 5f 2y'2 f(2) f(3) -- O.O0024f'J4f(E)f (3) -- 0.000001 l f J 4 f (3)2 -- 3.4 x 10 7f4f(4) -- O.O00015 f 2f'y'2 f(4) _ O.O000081f'2y'4 f(4) _ O.O0023 fy'4 f(2) f(4) -- 2.7 x 10 7yt6f(3) f(4) _ O.O000013 f 3 y,2 f(5) _ O.O000057 ff,y,4 f(5)
- O.O00023y'6f(2)f (5) - 6.7 x 10 7f2yt4f(6)
_
_
4.2 X lO-7ftyt6f (6)
-- 9.0 X l o - S f y ' 6 f (71 -- 3.2 X 10-9j8f(8)).
(14)
Notice that special care was taken so the m a x i m u m coefficient is only 3.6 x 10 -4. If we apply the resulting method (10) to the scalar test equation (2) we obtain (3) with H2
A(H)=
1,
B(H)=I--~(-t
H4
4!
H6
H8
6! -~ 8!
H 10
HI2
10! ÷ 12!
4142090H 14 390259648662642' (15) Based on Definitions 2 and 3 and on Theorem 1 it is easy to see that IC(H)I--IB(H)[ < 1
(16)
for all H E (0, 2.89). F r o m Definition 2 and Theorem 3 we have that t = (2.12 x 1 0 - 8 ) H 14
(17)
i.e., the method has a phase-lag of order 14. This is the higher order phase-lag that it can be obtained for this family of methods. Notice that both methods are also the highest algebraic order explicit methods in the literature. 4. Numerical
results - conclusion
In this section we illustrate the new eighth order method by applying it to the numerical solution of three problems. The first is an inhomogeneous problem, the second is the well known Bessel equation and the third is a nonlinear problem.
4.1. I n h o m o g e n e o u s
equation
We consider the following problem y"= 100y+99sinx,
y ( 0 ) = 1,
y ' ( 0 ) = 11,
whose theoretical solution is y(x) = cos (10x) + sin (10x) + sin (x).
(18)
Ch. Tsitouras, ZE. Simos / Appl. Math. Comput. 95 (1998) 15~6
24
Table 1 C o m p a r i s o n o f the e n d - p o i n t g l o b a l e r r o r s in the a p p r o x i m a t i o n s o b t a i n e d b y u s i n g M e t h o d s 1~, (F.Ev. = F u n c t i o n e v a l u a t i o n s ) F.Ev.
Method 1
Method 2
Method 3
Method 4
4000 6000 8000 10 000
2.3 8.3 8.0 1.3
3.7 1.4 1.4 2.3
2.2 8.3 8.3 1.4
2.0 1.5 4.5 3.8
x x × x
10 -7 10 -9 10 -I° 10 10
x x x x
10 10 10 10
7 ~ ~ lit
x x × x
10 v 10 9 10 -I° 10 -Ill
x x x ×
10 -7 10 9 10 -ll 10 -]2
Eq. (18) has been solved numerically for 0 ~
4.2. The Bessel equation We consider the following problem
y" = -
100 +
,
y(1) = -0.2459357644513483,
y'(1) = -0.5576953439142885
(19)
whose theoretical solution is: y(x)---- V / ~ J 0 ( 1 0 x ) .
(20)
Eq. (19) has been solved numerically in order to find the 10th root of the solution which is equal to 32.59406213134967. We use the methods mentioned in the previous example. In Table 2 we present the global error at the 100th root. The steps are chosen so that the function evaluations are of equal size.
4.3. A nonlinear example [4] We consider the following nonlinear problem studied by Chawla and Rao [4] y " + lOOy = siny,
y(O) =Y0 = O,
y'(O) =Y0 -- 1.
(21)
In this example, since there is not an analytical solution, we measure the global errors in computation of y(2Orc) = 0.000392823991.
Ch. Tsitouras, TE. Simos / Appl. Math. Comput. 95 (1998) 15-26
25
Table 2 Comparison of the end-point global errors in the approximations obtained by using Methods 1 4 (F.Ev. = Function evaluations) F.Ev.
Method 1
Method 2
Method 3
Method 4
6000 7000 8000 9000
1.6 4.8 1.7 6.5
2.6 7.1 2.2 8.0
2.2 6.4 2.2 8.5
4.4 7.1 1.5 3.8
x x × ×
10 9 10 Io 10 -t° 10 tl
x x × ×
10 10 10 10
9 Io Ill tl
x x × x
10 9 10 -I° 10 -t° 10-j~
x × × ×
10 l0 10 u 10 -I1 10 I:
Table 3 Comparison of the end-point global errors in the approximations obtained by using Methods 1 4 (F.Ev. = Function evaluations) F.Ev.
Method 1
Method 2
Method 3
Method 4
10 12 14 16
4.9 1.2 3.4 1.2
1.4 3.2 9.2 3.2
1.1 2.9 9.7 3.8
2.4 2.6 3.6 5.8
000 000 000 000
z x x x
10 o 10 9 10-1° 10 Io
x x x x
10 10 10 10
s ~ ~o Io
x x x x
10 10 10 10
s 9 io to
x x x ×
10 -9 10 -II 10 Ji 10 -12
Eq. (21) has been solved numerically for 0 ~
References [I] L.D. Landau, F.M. Lifshitz, Q u a n t u m Mechanics, Pergamon Press, New York, 1965. [2] R.L. Liboff, Introductory Q u a n t u m Mechanics, Addison Wesley, Reading, MA, 1980. [3] L. Brusa, L. Nigro, A one-step method for direct integration of structural dynamic equations, Int. J. Numer. Methods Eng. 15 (1980) 685-699. [4] M.M. Chawla, P.S. Rao, Numerov-type method with minimal phase-lag for the integration of second order periodic initial-value problems II. Explicit method, J. Comput. Appl. Math. 15 (1986) 329 337.
26
Ch. Tsitouras, T.E. Simos / Appl. Math. Comput. 95 (1998) 15-26
[5] M.M. Chawla, P.S. Rao, An explicit sixth-order method with phase-lag of order eight for y" -- f(t,y), J. Comput. Appl. Math. 17 (1987) 365. [6] T.E. Simos, An explicit almost P-stable two-step method with phase-lag or order infinity for the numerical integration of second order initial value problems, Appl. Math. Comput. 49 (1992) 261-268. [7] M.M. Chawla, P.S. Rao, Numerov-type method with minimal phase-lag for the integration of second order periodic initial-value problems, J. Comput. Appl. Math. 11 (1984) 277-281. [8] J.P. Coleman, Numerical methods for y " = f(x,y), via rational approximation for the cosine, IMA J. Numer. Anal. 9 (1989) 145-165. [9] P.J. van der Houwen, B.P. Sommeijer, Explicit Runge-Kutta (-Nystrom) methods with reduced phase errors for computing oscillating solutions, SIAM J. Numer. Anal. 24 (1987) 595 617. [10[ T.E. Simos, Explicit two-step methods with minimal phase-lag for the numerical integration of special second-order initial-value problems and their application to the one-dimensional Schrodinger equation, J. Comput. Appl. Math. 39 (1992) 89-94. [11] J.D. Lambert, i.A. Watson, Symmetric multistep methods for periodic initial-value problems, J. Inst. Math. Applic. 18 (1976) 189 202. [12] R.M. Thomas, Phase properties of high order almost P-stable formulae, BIT 24 (1984) 225238. [13] T.E. Simos, G. Mousadis, Some new Numerov-type methods with minimal phase-lag for the numerical integration of the radial Schrodinger equation, Comput. Math. Appl. 83 (1994) 1145 1153. [14] J.R. Dormand, M. El-Mikkawy, P.J. Prince, High order embedded Runge-Kutta-Nystrom formulae, IMA J. Numer. Anal. 7 (1987) 423~430. [15] E. Hairer, Unconditionally stable methods for the second order differential equations, Numer. Math. 37 (1984) 355. [16] Z. Jackiewicz, S. Tracogna, A general class of two step Runge-Kutta methods for ordinary differential equations, SIAM J. Numer. Anal. 32 (1995) 1390 1427. [17] E. Hairer, S.P. Norsett, G. Wanner, Solving Ordinary Differential Equations I (Nonstiff Problems), second edition, Springer, Berlin, 1993.