An international Journal
computers & mathematics
with applications
PERGAMON
Computers and Mathematics with Applications 44 (2002) 83 93 www.elsevier.com/locate/camwa
Stability Analysis of TwoStep RungeKutta Methods for Delay Differential Equations
Z. BARTOSZEWSKI Faculty of Applied Physics and M a t h e m a t i c s Gdafisk University of Technology Narutowicza 11/12, 80952 Gdafisk, Poland zbart @mif. pg. gda. pl Z. JACKIEWICZ* D e p a r t m e n t of M a t h e m a t i c s Arizona State University Tempe, AZ 85287, U.S.A. j ackiewiOmath, la. asu. edu
(Received October 2000; revised and accepted September 2001)
AbstractWe investigate stability properties of twostep RungeKutta methods with respect to the linear test equation
y'(t) = ay(t) + by(t  T), t > O,
u(t) = 9(t),
t C [~, 0],
where a and b are complex parameters. It is known that the solution y(t) to this equation tends to zero as t ~ oc if Ibl <  Re(a). We will show that under some conditions this property is inherited by any Astable twostep RungeKutta method applied on a constrained mesh to delay differential equations, i.e., that the corresponding method is Pstable. @ 2002 Elsevier Science Ltd. All rights reserved. KeywordsTwostep Pstability. RungeKutta methods, Delay differential equations, Absolute stability,
1.
INTRODUCTION
C o n s i d e r t h e i n i t i a l  v a l u e p r o b l e m for d e l a y differential e q u a t i o n ( D D E )
~'(t) = / ( t , y(t) = g(t),
y(t), y(t  ~))~
t e [to, T], t e [to  7, to],
(1.1)
T > 0~ w h e r e g is a specified initial f u n c t i o n a n d f satisfies c o n d i t i o n s w h i c h g u a r a n t e e t h e e x i s t e n c e o f t h e u n i q u e s o l u t i o n y t o (1.1). Such c o n d i t i o n s c a n b e f o u n d in, for e x a m p l e , [14]. The authors wish to express their gratitude to anonymous referees for their helpful remarks. *The work of this author was partially supported by the National Science Foundation under Grant NSF DMS9971164. 08981221/02/$  see front matter (~) 2002 Elsevier Science Ltd. All rights reserved. PII: S08981221 (02)001311 Typeset by ~ 4 ~  T E X
84
Z. BARTOSZEWSKI AND Z. JACKIEWICZ
Many numerical methods have been proposed for the numerical solution of problem (1.1); the recent surveys are given by Jackiewicz and Kwapisz [5] and Zennaro [6]. In this paper, we are concerned with twostep RungeKutta (TSRK) methods which can be defined by
If J+l = f (ti + cjh, UjYh(ti1) + (1  u j ) Y h ( t i )
\
/2
y h ( t i ~ ( h ) = ~ ] ( ( ) y h ( t i _ l ) t (1  ?](~) )Yh@i) [ h ~~ ( V s ( ~ ) K s 1 W s ( ~ ) / ~ ' S + l ) , s=l i = 0,1,...,n1, ~ E (0,1], nh = T  t o , ti = to + i h . Here, ~ is the number of stages, yh(t) is a continuous approximation to y(t), [(~+~ are approximations (possibly of low order) to y'(ti + cjh), and r/(~), vi((), and wi(~), i = 1,2 . . . . , , , are polynomials such that r](0) = 0, vi(0) = 0, and wi(0) = 0. These methods form a subclass of general linear methods introduced by Butcher [7] and could be possibly also referred to as twostep hybrid methods. They generalize kstep collocation methods (with k = 2) for ordinary differential equations (ODEs) studied by Lie and Norsett [8] and Lie [9], and TSRK methods for ODEs investigated by Byrne and Lambert [10], Renaut [11,12], Caira et al. [13], Jackiewicz et al. [14], Jackiewicz and Zennaro [15], and Jackiewicz et al. [16]. The discrete version of these methods (in somewhat different socalled Ynotation) was introduced by Jackiewicz and Tracogna [17] in the context of ODEs. They were filrther investigated by Jackiewicz and Tracogna [18], Butcher and Tracogna [19], Tracogna [20], Tracogna and \Velfert [21], Jackiewicz and Vermiglio [22], Hairer and Wanner [23], and Bartoszewski and Jackiewicz [24]. The variable stepsize continuous TSRK methods for ODEs were investigated by Jackiewicz and Tracogna [25]; they result in the formulation (1.2) when applied with a constant stepsize h to the DDE (1.1). The Astable TSRK methods have been constructed in [14,17,25]. Following [25], we will represent these methods by the following table of the coefficients:
Ul
all a~l
..
al~, a~
bll b~l
Wl(~)
??" . ??
''"
bl~ b,~,
Wlz(~ )
1](~)
vT(~ ")
wT(~)
=
It u
v(~)
vl(~)
v.(~)
We would like to stress that the methods (1.2) given by the above table of coefficients are more general than the methods studied in [8,9] (with k = 2) and in [1016] since for all these methods the coefficient matrix A is identically equal to zero? As illustrated in [17], the presence of extra parameters aij in (1?2) makes it possible to construct highorder methods with relatively few stages, and we refer to [17,24,25J for specific examples. Assuming that h = r / m for some positive integer m we can approximate the delayed arguments yh(ti + cjh  T) = yh(tim + cjh) by 12 Yh(tim + cjh) = rl(cj)Yh(ti ..... l) + (1  ~?(Cj))yh(tim) b h ~ (Vs(Cj)li]S_rn k ws(cj)Ki  r e + l ) , s ,
s=l
j = 1 , 2 , . . . , u. Substituting this relation into (1.2) with ~ = 1, we obtain
K~+I = f
t~ + cjh, u~y~_l + (1  ~j)y~ + h ~
8=1
( a . K : + bSq%~) ,
~jy~r~~ + (1  v~)y~m + h ~
"
(~jq~_.~ + 5 . K ; _ m + i
)) ,
(1.3)
s=l
8~1
Delay Differential Equations i = 1,2 . . . . . n  i, j = 1 , 2 , . . . , u , a . = ~,~(~).
85
where Yi = yh(ti), r] = r~(1), rb = V(cj), 7is = % ( c j ) , and
It is t h e purpose of this paper to investigate stability properties of (1.3) with respect to the linear test equation J ( t ) : ay(t) + by(t  7), t >_ O, (1.4) y(t) = g(t), t e lT, o1, where a and b are complex parameters. Stability properties of R u n g e  K u t t a m e t h o d s with respect to this test equation have been investigated by Koto [26,27], Zennaro [28], in 't Hout and Spijker [29], in 't Hout [301, Gugliehni [31], and Guglielmi and Hairer [32]. Stability properties of linear multistep m e t h o d s for DDEs with respect to (1.4) have been investigated by Cryer [33], Bickart [34,351, W i e d e r h o l t [36], and W a t a n a b e and Roth [37]. It was proved by Barwell [38] t h a t the solution y(t) to (1.4) tends to zero as t ~ oc if
Ibl <  Re(a),
(i.a)
and it was proved by Zennaro [28] t h a t if y(t) ~ 0 as t * ec, then Ibl _<  R e ( a ) . Hence, inequality (1.5) defines the interior of the region of asymptotic stability of equation (1.4). In this paper, we investigate under w h a t conditions the ~ y m p t o t i c stability properties of this test equation are inherited by the numerical a p p r o x i m a t i o n to (1.4) by the T S R K m e t h o d (1.3). T h e a p p r o a c h of this paper is based mainly on the technique proposed by Zennaro [281 in the context of R K m e t h o d s for DDEs. However, in spite of some similarities to the R K case there are also m a n y i m p o r t a n t differences, and the extension of the approach of [28] to T S R K m e t h o d s for D D E s is far from trivial. T h e m a i n result of the p a p e r is T h e o r e m 3.3 in Section 3, which states t h a t under sonic technical conditions the T S R K m e t h o d (1.2) for DDEs is P  s t a b l e if the underlying T S R K m e t h o d for O D E s is Astable. 2. PRELIMINARY RESULTS AND DEFINITIONS
A p p l y i n g (I.3) to the test equation (1.4) with h = 7/m, we obtain
t(~+i = a
(
u j y ~  i + (1  u j ) y i + h E
8=1
a j s K ~ + bj~KJ+i
))
,
+ b r b y i _ m _ t + (1  rJj)yi,~ + h
s=l
(Tj~K;?_ m + 5j~Ki ~ +~ ....
(2.i)
i = 0, 1 , . . . , j = 1 , 2 , . . . , u. P u t t i n g (~ = ha, /3 = hb, 35 = [7~1,... ,7],j]T, F = [~js]j,s=l, A = r5J s l j ,1,,l ~ K~ = [KJ,. " ' ~ K~] r , and e = [1, " ' " 1] T E R ' , equation (2.1) can be w r i t t e n as t s= hK~+l = o~(y~iu + y~(e  u) + h A K i + hBK,s,+i) + ~ (Yi ....  l u + y i  m (e  ~) + h F I ( i _ , , + h A K i . . . . 1), Yi+l = 7]yi1 + (1  ~q)Yi + v T hI(i + w T h K i + l , i = 0, 1 , . . . . Denote by {y.~.(m;ch~)}~=0 the solution to (2.2) with h = r / m . naro [28], we introduce the following definition. Following Zen
11
(2.2)
86
Z. BARTOSZEWSKI AND Z. JACKIEWICZ
DEFINITION. M e t h o d (1.3) is said to be stable for given (a,/3) i f the sequence { y i ( m ; a,/3)}i°°__0 tends to zero as t 4 oc for a n y integer m >_ 1. T h e region S p o f s t a b i l i t y o f (1.3) is the set o f all vaIues (a,/3) for which this m e t h o d is stable. M e t h o d (1.3) is said to be P  s t a b l e K i t s region of s t a b i l i t y S p includes the set {(a,/3) :1/31 <  R e ( a ) } . Observe t h a t it follows from (1.5) t h a t if the m e t h o d (1.3) is P  s t a b l e , then the numerical solution {y~(m; a , 3)}~0 defined by (2.2) tends to zero as i + oo whenever the analytical solution y(t) to (1.4) tends to zero as t + oo. Observe also t h a t a a n d / 3 d e p e n d on the stepsize h. A s s u m i n g t h a t the m a t r i x I  a B is nonsingular we can compute h K i + l from the first relation of (2.2), and s u b s t i t u t i n g the resulting expression into the second relation of (2.2) leads to the following vector recurrence equation: Y,+I = L Y / + MY/1 + N Y i  , ~ + I + RY/m q SY/ . . . . 1, i = 0, 1 , . . . , where Y / = [yi, h K T ] r and the matrices L, M , N , R, and S are defined by L= [1~+awT(Ic~B)l(eu) o~(I  oeB) l(e  u) L 00] v T +awq(IozB)lA] a(I  olB)lA (2.3)
] '
[ rl+OewT(IaB)lu M = L a(I  ctB)lu
N = /3(Z  a B )  I A
J '
[/3wT(I  a B )  l ( e  ~)
t~ = S = L /3(I  o~B)l(e  (,)
3wr(I
 aB)~r]
J '
/3(I  aB)lr "
I/3,.r(_z. ~ B )  I ~
L /3(1aB)i(,
0]
0
T h e characteristic equation of
(2.3)
is (2.4)
det (A'~+2I  Am+lL  M~M  A2N  AR  S) = 0.
Observe t h a t A is a root of this equation if and only if there exists x* E C "+1, x* ? 0, such t h a t Am+2X*  Am+ILx*  A " ~ M x *  A 2 N x *  ARx*  Sx* = 0. (2.5)
P u t x* = [p, xT] q, where p E C and x c C ~, and assume t h a t p ? 0. This seems to h a p p e n for most T S R K methods. Then (2.5) is equivalent to Ara+2  Am+l ((1  r] + c ~ w r ( I  c t B )  i ( e  u)) + (v T + a w T ( I

oLB)IA) x)


A m (rI + a w T ( I


aB)lu)
(2.6)
A2/3~,~(.r  a B )  ~ A x  .X (/3~r (Z  ,~B) ~ (e  ~)
+ / 3 w T ( I  c t B )  I F x )  / 3 w T ( I  OeB)I~ = O,
and Am+2x  A"~+' ( a ( I  a B )  l ( e  u) + a ( I  a S )  l A x )  Ama(I  aB)lu _ A2/3(I  o ~ B )  l A x
(2.7)
 ~ ( Z ( I  a B )  1 (e  ~) + Z ( I  a B )  l r x ) /3(I  aB)l~ = 0.
Consider also the following relation: A 2 w r x  A 2 + A (1  ~ + v r x ) +7] = O,
(2.8)
Delay Differential Equations
87
and the q u a d r a t i c function
(2.9)
where the rational functions
R l (oz, z) = 1  rl + v T x + wq(I  ctB  zA) 1 (o~(e  u) + z (e  *2) + oeAx + z F x )
and
R2(OZ, Z ) ~ $1 ~ w T ( /   (_~B  Z A )  1 (6t% ~ Z*2)
are well defined if the m a t r i x I  a B  z A is nonsingular. We have the following theorem. THEOREM 2.1. A s s u m e t h a t the m a t r i x I  a B  z A is singular if and only i f z is a pole o f R1 (a, z) or R 2 ( a , z ) . Then A ? 0 such that 3 / A m is not a pole o f R l ( a , z ) , and R 2 ( a , z ) satisfies (2.6) and (2.7) if and only i r a is a root o f (2.8) and (2.9). P R o o e OF NECESSITY. Assume t h a t A ? 0 such t h a t / 3 / A m is not a pole of R l ( a , z), and R 2 ( a , z) satisfies (2.6) and (2.7). Then it follows from the assumptions of the theorem t h a t the maZrix I  c~B  (/3/Am)A is nonsingular. Consider the relation
Am+2(I  eeB)x  A m + l ( a ( e  u) + ceAx)  AmCtU  A2/3AX  A (3 (e  *2) + ~ r x )  3*2 = 0,
which is equivalent to (2.7). Dividing this relation by A~ and then multiplying it by w T ( I a B  (/3/Am)A) 1 we o b t a i n
A2wTx.XwT (Io:B  AI3T x (oz(e u) + A/3T2,~ *2)+ oeAx+ A/3~ A)' (ePx)
(2.10)
To c o m p u t e A 2 w T x , we multiply (2.7) by W T and then compare the resulting relation with equation (2.6). This leads to
A2wmz = A2  k (1  r] + v T x )  r],
which is equivalent to (2.8). Substituting this equation into (2.10), we o b t a i n (2.9). This completes the proof of necessity.  PROOF OF SUFFICIENCY. From (2.9), it follows t h a t
,,,2 _ ,,, (1  ~ + vrx)  ~j
= ..,,S (,,  ~B  ~#dA)' (o4e  u) + ~~,~(~  .2) + o~Ax + Xf~ r.~ )
Using relation (2.8) it can be verified t h a t (2.7) is satisfied. Multiplying (2.8) by Am and (2.7) by w T and comparing the resulting relations, we then obtain (2.6). This completes the proof of sufficiency.  REMARK. A s s u m p t i o n s of the Theorem 2.1 are not very restrictive and seem to be satisfied for most T S R K m e t h o d s of practical interest. For example, t h e y are satisfied for the T S R K m e t h o d s constructed in Section 4.
88
Z, BARTOSZEWSKI AND Z. JACKIEWICZ
3. P  S T A B I L I T Y
PROPERTIES
OF TSRK
METHODS
I n t r o d u c i n g the n o t a t i o n by z = / ~ / A m, relation (2.9) can be rewritten in the form .~2 _ )~Rl(o~, z)  R2(c~, z) = 0. Now define the set P~ = {z : one of the roots of A2  ARI(OZ, z)  R2(ct, z) = 0 is on t h e unit circle and the other is inside or on the unit circle} and the q u a n t i t y era = min Izl.
z6F,~
W e have the following theorem. THEOREM
roots of
3. i. A s s u m e that z 0 is not a pole o/Rl(a, z) and R2(a, z). A s s u m e also that both
,~2 _ )~RI(OL,0) __ i~2(OL, 0 ) ? 0 (3.1)
are inside of the unit circle and that IZl < a s . Then ai1 roots of
are inside o f the unit circle for 311 integers m >_ 1. PROOF. Since the functions R,(c~,0) and R 2 ( a , 0 ) are well defined and the roots of inside of the unit circle, it follows t h a t the roots of
/~2 __/~R1 (o~, z)  Rl(O!, z) : 0
(3.1)
are
(3.3)
are also inside of the unit circle for all [z I < a a (see Figure 1 for a geometrical explanation). Im(z)
\F~
~,\
\
~
~
Re(z)
Figure 1. Geometrical interpretation of ac,,
Delay Differential Equations
89
A s s u m e now to the c o n t r a r y t h a t equation (3.2) has a root ,~ such t h a t IA[ >_ 1. T h e n
<
for all integers m >_ 1. This means t h a t for z = ~/Am we have Izl < ~ and one root of equation (3.3) has modulus greater t h a n or equal to one. This contradiction concludes the proof of t h e theorem.  T h e next t h e o r e m gives a characterization of ~r~ for some T S R K methods. THEOREM 3.2. A s s u m e that ~ = u, P = A, and A = B. Then a s = dist (c~, OSA) , for every c~ E SA, where SA is the stability region o f the T S R K m e t h o d for ODEs. PROOF. We have
Rl(Ct, z) = 1  ' / + v T z + w T ( I  (ct + z ) B ) 1 x ((ct + z ) ( e  u) + (ct + z ) A x ) = Rl(o~ + z, 0),
and
R2(c~, z) = ~ + (c~ + z ) w r ( I  (c~ + z ) B )  l u By the definition of SA, it follows t h a t one of the roots of
= R2(c~ + z, 0).
A2  ARI(OZ, Z)  R 2 ( ~ , z ) = )~2 _ ARI(O~ + z,0)  R2(o~ + z,0) = 0 has modulus equal to one and the second has modulus less t h a n or equal to one if and only if c~ + z E OSA. It follows from the definition of F~ t h a t z E F a if and only if c~ + z E OSA (see F i g u r e 2). Hence, o',~ = zinrf Izl = inf I(ct + z )  o ~ I =dist(c~,OSA),
,, c~+zEOSA
which completes the proof of the theorem. We are now ready to formulate and prove the main result of this paper. I m ( z ) ]
I
(
\
Re(z)
c~+z
//
/
/
Figure 2. Region S A and OSA.
90
Z. BARTOSZEWSKI AND Z..JACKIEWICZ
THEOREM 3.3. A s s u m e that the T S R K method for ODEs such that "~ = u, P = A, and A = B is Astable. Then the corresponding T S R K method (1.2) for DDEs is Pstable. PROOF. We have
Sp = {(c~,/3) : all roots of A2  AR1 (ct, ~~)  JR'2 (c~, k~7) = 0
are inside of the unit circle for m > 1 } It follows from T h e o r e m 3.1 t h a t S ) given by S~ = {(c~,/3) : b o t h roots of A2 /~/E~l(OZ, 0)
 /~2(O', 0) =
l
0
are inside of the unit circle and 1/31 < ~ } satisfies
S'p C Sp.
We have to show t h a t { ( a , / 3 ) : R e ( a ) < 0 and 191 <  R e ( s ) } c Sp
(3.4)
(this is tile definition of P  s t a b i l i t y ) . Take (c~,/3) such t h a t R e ( a ) < 0 and 1/31 <  Re(o,). Since the T S R K m e t h o d for O D E s is Astable we have
{a : Re(a) < 0} c Sa.
This means t h a t both roots of equation (3.1) are inside of the unit circle. It follows fl'om Theorem 3.2 and A  s t a b i l i t y t h a t ~ = dist (a, OSA) >  Re(or). Hence,
191 < Re(c~) < a s .
This means t h a t (c~,/3) ~ S},, and as a consequence of completes the proof.
(3.4)
it follows t h a t (ct,/3) C Sp.
This 
4. E X A M P L E S
OF
PSTABLE
TSRK
METHODS
FOR
DDES
In this section, we will illustrate by two examples how to construct T S R K methods for DDEs which are P  s t a b l e . We s t a r t with the T S R K method for ODEs given by 0.164905 0.210337 0.128015 0.198522 1.07121 0.284316 0.75 2.70983 1.12692 0 0.75 0.0293846
7]
vT
wX
=
0 0
which, as d e m o n s t r a t e d in [17], is Astable and has order p = 4 and stage order q = 4. We c o m p u t e next the continuous weights
such t h a t vi(1) = vi and F = A, and the continuous weights
~{(() : ( (~{,1 + w~,~(+ ~,{,~ ~ ) ,
i = 1,2,
Delay Differential Equations
91
such t h a t wi(1) = wi and A = B, where F and A are the matrices defined in Section 1. This leads to the linear systems of equations
vdcj) = aj~,
v i ( 1 ) = vi,
and
w i ( c j ) = bji,
w~(1) = wi,
i , j = 1, 2, for the coefficients vkt and wkl, k = 1, 2, l = 1, 2, 3. The solutions to the above systems correspond to the continuous weights vi(~) and wi(~) given by
Vl(~) = ~ (0.57142  0.559464 4 + 0.116058~2), v2(~) = ~ (  0 . 3 3 2 7 7 . + 0 . 1 5 1 5 2 5 4  0.10307~2), and
Wl(sc) = ~ (0.755371 + 0.49649s c  0.124943~2), 'w2(~) = ~ (0.00598078  0 . 0 8 8 5 5 t 3 ? + 0.111955,12).
It can be verified using T h e o r e m 3 in [25] t h a t the resulting T S R K m e t h o d for DDEs is convergent with uniform order p = 4; i.e., there is no superconvergence at the gridpoints. We t h i n k this is a very desirable p r o p e r t y of the m e t h o d since we can generate dense o u t p u t without any ~dditionat cost. T h e m e t h o d constructed above is also P  s t a b l e as can be easily verified using T h e o r e m 3.3. T h e o r e m 3.3 can also be used to construct T S R K methods of the t y p e considered in [8,9,14], i.e., with A  0. For example, s t a r t i n g with the Astable m e t h o d of order p = 4 for O D E s ul A B 0
= 0
/~/
vT
wT
0.462626
0 0 0.592719
0 0 0.457494
0.527766 0.0679367 0.0203561
1.06598 0.47028 0.392057
constructed in [14] and proceeding similarly as in the previous example, we o b t a i n the P  s t a b l e m e t h o d of uniform order p = 4 for DDEs. The coefficients rj(~), v~(4), and w~({) of this m e t h o d are given by r/(~ = ~ (0.835974 + 2.602294  1.3036942), c) vl(~) = ~ (1.07105 + 3.33407sc  1.6703~2) , v2(s = 4 (  0 . 8 2 6 7 + 2.57343 4  1.28923~2), c) and w l ( s = ~ (0.226373 + 0.073107~ + 0.173622sc2), c) w2(s = 4 (2.28815  3.378314 + 1.482214~). c) It can be verified by direct computations t h a t both methods constructed in this section satisfy assumptions of T h e o r e m 2.1. For example, the m a t r i x I  (c~ + z ) B corresponding to the l a t t e r m e t h o d is singular at z = 1.55644c~:t:0.834543 i which are also the poles of the function R1 (c~, z). A s y s t e m a t i c approach to the construction and implementation of highly stable T S R K m e t h o d s for DDEs will be the subject of future work.
92
Z. BARTOSZEWSKI AND Z. JACKIEWICZ
REFERENCES
i. J. Hale, Theory of F~mctional Differential Equations, SpringerVerlag, New York, (1977). 2. V.B. Kohnanovskii and V.R. Nosov, Stability of Functional DTffcrential Equations, Academic Press, London, (1986). 3. Y. K u a n g , Delay Differential Equations with Applications in Population Biology, Academic Press. Boston,
(1993).
4. Z. K a m o n t and M. Kwapisz, On ttre Cauchy problem for differentialdelay equations in a Banach space, Math. Nachr. 74, 173190, (1976). 5. Z. Jackiewicz and M. Kwapisz. T h e numerical solution of functional differential equations, a survey, ]flat. Stos. 33, 57 78, (1991). 6. M. Zennaro, Delay differential equations: Theory and numerics, In Theory and Numerics of Ordinary and Partial Differential Equations, Advances ill Numerical Analysis, Volume 4, (Edited by M. Ainsworth et al.), pp. 291 333, Clarendon Press, Oxford, (1995). 7. J.C. Butcher, The Numerical Analysis of Ordinary Differential Equations, John Wiley & Sons. Chichester, New York, (1987). 8. 1. Lie and S.P. Norsett, Superconvergence of multistep collocation, Math. Comput. 52. 65 79, (1989). 9. I. Lie, Local error estimation for multistep collocation methods, B I T 30, 126144. (1990). 10. G.D. Byrne and R.J. Lambert, P s e u d o  R u n g e  K u t t a m e t h o d s involving two points, J. Assoc. Comput. Mach. 13, 114 123, (1966). 11. R.A. Renaut\~Zilliamson, Numerical solution of hyperbolic partial differential equations, Ph.D. Thesis, Cambridge University, Cambridge, U.K., (1985). 12. R. Renaut, Twostep R u n g e  K u t t a m e t h o d s and hyperbolic partial differential equations, Math. Comput. 90, 563579, (1990). 13. R. Caira, C. Costabile and F. Costabile, A class of pseudo R u n g e  K u t t a methods, B I T 30. 642 649, (1990). 14. Z. Jackiewicz, R. R e n a u t and A. Feldstein, Twostep H u n g e  K u t t a methods, S I A M J. Numcr. Anal. 28, 11651182, (1991). 15. Z. Jackiewiez and M. Zemmro, Variablestepsize explicit twostep R u n g e  K u t t a methods, l~.[uth. Comput. 59, 421 438, (1992). 16. Z. Jackiewicz, R. R e n a n t and M. Zennaro, Explicit twostep R u n g e  K u t t a methods, Appl. Math. 4 0 , 4 3 3  4 5 6 ,
(1995).
17. Z. Jackiewicz and S. Tracogna, A general class of twostep R u n g e  K u t t a m e t h o d s for ordinary differential equations, S I A M J. Numer. Anal. 32, 1390 1427, (1995). 18. Z. Jackiewiez and S. Tracogna, A representation formula for twostep R u n g e  K u t t a methods, In Proc. of the Second Hellenic European Conference on Mathematics and Informatics, (Edited by E.A. Lipitakis). pp. 111 120, Hellenic Mathematical Society, Athens, (1994). 19. J.C. Butcher and S. Tracogna, Order conditions for twostep R u u g e  K u t t a methods. Appl. Numer. Math. 24, 351 364, (1997). 20. S. rlS'acogna, Implementation of twostep R u n g e  K u t t a m e t h o d s for ordinary differential equations, & Comp. Appl. Math. 76, 113 136, (1997). 21. S. ~lh'acogna and B. V~lelfert, Twostep R n n g e  K u t t a : Theory and practice, B I T 40. 7"75 799, (2000). 22. Z. Jaekiewicz and R. Vermiglio, General linear m e t h o d s with external stages of different orders, B I T 36, 688712, (1996). 23. E. Hairer and G. Wanner, Order conditions for general twostep R m l g e  K u t t a methods, S I A M J. Numer. Anal. 34, 20872089, (1997). 24. Z. Bartoszewski and Z. Jaekiewicz, Construction of twostep R u n g e  K u t t a m e t h o d s of high order for ordinary differential equations, Numerical Algorithms 18, 5170, (1998). 25. Z. Jackiewicz and S. Tracogna, Variable stepsize continuous twostep R u n g e  K u t t a m e t h o d s for ordinary differential equations, Numerical Algorithms 12, 347 368, (1996). 26. T. Koto, A stability property of Astable natural R u n g e  K u t t a m e t h o d s for s y s t e m s of delaydifferential equations, B I T 34, 262267, (1994). 27. T. Koto, A criterion for Pstability properties of R u n g e  K u t t a methods, B I T 38, 737 750, (1998). 28. M. Zennaro, Pstability properties of R u n g e  K u t t a m e t h o d s for delay differential equations, Nu.n*er. Math. 49, 305318, (1986). 29. K.J. in 't Hour and M.N. Spijker, Stability analysis of numerical m e t h o d s for delay differential equations. Numer. Math. 59, 807814, (1991). 30. K.J. in 't Itout, T h e stability of a clauss of R u n g e  K u t t a m e t h o d s for delay differential equations, Appl. Numer. Math. 9, 347355, (1992). 31. N. Guglielmi. On the asymptotic stability properties of R u n g e  K u t t a m e t h o d s for delay differential equations. Numer. Math. 77, 467485, (1997). 32. N. Guglielmi and E. tIairer, Order stars and stability for delay differential equations. Numer. Math. 83. 371 383, (1999). 33. C.W. Cryer, Itighly stable nmltistep m e t h o d s for retarded differential equations. S I A M .I. Numer. Anal. 11, 788797, (1974). 34. T.A. Bickart, P  s t a b l e and P[c~,13]stable integration/interpolation m e t h o d s in the sohltion of retarded differentialdifference equations, B I T 22, 464476, (1982).
Delay Differential Equations
93
35. T.A. Bickart, Fstable and F[c~,/3]stable integration/interpolation methods in the solution of retarded differentialdifferenceequations, J. Comp. Appl. Math. 8, 279 284, (1982). 36. L.F. Wiederholt, Stability of multistep formulas for delay differential equations, Math. Comp. 30, 283290, (1976). 37. D.S. Watanabe and M.C. Roth, The stability of difference formulas for delay differential equations, S I A M J. Numer. Anal. 22, 132145, (1985). 38. V.K. Barwell, On the asymptotic behavior of the solution of a differential difference equation, Utilitas Math. 6, 189194, (1974).