Elimination techniques and interpolation

Journal of Computational and Applied Mathematics 19 (1987) 125-132 North-Holland

125

Elimination techniques and interpolation M. GASCA * Departamento Ecuaciones Funcionales, Facultad de Ciencias, Universidad de Zaragoza, Zaragoza, Spain

E. L E B R O N Escuela Universitaria de Arquitectura T~cnica, Universidad Polit~cnica de Madrid, Madrid, Spain Received 19 December 1985

Abstract: A general recurrence interpolation formula due to Gasca and L6pez-Carmona (J. Approx. Th. 34 (1982)) is studied by means of elimination techniques in linear systems. The repeated application of this formula is investigated in order to show its equivalence with extrapolation methods.

1. Introduction Miihlbach's results [10,11] on the generalized Newton and Aitken-Neville interpolation formulae have been the origin of other papers which treat the same subject from several points of view [1-4,6-9,12-15]. Among them, the papers by Hhvie [6-9] and Miihlbach [12-15] have shown that elimination techniques in linear systems lead in a natural way to recurrence interpolation formulae. In the present note we shall show how these techniques can be applied to obtain the main formula of [3]. In [3] sufficient conditions were given to guarantee the validity of the formula, and, in this connection, we shall prove now that one of these conditions is also necessary if the remaining conditions are assumed to hold. Also, the existence and unicity of the solution of an interpolation problem are studied by considering several simpler problems. We apply repeatedly recurrence interpolation formulae to show in a simple way the equivalence of this process with elimination techniques as used in extrapolation methods. Finally a trigonometric Hermite interpolation problem is solved by these techniques. More details can be found in [5].

2. Elimination techniques and interpolation Let K be an arbitrary field. Consider the problem of finding a = aTx

* Partially supported by NATO Research Grant 027-81. 0377-0427/87/$3.50 © 1987, ElsevierSciencePublishers B.V. (North-Holland)

(2.1)

M. Gasca, E. Lebran / Elimination and interpolation

126

where a o = system

(a01 , a02,...

,

aon+m) y E K "+'', n, m ~ N and x ~ K n+m is the solution of the

Ax=b, A = (aij)l
(2.2)

b = (bl, b 2 , . . . , bn+m) T E K "+m. Assume that there exist s (1 < s < m + 1) regular square submatrices A h (h = 1, 2 , . . . , s) in the n first columns of A. (Note that we do not assume at this stage that A is regular. This regularity is studied in Theorem 2). If we denote by O~lh< O[2h< " " " < Olnh (h = 1, 2,,.., s) the indices of the rows of A h and set {~h=(alh,~2h,...,a~h}c(1,2,...,n+m}, s = {1,2,...,n+m}[.J Ih,

(2.3)

h=l

then the problem (2.1)-(2.2) is equivalent to the s simultaneous systems (h = 1, 2 , . . . , s) n+m

Y'. aijxj = bi,

i ~ Ih U J,

j=l n+m E aojXj -- a = O.

(2.4)

j=l

On applying blockwise elimination, with A h as pivot, we can eliminate Xl, x 2 , . . . , x, in the last equation of (2.4), to obtain n+m

~.,

a ~ j x j - a = b h,

(2.5)

j=n+l

where aohj= aoj-- (aol, ao2,..., a o , , ) A ; l ( a,~lh, a,~2h,..., aa,h) T, (2.6) bh= --(am, ao2,...,ao,,)A;l(b,~,~,

ba2h,...,banh) T.

Let V be a vector space over K, { vl, v2,..., Vn+m} a set of n + m linearly independent elements of V and L o, Li, i = 1, 2 , . . . , n + m linear forms on V such that Liv j = aij ,

0 <~i <~n + m,

1 ~
(2.7)

If we denote ( k = 1 , 2 , . . . , n + m) Vk = span{v1,..., oh},

(2.8)

n+m

Pn+m = E

XiUi E Vn+m,

(2.9)

i=1

then the problem (2.1) (2.2) (and analogously the system (2.4)) can be written in the form LiPn+ m = L i f ,

Lop,,+m - a = O,

l K i <~ n + rn,

(2.10)

M. Gasca, E. Lebr6n / Elimination and interpolation

127

where f ~ V satisfies Lif=

bi,

(2.11)

1 <~ i <<,n + m .

Now (2.5) is the well-known Newton generalized formula _ r,,n+j] , , h I - a = - L o p , , h,

x n + j L o ( vn+j

(2.12)

j=l

where P,,h (resp. ,',,h""+J~is the unique solution in V, of the interpolation problem LiPn, h = L i f ,

(resp.

(2.13)

i E Ih

L""+J-Liv,+jir,,h-

,

(2.14)

i~Ih).

By simplicity, we shall denote in (2.12) r~+J ,h -- -

-"+J "

(2.15)

Un+j -- I)n,h

Multiplication of (2.5) by )t h, for h = 1, 2,..., s and addition of the resulting equations lead to Y'~ Y'~ ) t h a ~ j x n + j - a = j=s

h=l

(2.16)

)thb h h=l

provided that Xl, X > . . . , )t, satisfy

{

~ 1 -t- )k 2 -1- " ' "

X l a l j + ..

+)k s =

1,

+)~sa~j=0,

(2.17)

j=1,2,...,s-1.

From (2.6) and (2.13)-(2.15) we can write { ~h

__ r

.n+j

a°J - ~°r"'h '

(2.18)

bh = - L o p ~ , h

and therefore (2.16) is the general recurrence interpolation formula obtained in [3]. By similar reasonings Mfihlbach [15] has just obtained another slightly different formula. (See [4] for the differences between both formulae or consider any simple example with s ~ 1 and s ~ m). Let us now assume s

card U I h = n + s - 1 , h=l

card(IhC3Ih+a)=n-1,

(2.19)

l <~h<~s-1,

det( ait)l~t~n,ic{o)u(lhnih+O

4= O,

l~h~s-1,

and denote 1

1

...

1

1

all

a21

...

a~x

L 0 r n,1 n+l

al,_l

a2s 1

...

a 0ss - 1

D= L

n+s-1

OFn ,1

• " "

L

r n+l 0 n,2

L 0Fnn,2+ s - 1

""

1 I .n+l * J 0 r n ,s

.

L

n+s-1

0rn,s

(2.2o) Then we have the following theorem.

M. Gasca, E. Lebrbn / Elimination and interpolation

128

Theorem 1. Under the assumptions (2.19),

D 4=0 ¢* det(aij)i~i,u...

Uls,l<~j<~n+s_ 1

4=0.

The proof is quite simple and can be found in [5]. The most interesting particular case corresponds to s = m + 1. Then D 4= 0 ¢* A regular. If we further assume that s = 2, m = 1, then A is regular iff

L 0 r n,1 n+l

4=

L 0 rn,2 n+l

or, equivalently, _,+1 Lolan, 1 4= L o P_,+1 n, 2

(see [14])

An interesting application of this result to multivariate interpolation can be found in [3]. The case s = 1 (D = 1) is trivially included in the theorem. For s < m + 1 the regularity of A is not a consequence of the requirement D 4= 0. By blockwise elimination in (2.2) we arrive directly at the following theorem. Theorem 2. Assume that D 4=0, s < m + 1 and (2.19) holds. Then A is regular if and only if n +j d e t ( L t r n + s _ l ) t ~ J , s < j < m 4= O,

where S

tOl

J={1,2,...,n+m}-

h=l

and -nn+j +s--1

~---

vn + j

-"+J

--1)n+s--1,

S

~
pn+j n+s--1 E V n + s _ l ,

t e U Ih.

T*-'tl'n+s-1 nn+J = LtVn+j '

(2.22)

h=l

When A is regular, x,+s, xn+s_l,..., Xn+m can be computed from m

E Xn+jLtrn+J-e =

Ltr,+~-l,

t ~J,

(2.23)

j=s

where r, +~_ 1 is defined similarly to r2++~_1 with v~+j replaced by f. Note that S

L .tr,+s_l n + j = ~ -A h t L- t r n ,n+j h ,

j=s, s+l,

...,

m

h=l s

Ltr,+,-1 = 2 XhtLtr,,h,

(2.24)

h=l

where, for each t

2 )~h,

h=l

1,

E x, t h t ~ _' ~' t"Z+n ,Jh= 0 ,

j=l,

2, • . . , s - - l ,

(2.25)

h=l

and L't" r ,..+J n,h , Ltr~,h can be obtained in (2.4) by blockwise elimination similarly to (2.5) (see (2.18)).

M. Gasca, E. Lebrbn / Elimination and interpolation

129

3. Repeated elimination The above elimination technique can be repeatedly used if adequate conditions hold, and seen to be equivalent to extrapolation methods (see [1,2,6-8,13,14]. Consider the system (2.1)-(2.2) where, by simplicity, n + m is replaced by n, and assume that all the square matrices with the first na, /72,..., nq columns of A are regular (n I < n 2 < nq •

•

•

<

=/7).

Let us take n - n I + 1 subsets I,,,h, h = 1, 2 , . . . , n with cardinal/71 and n--nl+l U I,~h = {1, 2,---, n} • h=l

As in the previous section, we can construct n - n 1 + 1 systems (2.4) and apply blockwise Gauss elimination to obtain xjLo<,h-LoPn

= --LoPnlh,

h=l,2,...,n-n

I + 1.

(3.1)

j=nl+l

If we find n - n 2 + 1 subsets 1,2,h c {1, 2 , . . . , n}, h = 1, 2 , . . . , n - n 2 + 1 such that n2-na+l

ln2h =

U

Inlk,h,

k=l

lnlk,hE(Inll, lnl2,...,Inln_nl+l }

for each h,

(3.2)

and that (cf. (2.19)) card In2h = n2, card(In,k, h 0 In1k+ 1,h ) ~- /71 --

1,

k = 1, 2 , . . . , n 2 - nl,

(3.3)

det( Ltvj)lt~(l.l~,hnl.,k+,,Ou {o}--/=0, k = l , 2 , . . . , n 2 - n l ; then, as in (2.16), we can write n2--nl+l

nz--nl+]

~,

Xj

E

j=n2+l

~knlk,htOrJk,h

-- LoPn = -

E

•nlk,hLoPnlk,h,

(3.4)

k=l

k=l

that is

(3.5)

xjLrJ2,h _ Lopn = _ Lpn2,h. j=n2+l

In this way a situation analogous to (3.1) is reached and we can attempt to iterate the elimination process provided that conditions similar to (3.3) hold. At the last but one step we would obtain n

E

xjLorJq_l,h-Lop,=-Lop, l_l,h,

h=l,2,...,n-nq_l+l,

(3.6)

j=nq_l+l and finally n--nq_l+l -- ZoPn --- -

E k=l

~knq_ak,ltOPnq_~k,l"

(3.7)

M. Gasca, E. Lebrbn / Elimination and interpolation

130

The table Lopnj plays the same role as Aitken or Neville's interpolation tables. If we only assume

det(aij)l<,i,i<,, 4:0,

t= l , 2 , . . . , q .

(3.8)

then we can use analogously the generalized N e w t o n formula in the elimination process, giving rise to the blockwise Gauss elimination method.

4. Example: trigonometric Hermite interpolation Let P>~+I be the element of the space Pm= span{ 1, cos x, sin x , . . . , cos mx, sin mx } determined by the conditions

P2m+l(ti)=f(ti), i = 1 , 2 , . . . , m + 1 , p~m+l(ti)=f'(ti), i = 2 , 3 , . . . , m +

1,

(4.1)

where t i e [a, a + 2,~), a ~ R , i = 1, 2 , . . . , m + 1 and ti4: tj if i:/:j. For t ~ [a, a + 2v), t 4=t i Vi, p2m+l(t) can be obtained by the elimination process (3.1)-(3.7) with n = 2 m + 1, q = m , n j = 2 j + 1 ( j = 1, 2 , . . . , m). Denote

(Llf=f(tl),

L2J=f(tj+l),

j:

I L z j + l f = f (tj+a),

1, 2 , . . . , m ,

(4.2)

j=l,2,...,m,

~Lof = f ( t ) , and, for j = 1, 2 , . . . , m

I2j+1a={1, Z , . . . , Z j + 1}, I2j+1,2i={2i, 2i+ 1 , . . . , 2 i + 2 j } ,

i=l,2,__.,m-j,

I2j+1.2~+1={2i, 2 i + 2 , 2 i + 3 , . . . , 2 i + 2 j + 1 } ,

(4.3)

i=l, 2,...,m--j.

For h -- 1, 2 , . . . , 2m - 2k + 1 we have

I2j+l,h= I2j_l,h U I2j_l,h+l l.J Z2j_l,h+2, and therefore (3.3) and similar conditions hold for I2j+l,k,h = I2j+l,h+k_l, k = 1, 2, 3.

j = l, 2,..., m.,

(4.4)

The system to be considered is X 1 nt- X 2 COS t 1 + X 3 s i n

11 -t- • • •

-}-X2m

COS

mt 1 +

X2m+l

sin

mt I

=f(tl)

,

x I + x 2 cos t2 + x3 sin t 2 + • • • +X2m COS mt2 + X2m+l sin mt 2 = f ( t 2 ) , --x 2 sin t 2 + x 3 cos t 2 + . . . .

X2mm sin mt 2 + X2m+lm COS mt 2 = f ' ( t 2 ) , (4.6)

x I + x 2 cos tin+1 + "'" +X2m COS mtm+ 1 + X2m+l sin m t m + l : f ( t m l a ) , --X 2 sin tin+ 1 + . . . . Xzmm sin mtm+ 1 + X2m+lm COS mira+ 1 = f ' ( t m + l ) , Xa + x2 cos t + x 3 sin t + .. • +X2m COS mt + X2m+l sin mt --P2m+l(t) = O,

M. Gasca, E. Lebrbn / Elimination and interpolation

131

and the elimination process leading to (3.7) is easily performed by mean of the following scheme: For i = 0, 1 , . . . , m - 1, do F / = 1 - cos(ti+ 1 - ti+2) , A i = (1 - cos(t - t i + 2 ) ) / F i , B i = (cos(ti+

1 -

ti+2)

-

cos(/-

ti+2))/Fi,

C i = (sin(ti+2 - t) + sin(t D, = (1 - cos(t - t i + l ) ) / F i , E i ~- ( c o s ( / i +

1 -

ti+2)

-

ti+l)

cos(t

+

-

sin(t~+l - t~+2))/Fi,

ti+l))/Fi,

LoP3,2i+1 = A i f ( ti+l) - B i f ( ti+ 2) - C i f ' ( ti+ 2) , If i > 1 then LoP3,2i = - E i f ( t i + l ) + C i f ' ( t i + l ) + D i f ( t i + 2 ) ; For s = 2, 3 , . . . , m do L 0 r 32s ,2i+1 = COS st -- A i cos s t i + 1 -Jr- B i c o s s t i + 2 - C i s s i n s t i + 2 , L 0 r32s+1 , 2 i + 1 = sin st - A i sin sti+ 1 + B~ sin sti+ 2 + Cis cos sti+2, If i >~ 1 then Lor2~i = cos st + E i cos sty+ 1 + Cis sin sti+ 1 - D i cos sli+2, L0r3~,~/+1 = sin st + E i sin sti+ 1 - Ci s cos sti+ 1 - D i sin sti+ 2 enddo enddo.

( N o w the system 2m-2 E

r

j+3

X j + 3LOr~,h

-- L o P 2 m + l = -

LoP3,h, h = l , 2 , . . . , 2 m - 1

j=l

can be written.) For k = 2, 3 , . . . , m do For h = 1, 2 , . . . , 2 m + 1 - 2k, do R

- 1

r 2k

1

2k+1

r

2k

r

2k

•

~.2k+l

h -- ~-'O'2k-l,h+l~"OF2k-l.h+2 -- LOF'2k-1 h + 2 J t - ' 0 1 2 k - 1

2k-1 '

-I- I

2k

T

2k+1

'

1

.2k+1

h+l '

-- ~ O r 2 k - 1 h + 2 L o r 2 k - 1 h -- L o r 2 k - 1 h*-'Or2k-1 h + 2

-l-r 2k ' r 2k+l ' r 2k ' r 2~+1 -- Z o r 2 k - 1 h L o r ' ) k - 1 h + l - - L o r ' 2 k - 1 h+ l L o r 2 k - I h, __ . r 2k ' r 2k2~1 r 2k' -2k+1 \.n L r0 2 k - 1 h + 2 L r0 2 k - 1 h + l ) / I t 2 k - l , h , 2k-ll h - - ~ L r0 2 k - 1 h + l L r0 2 k - l , h + 2 - ~k ' ' __it 2k ' r 2k+l r 2k -~" 2k+l ~'~n 2k-12h -- (L0r2k-I h+ 2 L o r 2 k - 1 h -- L O r 2 k - 1 h L o r 2 k - 1 h+ 2 ) / l ~ 2 k - 1 h, ' ' 2k ' 2k+1 ' 2k ' 2~+1 ' )t 2 k - l , 3 , h - - ~ ( L r0 2 k - l , h L r0 2 k - l , h + l L r0 2 k - l , h + l L r0 2 k - l , h ) / R 2k-l,h ,

~.

L o P 2 k + l,h = ~ 2k_ l,l,h L o P 2 k _ l,h q'- )k 2k_ l,2,h L o P 2 k _ l,h + 1 q- ~ . 2 k _ l , 3 , h L o P 2 k _ l , h + 2

If k < m then For j = 1, 2 , . . . , 2 m r

2k+l+j_

Lor2k+l,h

2k, do

)t

-- ,, 2 k - 1 1 h _~_).

' ,

r

2k+l+j--

L r0 2 k + l . h f

~

r2k+l+j 1,h+2

-- '~2k-l,3,h~-'O'2k

enddo enddo enddo.

tA

r 2k+l+j L r0 2 k - l h + 1 2k-12h ' , ,

132

M. Gasca, E. Lebrbn / Elimination and interpolation

(At this point, for each k, the system 2rn- 2k

E

r

2k+l+j

X 2 k + l + j L O F ~ k + l,h

--

LoP2m+ l = _ L o P 2 ~ + l,h ,

h = 1, 2 , . . . , 2 m - 2 k + 1

j=l

can be written.) Finally we obtain L o P 2 m + l = LoP2m+l,1.

References [1] C. Brezinski, The Miihlbach-Neville-Aitken algorithm and some extensions, B I T 20 (1980) 444-451. [2] C. Brezinski, A general extrapolation algorithm, Numer. Math. 25 (1980) 175-187. [3] M. Gasca and A. L6pez-Carmona, A general recurrence interpolation formula and its applications to multivariate interpolation, J. Approx. Theor. 34 (1982) 361-374. [4] M. Gasca, A. L6pez-Carmona and V. Ramirez, A generalized Sylvester's identity on determinants and its applications to interpolation problems, in: Schemp and Zeller, Eds., Multivariate Approximation Theory H (Birkhauser, Basel, 1982) 171-184. [5] M. Gasca and E. Lebron, Elimination techniques and interpolation, Preprint Universidad Zaragoza, 1985. [6] T. H~vie, Generalized Neville type extrapolation schemes, B I T 19 (1979) 204-213. [7] T. H~vie, Remarks on the Mfihlbach-Neville-Aitken algorithm, Math. and Comp. no. 2/80, Dept. of Numer. Math. Univ. Trondheim, Norway. [8] T. Hhvie, Remarks on a unified theory for classical and generalized interpolation and extrapolation, B I T (1981) 465-474. [9] T. H~vie, An introduction to generalized interpolation/extrapolation with particular emphasis on Numerical quadrature. Part I, Math. and Comp. No. 1/83. [10] G. Mfihlbach, The general Neville-Aitken algorithm and some applications, Numer. Math. 31 (1978) 97-110. [11] G. Mi~hlbach, The general recurrence relation for divided differences and the general Newton interpolation algorithm with applications to trigonometric interpolation, Numer. Math. 32 (1979) 393-408. [12] G. Miihlbach, An algorithm approach to finite linear interpolation, in: E.W. Cheney, Ed., Approximation Theory 111 (Academic Press, New York, 1980) 655-660. [13] G. Miihlbach, On two general algorithms for extrapolation with applications to numerical differentiation and integration, in: De Bruin and Van Rossem, Eds., Pad~ Approximations and lts Applications, Amsterdam 1980, Lecture Notes Math. 888 (Springer, Berlin, 1981). [14] G. Miihlbach, Extrapolation algorithms as elimination techniques with applications to systems of linear equations, Rep. 152, Inst. fiir Math., Univ. Hannover, 19822 [15] G. Miihlbach, Two composition methods for solving certain systems of linear equations, Numer. Math. 46 (3) (1985) 339-350.