当前位置:首页 >> >>

Explicit Runge-Kutta Methods for the Numerical Solution of Singular Delay Differential Equa


UNIVERSITY OF MANCHESTER

Explicit Runge-Kutta Methods for the Numerical Solution of Singular Delay Di erential Equations C.A.H. Paul & C.T.H. Baker Numerical Analysis Report No. 212 (revised) April 1992

University of Manchester/UMIST Manchester Centre for Computational Mathematics Numerical Analysis Reports

DEPARTMENT OF MATHEMATICS
Reports available from: Department of Mathematics University of Manchester Manchester M13 9PL England And by anonymous ftp from:
vtx.ma.man.ac.uk (130.88.16.2)

in pub/nareps

Explicit Runge-Kutta Methods for the Numerical Solution of Singular Delay Di erential Equations
C.A.H. Paul and C.T.H. Baker April 27, 1992

In this paper we are concerned with the development of an explicit Runge-Kutta scheme for the numerical solution of delay di erential equations (DDEs) where one or more delay lies in the current Runge-Kutta interval. The scheme presented is also applicable to the numerical solution of Volterra functional equations (VFEs), although the theory is not covered in this paper. We also derive the stability equations of the scheme for the ODE
y0 (t) = y(t);

Abstract

and the DDE

); where the delay and the Runge-Kutta stepsize H are both constant. In the case of the DDE, we consider the two distinct cases: (i) H , and (ii) < H . We evaluate the performance of the scheme by solving several types of singular DDE and a VFE.

y0 (t) = y(t) + y(t ?

Key words. Explicit Runge-Kutta, singular delay di erential equations AMS subject classi cations. primary 65Q05

1 Introduction
We are concerned with extending explicit Runge-Kutta ordinary di erential equation (ODE) solvers to solve delay di erential equations (DDEs). Any Runge-Kutta method may be expressed in the form of a Butcher tableau. Thus, a -stage Runge-Kutta method to solve an ODE from an initial point t0 with initial value y (t0) and derivative function f (t; y (t)) may be associated with the Butcher tableau

c A b

c1 a11
. . . . . .

a1
. . .

c a1 b1
1

a . b

(1)

Mathematics Department, University of Manchester

e The Runge-Kutta method represented by (1) yields an approximation y (t) to the true solution y(t) of the di erential equation, which is given by e e y(tn + H ) = y(tn ) + H
where

X
i=1

bi fni ;

e y(t0 ) tni fni e y(tni)

= = = =

y(t0 ); tn + ci H; e f (tni ; y(tni )); e y(tn ) + H P aij fnj :
j =1

(2)

For coincident fci g, these formulae contain an ambiguity which is resolved in an obvious way. Clearly the Runge-Kutta formulae (2) may only be solved directly if aij = 0 for all j i. In this case the Runge-Kutta method is called explicit , and the corresponding matrix A is strictly lower triangular. The term \order of a Runge-Kutta method" refers to the global order (Defn. 1.3). e The order of a Runge-Kutta method governs how close the approximate solution y (t) is to the true solution y (t) for small stepsizes H , given smoothness. The order equations for a Runge-Kutta method upto order three are:

Pb = 1 i i Pb c = 1 i i 2 Pib c2 = 1 9 > i i 3 = P bi a c = 1 > i ij j 6 ; i;j

order 1; order 2; order 3:

The attainable order of a method is determined by the order equations, of which there are two sub-types that are of interest to us here.

De nition 1.1 The p-th order quadrature equation is
X
i

1 bicp?1 = p : i

(3)

The quadrature order equations1 are associated with the integration of pure quadrature terms { i elementary di erentials of the form @ fi . If the derivative function also depends on y (t), the non@t quadrature order equations must be considered. In 2] the totality of order conditions is presented in terms of simplifying conditions.

De nition 1.2 A Runge-Kutta method is local order p+1 for su ciently smooth derivative function
f (t; y(t)) and solution y(t), if all the order equations upto and including the p-th order equations
are satis ed.
1

In the shorthand notation of Butcher 2], the order equations (3) upto order m may be written as B (m).

2

if for xed T > t0

De nition 1.3 A Runge-Kutta method is of global order p (or a p-th order Runge-Kutta method)
e jy(t) ? y(t)j = O(H p)

for all t 2 t0; T ], as H ! 0.

Theorem 1.4 A Runge-Kutta method is global order p if it is local order p + 1. e De nition 1.5 A continuous extension y(tn + H ) of a -stage Runge-Kutta method is an approximant to y (t) over the interval tn ; tn + H ], which may be written as

e e y(tn + H ) = y(tn ) + H

X
i=1

bi ( )fni for 0

1:

(4)

There exist other continuous extensions in the literature, such as highly continuous extensions 4] and multistage Runge-Kutta continuous extensions 5].

e De nition 1.6 A continuous extension y(tn + H ) has continuous quadrature order q if X
i

De nition 1.7 A continuous extension (tn + H ) of a su ciently smooth function (t) has local
order p if
0

bi( )cs = s + 1 for s = 0; . . . ; q ? 1: i

s+1

(5)

max1 (tn + H ) ?

(tn ) + H

X
i

bi( ) 0 (tn + ci H )

!

= O(H p):

e The order of a continuous Runge-Kutta method (that de nes an extension y (tn + H )) is detere e mined by replacing H by H , fbig by fbi ( )g and y (tn + H ) by y (tn + H ) in the order equations 2]. The continuous order equations upto order three are thus Pb ( ) = i i P b ( )c = i i i P b ( )c2 =
i i
2 3 6
2

3

P bi ( )a c = i ij j
i;j

3

9 > = > order 3: ;

order 1; order 2;

(6)

A Runge-Kutta method can be combined with a continuous extension to produce a continuous Runge-Kutta method, which may be associated with the continuous Butcher tableau

c A b( )

c1 a11 c
. . .

Generally the values fbi( )j =1 g coincide with the values fbig.

a1 b1( )

. . .

a1 a . b()
. . . (7)

De nition 1.8 The full-step order of a continuous Runge-Kutta method is the order of the method
at = 1.

3

2 Classes of DDE
De nition 2.1 A delay di erential equation (DDE) is a di erential equation in which the derivative function depends on a nite number of previous values of the solution y (t).

The simplest DDE is a pure delay di erential equation with constant delay (t) = , for example

y0 (t) = y(t ? ) (t t0 ); (8) (t) = f (t); where the initial function (t) is de ned on the interval t0 ? ; t0]. Over the interval t0 ; t0 + ],
(8) is a pure integration problem. (t; y (t)).

De nition 2.2 A DDE is state-dependent if the delay (t) is dependent on the solution y(t), thus
such point t

De nition 2.3 A DDE is singular at the point t if the delay satis es (t ) = 0. If there is no
t0 , then the DDE is non-singular.
There exist several classes of singular DDE; the simplest being initial value delay di erential equations (IVDDEs). An IVDDE has an initial delay (t0 ) = 0, for example

y0(t) = y(t2 ) (0 t 1); y(t0 ) = (0):

(9)

In this case (t) = t ? t2 ; so this is also a varying-delay DDE. The singular points of (9) are the real roots of t ? t2 , thus ft g = f0; 1g. Another example of a singular varying-delay DDE is

y 0(t) = y(t ? jt ? 1j) (t 0):
Clearly as t ! 1 the delay jt ? 1j ! 0, so that ft g = f1g. Finally there are state-dependent DDEs, such as

(10)

y 0 (t) = y(t ? 1(t; y(t ? 2 (t)))):

(11)

This state-dependent DDE is singular if either the primary delay 1 (t; y (t ? 2(t))) = 0 or the secondary delay 2 (t) = 0 for some t t0 . However since 1 is dependent on y (t), whether (11) is singular is also determined in part by the initial conditions. Another class of singular equations is the class of Volterra functional equations (VFEs), such as

y 0(t) =

t? (t) Z t? (t)

K (s; t; y(s)) ds (t t0 );

(12)

where (t) (t) and (t) 0 for all t t0 . If (t) < H for some t t0 , the integral term requires solution values from the current interval. If (t ) = 0, the VFE is singular at t . Although VFEs are not DDEs, the scheme outlined in this paper for solving singular DDEs is also suitable for solving (singular) VFEs. 4

3 Adapted Runge-Kutta ODE Methods for DDEs
A general DDE with one delay may be written as

y 0(t) = f (t; y(t); y(t ? (t)));
in which case the formula for fni in (2) becomes

e e fni = f (tni; y(tni ); y(tni ? (tni ))):
Thus in order to extend an ODE solver to solve DDEs, it is necessary to include an approximant to e y(t), so that y(t) may be evaluated at non-meshpoints tni ? (tni). The simplest DDE solver is a continuous Runge-Kutta method (7), in which the continuous extension is used to evaluate previous solution values. This method takes no account of the derivative discontinuities that can arise from the initial point t0 or from the initial function (t). E cient DDE methods track the position and order of these discontinuities and include them as meshpoints 10, 11]. This is re ected by a result in 1]:
discontinuities upto order r only occur at meshpoints, the global order of convergence of the resulting DDE method is min(p; q; r).

Theorem 3.1 Given a Runge-Kutta method of order p with an approximation scheme of order q, if

An explicit Runge-Kutta DDE method may be associated with a continuous Butcher tableau (7), where aij = 0 for all j i. A natural formula for the continuous extension is

e e y(tn + H ) = y(tn ) + H
where

X
i

bi( )fni for 0

1;

(13)

e e fni = f (tni ; y(tni ); y(tni ? (tni))); iP ?1 e e y(tni ) = y(tn ) + H aij fnj :
j =1

(14)

Problems can arise when using explicit Runge-Kutta methods to solve DDEs if the lag term e tni ? (tni ) > tn for some tni , since no approximant for y(tni ? (tni )) is available. Thus the explicit Runge-Kutta formulae (14) become implicit for certain DDEs. The requirement of maintaining an explicit method imposes a restriction on the stepsize, which is satis ed when

t ? (t) tn for all t 2 ftni g:

(15)

For some DDEs, a strategy to avoid the Runge-Kutta formulae (14) becoming implicit is to reduce the stepsize H until the stepsize restriction (15) is satis ed. However this strategy (a) reduces the e ciency of the DDE solver, and (b) fails in the case of singular DDEs. 5

We consider two possible methods for solving singular DDEs { both of them iterative. One strategy is to solve the resulting implicit Runge-Kutta formulae using a Newton iteration. However this approach requires us to nd the partial derivatives of the derivative function, which in the case of DDEs (especially state-dependent DDEs) is non-trivial. Although this problem may be overcome by using numerical approximations to the partial derivatives, this is a computationally expensive approach. The second strategy is to use a natural iteration in which the approximant for the current interval is improved after each iteration. We show that, using this approach, the m-th iterate of a -stage Runge-Kutta method is at most a m -stage Runge-Kutta method. Thus it is possible to perform a rigorous analysis of the resulting method.

e De nition 3.2 An m-th iterated continuous extension (ICE) ym(tn + H ) of a continuous Rungee e Kutta DDE method (13), with starting iterate y0 (tn + H ) := y (tn ), is de ned as e e ym(tn + H ) = y(tn) + H e e y0(tni ) = y(tn) + H
r znk

X

i X

m e bi ( )f (tni; ym?1 (tni ); zni);
0 e aij f (tnj ; y0 (tnj ); znj );

(16)

j 8 > y(tnk ? (tnk )) if tnk ? (tnk ) tn ; >e < = > yr?1 (tnk ? (tnk )) if tnk ? (tnk ) > tn and r > 1; e >e : y(tn ) if tnk ? (tnk ) > tn and r 1;

where 0

1 and m 1.

e e De nition 3.3 The starting iterate y0(tn + H ) for an ICE ym(tn + H ) is an initial approximant
to y (t) over the interval tn ; tn + H ].

Two examples of starting iterates are

e y(tn), and e the approximant y (tn?1 + H ) where > 1.

Theorem 3.4 The m-th ICE of a -stage Runge-Kutta method is the rst ICE of a m -stage
Runge-Kutta method with continuous Butcher tableau

c cm A Bm 0m b( )

0 0 B 0 , 0 0 b( ) where B = bj (ci)], Bm = diag (B; .{z. ; B), cm = (c; .{z. ; c), 0m = (0; .{z. ; 0). | . } | . } | . } c
m?1 m m?1

. . . . . . . . .

A B 0
. . .

... ...

0 0

0
... ... ... . . . . . . . . .

(17)

6

Theorem 3.5 The continuous order equations (6) upto order p of A and c B are equivalent for m p ? 1: (18) cm Bm 0m b( ) b( ) Example 3.6 Consider a typical fourth-order continuous order equation for a -stage continuous
Runge-Kutta method. The second ICE has the Butcher tableau

b c
so that

b A b b( )
2 Xb

c c
X
j;k=1

A 0 B 0 , 0 b( )

i;j;k=1

2 Xb

bi ( )bij ajk ck = a b b
=

i= +1 2 X i= +1

bi( )

b a b b aij bjk ck since bi( ) = 0 for i = 1; . . . ; ;
bj (ci? )ajk ck :

bi? ( )

X

Thus if m < p ? 1, the continuous order equations still depend on the original Runge-Kutta matrix A. The third ICE has the Butcher tableau

j;k=1

b c

b A b b( )
2 X

c A 0 0 c B 0 0 c 0 B 0 0 0 b( )

so that (for the same continuous order equation),
i;j;k=1
3 Xb

bi ( )bij bjk ck = a a b
= =

i=2 +1 3 X

j = +1 k=1 2 X bi?2 ( ) bj? i=2 +1 j = +1

3 Xb

bi ( )

bij a

X

b b b ajk ck since bi ( ) = 0 for i = 1; . . . ; 2 ;
(ci?2 )

X
k=1

bk (cj? )ck ;

X

i;j;k=1

bi( )bj (ci)bk(cj )ck by a shift of subscripts.

This nal summation corresponds to the same fourth-order continuous order equation, but for a -stage Runge-Kutta method with Butcher tableau

c

Thus if m p ? 1, the continuous order equations depend only on the matrix B.

B . b( )

Theorem 3.7 Given a Runge-Kutta method of arbitrary order with a q-th order continuous quadrature extension, the m-th ICE is order min(m + 1; q ) if t ? (t) tn for all t 2 ftni g.

7

order continuous quadrature extension, so that

Example 3.8 Consider a typical eighth-order continuous order equation, combined with an eighthX
i

bi( )cs = s + 1 for s = 0; . . . ; 7: i

s+1

(19)

For the m-th ICE, where m 7,
i;j;k;l m Xb

b a bj b a b bi ( )cibij c2ajk bkl cl =
= = = =

X
i;j;k;l

bi( )ci bj (ci)c2bk (cj )bl(ck )cl by (18); j bi( )ci bj (ci)c2bk (cj ) j

X

X
l

i;j;k

bl(ck )cl ;

!

X

i;j;k

1 X bi ( 2 i;j

2 bi ( )cibj (ci)c2 bk (cj ) c2k by (19); j

)cibj (ci )c2 j
6 )ci c6i ;

c3 j

= 288 : In regards to the general order equation for an elementary weight (35): (0) = 0 ; (0) = (1) = i ; (1) = (2) = j ; (2) = (3) = k ; (3) = l; s0( (0)) = 1 ; s0 ( (1)) = 2 ; s0 ( (2)) = 0 ; s0 ( (3)) = 1:

1 X bi ( 6 i 8

3;

p

Theorem 3.9 Given a Runge-Kutta method with a p-th order continuous extension, which is also
a q -th order continuous quadrature extension, the order of the m-th ICE is min(m + p ? 1; q ) if t ? (t) tn for all t 2 ftni g.

For the second ICE, the corresponding sixth-order continuous order equation is
h;i;j;k h;i;j;k

is also a seventh-order continuous quadrature extension. A typical fth-order continuous order equation is 5 X (20) bi ( )aij aik cj ck = 20 : i;j;k
2 Xb X bh ( )bi (ch)aij aik cj ck see (34); bh ( )bhi aij bik cj ck = a b a bb

Example 3.10 Consider the second and third ICEs of a fth-order continuous extension, which

c5 h bh ( ) 20 by (20); h 6 p by (5): = 120
= 8

X

For the third ICE, the corresponding seventh-order continuous order equation is
g;h;i;j;k
3 X b bg ( )bgh bhiaij aik cj bk = a a b b bc

X
g;h;i;j;k

bg ( )bh (cg )bi(ch)aij aik cj ck see (34);

=

c6 g bg ( ) 120 by (5) and (20); g 7 p = 840 by (5):
2 Xb bh ( )bhi ahj bik ck cj by assumption; a b a bb

X

Now consider a continuous order equation for the second ICE where the elementary weight is not a direct extension of an elementary weight of the original Runge-Kutta method.
h;i;j;k;l;m
2 X b bh ( )bhiahj bik akm bjl a b a b a

= = = =

h;i;j;k

X

h;i;j;k

bh ( ) bi(ch)aik ck bj (ch )cj see (34); c2 by (20); h 2

X

ture extension of a Runge-Kutta method if

Theorem 3.11 Given distinct abscissae fc1; . . . ; c g, there exists a q-th order continuous quadraq
, and

= 72 by (5):

h 1 X bh ( 12 h 6

3 bh ( ) c6h

)c5 ; h

p

the values fbi( )j =1 g are used to replace the values fbi g.

Proof Given distinct abscissae fc1; . . . ; c g, there exist polynomials fbi( )g such that
X
i=1

bi ( ) (ci) =

Z
0

(s) ds

for all polynomials of degree q < . The polynomials fbi ( )g are constructed from the cardinal interpolation polynomials associated with any q of the abscissae,

li (s) =
The polynomials fbi ( )g are then de ned as

q Y (s ? cj ) : j (ci ? cj )
j6=i
=1

8R > li(s) ds for each of the q corresponding ci, < bi ( ) = > 0 :0 otherwise:
9

continuous extension, if t ? (t) > tn for some t 2 ftni g, the order of the DDE method is p after p iterations if the current interval contains no jump discontinuities of order r p,

Theorem 3.12 Given an iterated Runge-Kutta scheme (21) for a DDE based on a q-th order

e e y0(tn + H ) := y(tn) for 0
p q.

1,

order of any DDE method cannot exceed r ? 1. Thus for a DDE method of order p based on an iterated Runge-Kutta scheme (Defn. 3.2), jump discontinuities of order r p must only occur at meshpoints. Also by Theorem 3.1, the order of the DDE method cannot exceed the order of the approximation scheme used; so the order of the DDE method cannot exceed q . A q -th order continuous extension is at least a q -th order continuous quadrature extension. Therefore, by Theorem 3.7, the order of the p-th ICE is at least min(p + 1; q ) if t ? (t) tn for all t 2 ftni g. Now consider the iterated Runge-Kutta scheme when t ? (t) > tn for some t 2 ftni g:
m e e ym (tn + H ) = y(tn ) + H P bi( )fni ; i 0 e e y0 (tni ) = y(tn ) + H P aij fnj ; j r fnk

Proof By Theorem 3.1, if there exists a discontinuity of order r in the current interval, then the

where

8 < f (t ; y (t ); zr ) e = : nk r?1 nk 0 nk e f (tnk ; y0 (tnk ); znk )

if r 1; if r = 0:

(21)

e e For the rst ICE, for each tni ? (tni ) > tn , the initial estimate for y (tni ? (tni )) is y (tn ) by 1 is at least order zero, so the approximant assumption. Therefore the derivative function estimate fni after one iteration is O(H ) accurate. After each iteration the accuracy of the approximant increases by order H , as shown by
m e e ym (tn + H ) = y(tn) + H P bi ( )fni ; i P e e e = y (tn ) + H bi ( )f (tni ; y(tni ) + O(H m ); y(tni ? (tni )) + O(H m?1 )) by (21); i P e = y (tn ) + H bi ( )fni + O(H m); i e = y (tn + H ) + O(H m ):

Thus after p iterations the order of the approximant is O(H p), if p q and no discontinuities upto order p exist in the current interval.

Theorem 3.13 Given a p-th order Runge-Kutta DDE method with a continuous extension of order

q, and an iterated Runge-Kutta scheme (21) for solving the DDE when t ? (t) > tn for some t 2 ftni g, the order of the resulting Runge-Kutta DDE method is min(p; q; r) if jump discontinuities upto order r are included in the meshpoints.
10

by Theorem 3.1. If t ? (t) > tn for some t 2 ftni g, the order of the Runge-Kutta DDE method is min(p; q; r) by Theorem 3.12, if the Runge-Kutta method is iterated at least min(p; q; r) times.

Proof If t ? (t) tn for all t 2 ftnig, the order of the Runge-Kutta DDE method is min(p; q; r)

3.1 Parallel Runge-Kutta methods
There is a connection between the iterated Runge-Kutta scheme represented by the continuous Butcher tableau (17) and parallel Runge-Kutta methods. The implication for solving the RungeKutta equations on a parallel computer gives rise to our terminology.

De nition 3.14 An (m + 1)-stage explicit parallelisable Runge-Kutta method is a Runge-Kutta

method in which the values ffni g may be grouped into an ordered sequence of m + 1 batches (one of which contains more than one element), in which the members of each batch depend only on members of preceding batches. continuous Runge-Kutta ODE method of order p based on these abscissae. If c1 = 0 and cp = 1, the rst-same-as-last strategy may be used (which sets fn;1 := fn?1;p ) and the method may be reduced to a (p ? 1)-stage parallelisable explicit continuous Runge-Kutta ODE method.

Theorem 3.15 Given p distinct abscissae fc1; . . . ; cpg, there exists a p-stage parallelisable explicit

Proof A proof is based on constructing a suitable tableau (17) (where m = p ? 1). Given p distinct

abscissae fc1 ; . . . ; cpg, there exists a p-th order continuous quadrature extension by Theorem 3.11. Therefore the (p2 ? p)-stage Runge-Kutta method (17) has at least a p-th order continuous extension by Theorem 3.7. However, since the rst abscissae c1 = 0 (for all explicit Runge-Kutta methods), the rst stage need only be evaluated once. Also, from the de nition of the p-th ICE (16) for p 2, it is clear that the p stages of the p-th ICE are independent of each other. Thus they can all be calculated in one parallel stage. For the rst p stages of the Runge-Kutta method, if ai1 = ci for all i and aij = 0 otherwise, the last p ? 1 stages only depend on the rst stage. Therefore they can all be calculated in one parallel stage. So after the rst iteration, the Runge-Kutta method is order two and is equivalent to two parallel stages. Each further iteration increases the order of the method by one (upto a maximum order of min(m + 1; p)) and is equivalent to one parallel stage. Thus the resulting p-th order parallel Runge-Kutta method has only p parallel stages. If the rst-same-as-last strategy is used, the resulting parallel Runge-Kutta method is clearly only p ? 1 parallel stages.
explicit continuous Runge-Kutta singular DDE method of order p based on these abscissae. If the rst-same-as-last strategy is used, this method may be reduced to a p-stage parallelisable explicit continuous Runge-Kutta singular DDE method.

Theorem 3.16 Given p distinct abscissae fc1; . . . ; cpg, there exists a (p + 1)-stage parallelisable

Proof The proof is the same as for Theorem 3.15, except that the Runge-Kutta method is iterated
meshpoints, the expected order of the iterated Runge-Kutta scheme is maintained. 11

p times (as required by Theorem 3.12). If derivative discontinuities upto order p are included in the

4 Stability
In this section we examine the stability of the iterated Runge-Kutta scheme (16). There are three possible types of stability analysis, the type being determined by the nature of the stability equation. We consider an example stability equation of each type: (i) y 0(t) = y (t) for an ODE stability analysis; (ii) y 0(t) = y (t) + y (t ? ) where

H , for analysis of a non-singular DDE;

(iii) y 0(t) = y (t) + y (t ? ) where < H , for analysis of a `singular' DDE.

4.1 ODE stability
Using the notation of (34), the stability polynomial for the m-th iterate is

b b e R( H ) = 1 + H bT (Im ? H A)?1 b;
i=0 where e and b are vectors 1; 1; . . . ; 1]T of the appropriate e lower triangular, (I ? H A)?1 may be expanded as

= 1 + H bT

m?2 X

( H )iBi e + ( H )mbT Bm?1 (I ? H A)?1e; length. Using the fact that A is strictly

(I ? H A)?1 = Thus

?1 X
i=0

( H )iAi :
?1 X
i=1

R( H ) = 1 + bT

m?1 X

If ( H B) < 1, since (I ?

i=0 ?1 exists, H A)

( H )i+1Bi e + bT
1 X

( H )m+i Bm?1 Ai e

R( H ) ! 1 + bT

! 1+

i=0 H bT (I

( H )i+1Bi e as m ! 1;

? H B)?1e:

Now suppose that D denotes the disk in the complex H -plane such that ( H B) < 1, then the intersection of D with the stability region of the iterated Runge-Kutta scheme tends to the intersection of D with the stability region of the implicit Runge-Kutta method with tableau (cf. Theorem 3.5) with = 1 c B . (22) b( ) How close the stability regions of the iterated methods can become to the stability region of (22) is determined by the restriction that ( H B) < 1. For the fourth-order continuous embedded and the fth-order hermite approximants (after 250 iterations), the stability regions are 12

Iterated CEF stability region : y'(t) = y(t)
6 6

Iterated hermite stability region : y'(t) = y(t)

4

4

2

2

( H) space

0

Stability region

( H) space

0

Stability region

-2

-2

-4

-4

-6 -6

-4

-2

0

2

4

6

-6 -6

-4

-2

0

2

4

6

( H) space

( H) space

The stability regions of (22) for the fourth-order continuous embedded and fth-order hermite approximants are
B matrix stability region : y'(t) = y(t)
20 15 10 5 0 -5 -10 -15 -20 -4

( H) space

B matrix stability region

-3

-2

-1

0

1

2

( H) space

For the continuous extensions considered in Section 5, the B matrices have the following spectral radius (B): Low order hermite 0:2887; Continuous embedded formulae 0:2308; High order hermite 0:1950; High order quadrature formulae 0:1751: Note that D is a circular disk and in the gures the stability regions are semicircles of radius 1= (B), the radius of D. We have found, in the cases investigated, that the stability region of the iterated Runge-Kutta scheme (16) is well approximated by the intersection of the stability region of (22) and the semicircle of radius 1= (B). (This also happens to be the stability region of the truncated stability polynomial for (22).) 13

4.2 Non-singular DDE stability
In order to analyse the stability of

y 0(t) = y(t) + y(t ? );

(23)

we impose restrictions on the Runge-Kutta stepsize H and the delay . We require both H and to be constant, and in addition that H . In the case that = NH , it is possible to solve (23) without an approximant. However in general, and for the case that = (N + #)H where 0 # < 1, e the stability also depends on the continuous extension y (t). The most general form of (23) allows the parameters and to be complex; however we restrict ourselves to considering and real. Other authors have considered (23) with = const and complex; and in some cases with = 0. For a rigorous analysis of (23) with constant delay and stepsize, we refer the reader to Paul and Baker 6]. In order to apply the stability analysis of 6], we use Theorem 3.4 and express the m-th iterate of the iterated Runge-Kutta scheme as a m -stage Runge-Kutta method. We now use two of the main results from 6], applying them to the m -stage Butcher tableau: For the linear stability equation (23), if = NH the iterated Runge-Kutta scheme is stable if all the zeros i of SH ( ; ; ) := det N +1 I ? N W ? X ? Z]; satisfy j i j 1 for all i and if j i j = 1, then i is semi-simple. This condition on the roots i is b known as the root condition. Using the notation of (34), if we write S S( H ) = (Im ? H A)?1 b e and r( H ) = 1 + H bT Sb, then

0 Br ( H ) B B W=B B H Sb B e @

0 C C

0 B0 B C; X = B B C B0 C B 0 C @ A 1

b b H bT SB C C

0 T b e B H b Sb B C; Z = B B C B H Sb C b B H SB C e @ A 1

0 C C

1

C: C C 0 C A

b We now address the case = (N + #)H for 0 # < 1. If # ci we de ne i = 1, i = 0 and i = b b b ci ? #. Alternatively, if # > ci we de ne i = 0, i = 1 and i = ci ? # +1. Writing = 1 ; . . . ; m ]T , b = 1; . . . ; m ]T , = diag( 1; . . . ; m ), ? = diag( 1; . . . ; m ) and B(#) = bj ( i )], the iterated b Runge-Kutta scheme is stable if all the zeros i of
SH ( ; ; ) := det
satisfy the root condition, where
N +2 I ? N +1 U ? 2 V ?

W ? X ? Z] 0 C C
1

0 Br ( H ) B B U=B B H Sb B e @

0 C C

0 bT b B0 H b S B(#) B C; V = B B C B0 H S B(#) C b B 0 C @ A 1
14

0 T 1 b B Hb S C C B C; W = B B C B HS C B C @ A

C; C C 0 C A

0 B0 B B X=B B0 B @

b b H bT S?B(#) C C b H S?B(#)

1

0 T b B Hb S B C C; Z = B B C B HS C B A @

0 0

1 C C C C: C C A

The following gures represent the stability regions for the original Dormand and Prince RungeKutta method and the iterated Runge-Kutta method applied to (23). To aid comparison, they are all plotted on the same axis.

15

Dormand & Prince CEF method : y'(t) = y(t) + y(t- )
4 3 2 1

=H = 3H = 5H

H space

0 -1 -2 -3 -4 -4

Stability region

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

H space Iterated Dormand & Prince CEF method : y'(t) = y(t) + y(t- )
4 3 2 1

=H = 3H = 5H

H space

0 -1 -2 -3 -4 -4

Stability region

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

H space

16

Dormand & Prince CEF method : y'(t) = y(t) + y(t- )
4 3 2 1

= 1.25H = 3.25H = 5.25H

H space

0 -1 -2 -3 -4 -4

Stability region

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

H space Iterated Dormand & Prince CEF method : y'(t) = y(t) + y(t- )
4 3 2 1

= 1.25H = 3.25H = 5.25H

H space

0 -1 -2 -3 -4 -4

Stability region

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

H space

17

Dormand & Prince CEF method : y'(t) = y(t) + y(t- )
4 3 2 1

= 1.5H = 3.5H = 5.5H

H space

0 -1 -2 -3 -4 -4

Stability region

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

H space Iterated Dormand & Prince CEF method : y'(t) = y(t) + y(t- )
4 3 2 1

= 1.5H = 3.5H = 5.5H

H space

0 -1 -2 -3 -4 -4

Stability region

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

H space

18

Dormand & Prince CEF method : y'(t) = y(t) + y(t- )
4 3 2 1

= 1.75H = 3.75H = 5.75H

H space

0 -1 -2 -3 -4 -4

Stability region

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

H space Iterated Dormand & Prince CEF method : y'(t) = y(t) + y(t- )
4 3 2 1

= 1.75H = 3.75H = 5.75H

H space

0 -1 -2 -3 -4 -4

Stability region

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

H space

19

Dormand & Prince hermite method : y'(t) = y(t) + y(t- )
4 3 2 1

=H = 3H = 5H

H space

0 -1 -2 -3 -4 -4

Stability region

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

H space Iterated Dormand & Prince hermite method : y'(t) = y(t) + y(t- )
4 3 2 1

=H = 3H = 5H

H space

0 -1 -2 -3 -4 -4

Stability region

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

H space

20

Dormand & Prince hermite method : y'(t) = y(t) + y(t- )
4 3 2 1

= 1.25H = 3.25H = 5.25H

H space

0 -1 -2 -3 -4 -4

Stability region

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

H space Iterated Dormand & Prince hermite method : y'(t) = y(t) + y(t- )
4 3 2 1

= 1.25H = 3.25H = 5.25H

H space

0 -1 -2 -3 -4 -4

Stability region

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

H space

21

Dormand & Prince hermite method : y'(t) = y(t) + y(t- )
4 3 2 1

= 1.5H = 3.5H = 5.5H

H space

0 -1 -2 -3 -4 -4

Stability region

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

H space Iterated Dormand & Prince hermite method : y'(t) = y(t) + y(t- )
4 3 2 1

= 1.5H = 3.5H = 5.5H

H space

0 -1 -2 -3 -4 -4

Stability region

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

H space

22

Dormand & Prince hermite method : y'(t) = y(t) + y(t- )
4 3 2 1

= 1.75H = 3.75H = 5.75H

H space

0 -1 -2 -3 -4 -4

Stability region

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

H space Iterated Dormand & Prince hermite method : y'(t) = y(t) + y(t- )
4 3 2 1

= 1.75H = 3.75H = 5.75H

H space

0 -1 -2 -3 -4 -4

Stability region

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

H space

23

4.3 Singular DDE stability
Finally we analyze the iterated Runge-Kutta scheme when applied to a DDE which cannot be solved by the straightforward application of an explicit Runge-Kutta DDE method. We avoid the obvious choice of an initial value DDE for the stability equation, because this would require a varying delay. Instead we remain with the linear stability equation (23), but we consider a delay < H . Using the notation of section 4.2 for f i g, , , , ? and B(#), after the rst iteration for the n-th interval

fn1] = e yn + H Af n1] + yn?1 + H ?B(#)fn?1 + yn ; = (I ? H A)?1 e yn + yn?1 + H ?B(#)fn?1 + yn ] :
For the i-th iterate, where i > 1,

may be expressed as a vector recurrence

fni] = e yn + H Bf ni?1] + yn?1 + H ?B(#)fn?1 + yn + H B(#)fni?1]: ] Given that yn+1 = yn + H bT fns] and fn?1 fns?1 , where s is the last iteration, the iterated scheme W n=X
n?1 + Y n?2 ;

1] s where m = ym+1 ; H fm ; . . . ; H fm]]T . This vector recurrence is stable if all the roots i of

SH ( ; ; ) := det 2W ? X ? Y]
satisfy the root condition, where

0 B B B B B W=B B B B B @

0

1 0 . . . . . .

? H B ? H B(#)

0 I 0

...

0 I 0

0
... ...

0 1 B B H Se + H S B B B X = B He + H B B . B . B . @ 0 B 0 B HS B B B Y = B ... B B . B .. B @
H He + H

? H B ? H B(#)

1 ?bT C C 0 C . C; . C . C C C C 0 C A 1 C C A

I

0
. . . . . . . . .

0
. . . . . . . . .

0 C C H S?B(#) C C C H ?B(#) C ; C C
H ?B(#)
. . .

0 0
. . . . . . . . .

0

0C
. . . . . . . . .

0
24

0

1 C C C C: C C C C C A

5 Numerical Methods
For a practical application of the iterated Runge-Kutta scheme (21), we consider the fth-order Dormand and Prince explicit Runge-Kutta ODE method 3]. The following continuous extensions are investigated: low order hermite, continuous embedded formula, high order quadrature formula, high order hermite. The four approximation schemes have the following continuous extensions and B matrices: Third degree Hermite formula

0 0 52 0 1113 36 0 371 448 0 1113 0 352000 811377 500 0 1113 35 500 384 0 1113 Fourth degree continuous embedded formula

0 0 B 6599 B 48000 B B 2667 B B 16000 B B 341 B = B 3000 B B 433 B B 4374 B B 35 B 384 @
b1( b2( b3( b4( b5( b6( b7(
) ) ) ) ) ) )

b1( b2( b3( b4( b5( b6( b7(

) ) ) ) ) ) )

= = = = = = =

500 1000 ? 1113 3 + 371 2 ; ? 125 3 + 125 2; 96 64 2187 3 ? 6561 2 ; 3392 6784 11 ? 11 3 + 28 2 ; 42 3 ? 2: 13 192 9 64 7 12 1375 2187 125 192 125 192

0;

157 3 ? 221 2 + 192 128

;

0

848000 ?59049 848000 ?15309 53000 ?33 106 ?2187 6784 ?2187 6784

?28431

0

143 10500 99 3500 44 375 1936 15309 11 84 11 84

0

?4 C C 125 C C ?63 C 1000 C C ?16 C C 125 C C ?64 C 729 C C

0 1

C 0 C A 0

= = = = = = =

1163 ? 1152 4 + 1039 3 ? 1337 2 + ; 360 480 0; 7580 4 18728 3 4216 2 3339 ? 3339 + 1113 ; 415 9 ? 192 4 + 2 3 ? 27 2; 16 8991 4 + 2673 3 ? 2187 2 ; ? 6784 2120 8480 187 4 ? 319 3 + 33 2 ; 84 105 35 0:

(24)

25

0 0 B 79241 B 720000 B B 152397 B B 1280000 B B 917 B = B 11250 B B 25042 B B 295245 B B 35 B 384 @
35 384

Sixth degree quadrature formula The polynomials fbi( )g satisfy the following condition

0 0 0 0 0 0 0

46028 417375 38559 185500 201344 417375 10316288 21907179 500 1113 500 1113

0

?839 24000 ?6129 128000 127 375 9404 19683 125 192 125 192

0

?9963 4240000 6561 67840000 ?4131 66250 ?1042 7155 ?2187 6784 ?2187 6784

0

297 17500 5841 280000 ?176 4375 704 688905 11 84 11 84

0

01 C 0C C C C 0C C C 0C C C 0C C C C 0C A 0

X
i=1

bi( )cs = i

s+1

s+1

for s = 0; . . . ; 5:

However the polynomials fbi( )g do not satisfy fbi( )j =1 g = fbig for the original Runge-Kutta method. This can be seen from table (26) in the `order with respect to B' column.

b1( b2( b3( b4( b5( b6( b7(

) ) ) ) ) ) )

= = = = = = =

? 125 6 + 1435 5 ? 535 4 + 4663 3 ? 281 2 + ; 32 96 24 288 48 3125 6 ? 33625 5 + 60125 4 ? 20125 3 + 500 2 ; 124 372 496 279 31 10000 6 + 104000 5 ? 132200 4 + 214400 3 ? 12800 2 ; ? 371 1113 1113 3339 1113 625 6 ? 5375 5 + 1775 4 ? 6275 3 + 25 2 ; 16 48 16 144 4 2460375 6 + 6790635 5 ? 3247695 4 + 2499741 3 ? 177147 2 ; ? 52576 52576 26288 52576 26288 375 6 ? 985 5 + 3635 4 ? 257 3 + 12 2 ; 28 28 112 21 7 0:
26197 111600 7059 24800 208 6975 512000 16474671 125 4464 125 4464

0 0 B 11561 B 180000 B B 80889 B B 1280000 B B 487 B = B 5625 B B 45932 B B 531441 B B 25 B 288 @
25 288

0

Fifth degree Hermite formula In order to obtain a fth-order hermite approximant it is necessary to have three fth-order solution points { the natural choice of support points is t = 0; 1 ; 1. In order to obtain a fth-order 2 1 solution at the point t = 2 , it is necessary to add two extra stages to the Runge-Kutta method (in our case due to Shampine 8]).
1 2 1 2

83475 ?297 5300 34816 83475 81920000 197164611 200 477 200 477

?9032

0

41 900 1047 25600 152 225 382000 531441 25 36 25 36

0

32860000 ?91939293 2103040000 ?531441 1026875 ?20788 44361 ?19683 52576 ?19683 52576

?1594323

0

2561 210000 441 40000 1424 13125 131072 1240029 7 48 7 48

0

01 C 0C C C C 0C C C 0C C C 0C C C C 0C A 0

?33728713 2 ?30167461 7739027 ?19162737 104693760 21674880 17448960 123305984 7157 10825 ?220887 0 70925
75776 164724 113664 4016128

80069 3530688

0

?26949 363520 ?107
5254

?5
74

26

The resulting approximant may be expressed as a continuous extension with polynomials fbi( )g:

0 0 0 0 0 0 0 0 1 C 3317 13947637 220196 0 1029525 35520 ?38790819 689587500 ?163052 ?128 ?384 C 784400000 8209375 4625 3125 C C ?790016571 94914171 92643 0 274540 109827 12550400000 3677800000 ?200718 ?441 ?441 C 947200 8209375 9250 3125 C C 16640 1024 ?1139427 13044592 ?446768 ?128 384 C C 0 41181 2775 6128125 172396875 8209375 4625 3125 C C 344550400 ?4684 117119552 ?8731328 ?2560 3584 C 0 810565623 1165600 2184813 17649 1085852061 155121723 242757 59049 C C C 125 ?2187 11 500 0 1113 0 0 0 C C 192 6784 84 C 125 ?2187 11 500 0 1113 0 0 0 C C 192 6784 84 C 10825 ?220887 80069 ?107 ?5 70925 C 0 164724 113664 0 C 4016128 3530688 5254 74 A 7157 70925 10825 ?220887 80069 ?107 ?5 0 164724 113664 0 75776 4016128 3530688 5254 74 The orders of the continuous extensions di er according to which matrix A or B they are associated with. In table (26) the orders are given for

0 0 B 1367209 B 14800000 B B 22966881 B B 236800000 B B 11147 B 115625 B B 205348 B = B 2184813 B B B 35 B 384 B B 35 B 384 B B B 7157 B 75776 @

b1( b2( b3( b4( b5( b6( b7( b8( b9(

) ) ) ) ) ) ) ) )

= = = = = = = = =

? 4000 5 + 1245700 4 ? 398800 3 + 413200 2; 371 41181 13727 41181 ? 125 5 + 83775 4 ? 44725 3 + 225 2 ; 8 2368 1776 37 6561 5 ? 4428675 4 + 798255 3 ? 98415 2 ; 848 251008 62752 31376 23529 ? 22 5 + 527571 4 ? 285659 3 + 18389 2 ; 7 73556 55167 5 ? 21872 4 + 14847 3 ? 3483 2 ; 4 2627 2627 2627 ? 40 4 + 80 3 ? 40 2; 37 37 37 5 ? 40 4 + 32 3 ? 8 2 : 16

0;

29 5 81685 4 24433 3 6839 2 16 ? 14208 + 3552 ? 1776 +

;

(25)

c A and c B . b( ) b( )
Continuous Extension Low order Hermite Embedded formula Quadrature formula High order Hermite Full-step Order Continuous Order wrt A wrt B wrt A wrt B 5 4 3 3 5 5 4 4 2 6 2 6 5 5 5 5

(26)

27

To compare the two full-step fth-order iterated Runge-Kutta methods with the original Dormand and Prince method, we consider the sixth-order defects:
Elementary Weight Original Embedded Hermite Elementary Weight Original Embedded Hermite ?1 ?1 ?1 bi aij c2 ajkck 1 1 ?1 bic5 i j 5400 5400 5400 10800 10800 54000 ?1 ?1 ?1 bi aij ajk ck ajm cm 1 1 ?1 bic3 aij cj i 10800 10800 10800 21600 21600 108000 ?1 ?1 ?1 bi ci aij ajk c2 ?1 ?1 bici aij cj aik ck 0 k 21600 21600 21600 64800 64800 ?1 ?1 ?1 bi aij cj ajk c2 1 1 ?1 bic2 aij c2 i j k 16200 16200 16200 16200 16200 81000 ?1 ?1 ?1 bi aij ajk c3 1 ?1 biaij c2 aik ck 0 j k 32400 32400 32400 21600 108000 ?1 ?1 bi ci aij ajk akmcm ?1 ?1 ?1 bici aij c3 0 j 21600 21600 3600 129600 129600 1 1 ?1 bi aij cj ajk akm cm 1 ?1 biaij c4 0 j 5400 5400 27000 32400 162000 ?1 ?1 bi aij ajk ck akm cm 1 ?1 bic2 aij ajk ck 0 0 i 32400 32400 43200 216000 ?1 ?1 bi aij ajk akm c2 1 ?1 biaij aik ck ajm cm 0 0 m 64800 64800 64800 324000 ?1 ?1 bi aij ajk akm amn cn 1 1 ?1 bici aij cj ajkck 0 43200 43200 3600 129600 648000

Each defect shown is the amount by which the quantity stated algebraically di ers from the value necessary to satisfy the sixth-order condition. The defects of the iterated schemes correspond to the sixth-order defects of the Runge-Kutta method represented by the tableau (17) with m = 5 (as required by Theorem 3.12). The 2-norms of the defects are 5:03 10?4, 3:27 10?4 and 2:37 10?4 respectively.

6 Numerical Results
In this section we show that the iterated Runge-Kutta scheme (21) developed in Section 3 can indeed solve the types of DDE suggested in Section 2. The two continuous extensions used are the fourth-order continuous embedded formula (24) and the fth-order hermite formula (25). If a constant stepsize Runge-Kutta method is employed, it should be possible to demonstrate the order of convergence of the numerical methods2 .

De nition 6.1 The scaled relative error is de ned as
where p is the (expected) order of the Runge-Kutta DDE method.

e t) ( E (t) := y(H p?(y)t) yt

(y (t) 6= 0);

Example 6.2 For the rst example consider (10)
y 0(t) = y(t ? jt ? 1j) (t 0); (27) (t) = 1 (t 0); which has the solution (28) with a discontinuity in the (n + 1)-st derivative at t = 2n ?1 for n 0. 2n
2

The graphs presented di er from those that appear in earlier versions of this report, in that the numerical testing was carried out using a later version of the DDE code. However, the qualitative and quantitative results are unchanged.

28

y(t) =
where
n
0

n= ln (1?2t) +1 ? ?ln X ?(i+1)i 2i t ? 2i + 1 i i=0 n X i=1

2

2

i!

n?i ;

(28)
2

=

i

i i?1 ? 1 ? 2i?n 1 ? 2i?n 2

2

?(i?1)i

= 1; z] = integer part of z: Table (29) gives the relative and scaled relative error (where p = 5) in the numerical solution. Embedded Formula Hermite Formula Stepsize Relative Scaled Relative Scaled Error Error Error Error 1 t= 2 0:05 1:48E? 16 4:74E ? 10 2:96E? 16 9:47E? 10 (29) 0:025 1:33E? 15 1:36E? 7 1:33E? 15 1:36E? 7 0:0125 1:48E? 15 4:85E? 6 1:48E? 15 4:85E? 6 t=1 0:05 1:85E? 8 5:92E? 2 1:87E? 8 5:98E? 2 0:025 3:25E? 11 3:33E? 3 3:21E? 11 3:29E? 3 0:0125 1:10E? 13 3:59E? 4 8:43E? 14 2:76E? 4
10 -1 10 -2 10 -3
Scaled absolute relative error

i!

n?i ;

10 -1

CEF 0.05 CEF 0.025
Scaled absolute relative error

10 -2 10 -3 10 -4 10 -5 10 -6 10 -7 10 -8 10 -9 10 -10

HER 0.05 HER 0.025 HER 0.0125

10 -4 10 -5 10 -6 10 -7 10 -8 10 -9 10 -10

CEF 0.0125

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

t

t

Example 6.3 The next example is an IVDDE (9), which has a smooth solution (31).
y 0(t) = y(t2) (0 t 1); (0) = 1:
29 (30)

y(t) =

1 X
n=0

m=1

t(2n ?1) n Q (2m ? 1) (0 t 1):
Hermite Formula Relative Scaled Error Error 3:07E ? 5 3:46E ? 5 0:00 1:98E ? 3 2:70E ? 3 2:86E ? 3

(31)

Stepsize 0:02 0:01 0:005 0:02 0:01 0:005
10 -2

Embedded Formula Relative Scaled Error Error 7:50E? 14 3:17E? 15 0:00 8:96E? 12 3:57E? 13 1:25E? 14

2:34E? 5 3:17E? 5 0:00

9:82E? 14 3:46E? 15 0:00 t=1 2:80E? 3 6:34E? 12 3:57E? 3 2:70E? 13 3:99E? 3 8:94E? 15
10 -2

1 t= 2

10 -3
Scaled absolute relative error Scaled absolute relative error

10 -3

10 -4

10 -4

10 -5

10 -5

CEF 0.02
10 -6

HER 0.02
10 -6

CEF 0.01 CEF 0.005

HER 0.01 HER 0.005

10 -7

0

0.2

0.4

0.6

0.8

1

10 -7

0

0.2

0.4

0.6

0.8

1

t

t

and secondary delay 2 = 0 for t = 0. is the root of 1 ? t + 2t2 ? t3 which lies in the interval 1; 2].

Example 6.4 This example is a state-dependent DDE (11) with primary delay 1 = 0 for t 2 0; ]
y0 (t) = y(t ? y(t ? t2 )) (t 0); (t) = t2 (t 0):
(32)

8 < y(t) = : 0 f (t) ? f ( )
30

0 t

; t < 1;

6 1 where f (x) = 1 x9 ? 1 x8 + 7 x7 ? x6 + x5 ? 1 x4 + 3 x3 . 9 2 2

Stepsize 0:02 0:01 0:005 0:02 0:01 0:005
10 6 10 5 10 4
Scaled absolute relative error

Embedded Formula Relative Scaled Error Error 1:97E? 8 1:63E? 9 1:32E? 11 5:82E? 12 4:99E? 13 2:45E? 14 6:17 16:3 4:22

1:97E? 8 1:63E? 9 1:32E? 11 t=5 1:82E? 3 5:82E? 12 4:99E? 3 4:99E? 13 7:85E? 3 2:45E? 14
10 6 10 5 10 4
Scaled absolute relative error

t = 21 2

Hermite Formula Relative Scaled Error Error 6:17 16:3 4:22 1:82E ? 3 4:99E ? 3 7:85E ? 3

10 3 10 2 10 1 10 0 10 -1

10 3 10 2 10 1 10 0 10 -1

CEF 0.02
10 -2 10 -3

HER 0.02
10 -2 10 -3

CEF 0.01 CEF 0.005
0
1 2 3 4

HER 0.01 HER 0.005
0
1 2 3 4

5

5

t

t

y(0) = 1 and y(t) = et for t 0. Thus

Example 6.5 The nal example is a VFE (12) with (t) = 1, (t) = 0 and initial conditions
Zt 0 (t) = y (t ? 1) + y (s) ds (t y t?1
has the analytical solution 0); (33)

y(t) = et for t 0:

The integral term was evaluated using an adaptive quadrature scheme 7] applied to the continuous extension, and the following errors resulted. 31

Stepsize 0:05 0:025 0:0125 0:05 0:025 0:0125
10 -4

Embedded Formula Relative Scaled Error Error 6:14E? 12 1:79E? 13 1:88E? 14 1:31E? 11 3:43E? 13 3:41E? 14

Hermite Formula Relative Scaled Error Error t=5 1:97E? 5 1:67E ? 12 5:35E ? 6 1:84E? 5 5:46E ? 14 5:59E ? 6 6:16E? 5 2:57E ? 14 8:42E ? 5 t = 10 4:19E? 5 3:34E ? 12 1:07E ? 5 3:51E? 5 1:47E ? 13 1:50E ? 5 1:12E? 4 2:43E ? 14 7:98E ? 5
10 -4

Scaled absolute relative error

10 -6

Scaled absolute relative error

10 -5

10 -5

10 -6

CEF 0.05
10 -7

HER 0.05
10 -7

CEF 0.025 CEF 0.0125
10 -8 0 5 t 10 10 -8 0

HER 0.025 HER 0.0125
5 t 10

7 Comments
It must be noted that order results are asymptotic results. Thus if a constant stepsize is too large, the order of the method is not apparent in the error in the numerical solution. However if the stepsize is too small, the error in the numerical solution is dominated by round-o error. First, comparing the results for the continuous embedded and hermite approximants, it is clear that the higher order continuous hermite approximant (cf. (26)) performs marginally better. The graphs for (27) re ect that initially the error in the numerical solution is due to round-o error. This can be veri ed by comparing the numerical results against the unit round-o for double precision FORTRAN which is 10?15 . The sharp jump in the error is due to a derivative discontinuity occurring within a Runge-Kutta interval, as opposed to occurring at a meshpoint. Due to the di ering 32

stepsizes, the rst derivative discontinuity within a Runge-Kutta interval occurs at t = 0:875, 0:9375 and 0:96875 for H = 0:05, 0:025 and 0:0125 respectively. The graphs for (30), (32) and (33) all indicate that the order of the Runge-Kutta DDE method is ve, by the fact that the lines virtually coincide. However for (30), the error in the numerical solution corresponding to the smallest stepsize is dominated by round-o error. For (33), again for the smallest stepsize, the dips in the graph correspond to the sign of the relative error changing. Also the error in the numerical solution is again dominated by round-o error. However, considering that the di erential equation contains a Volterra term, the agreement in the numerical solutions is surprising.

8 Conclusion
It is clear just from the small number of illustrative examples that the iterated Runge-Kutta scheme can successfully solve singular DDEs explicitly. However the p2 sequential Runge-Kutta stages may be prohibitive, especially when a more straightforward method is available { such as extrapolation. But extrapolation cannot be used for IVDDEs or the state-dependent DDE (32) included in our examples. Also, we consider the persistent extrapolation that would be required to solve the VFE (33) to be unreliable. Further, extrapolation is only reasonable if no low order derivative discontinuities are crossed and, by their nature, a derivative discontinuity always occurs at the singular points ft g of a DDE. There is also a similar algorithm for solving singular DDEs due to Tavernini 9], but it is harder to implement. e The iterated Runge-Kutta scheme (16) is de ned so that the current estimate for ym (tni ) is

e e ym (tni ) = y(tn) + Hn

X
j

m e bj (ci)f (tnj ; ym?1 (tnj ); znj ):

e However it is possible to use a current estimate for ym (tni ), so that e e ym (tni ) = y(tn) + Hn X
j m e aij f (tnj ; ym (tnj ); znj ):

(Using bj (ci ) in place of aij simpli es the analysis, but the resulting Runge-Kutta method is implicit.) The disadvantages of this alternative scheme are: In general there is no scope for parallelization. The analysis of this iterated scheme is considerably harder than the one presented for (16). From numerical experimentation, it appears that (16) is more e cient. Finally, the parallel continuous explicit Runge-Kutta ODE and DDE methods suggested in section 3.1 were unexpected. Despite the fact that only the continuous quadrature order equations need to be solved (because the non-quadrature order equations are automatically satis ed), so that constructing high order parallel continuous Runge-Kutta methods is straightforward, further investigation was postponed. 33

Appendix
Theorem 8.1 The m-th ICE of a -stage Runge-Kutta method is the rst ICE of a m -stage
(4) of the -stage Runge-Kutta method. Consider the case m = 2; the i-th projected gradient of the second ICE is

Proof Consider the case m = 1; this is trivial since the rst ICE is simply the continuous extension
1 0 X e e e f (tni ; y1 (tni); z1( ni)) = f @tni ; y(tn) + H bj (ci )f (tnj ; y0 (tnj ); z0( nj )); z1 ( ni)A ; j =1 8 < y(s) if s tn; e zi (s) = : e yi (s) if s > tn ; where ni = tni ? (tni ). Thus each projected gradient of the second ICE may be expressed in terms

Runge-Kutta method.

of an extra stage of the original Runge-Kutta method, where bj (ci) corresponds to an akl . Therefore e y2(tn + H ) is equivalent to the rst ICE of a 2 -stage Runge-Kutta method with Butcher tableau 0 b (c ) 1 c A 0 b (c1) 1 1 . C: . C B c B 0 , where B = B ... . A @ b1(c ) b (c ) 0 b( ) We now assume as an inductive hypothesis that the m-th ICE of the original Runge-Kutta method is equivalent to the rst ICE of an m -stage Runge-Kutta method with Butcher tableau

cm

Considering the (m + 1)-st ICE of the original Runge-Kutta method, from (16) it is clear that the projected gradients of the (m + 1)-st ICE depend only on the projected gradients of the m-th ICE. Also, since the same recursive de nition is used, it is also the same linear combination of projected gradients that form the (m + 1)-st ICE as the second ICE. Thus the (m + 1)-st ICE may be written as the m-th ICE with an extra stages added, which is equivalent to an (m + 1) -stage RungeKutta method. Since the same linear combination of projected gradients is used, the same matrix B applies. Therefore the continuous Butcher tableau of the resulting Runge-Kutta method is

A Bm . 0m b( )

cm+1

Theorem 8.2 The continuous order equations (6) upto order p of A and c B are equivalent for m p ? 1: cm Bm 0m b( ) b( )
34

A Bm+1 . 0m+1 b( )

Proof Let the operator mean \has the same order as" and consider the m-th ICE. The resulting
Runge-Kutta method has the following Butcher tableau by Theorem 3.4.

b c

b A b b( ) 8 > > >c > m > > > > > < > > >c > m > > > > > :

cm

A B 0
. . .

0 A Bm 0m b( ) 0T ?p+3 m

... ... ... ... ...

0 0

0
. . . . . . . . . (34)

0m

0 B

b( )
for p > m + 1:

0

0m?p+2 Bp?1 0 0m b( ) c B for p m + 1 by a shift subscripts. b( )

...

0T ?p+3 m
. . . for p m + 1;

Theorem 8.3 Given a Runge-Kutta method of arbitrary order, with a q-th order continuous quadrature extension, the m-th ICE is order min(m + 1; q ) if t ? (t) tn for all t 2 ftni g.

Outline Proof By Theorem 3.5, the continuous order equations upto order m + 1 of the m-th
ICE are equivalent to the order equations of

c

B . b( )

The continuous quadrature order equations (5) are automatically satis ed because the coe cients fbi( )g and fcig remain unchanged in the iterated Runge-Kutta method (16). Now consider the nonquadrature continuous order equations upto order r = min(m + 1; q ). The general order equation for an elementary weight of order r may be expressed as (cf. Example 3.8):

X
i

bi( )

1? X X r?Y S a | {z } j=0

(j )

r?1?S

s ( (j )) (j )c (j ) =
0

r

for some constant ;

(35)

where the inner summation is over all (j ) for j 1 and

X
j

aij = ci ; S =

X
j

s0( (j )); a0

(j ) := 1 for any j

and a (j ) (j ) = b (j )(c (j )):

(36)

35

Therefore (35) may be reduced to

1 P b ( ) X X r?Q?S b (c )cs ( (j)) by (36) i i | {z } j=0 (j) (j) (j) r?1?S ! 2 P b ( ) X X r?Q?S b (c )cs ( (j)) P b s ( (r?1?S )) = i (c )c i | {z } j=0 (j) (j) (j) (r?1?S)=1 (r?1?S) (r?1?S) (r?1?S) r?2?S ! ? 2 P b ( ) X X r?Q?S b (c )cs ( (j)) cs r? r?S?S if s0 ( (r ? 1 ? S )) < q: = i s ( (r?1?S ))+1 i | {z } j=0 (j) (j) (j)
0 0 0 0 0( ( ( 1 1 ))+1 ) 0

r?2?S

8 < s ( (l)) + sk?1 ( (l)) + 1 if l = r ? k ? S; sk ( (l)) := : k?1 sk?1 ( (l)) otherwise and repeating this process of rearrangement r ? 2 ? S times yields 1? r X X X r?Y S a (j) (j)cs ((j) (j)) = r?Q S ; bi ( ) 1? | {z } j=0 i si ( (r ? 1 ? S ? i)) + 1 r?1?S
0

Now setting up a sequence of the form

i=0

which is the correct continuous order equation. Since this result holds for all continuous order equations upto order min(m + 1; q ), the order of the m-th ICE is min(m + 1; q ).

Theorem 8.4 Given a Runge-Kutta method with a p-th order continuous extension, which is also
a q -th order continuous quadrature extension, the order of the m-th ICE is min(m + p ? 1; q ) if t ? (t) tn for all t 2 ftni g.

Outline Proof A q-th order continuous quadrature extension satis es
X
i

bi( )cs = s + 1 for s = 0; . . . ; q ? 1: i

s+1

(37)

For m = 1; this is trivial since the rst ICE is the original p-th order continuous extension. The s-th order continuous order equations of the original Runge-Kutta method are of the form

X
(1)

b

(1)( )

? X X sY1 a | {z } j=1
s?1

(j ) (j ) =

s;

(38)

for s p and some xed constant . For m 2, consider the continuous order equations for the m-th ICE upto order m + p ? 1. From (35) and (38), the r-th order continuous order equations for the m -stage Runge-Kutta are equivalent to rY ?m Y X X X m?1 b (j)(c (j)) a (k+m?1) (k+m?1) ; b (1)( )
(1)

| {z } j=1
r?1

k=1

by a shift of subscripts.

36

If the elementary weight is an extended elementary weight of the original Runge-Kutta method, we have (for example):

X
(1)

b

?m X X m?2 Y X X X rY a b (j)(c (j)) b (m?1)(c (m?1)) | {z } j=1 | {z } k=1 (k+m?1) (k+m?1) (m?1) r?m m?2 m?2 X X XY m+1 = b (c ) cr?m?1) for r ? m p ? 1 by (38); (39) b (1)( ) ( | {z } j=1 (j) (j) (1) m?2 0 1 X X X @m?3 Y X m+1 = b (1)( ) b (j)(c (j)) b (m?2)(c (m?2))cr?m?1)A ; (40) ( | {z } j=1 (1) (m?2) m?3 m+2 cr?m?2) X X X m?3 Y ( = b (1)( ) if r q by (37): (41) b (c ) | {z } j=1 (j) (j) r ? m + 2 (1)
(1)( )

m?3

Repeating (40) and (41) m ? 2 times in total yields:

X
(1)

b

(1)(

)

rY ?m Y X X m?1 b (j)(c (j)) a | {z } j=1 k=1 r?1

(k+m?1) (k+m?1) = m?2 Q

r

i=1

(r ? m + i + 1)

;

which is correct. For elementary weights which are not direct extensions of elementary weights of the original Runge-Kutta method, the same proof applies. However stage (39) is repeated more than once (cf. Example 3.10). Thus for r min(m + p ? 1; q ), all the continuous order equations are satis ed. So the resulting continuous Runge-Kutta method is order min(m + p ? 1; q ).

References
1] Arndt, H., Van der Houwen, P.J. and Sommeijer, B.P., Numerical Interpolation of Retarded Di erential Equations, Delay Equations: Approximation, Theory and Applications, International Series of Numerical Mathematics 74 (1985), Birkhauser, Boston, pp. 41{51. 2] Butcher, J.C., The Numerical Analysis of Ordinary Di erential Equations, Wiley, Chichester (UK), 1987. 3] Dormand, J.R. and Prince, P.J., A Family of Embedded Runge-Kutta Formulae, J. Comp. Appl. Math. (1980) 6, pp. 19{26. 4] Higham, D.J., Highly Continuous Runge-Kutta Interpolants, Technical Report 220/89, Dept of Computer Science, University of Toronto (1989). 5] in't Hout, K.J., A new interpolation procedure for adapting Runge-Kutta methods to delay di erential equations, IMACS '91 Proceedings (eds. R. Vichnevetsky & J.J.H. Miller), (1991) 1 p. 309 INCA, Dublin; BIT (1992) 32 pp. 634{649. 6] Paul, C.A.H. and Baker, C.T.H., Stability boundaries revisited { Runge-Kutta methods for delay di erential equations, Technical Report No. 205, Univ. of Manchester (1991). 7] Paul, C.A.H., A fast, e cient, very low storage, adaptive quadrature scheme, Technical Report No. 213, Univ. of Manchester (1992). 8] Shampine, L.F., Some practical Runge-Kutta formulas, J. Math. Comp. (1986) 46, pp. 135{150.

37

9] Tavernini, L., One-Step Methods for the Numerical Solution of Volterra Functional Di erential Equations, SIAM J. of Num. Anal. 8 (1971), pp. 786{795. 10] Wille, D.R. and Baker, C.T.H., The Tracking of Derivative Discontinuities in Systems of Delay di erential equations, Technical Report No. 185, Univ. of Manchester (1990). 11] Wille, D.R. and Baker, C.T.H., A short note on the propagation of derivative discontinuities in Volterra-delay integro-di erential equations, Technical Report No. 187, Univ. of Manchester (1990).

38


赞助商链接
相关文章:
更多相关文章: