Polynomial and rational interpolation in the numerical solution of stiff systems

270

POL~q(~qlALANDP,ATIONAL INTEBPOLATION IN THE NUMERICAL SOLb'TION OF STIFF STYSTEMS S.D. C o l q u l t t and J a c k W i l l i a m s Department o f N a t h e m a t i c s , J n i v e r s l t y o f Manchester Manchester, 1413 9PL, E n g l a n d AbStract

We a n a l y s e t h e a c c u r a c y o f t h e p o l y n o m i a l I n t e r p o l a t i o n scheme

used i n standard v a r i a b l e o r d e r , v a r i a b l e s t e p BDF codes f o r s t i f f

systems.

i n p a r t i c u l a r the s t r a t e g y used i n the NAG cooe i s considered and shown t o be r e l i a b l e and robust on a l a r g e number o f t e s t problems.

For the case

o f monotonic s o | u t i o n components, a comparison i s made w i t h a p a r t i c u l a r l y simple i n t e r p o l a t i o n scheme based on a r a t i o n a l q u a d r a t i c . 1.

Introduction

In the numerical s o l u t i o n of i n i t i a l

v a l u e problems

for y

|

- f(x,y),

y ~ I~m,

x ~[a,b],

(1) yCa) " Y 0 , i t i s now con~con p r a c t i c e t o improve t h e e f f i c i e n c y of codes by using interpolants.

This a p p l i e s t o l i n e a r m u l t l s t e p methods f o r s t i f f

and

n o n s t i f f problems and more r e c e n t l y t o Runge K u t t a methods f o r n o n s t i f f problems, Gladwell et a l [ 8 ] , E n r l g h t et a l

[6].

T y p i c a l l y , the main

purpose o f f n t e r p o l a n t s is t o a l l o w a s t e p s f z e s e l e c t i o n which i s u n r e s t r i c t e d by the u s e r ' s p r l n t p o i n t requirements ( p o s s i b l y beyond the end o f the s p e c i f i e d range [ a , b ] ) ,

approximations at the u s e r ' s p r i n t po#n¢~

be i n g computed by a s , ; l t a b l e i n t e r p o l a n t . This is p a r t i c u l a r l y e s s e n t i a l when dense output is required.

virtually

u s e f u l and

Polynomial i n t e r p o l a n t s

are used i n t h i s f a s h i o n i n many implementations o f the backward

differentiation library

formulae (BDF) f o r s t i f f

systems e . g . i}O2 QBF i n the NAG

[ 3 ] , DCO3 i n the Harwell code [4] and i n LSODE [10] and i t s

Variants. Other uses may i n c l u d e the c a l c u l a t i o n o f continuous numertca~ solutions

(~)ElsevierSciencePublishingCo., Inc., 1989 655 AVenueof the Americas,New York, NY 10910

0096-3003/89/$03.50

271

for producing "defect"

visually

Finally,

finding

in

x

pleasing

graphical

tnterpolants

are

y(x)

y'(x)

in general

variable

order

current

accepted

given

locally

on codes for

stiff

by interpolation. systems;

specifically

terms ~hose codes which implement the variable

BDF r a n g i n g step

from orders

from

Xn_ 1

the Nordsieck array

to

one to five. x n - Xn_ 1 + h n ,

is available

at

Typically,

we

step,

for a

with an order

the end of the step,

p

and is

by

(Cn0 , C n l . . . . . where

Yn

and

Cnp ) !

h

(Yn,

Is the approximation t

xn

out root

- O.

a r e now a p p r o x i m a t e d

p a p e r we c o n c e n t r a t e

consider

formula,

in carrying

the

in problems of the form

and

In this

or for monitoring

indispensible

C(x,y(x),y'(x))

where

output

to the true

hPn,(p)) p!~n ,

solution

y(x)

of

(1)

at

H

Yn, Yn . . . . .

y(r)(xn),

' h2 u nYn, ::~-Yn . . . . .

can be interpreted

r - 1,2 .... p.

intepolating

From t h e s e

polynomial

of degree

as approximations

scaled p

deriYative

to the derivatives

approximations

an

can be formed with the Taylor

representation,

and for higher that

Xn_ 1 < x < x n derlvat'ves).

~p,n(X)

interpolates

provides

From t h e u n d e r l y i n g

can also certain

be defined

fixed

values

or variable

discussion.

of

Yn--t

to

derivation

(and possibly

o f t h e BDF, I t

and back information

d e p e n ~ ~n w h e t h e r

i - 1,2 ..... ~he BDF i s

form, see Jackson

c a s e we a l s o

y(x)

as the polynomial

~p,n ( X n - l ) " Y n - l ,

coefficient

In either

uniquely

comvuted current

~ p , n ( X n ) " Yn - Yn, The a c t u a l

approximations

which

according p

to, (3)

implemented

and Sac~s-Davls

ha, e the approximate

I s known

in i~s

{11] f o r

derivative

272

condition I

!

h n ~ p , n ( X n ) - On, I - hnY n ~ h n f ( x n , Y n ) w h e r e now

f(xn,Y n)

denotes

computed approximation Cn, 1

is usually

correc¢or,

to

obtained

Cn, 1

t h e ~xa0¢ v a l u e

Y(Xn).

i - 1,2 .....

or

b)

p > 1

predict

implicit

EPISODE ( B y r n e

form,

hn_ i - hn,

and c o n s t a n t o r d e r t h r o u g h o u t t h e since the associated

a value at

Xn)

Pn(x)

~th

p

for example

p o l y n o m i a l (used t o

i s b~sed on i n t e r p o l a t i n g The a c t u a l

I - 1,2 .....

steps.

order predictor

interpolated

correcter

polynomial

then satisfies t

i

Cn(xn) " Yn,

C n ( x n - i h n ) - PnCx n - ihn~ , Therefore,

at

xn

which are stored however,

differently

it

is th-

scaled

in the Nordsieck

faylor

array

i - 1,2 ..... coefficients

at

Xn,

a n o r d e r c h a n g e h a s b e e n made t h e a n d now t h e s i t u a t i o n

relationship

between the

stepsize/order Finally,

coefficient

subject to constant stepstzes

Cn(xn) = Yn,

If,

In the

is only true when:

p a s t v a l u e s a t e a u a l l v spaced p o i n t s . Cn(x)

Yn

since

1

p-l, This arises

relationship

p.

implemented fixed

DO2QBF i n NAG a n d LSODE, t h i s p -

with

Is the

we a l w a y s h a v e

i n t h e more w i d e l y

a)

that

Yn

f o r m o f t h e BDF, f o r e x a m p l e

Yn-i - Yn-i, However,

and

residual.

coefficient

and Hindmarsh [I])

f(x,.)

is an approximate

by the condition

~,t~lds a zero

In the variable

This

of

in all

an approximation

changes will situa¢tons to

{Yn-i}

be d i s c u s s e d

and

Y(Xn_i),

Pn(X)

polynomial

~g,n(x)=--Cn(x).

may b e o b t a i n e d

dependent.

{Yn-t} subject

further

a n d c o d e s w~ r e g a r d

of this

and hence

is implementation values

p.

The to

in the next section. each

wh!ch is only available

Yn-i

value

as being

a s a ~n~£,;~ed

273

q u a n t i t y when, b ec a u s e o f t h e p a r t i c u l a r

order/stepsize

values,

Yn-i " Y n - i -

~Yn-p

Yn-1

Yn-2

Yn

fp,n(x) Xn-p

Xn-2

Xn-1

Xn

Fi=.,ure :

In v i e w o f t h e i m p o r ta n t r o l e p l a y e d by software,

it

reliability

is of considerable interest o f t h e above i n t e r p o l a t i o n

performance with a p a r t i c u l a r l y rational

~p,n(X)

in p r e s e n t day

t o c o n s i d e r t h e a c c u r a c y and

scheme.

I n a d d i t i o n we compare i t s

simple i n t e r p o l a t i o n

scheme b a s e d on a

q u a d r a t i c which, when a p p l i e d t o monotonic d a t a ,

interpolating

y i e l d s an

f u n c t i o n which i s g l o b a l l y c o n t i n u o u s l y d i f f e r e n t i a b l e

and

monotonic.

The n u m e r i c a l t e s t i n g

has been c a r r i e d out u s i n g t h e c o n s t a n t

coefficient

BBF implemented i n t h e NAG code DO~QBF, and run on t h e CDC7600

(a wo r k i v ~ p t o c : ~ i o v o f 1~ decimal d i g i t s ) .

conclusions will '_ ~ ' ~ l ~ 2.

objective

is to relate

the e r r o r

Yn-i,

in (2),

in the q u a n t i t i e s

Lp, n ( y ( x ) )

d e n o t e t h e Lagrange i n t e r p o l a t i n g

defined ~

Ln n (Y(Xn-r)) " Y(Xn-r), r - 0 ,1 . . . . . p.

~p,n(X) - y(x) - (Xp,n(X) -

values

t - 0 , 1 , 2 . . . . . p,

and pure i n t e r p o l a t i o n

xe(xn_ 1, Xn),

via (3).

error.

Let

polynomial of degree

p

Then - y(x)),

e r r o r s due t o i n e x a c t d a t a Expressing

Lagrange form and assuming a p p r o p r i a t e d i f f e r e n t i a b i l l t y , y i e l d s zo r

~p,n ix) - y(x)

L p , n ( Y ( x ) ) ) + (Lp,n(Y(X))

t h e e r r o r c o n s i s t s o f two p a r t s ;

(Yn-i}

the o v e r a l l

- approximation properties

to the error

showing t h a t

is that

Rp? ,~able to other similiar ~mplementations,

Polynomial i n t e r p o l a t i o n Our f i r s t

Our b e l i e f

~,n(X) the error

i~ i t s formula

274

"

IITp.n(X)

y(X)ll

II i ~-o

(

1

÷ ~ where

y(P+l)([x)

et(x) [~n-t - y(Xn-t)ltt II ~ (X-Xn_ i ) y ( p + l ) ( ~ x ) l l l-O

denotes a vector of

component i s e v a l u a t e d at some p o i n t

(p+l)st

{Xn_i}y. O.

( Xn_lm,~
+ ~ 1

in which the r t h

~r(X) e (Xn_p, Xn); (ill(y)}

the Lagrange b ~ s i s f u n c t i o n s f o r the ~et IITp,n(X) - y ( X ) l l

derivatives

are

Therefore,

O~g~p I I ~ n _ i - Y ( X n _ t ) l l

Xn-m~]~'X
Now l e t

Ap -

x Xn-l~X(Xn

the s e t {Xn_i)i.O;

i ~-0 I f l l ( x ) l '

the m o d i f i e d Lebesgue c o n s t a n t f o r

i t s v a l u e depends c r i t i c a l l y

p o l a t i n g p o i n t s and on

p.

Consfde- now the next term in (4), and l e t

For

p - 1,

x e (Xn_l. ~'n),

P(x) (

2

hn/4.

P(x) - IX - Xnll x - Xn_l! /

( Let

h - O (max i ( p - 1 hn_ 1,

on the p o s i t i o n o f the I n t e r

For

P(x) - I

p > 1,

~ (X-Xn_i)l. i-O

x e(Xn_ 1, Xn),

~ IX - Xn_ll 1-2

2

~-

(hn+hn- 1) (hn+hn_l+hn_2)... (hn+hn_l+..+hn_p+ 1) then for

x e (Xn_ 1, Xn)

P(x) ( ~1

p l h p+I.

Combining these r e s u l t s ~n (4) y i e l d s ,

+

1

hP+1 I ly(p+l.) (~x) I I.

(5)

275

As f a r as t h e second term i s c o n c e r n e d we n o t e t h a t t h i s

i s o f t h e same form

as t h e p r i n c i p a l

~e)

[Xn_l,

step

local t r u n c a t i o n error (or local e r r o r

for the

Xn] , ee " ~ - T1

where t h e d e r i v a t i v e s

hnP+lll y ( p + l ) l l ,

are d e f i n e d at p o i n t s i n

( ~ n - p , Xn)-

The codes

a t t e m p t a t each s t e p t o keep t h i s v a l u e (or some r e l a t e d q u a n t i t y ) than the u s e r ' s control,

tolerance

less

T O L . Hence f o r s i m p l e a b s o l u t e l o c a l e r r o r

(5) t a k e s the form I I ~ p , n ( X ) - y(X) l l

( Ap

max 0(i(p

IlYn_i-Y(Xn_i)li

+ kp TOL, where i f

hP+llly(p+l)(~x)ll

(6)

~ hnP+llly(p+l)ll,

then

kp ~ 1 / 4 .

From t h e

p o i n t o f view o f g a i n i n g i n s i g h t i n t o the i n t e r p o l a t i o n p r o c e s s t h e form o f

(6) seems a p p r o p r i a t e - most codes c o n t r o l i n d i r e c t l y v a l u e o f TOL i n the c o n t r o l o f t h e y may not s a t i s f y Stetter

[13],

ee) the v a l u e s o f

llYn_i-Y(Xn_i)ll

the d e s i r a b l e t o l e r a n c e p r o p o r t i o n a l i t y

although

condition,

Seward [ 1 2 ] .

We now d i s c u s s f u r t h e r the r o l e o f m ~ h a n g e d by mapping

Ap.

The v a l u e o f

[Xn_p, Xn] ------> [ - 1 , 1 ]

x -- ( xn - Xn-D 2

so w i t h

(through the

Ap

is

by t h e t r a n s f o r m a t i o n

) t + Xn-° + xn "2 '

Xn_p+ i ~ t i ,

Ap -

max t p _ 1 ( t (1

i~ 0 l e i ( t ) , , -

We may r e g a r d the v a l u e o f interpolation process.

Ap

fli(t)

J ; ? (t i - t J)

as a measure o f t h e s t a b l i i t y

o f i=h~

I t i n d i c a t e s the amount by which e r r o r s

in the

c u r r e n t and back v a l u e s can be d i r e c t l y p r o p a g a t e d i n t o the i n t e r p o l a t e d value.

In all

implementations

Ap - 1,

for

p - 1.

Otherwise its

actual

276

value will

be code ~ependent,

t o = -1 < t 1 < t 2 < ... relative

p+l

Similarly,

makin 8 one step there

at

the distribution

is restricted

in stepsizes

c h a n g e i s made, a l l

steps.

that


increase/decrease

a stepsize

since

a c c o r d l n 8 t o maximum/minimum

and order

changes.

Typically,

c o d e s k e e p t h e same s t e p s l z e

followin 8 an order

t h e new o r d e r

a r e no r e j e c t e d

of the points

change;

and this

least

t h e NAG c o d e DO2QBF

with the old stepsize.

steps

for at

when

We a s s u m e now

is the situation

we a n a l y s e

below. For

p ~ 5 ( o r 6)

for the various

it

is therefore

stepsize

commonly o c c u r r i n g

possible

distributions

constant

to compute the values

which the codes allow.

order,

constant

3

4

step,

the values

of

Ap

For the are

approximately:

2

P Ao

1.25

Now we c a n t r e a t DO2QBF.

in detail

at

Xn_l,

if

Yn-1

Yn " Yn

as a result

p - 1

of an order

has been obtained and

Yn-i change

and

Yn-I

in step

at

Referring

order

to Figure

1,

by

if

Yn

1)

Yn " Yn- A l t e r n a t i v e l y

of an order

change at

all

Xn_2, t h e n

increases

tti¢ v a l u e s

as the

satisfy.

i - 0,1,2,...,P,

and stepsize

throughout.

change is to change

order

throushout

the

values

is similar

to above,

Xr_l,

i n t h e NAG c o d e

change (increase/decrease

down t h e b l o c k u n t i l

the purpose of an order of fixed

can arise

The n u m b e r o f m a t c h i n g v a l u e s

further

correspondln 8 to equal

4.55

that

t o 5.

as a result

~n-1 - Yn-1.

change occurs

the case

3.11

t h e n we o n l y h a v e t h e m a t c h l n 8

Yn-i - Yn-l,

since

2.21

6

the situations

The c o d e u s e s o r d e r

Is obtained

order

1.63

5

hn ~ hn_l,

p

then

steps, if

This situation

(Increase)

the step.

the situation Yn

Yn = Yn

In

of matching

is obtained an6

~rises

via a

Yn-1 - Yn-1

only.

2"t7

The m a t c h i n g Clearly

sn~s=~es

as the step

in any given block,

there

change occurs exist

further

down t h e b l o c k .

a t m o s t two s t e p s i z e s

a n d two

orders. The s i t u a t i o n all

possible

c a n now b e i l l u s t r a t e d

stepsize

distributions

The c a s e

of equal steps

includes

the case of an order

by

h

schematically

as follows,

can be Indentified

has been considered change at

for a given block,

- w h e r e we n o t e

Xn_ 1 .

in which

The s t e p

that

sizes

this

are denoted

and

h

p-2

h

p-3

p-4

p-5

Values of

g

Ap

h

h

h

h

h

h

~

h

F~

~

h

h

h

h

h

h

h

h

h

h

~

E

h

E

F~

r~

~

can now be calculated for each distribution,

where we

278

consider the cases of the l a r g e s t v a l u e In BO2QBF, e x c l u d i n g s t a r t situation

1 h - ~ h.

Increase up),

~ - 10h

h - 5h,

(this

Is the d e f a u l t

and a t y p i c a l d e c r e a s e

The a p p r o p r i a t e v a l u e s a r e shown i n T a b l e s 1, 2

and 3, and f o r each p

P

correspond to the d i s t r i b u t i o n s

shown above.

Ap, h - 1Oh

2

5.55

3

36.2,

5.09

4

212,

62.8,

5.32

5

1102,

729,

83.8,

6.08

Table 1

P

Ap, ~ - 5h

2

3.08

3

11.0,

3.17

4

37.7,

17.0,

3.64

5

121,

99.8,

22.1,

4.53

Table 2

r

^p, 5 -

1 T+h

2

1.03

3

1.05,

1.34

4

1.09,

1.41,

1.83

5

1.13,

1.47,

1.98,

Table 3

2.58

279

i n o r d e r t o make a s i m p l e p r e l i m i n a r y c o m p a r i s o n b et w een p o l y n o m i a l interpolation

and t h e p r o p o s e d r a t i o n a l

interpolant

of the next section

seems necessary to c o n s i d e r th® asymptotic form o f the e r r o r i n ( 5 ) . this

rigorously

i s not s t r a i g h t f o r w a r d ,

s t e p and o r d e r t o o b t a i n

Yn-i,

Yn-i

is s a t i s f i e d

for all

maximum step s i z e . 0(hP)

" Y(Xn-l)

h ( h O,

for so~

i - 0,1 ....

h o > 0;

P

where a g a in

h

i s the

Hence, f o r m u l a (5) shows t h a t the i n t e r p o l a n t s are i s d o m i n a t e d by t h e e r r o r

in the

values.

The above r e s u l t s interpolation.

c o n f i r m t h e e x p e c t e d form o f t h e e r r o r

In a d d i t i o n

the i n t ~ r p o l a n t s y(x)

However, some i n s i g h t may

and t h a t

+ O(hP),

a c c u r a t e and t h a t t h e e r r o r

{Yn-ll

p

To do

s i n c e t h e c o d e s combine v a r i a b l e

i - O,l,...p.

he g a i n e d by assuming f i x e d o r d e r

i t may be d e s i r a b l e

t o c o n s i d e r t o what e x t e n t

are able to share other properties

as r e p r e s e n t e d by computed i n f o r m a t i o n .

in polynomial

of the t r u e s o l u t i o n

Many s t i f f

systems, have

s o l u t i o n components which a r e w h o l l y monotonic i n c r e a s i n g / d e c r e a s i n g a range [ a , b ] , tnterpolant

or have o n l y a s m a l l n~'daber o f t u r n i n g p o i n t s .

which i s a l s o monotonic would t h e r e f o r e

particularly

if visually

In g e n e r a l , globally

i - 0,1 . . . . p

be v e r y d e s i r a b l e

Yn-p < Yn-p+l < - "

which i s

{Yn-i, i n - i } ,

monotonic i n c r e a s i n g sequence, < Yn

For c o ~ ' c n i e n c e t h e v a l u e s a r e now t r e a t e d

c o r r e s p o n d t o a s o l u t i o n component which i s m o n o t o n i c . strictly

-

[8] f o r n o n s t ~ f ~ p r o b l e m s .

Suppose now t n e computed v a l u e s

corrospond to • strictly

f n - i > 0.

An

t h e above p o l y n o m i a l scheme p r o d u c e s an i n t e r p o l a n t

C0[a,b].

over

a c c e p t a b l e p l o t s a r e r e q u i r e d i n t h e c o n t e x t o f low

t o medium a c c u r a c y ; s e e C l a d w e l l e t a l

with

it

d e c r e a s i n g s e q u e n c e can be c o n s i d e r e d s i m i l a r l y ,

as s c a l a r s

The ca~e o f a if

p - 1,

and

280

then it

is clear that

tp,n(X)

t

satisfies

Tp,n(x)

> 0,

x e

[xn_ X, Xn].

!

if

p - 2,

then since

Wp,n(Xn) " f n

a simple argument us in g the Nean

t

Value theorem shews t h a t

r p , n (x) > 0

x e [Xn_ 1, Xn].

and 2 E~DFformulae y i e l d monotonic i n t e r p o l a n t s . m o n t o n o c i t y cannot he guaranteed, a l t h o u g h i t a s y m p t o t i c a l l y as a later 3.

h ~ 0.

The c a s e s

p > 2

Hence the o r d e r s 1

However, f o r

p > 2

i s t o be expected are i n v e s t i g a t e d n u m e r ' c a l l y in

section.

C1 monotonic r a t i o n a l i n t e r p o l a t i o n In [ 9 ] , G r e g o r y and D e l b o u r g o d e s c r i b e a p a r t i c u l a r l y

interpolation

scheme which i s g l o b a l l y

C1

simple

and monotonic when a p p l i e d t e

mon o t o n i c d a t a . Suppose t h a t t h e code y i e l d s t h e computed s e q u e n c e {Yn,

in}.

n - N, N+I . . . . . N

where

YM < YN+I < - . - < YN; f n > O. The v a l u e s a r e a g a i n t r e a t e d a s s o l u t i o n components, t h a t strictly generally

is scalars;

These v a l u e s w i l l

d e c r e a s i n g s e q u e n c e c a n be t r e a t e d s i m i l a r l y .

i n v o l v e v a r i a b l e stepslze and o r d e r i n t h e i r c a l c u l a t i o n .

a pJecewise r~tional quadratic function which i s m o n o t o n i c i n c r e a s i n g on R(xn) " Yn, The i n t e r p o l a n t

R'(Xn) " i n .

0 - (x - X n _ l ) / h n ,

e e [0,I],

~.D[~_1o2 " Yn-Z + ~ n - I

~ (in +

Yn " Y(Xn),

fn-i

J~+2. . . . . N,

In [9]

is constructed and s a t i s f i e s

N, H+l . . . . . N. as f o l l o w s .

then for

Let

x e[Xn_ 1, Xn],

with

R(x) - Yn-1 + P n - 1 ( e ) / Q n - l ( 0 )

+ fn_zo(z ~ze)] - 2~n-1)

0

e)

(7)

t h e t r u e s o l u t i o n v a l u e s o f (1) a s s u r a e d m o n ¢ o n i c

and f o u r t i m e s c o n t i n u o u s l y d i f f e r e n t l a b l ¢ , x e [Xn_l, Xn]

R ' ( x ) > O,

n-

can be c o n s t r u c t e d d i r e c t l y

An-1 " (Yn - Y n - 1 ) / h n , n - N + I ,

Suppose now

R ( CI[XM, x N]

[xN, ~N],

a

t h e n from [ 9 J ,

for

281

l y ( x ) - R(x)l ( ~

~c

Ily'lt

max ( l y ' ( X n _ l ) - f n . l l ,

Ily°II +

(Ily(4)II

l y ' ( x n) - f n l )

Ily(3)112 + ? I l y ( 2 ) I I

hn

ily(3)II) (8)

where

c

is a constant

independent of

1 c ) ~

it.it

and

, min

,

n

satisfylrg

y'(x),

IX N, x ~ j

iS t h e u n i f o r m norm on

(9)

[xM, XN].

i t i s c l e a r t h a t the formula in (7) i s p a r t i c u l a r l y In all

f o r m s o f BDF codes.

we i n t e r p r e t

From ( 8 ) ,

for exact

values,

t h e a c c u r a c y as b e i n g a t b e s t O ( h 4 ) ,

h ~ O; t h e p r e c i s e

h - max hn, as n Form d e p e n d l , , g ~n t h e a c c u r a c y o f t h e d e r i v a t i v e

For a c t u a l comvuted v a l u e s

{yn } ,

and denote by

r~x)

We t r e a t

r ( X n . 1 ) - Y(Xn_l),

r(xn)

- Y(Xn)

r'(Xn-1)

r'(xn)

" fn,

-

the i n t e r v a l

R(x~

[Xn_ I , Xn]

the above r a t i o n a l I n t e r p o l a n t which s a t l s f i e s ,

" in-l,

x ~ [xn_ 1, X n ] , R(x)

values

the e r r o r in the i n t e r p o l a n t

can be t r e a t e d a s y m p t o t i c a l l y as f o l l o w s .

then for

{Yn}

simple to implement

y(x) -

p r o c e e d i n g as i n t h e polynomial ca~e,

(R(x) - r(x))

hence e m p l o y i n g the u n i f o r m norm on

ilR(x) - y(x)|!

(IIR(X)

+ (r(x)

- y(x)),

[Xn_ 1, X n ] ,

- r(x)l!

+

IIr(x)

Now, t h e second t e r m I s bounded a c c o r d i n g t o ( 8 ) .

- y(x)ll.

To t r e a t

we a s s u m e Pl

Yn " Y(Xn) + O ( h ~

Yn-1 " Y ( X n - l )

+ O(hn-1)

then from ( 7 ) , a s t r a i g h t f o r w a r d s u b s t i t u t i o n y i e l d s , r ( x ) - y(xn_ 1) + P n . l ( 0 ) / Q n - l ( 0 )

(10)

the first

term

282

e(x) - Yn-1 + P . . l ( e ) + 0(1~) en-l(e) + 0(~) h-

where

min (h n, hn_l),

and

q-

Now from t h e p r o o f o f t h e r e s u l t is unlformlyboundedbelow via

min (p, P l ) ,

i n (8) i n [ 9 ] , Qn_l(0) ) c ,

:" I s shown t h a t

where

c

Qn_l(9)

satisfies

(9).

Hence, t h e r e f o l l o w s , R ( x ) - Yn-1 + P n - l ( O ) / Q n - l ( O )

and f i n a l l y , lie(x)

+ O(hq),

(10) g i v e s ,

- yfx)ll

¢ o([q)

+ hnAn(f)max(ly'(Xn_l)

- fn_l l,

Iv'tXn)-fnl)

4 ÷ hn Bn(f),

where t h e c o n s t a n t s

An(f), Bn(f)

accuracy of the interpolant (Yn, i n )

and

(Yn-l, in-l)

consists of errors (here it

are defined appropriately.

is naturally values.

d e t e r m i n e d by t h e a c c u r a c y o f t h e

As in the polynomial case,

i s not p o s s i b l e t o r e l a t e

error

this contribution to the local error, process).

le,

At b e s t t h e o v e r a l l a c c u r a c y

O(h4). For t h z s p e c i a l c a s e o f c o n s t a n t o r d e r and s t e p s i z e

order p

BDF y i e l d s

f n ~ y ' ( x n ) + 0(hP - 1 )

lie(x) - y(X)ll as e x p e c t e d . satis:~ctory

h

throughout,

the

hence

~ 0(hP) + 0 ( h 4 ) ,

The above i n t e r p o l a n t

can t h e r e f o r e be e x p e c t e d t o be

f o r low and moderate a c c u r a c y when implemented i n BDF c o d e s .

in a d d i t i o n i t has t h e a d v a n t a g e o f p r o v i d i n g a g l o b a l l y compared t o t h e

4.

the error

i n t h e d a t a v a l u e s , a l o n g w i t h pure i n t e r p o l a t i o n

owing t o t h e n o n l i n e a r i n t e r p o l a t i o n is

Hence t h e

CO

C1

lnterpolant

polynomial i n t e r p o l a n t .

Numerical R e s o l t s In t h i s s e c t i o n w ~ p r e s e r l t t h e r e s u l t s

of numerical t e s t s

carried

out w i : h t h e code D02QBF, i n which we compete t h e p o l y n o m i a l i n t e r p o l a t i o n

283

scheme b a s e d o n carried

Tp,n(X),

with

(26

o u t o n 29 systems

DETEST), E n r i g h t modified

the ratio~i

From t h e s t i f f

and P r y c e [ 5 ] ) .

with

system size

test

with

respect

defined

a range of

tolerance

[Xn_l,Xn]

m

50.

-

to the

~2 - norm;

Yn

Ynj

error

Ilynj

was c a r r i e d

y(x)

was a v a i l a b l e

real error

every step

different

high accuracy Our e x p e r i e n c e values

if,

was t h a t

this

as a good general

In the Following

tables

~

errors

norm,

~nJ-

The t r u e

For the Krogh problems the

Y(~nJ)

values

m -

were

an interpolation

throughout

and

Yn

the

In

value

values

- leading

In

unexpected

to a high value

the

to the accuracy.

of results

those

In .

s o o v e ~ r a l l we r e g a r d

we g i v e

For e a c h v a l u e

o f TOL t h e

max I n o v e r t h e r a n g e o f i n t e g r a t i o n For the polynomial n and rational (RAT) i n t e r p o l a t i o n schemes. In the rational c a s e we have to exclude

is

is not the only

t h e c o d e may h a v c o b t a i n e d

happened very rarely, guide

would be

the renge o£ integration,

However, since

process

of

necessarily

in

t h e c o d e DCO3 [4] w i t h t h e s m a l l e s t

From u n i t y .

Yn-1

to the global

value at the point

Essentially,

the process,

in the

9

- y(xn_l) ll , ily n - y(Xn) ll )

Otherwise

using

o f TOL.

as satisfactory

guide to assessing

to the

analytically

m - 50.

numerically

recommended value

not greatly

respect

- y(Enj)ll/maX[llyn_l

modification

regarded

out within

J - 1,2 .....

denotes the interpolated

approximated

v~lues

a large problem the

the problems were run over

was measured r e l a t i v e

by computing with

max 1(j(9

solution

In

s e t STlYrST ( s t i f f

at the points

and

and its

all

were

b y t h e v a l u e TOL.

interpolation

and t h e i n t e r p o l a t i o n

where

The t e s t s

The code u s e d a m i x e d l o c a l

~nJ " Xn-1 + J h n / l O ,

-

system test

In order to include

specified

F o r each p r o b l e m ,

In

R(x).

K r o g h p r o b l e m due t o Gear and Saad ~7] was used ( n o n l i n e a r ,

etgenvelues)

Yn-1

scheme

steps

[ X n _ l , Xn]

For which the

(POL)

284

monotonlclty conditions:

Yn-1 < Y n ,

f n - 1 > 0, f n > 0,

d e c r e a s i n g d a t a v a l u e s ) a r e not s a t i s f i e d . involved, i f at a l l ,

We a l s o count t h e

in the r a t i o n a l

interpolation

scheme was l e s s t h a n o r e q u a l t o t h e c o r r e s p o n d i n g p o l y n o m i a l e r r o r . i s t h e n e x p r e s s e d as a p e r c e n t a g e o f t h e t o t a l (RAT ~ POL).

number o f i n t e g r a t i o n

This steps

In a d d i t i o n we i n c l u d e t h e maximum o r d e r o f t h e 6DF f o r m u l a .

We ; ' e g a r d t h e f o l l o w i n g two s e t s o f r e s u l t s all

for

For most o f t h e pro~lems t h i s

o n l y a v e r y small number o f s t e p s .

number o f s t e p s f o r which t h e maximum e r r o r

(similarly

t h e tesL problems can be found i n C o l q u i t t

as t y p i c a l .

Results for

[2].

P r o b l e m , K r o g h (m-~)

TOL

POLY .

max i n . .

.

.

RAT

~ steps

Naximum

max ! n . . .

RAT ( POL

order over ran~e

!0 -3

1.69

1.17

81

3

10 - 4

1.34

1.17

79

4

10 - 5

1.44

1.29

76

5

10 - 6

1.3!

1.86

72

5

10 - 7

1.18

3.35

72

5

10 - 8

1.38

10.6

64

5

10 - 9

1.91

11.4

52

5

10 - 1 0

1.35

67

48

5

10 -11

3.12

123

/4/+

5

285

Large Krogh P r o b l em (m - 50)

TOL

POL¥ max I n

RAT max ! n

10 -1

1.22

1.07

91

1

! 9 "'2

1.69

1.04

93

3

10 -3

1.12

1.07

87

3

10-4

1,09

1.07

85

4

10 -5

1.13

1.08

82

5

10 - 6

1.08

1.74

69

5

10 -7

1.04

5.98

58

5

10 -8

1.06

7.38

43

5

10 -9

1.25

8.38

34

5

On t h e b a s i s o f a l l

% steps RAT ( POL

the r e s u | t s

we f e e l t h a t

Naximum order over range

the i n t e r p o l a n t s

t h e way p r e d i c t e d by t h e t h e o r y i n t h e p r e v i o u s s e c t i o n s . order polynomial lnterpolant

adjusts

requirements, w h i l s t the r a t i o n a l the t a b l e s

show, t h e r a t i o n a l

itself

behave i n

The v a r i a b l e

to the loc~! accuracy

interpelant

is at best

fourth order.

An

scheme i s w e l l s u i t e d t o low t o moderate

a c c u r a c y , where t h e o r d e r o f t h e polynomial l n t e r p o l a n t Higher a c c u r a c y r e n d e r s t h e r a t i o n a l

is restricted.

s~heme u n s a t i s f a c t o r y .

A l s o , f o r a i ! t h e p r o b l e m s , we i n v e s t i g a t e d t h e m o n o t o n i c i t y b e h a v i o u r o f ~p,n(X) w

~p,n ( x ) ,

for all

appropriate steps

we found t h a t

By c o n s i d e r i n g

t h e i n t e r p o l a n t was monotonic ~or a l l

smal! p ~ r c e n t a g e o f s t e p s following

[Xn_l, Xn].

a turning point,

( i - 2%). ~p,n(X)

t h e n t h e d a t a used in i t s c a l c u l a t i o n

For one o r two s t e p s

but a v e r y

immediately

c o u l d o f t e n l a c k m o n o t o n i c i t y - but was not g e n e r a l l y m o n o t o n i c .

286

5.

Conclusions We regard the polynomial i n t e r p o l a t i o n scheme as used i n f i x e d

c o e f f i c i e n t BDF t m p l e ~ n c z ~ t o n as a v e r y s a t i s f a c t o r y scheme,

it

is

r e l i a b l e and robust. The o n l y c r i t i c i s m

Is t h e p o s s i b i l i t y

of producing a slight

a c c u r a c y when b e g i n n i n g t o make l a r g e r e l a t i v e e m e r g i n g from t h e . f a s t possibility

C1 l n t e r p o i a n t s ,

a c c u r a c y i s not e s s e n t i a l .

6.

stage ( u s u a l l y at

low o r d e r ) . Ap

that

is,

when

This

~hen

I f smoothness i s o f importance t h e r a t i o n a l

incorporated to give

available

increases in step sizes,

b e i n g s u g g e s t e d by t h e l a r g e v a l u e s o f

- 5h, 10h. easily

transient

d e c r e a s e in

scheme can be v e r y

providing high

However, t h e polynomial scheme must s t i l l

be

i n o r d e r t o deal w i t h t u r n i n g p o i n t s i n t h e s o l u t i o n .

Acknowledgement The a u t h o r s a r e g r a t e f u l

t o P r o f e s s o r W. E n r i g h t f o r h i s h e l p f u l d i s c u s s i

on d u r i n g t h e p r e p a r a t i o n o f t h i s p a p e r . REFERENCES 1.

BYRIqE, C.D. and HIBIIRARSH, A.C. (1975). A p o l y a i g o r i t ~ numerical solution of ordinary differential

f o r the

equations.

AOq T r a n s . on N a t h - S o f t w a r e . 1, 71-96. 2.

S.D. COLQUITT. (1986).

Polynomial and R a t i o n a l

numerical solution of stiff

Interpolation

s y s t e m s , N. Sc. t h e s i s ,

in the

University of

ICmchester. 3.

BO2QBF (1985),

4.

DCO3 (1983),

NAG l i b r a r y ,

Oxford, England.

Harwell s u b r o u t i n e l i b r a r y ,

mark 11.

Harwell, Oxfordshlre,

England. 5.

ENRICHT, W H. and PR¥CE, J.D. assessing Initial

(1983).

v a l u e methods.

Two FORTR/~packages f o r Dept. Comp. So. Tech. Report

No. 1 6 J / 8 3 , Univ o f T o r o n t o , Canada.

287

6.

ENRIGHT, W.H., JACKSON, K.R., NORSE]r, S.P. and THONSEN, P.G. (1985). I n t e r p o l a n t s f o r Runge-Kutta formulas.

Dept. Comp. So. Teeh.

R e p o r t . No. 180/85, Univ. o f Toronto, Carmda. 7.

CEAR, C.W. and SAAD, Y. (1983).

lterat!ve solution of linear equations

;~ O~F codes, S ! A N J . S e t . S t a r . Comput., 4, pp.583-601. 8.

CLAD~J.~Y, i . ,

$HAN~INE, L . F . , BACA, L.S. =..,~_ D ~ ,

R.W. (1985).

P r a c t i c a l a s p e c t s o f i n t e r p o l a t i o n in Runge-Kutta codes, Numerical A n a l y s i s Report No. 102,

U n i v e r s i t y of Manchester,

England. 9.

GREGORY, J.A. and DELBOIJRGO, R. (1982).

Ptecewlse r a t i o n a l q u a d r a t i c

I n t e r p o l a t i o n to monotonic d a t a , IMA J. o f Num. Anal. 2, pp. 123-130. 10. HINI~4ARSH, A.C. (1980).

LSODE and LSODI, two new i n i t i a l

ordinary differential

value

e q u a t i o n s o l v e r s , AOI-SICNIM News l e t t e r ,

15, No.4, pp.10-11. 11. JACKSON, K.R. and SACKS-DAVIS, R. (1978),

An a l t e r n a t i v e

o f v a r i a b l e s t e p s i z e m u l t l s t e p formulas f o r s t i f f

Implemenatiou

ODEs.

Dept. Comp. $c. Tech. Report. Mo.121, Univ. o f Toronto, Canada. 12. S~ARD, W.L. (1985). solving stiff

Defect and l o c a l e r r o r c o n t r o l in codes f o r

i n l ~ l a l - v a l u e problems.

Dept. Comp. Sc. Teeh.

R e p o r t , No.184/85, Univ. o f Toronto, Canada, 13. STETTER, H.J. (1978). ODE-solvers.

C o n s i d e r a t i o n s concerning a t h e o r y For

Numerlcal Treatment oF D i f f e r e n t i a l

Equations,

R. B u l i r s c h , R.D. G r i g o r l e f f and J. $chroder, e d s . Lecture Notes in Maths. 631, S p r l n g e r - V e r l a g , pp.188-200.