<<

. 60
( 83 .)



>>

n’1
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)}=
( (

<<

. 60
( 83 .)



>>