<<

. 47
( 83 .)



>>

. .
. .
. .
Lj’1 ; u1 : : : : : : : : : Lj’1 ; uN
1
u1 ::: ::: ::: uN
˜j = (31)
det V Lj+1 ; u1 : : : : : : : : : Lj+1 ; uN
. .
. .
. .
LN ; u1 : : : : : : : : : LN ; uN
the adjoint Vadj of V equals

Vadj = (det V )C (32)

where C has entries cj; t deÿned by (29). Hence

V ’1 = C : (33)

It is remarkable [6] that in case of q simple ÿnite poles and a pole at inÿnity of multiplicity
n0 ; q + n0 = N; and N simple nodes the inverse of V can be factorized as
D1 0
V ’1 = V D2 ; (34)
0 H (s)
where D1 ; D2 are diagonal matrices of dimensions q and N , respectively, and where H (s) is a
triangular Hankel matrix of the form
« 
s1 s2 s3 : : : snq
¬s ·
s3 : : : ·
¬2 ·
¬ ·
¬s ::: · ·
H (s) = ¬ 3 ·:
¬· · ·
¬ ·
· 
snq

2.2. The Neville“Aitken interpolation formula

In [7] a Neville-Algorithm is given which computes the whole triangular ÿeld
u1 ; : : : ; uk
k=1; :::; N
(pik f)i=1; :::; N ’k+1 ; pik f = p f ∈ Uk (35)
ai ; : : : ; ai+k’1

of interpolants recursively where pik f agrees with the function f on {ai ; : : : ; ai+k’1 }. In [7] this
algorithm was derived from the general Neville“Aitken algorithm [11,3] via explicit formulas for
the Cauchy“Vandermonde determinants [17]. In [5] we were going the other way around and have
given a di erent proof of the Neville“Aitken algorithm which is purely algebraic. In this survey we
will derive the algorithm by a simple direct argument.
G. Muhlbach / Journal of Computational and Applied Mathematics 122 (2000) 203“222 209


Proposition 3. Let k ∈ N;
∈ Ck+1 ;
Ak+1 = (a1 ; : : : ; ak ; ak+1 ) = ( 1 ; : : : ; 1; 2; : : : ; p’1 ; ; p; : : : ; p)
m1 mp

with 1 ; : : : ; p pairwise distinct and m1 + · · · + mp = k + 1; Ak = (a1 ; : : : ; ak ); Ak = (a2 ; : : : ; ak+1 )
and a1 = ak+1 . Let Uk+1 = (u1 ; : : : ; uk+1 ) be a CV-system associated with the pole system Bk+1 =
(b1 ; : : : ; bk+1 ). Suppose Ak+1 © Bk+1 = …. Let p1 ∈ Uk interpolate f at Ak and p2 ∈ Uk interpolate
f at Ak and let p3 ∈ Uk+1 interpolate f at Ak+1 .
(i) If bk+1 ∈ C then
p1 (z)(z ’ ak+1 )(bk+1 ’ a1 ) ’ p2 (z)(z ’ a1 )(bk+1 ’ ak+1 )
p3 (z) = : (36)
(ak+1 ’ a1 )(z ’ bk+1 )
(ii) If bk+1 = ∞ then
p1 (z)(z ’ ak+1 ) ’ p2 (z)(z ’ a1 )
p3 (z) = : (37)
a1 ’ ak+1

Proof. (i) Call the right-hand side of (36) p3 . It belongs to Uk+1 in view of
˜
z ’ ak+1 bk+1 ’ ak+1 z ’ a1 bk+1 ’ a1
=1+ and =1+
z ’ bk+1 z ’ bk+1 z ’ bk+1 z ’ bk+1
by the partial fraction decomposition theorem. Obviously, p3 interpolates f at Ak+1 if all nodes are
˜
simple since the weights add to one and each of the unknown values p1 (ak+1 ); p2 (a1 ) has factor 0.
It is a consequence of Leibniz™ rule that this holds true also in case of multiple nodes. In fact,
d
p3 |z= i
˜
dz
bk+1 ’ a1 d z ’ ak+1 bk+1 ’ ak+1 d z ’ a1
= p1 (z) ’ p2 (z)
ak+1 ’ a1 d z z ’ bk+1 z= i ak+1 ’ a1 dz z ’ bk+1 z= i

bk+1 ’ a1 d d z ’ ak+1
= p1 (z)|z= |z=
ak+1 ’ a1 dz dz z ’ bk+1
i i
=0

bk+1 ’ ak+1 d d z ’ a1
’ p2 (z)|z= |z=
ak+1 ’ a1 dz dz z ’ bk+1
i i
=0

bk+1 ’ a1 d d z ’ ak+1
= f(z)|z= |z=
ak+1 ’ a1 dz dz z ’ bk+1
i i
=0

bk+1 ’ ak+1 d d z ’ a1
’ f(z)|z= |z=
ak+1 ’ a1 dz dz z ’ bk+1
i i
=0
d bk+1 ’ a1 z ’ ak+1 bk+1 ’ ak+1 z ’ a1
= f(z) ’ f(z)
dz ak+1 ’ a1 z ’ bk+1 ak+1 ’ a1 z ’ bk+1 z= i

d
= f(z)|z= i
dz
since the weights add to one. This is evident for all i = 1; : : : ; p and = 0; : : : ; mp ’ 1 except for
i = 1 and = m1 ’ 1 or i = p and = mp ’ 1. If i = 1 and = m1 ’ 1 then the unknown derivative
210 G. Muhlbach / Journal of Computational and Applied Mathematics 122 (2000) 203“222


(d=d z)m1 ’1 p2 (z)|z= 1 has the factor (z ’ a1 )=(z ’ bk+1 )|z= 1 which vanishes. Similarly, for i = p and
= mp ’ 1 the unknown derivative (d=d z)mp ’1 p1 (z)|z= p has the factor (z ’ ak+1 )=(z ’ bk+1 )|z= p
which vanishes.
(ii) Obviously, the right-hand side of (37) belongs to Uk+1 and satisÿes the interpolation conditions
as is shown by the same argument used in the proof of (i).

Remarks.

1. Letting bk+1 ’ ∞ in (36) gives an alternative proof of (37). Observe that (37) is the classical
Neville“Aitken recursion for interpolation by polynomials. It has to be used anytime a pole ∞
is inserted.
2. p3 (z) is a convex combination of p1 (z) and p2 (z) if a1 ; ak+1 ; z; bk+1 are real and a1 6z6ak+1 ¡ bk+1
or bk+1 ¡ a1 6z6ak+1 .
3. Proposition 2 constitutes an algorithm of arithmetical complexity O(N 2 ) to compute the triangular
ÿeld (35) recursively from initializations pi˜ f ∈ U˜ (˜¿1) which are generalized Taylor inter-
polants: pi˜ f agrees with f at (ai ; : : : ; ai+˜’1 ) where all nodes are identical ai = · · · = ai+˜’1 =: a:
From (26) and (27) immediately
l’1
d
pi˜ (f) = f(a)!1 (38)
dz
=0
with
1 ˜’ ’1 (d=d z) B˜ (a) +
!1 = (z ’ a) ∈ U˜ (39)
B˜ (z) =0 !
are derived.

2.3. Newton™s formula and the interpolation error

Given a complex function f which is deÿned and su ciently often di erentiable at the multiple
nodes of system (8) then (9) constitutes a system of linear equations for the coe cients cj of the
generalized polynomial
N
p = pf =: cj uj ∈ UN ; (40)
j=1

where the coe cient
u1 ; : : : ; uN f
N
cj = c1; j (f) =
a1 ; : : : ; aN j
will be referred to as the jth divided di erence of f with respect to the systems UN and AN . In
[6,9] for consistently ordered poles which are assumed to be simple if ÿnite and for simple nodes
algorithms for solving system (9) recursively are derived whose arithmetical complexity is O(N 2 ).
In this section we are going to derive Newton™s formula for the interpolant obtained in [14] and
a procedure to compute the divided di erences
u1 ; : : : ; uk+1 f
ai ; : : : ; ai+k k + 1
G. Muhlbach / Journal of Computational and Applied Mathematics 122 (2000) 203“222 211


recursively [14], see also [10,16]. For the latter a new short proof is given. This way we will
establish an algorithm to compute the interpolant (40) in the general case of multiple nodes and
multiple poles in Newton™s form recursively whose arithmetical complexity again is O(N 2 ). Later,
in Proposition 6 of Section 4, we will derive an algorithm solving the linear system (9) recursively
in the general case of multiple poles and multiple nodes avoiding the additional partial fraction
decomposition.

Proposition 4. If p1 f = : k c1; j (f)uj ∈ Uk interpolates f at Ak = (a1 ; : : : ; ak ) and p1 f =
k+1
k k
j=1
: k+1 c1; j (f)uj ∈ Uk+1 interpolates f at Ak+1 = (a1 ; : : : ; ak ; ak+1 ); then
k+1
j=1
k+1 k k+1 k
p1 f = p1 f + c1; k+1 (f)r1 uk+1 ; (41)
where
Ak (z) Bk (bk+1 )
k k
r1 uk+1 (z) = uk+1 (z) ’ p1 uk+1 (z) = (42)
Bk+1 (z) Ak (bk+1 )
with Ak the node polynomial associated with Ak ; Ak (bk+1 ):= k — (bk+1 ’ aj ) and with Bk ; Bk+1
j=1
the pole polynomials associated with the pole systems Bk and Bk+1 ; respectively. Furthermore;
u ; : : : ; uk+1 f
c1; k+1 (f) = 1
k+1
a1 ; : : : ; ak+1 k + 1
with ±
u1 ; : : : ; uk f u1 ; : : : ; uk f

 ’

 a2 ; : : : ; ak+1 k a1 ; : : : ; ak k


 if a1 = ak+1 ;

u1 ; : : : ; uk+1 f Ak’1 (bk )
Bk (bk+1 )
ak+1 ’a1
· ·
= (43)
(ak+1 ’bk+1 )—
a1 ; : : : ; ak+1 k + 1 Ak (bk+1 ) Bk’1 (bk )

 

 u1 ; : : : ; uk’1 ; f u1 ; : : : ; uk’1 ; uk
 det V
 det V

a; : : : ; a; a a; : : : ; a; a
k —
if in the second case all nodes are identical a1 = · · · = ak+1 = a. Here Ak’1 (bk+1 ):= j=2 (bk+1 ’ aj ).
For any z ∈ C \ Bk
Ak (z)
k k
r1 f(z) = f(z) ’ p1 f(z) = [a1 ; : : : ; ak ; z](Bk f) (44)
Bk (z)
with
k
[a1 ; : : : ; ak ; z](Bk f) = [a1 ; : : : ; ai ]Bk [ai ; : : : ; ak ; z]f + [a1 ; : : : ; ak ; z]Bk f(z) (45)
i=1

denoting the ordinary divided di erence where [a1 ; : : : ; ak ; z]Bk = 0 i at least one pole is at inÿnity.
Moreover; if bk+1 ∈ C is chosen arbitrarily;
Ak (z)Bk (bk+1 )
u ; : : : ; uk+1 f
r1 f(z) = 1
k
: (46)
a1 ; : : : ; ak ; z k + 1 Bk+1 (z)Ak (bk+1 )
k+1
Proof. Consider linear system (9) for the coe cients of p1 f:
«  « 
k+1
c1; 1 (f) L1 ; f
u1 ; : : : ; uk+1 . .
¬ ·¬ ·
. .
V =
 
. .
L1 ; : : : ; Lk+1 k+1
Lk+1 ; f
c1; k+1 (f)
212 G. Muhlbach / Journal of Computational and Applied Mathematics 122 (2000) 203“222


bordered by the equation
k+1

<<

. 47
( 83 .)



>>