A( j) ÿi tli ;

= a(tl ) + ™(tl ) j6l6j + n: (2.4)

n

i=0

A closed-form expression for A( j) can be obtained by using divided di erences. In the sequel we

n

(s)

denote by Dk the divided di erence operator of order k over the set of points ts ; ts+1 ; : : : ; ts+k . Thus,

for any function g(t) deÿned at these points we have

«

¬ s+k 1 ·

s+k k

¬ ·

(s) (s)

Dk {g(t)} = g[ts ; ts+1 ; : : : ; ts+k ] = · g(tl ) ≡ cki g(ts+i ): (2.5)

¬

i=s tl ’ ti

i=0

l=s

i=l

Then A( j) is given by

n

Dnj) {a(t)=™(t)}

(

A( j) = : (2.6)

n

Dnj) {1=™(t)}

(

As is clear from (2.6), A( j) can be expressed also in the form

n

n

( j)

A( j) = ni a(tj+i ); (2.7)

n

i=0

( j)

where ni are constants that are independent of a(t) and that depend solely on the tl and ™(tl ) and

( j)

n ( j)

satisfy i=0 ni = 1. The quantity deÿned by

n

n

( j)

( j)

= | ni | (2.8)

n

i=0

(note that nj) ¿1) plays an important role in assessing the stability properties of the approximation

(

A( j) with respect to errors (roundo or other) in the a(tl ). As has been noted in

n

various places, if l is the (absolute) error committed in the computation of a(tl ); l = 0; 1; : : : ;

( j) ( j)

then |A( j) ’ An |6 nj) (maxj6l6j+n | l |), where An is the computed (as opposed to exact) value of

(

n

A( j) . Concerning nj) we have a result analogous to (2.6), namely,

(

n

n

|Dnj) {u(t)}|

(

( j)

( j)

= | ni | = ( j) ; (2.9)

n

|Dn {1=™(t)}|

i=0

where u(t) is arbitrarily deÿned for all t except for t0 ; t1 ; : : : ; where it is deÿned by

u(tl ) = (’1)l =|™(tl )|; l = 0; 1; : : : : (2.10)

This is a result of the following lemma that will be used again later in this paper.

A. Sidi / Journal of Computational and Applied Mathematics 122 (2000) 251“273 259

(s)

Lemma 2.1. With Dk {g(t)} as in (2:5); we have

k

(s) (s)

|cki |hs+i = (’1)s Dk {u(t)}; (2.11)

i=0

where hl are arbitrary scalars and

u(tl ) = (’1)l hl ; l = 0; 1; : : : ; (2.12)

but u(t) is arbitrary otherwise.

(s) (s)

Proof. The validity of (2.11) follows from (2.5) and from the fact that cki =(’1)i |cki |; i=0; 1; : : : ; k.

The results in (2.6) and (2.9) form the basis of the W-algorithm that is used in computing both

the A( j) and the nj) in a very e cient way. For this we deÿne for all j and n

(

n

Mn( j) = Dnj) {a(t)=™(t)}; Nn( j) = Dnj) {1=™(t)}

( (

Hn( j) = Dnj) {u(t)}

(

and (2.13)

with u(tl ) as in (2.10), and recall the well-known recursion relation for divided di erences, namely,

( j+1) ( j)

Dn’1 {g(t)} ’ Dn’1 {g(t)}

Dnj) {g(t)}

(

= : (2.14)

tj+n ’ tj

(See, e.g., [15, p. 45].) Here are the steps of the W-Algorithm:

1. For j = 0; 1; : : : ; set

M0( j) = a(tj )=™(tj ); N0( j) = 1=™(tj ) H0( j) = (’1)j =|™(tj )|:

and (2.15)

2. For j = 0; 1; : : : ; and n = 1; 2; : : : ; compute Mn( j) ; Nn( j) ; and Hn( j) recursively from

( j+1) ( j)

Qn’1 ’ Qn’1

Qnj)

(

= (2.16)

tj+n ’ tj

with Qnj) equal to Mn( j) ; Nn( j) ; and Hn( j) .

(

3. For all j and n set

Mn( j) |Hn( j) |

A( j) = ( j)

and = : (2.17)

n n

Nn( j) ( j)

|Nn |

Note that the W-Algorithm for A( j) was originally developed in [8]. The recursion for nj) was

(

n

given recently in [10]. Stability and convergence studies for GREP(1) can be found in [10], and more

recently in [13].

260 A. Sidi / Journal of Computational and Applied Mathematics 122 (2000) 251“273

Let us now assume that A(y) and A depend on a real or complex parameter and that we would

™ ™

like to compute (d=d )A ≡ A assuming that A is the limit or antilimit of (d=d )A(y) as y ’ 0+.

™

We also assume that (y) and ÿi in (2.1) are di erentiable functions of and that A(y) has an

asymptotic expansion as y ’ 0+ obtained by di erentiating that of A(y) given in (2.1) and (2.2)

term by term. Thus

∞ ∞

A(y) ∼ A + ™ (y) ™

™ ™ ir

ÿi yir

ÿi y + (y) as y ’ 0 + : (2.18)

i=0 i=0

Here ™ (y) ≡ (d=d ) (y) and ÿi ≡ (d=d )ÿi in keeping with the convention of the previous section.

™

™

We can now approximate A by applying the extrapolation process GREP(2) to (2.18). The approx-

™

imations Bnj) to B ≡ A that result from this are deÿned via the linear systems

(

(n’1)=2 n=2 ’1

ÿ1i yl + ™ (yl )

B(yl ) = Bnj) + (yl )

( ir ir

ÿ2i yl ; j6l6j + n; (2.19)

i=0 i=0

™

where B(y) ≡ A(y) as before. (Compare (2.18) and (2.19) with (1.7) and (1.8), respectively.) Now

™ ™

the Bnj) converge to A, but their rate of convergence to A is inferior to that of the corresponding A( j)

(

n

to A. We, therefore, would like to employ the approach of [14] hoping that it will produce better

results also with GREP(1) .

2.2. (d=d )GREP(1) and its implementation

Let us di erentiate (2.7) with respect to . We obtain

n n

™( j) ( j)

™( j) a(tj+i );

An = ni a(tj+i )

™ + (2.20)

ni

i=0 i=0

™

where ™( j) ≡ (d=d ) ( j) and a(t) ≡ (d=d )a(t) ≡ A(y).

™

ni

ni

™( j)

It is clear that, unlike Bnj) in (2.19) that depends only on a(t); An depends on both a(t) and a(t).

(

™ ™

( j)

™

Also the stability of An is a ected by errors both in a(tl ) and a(tl ). In particular, if l and Ál are the

™

( j)

™( j) ™

(absolute) errors in a(tl ) and a(tl ), respectively, then |An ’ An |6 nj) [maxj6l6j+n max(| l |; |Ál |)],

(

™

( j)

™( j)

™

where An is the computed (as opposed to exact) value of An , and

n n

( j)

| ™( j) |:

( j)

= | ni | + (2.21)

n ni

i=0 i=0

We shall call this extension of GREP(1) simply (d=d )GREP(1) .

™( j)

2.2.1. Computation of An

™ ( j)

Let us start by di erentiating A( j) = Mn( j) =Nn( j) : Upon denoting (d=d )Mn( j) = M n and (d=d )Nn( j) =

n

™ ( j)

N n , we have

™ ( j) Mn( j) N ( j)

™n

Mn

™( j)

An = ( j) ’ : (2.22)

[Nn( j) ]2

Nn

A. Sidi / Journal of Computational and Applied Mathematics 122 (2000) 251“273 261

™ ( j)

Now Mn( j) and Nn( j) are already available from the W-algorithm. We need only compute M n and

™ ( j)

N n , and these can be computed by direct di erentiation of (2.16) along with the appropriate initial

conditions in (2.15).

2.2.2. Computation of an upper bound on nj) (

™( j)

The assessment of stability of An turns out to be much more involved than that of A( j) , and it

n

™ ( j) .

requires a good understanding of the nature of M n

First, we note that, as the tl are independent of , Dnj) and (d=d ) commute, i.e., (d=d )Dnj) {g(t)}=

( (