- Email: [email protected]
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.
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.