On the spline-based wavelet differentiation matrix

MATHEMATICS Applied Numerical Mathematics 17 (1995) 33-45

ELSEVIER

On the spline-based wavelet differentiation matrix Leland Jameson Institute for Computer Applications in Science and Engineering, NASA Langley Research Center, Hampton, VA 23681, USA

Abstract

The differentiation matrix for a spline-based wavelet basis without decomposition at multiple scales is constructed. It is shown that regardless of the properties that one imposes on a spline basis, the differentiation matrix is always the one encountered when the basis is simply the B-spline basis. The order of accuracy of this differentiation matrix is shown to display superconvergence. Furthermore, it is shown that spline-based bases generate a class of compact finite difference schemes. With this knowledge boundary conditions of the proper order of accuracy can be imposed.

Keywords: Wavelet approximation; Differentiation matrix; Superconvergence; Compact finite difference schemes

I. Introduction

The use of wavelets as a basis set for the numerical solution of partial differential equations (PDEs) is a very active area of research. Wavelets, Daubechies or spline-based, provide a very convenient structure for dividing information according to the frequency of the information at a given location. This splitting of the data by projecting onto a wavelet basis seems ideal for the solution of nonlinear PDEs where shocks or high frequencies might exist in a small portion of the domain, whereas in the remainder of the domain the solution might be comprised of very low frequency information. In this paper the goal is to understand the character of a spline-based wavelet numerical method. Numerical results for such a method can be found in [14]. This understanding will be found through the differentiation matrix. In [5] the differentiation matrix was found for Daubechies wavelets with periodic boundary conditions, and a very high degree of differentiation accuracy known as superconvergence was proven to exist. Furthermore, it was seen that a Daubechies-based wavelet method corresponds to a finite difference method with grid refinement in regions where small scale data is present. In the more complicated scenario where periodic boundary conditions are no longer assumed, it was seen in [6] that the superconvergence is lost. 0168-9274/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 0168-9274(95)00005-4

L. Jameson/Applied Numerical Mathematics 17 (1995) 33-45

34

In this paper the differentiation matrix for a spline-based wavelet basis with periodic and nonperiodic boundary conditions will be constructed. It will be shown that a spline-based wavelet basis, also, displays superconvergence when periodic boundary conditions are assumed. Furthermore, the superconvergence is lost when periodicity is no longer assumed. In addition, it will be seen that spline-based wavelet methods generate a class of compact finite difference schemes.

2. D e f i n i t i o n s a n d f r a m e w o r k

Please note that the framework and even much of the notation in this paper come from a paper by M. Unser and A. Aldroubi, see [12].

2.1. B-splines All B-splines in this paper are central B-splines: centered at 0. The B-spline of order n is denoted by Bin(x). All splines of order n, denoted by sn(x), are linear combinations of B-splines of order n, see [8,9]:

sn(x)= ECkfln(X--k),

(1)

k

and the space spanned by splines of order n is denoted by S n. The B-spline fl°(x) is the box function, which is 1 for x ~ [ - 7, 1 ~-] 1 and zero otherwise. All B-splines of higher order are generated from /3°(x): ~n(x)

= 8 0 * /~n-l(x),

(2)

where " * " denotes the convolution operator. Other useful forms of this equation are

~2n+l(X)=fln* fln(x),

(3)

fl2n+l(x) = j_ ~jn(y)fln(x--y) dy,

(4)

and

and finally, ~2n+ l( x ) = ( ~ n , [~n( x -- " ) ) ,

where for

(5)

f,g ~ L2(~),

(f, g) = f _ ~ f ( x ) g ( x ) dx.

(6)

We will, also, need the samples of B-splines at the integers:

b"( k ) = fln( X ) l x=k.

(7)

In order to use B-splines to construct a multiresolution analysis it is necessary that contractions and expansions of nth-order B-splines also be contained in S n. This is only true when n is odd, which is, therefore, assumed for this paper.

L. Jameson/Applied Numerical Mathematics 17 (1995) 33-45

35

2.2. Fourier transforms o f B-splines

From (2) we see that the Fourier transforms of B-splines are generated recursively. That is, the Fourier transform of a convolution is the product of the Fourier transforms of the convolved functions. Therefore, the Fourier transform of a B-spline is /~n(~) = s i n c n + X ( ~ ) ,

(8)

where sinc(~) - sin(~)/~ and is the Fourier transform of/3°(x). 2.3. Spline-based wavelet bases

As mentioned above, the subspace of nth-order splines, S n, is spanned by nth-order B-splines. The following theorem will define the scaling functions for a spline-based wavelet basis, see [13]: Theorem. The set o f functions { ~bn( x - k): k e Y} with

(~n(x) =

E

pk~n( x-k

)

(9)

k= - ~

is a basis o f S n provided that the sequence {p} is an invertible convolution operator f r o m l 2 into itself.

The requirement that {p} be invertible allows one to find the B-spline expansion from a given scaling function expansion, i.e., the scaling function and B-spline expansions contain exactly the same information and {p} simply maps between the two. The sequence {p} is chosen such that the scaling function satisfies whatever properties are specified. That is, one can specify that the scaling functions and wavelets be orthogonal which places restrictions on {p}. There are many ways to define the scaling function, but in this paper it will be seen that in the calculation of the differentiation matrix that the sequence {p} will "divide out" producing a result which holds for all spline-based wavelet bases produced using the current framework. Throughout the paper I will use the notation ~b~(x) to denote translation of ~bn(x) by k: ch~,(x ) = c~n( x - k ).

(10)

Define the dual of ~bn, 4~", as ~n(x)=

Y'~ r k [ 3 " ( x - - k ),

(11)

k ~ -o~

where the sequence {r] is chosen such that 6 k = (~b0,'" 4,7,).

(12)

Later the relationship between {r} and {p} will be made precise. However, as mentioned above, {p} and {r} will have no effect on the differentiation matrix.

L. Jameson/Applied NumericalMathematics 17 (1995) 33-45

36

The wavelet is defined as oo

On(x) = E k ~

qkfl"(x-k),

(13)

-oo

and its dual is oo

E k=

(14)

--o¢

where fl is the "B-spline wavelet". For this paper this B-spline wavelet is not needed and will not be dealt with, see [12]. Furthermore, as with the scaling function, the sequences {q} and {s} defining the wavelet and its dual will be of no consequence in the calculation of the differentiation matrix.

3. Multiple scale wavelet decompositions One of the underlying strengths of a wavelet decomposition is that it is possible to easily decompose a function, or signal, into its components at various scales. To illustrate, suppose that we are restricted to a finite number of dimensions, say d. Denote the finest scale approximation subspace by V0, and the coefficients of the expansion in V0 by ~. A usual wavelet decomposition would appear as, V o = Wl (~)W2~ " ' " ¢ Wj(~) Vj,

(15)

where, as is customary, Wj denotes the wavelet subspace at scale j. A function projected onto the above subspaces would, therefore, be represented by d coefficients, say b". U n d e r the assumption of finite dimensionality, there will be an invertible matrix mapping from ~-"to b~: ~=P~. That is, if g(x) ~ Vo with the expansion, d-1

g(x) = E s°ch°(x),

(16)

i=0 and one projects into V1 and Wa by d/2-1 d-1

Prig(x) = ~ j=O

Y'. s°(ch °, ~])dp)(x)

(17)

i=0

and d/2-1 d-1

Pw, g(x) = E j=O

E Si° (~)i,o ~ I ) ) O ) ( X ) ,

(18)

i=0

then the matrix form of the above linear transformations is denoted P:

pj.i= -1 ), , ( bo, qSj

(19)

L. Jameson/Applied Numerical Mathematics 17 (1995) 33-45

37

and

ed/2+/,i = ( 4)/0, q;1),

(20)

for i = 0 , . . . , d - 1 and j = 0 , . . . , ½ d - 1. Suppose, now, that we have found the matrix D mapping from ~' to s, where s represents the coefficients of the approximation in V0 of the derivative of a function. Of course, no one would choose to work with wavelets and stop the decomposition in the subspace V0. The matrix D does, however, characterize differentiation in any combination of wavelet subspaces, and the coefficients ~ can be found from

c =PDP-I~.

(21)

To restate, D dictates the character of wavelet differentiation. The primary concern in this paper is to understand D. Therefore, the remainder of this paper will not mention wavelet decompositions at multiple scales.

4. Differentiation in a wavelet subspace

In this section and in the next section all domains will extend from - ~ to + ~. The goal here is to understand differentiation in a wavelet subspace. We begin with a function f ( x ) ~ L2(~). We project f(x) into a wavelet subspace, V0 cL2(~). The projection is as follows:

Pvof(x) = E ( f , 6k)C~(X),

(22)

( ~ k ' ~1 ) = ~k,l"

(23)

k where 4~ is the dual scaling function defined by

If we first differentiate (22) to get d

dxPvof(x) = Y'~( f, ~k)~k(X), k

and then project back into V0 we get d Pvo-~xPvof(x) = E ( f , 4~k)(4;k, @)~b~,(x). l,k

(24)

(25)

The inner product (4~k, @) is the key to the understanding of spline-based wavelet differentiation. When we restrict ourselves to finite dimensions, later in the paper, the l,kth element of the derivative projection matrix is exactly (4;k, @).

4.1. The scaling function and its dual As mentioned above, the following relationship exists between the scaling function and its dual:

(~k = ((~0,'n ~b~),

(26)

38

L. Jameson/Applied Numerical Mathematics 17 (1995) 33-45

where for spline wavelets ~b(x) and ~(x) were defined by ~bn(x) = E p d 3 n ( x - k ) k

(27)

~f)n(x)=

(28)

and EFkfn(x--k),

where the superscript "n" denotes the order of the spline and fin(x) denotes the nth-order B-spline. We will now see the result of inserting the sequences {r} and {p} into (26). In terms of the definitions of the scaling function and its dual, Eq. (26) is

(c~k, ~o) = f Y'.rmfl~( x - m) Eplfln(x --l-- k) dx. m

(29)

1

Note that, since all summations and integrals are from - ~ to 0% shifts in the dummy variable can be made without altering the summation and integration limits. Now, interchange the order of summation and integration as well as combine the summations to get

6k= Y ' ~ r m p t f f l " ( x - m ) f l " ( x - l - k

) dx.

(30)

l,m

Shift the integration variable and rename as

b2n+l I+k - m = f

fn (x)fln(x+m_l_k)dx

(31)

to get

(~k ~- E

(32)

F m P l O--2n+l l+k_ m .

l,rn

Reverse the time on the sequence {p}, P 't-- P - t , to get (~k

(33)

E P k - I E - k 2 n i+ml U l _ m l m

This double summation is a double convolution 6k= p' * r * b2n+l(k),

(34)

and the convolution theorem yields the product of the discrete-time Fourier transforms,

(35)

1 =~t(~)~(~)~2n + 1(~). Solve this equation for P(~:) to get,

^

F(~) = (pt(~)b2n+X(~))

-1

.

This relationship between P(~) and/~(~) will be needed in the next section.

(36)

L. Jameson/Applied Numerical Mathematics 17 (1995) 33-45

39

5. Theoretical accuracy

The accuracy of the spline-based differentiation matrix can be seen by observing the behavior of the semi-discrete Fourier transform of p

dq = ( c~o, C~q)

(37)

near ~ = 0. It will be seen that the Taylor series of a~(~) about ~: = 0 is of the form d(~) = isc + O(sc2n+3) + " " ,

(38)

and this implies accuracy of order 2n + 2. 5.1. The semi-discrete Fourier transform of {d} As above, define

dq a s

dq = ( ~o, ~q),

(39)

and recall from the previous section that this is dq= Ep,rmf[3n(x-m)~n(x-l-q)

dx.

(40)

l,m

Again, note that as of now all summations and integrals are from - ~ to ~. Define al+q_ m :

y [~n(x)[~n(x - "1-m -- l -- q) d x

(41)

to get dq = EPtrmat+q_m .

(42)

l,m

Again, define p't = P - t to get dq = Y'~p'q_tErmat_m. l m

(43)

We see that the right-hand side is a double convolution. We can, therefore, use the convolution theorem to get a~(~:) =/~'(~:)P(~:)~(~:).

(44)

Recall from the previous section that the following expression for P(~) was found: F(~) =

(/~'(~)b2n+l(~)) -1

(45)

Using this expression for ?(~:) in the expression for d(sc) we get a~(~) = ~2n+l(sc ) .

(46)

From this expression we have arrived at a point where the accuracy of a spline-based wavelet differentiation matrix can be stated in terms of the accuracy of a B-spline basis--a result more

L. Jameson/Appfied Numerical Mathematics 17 (1995) 33-45

40

than 20 years o l d - - a n d can be summarized as follows: Using a Galerkin projection, the differentiation matrix constructed from n-order B-splines is accurate of order 2n + 2, see [10]. This roughly doubling of the accuracy is known as superconvergence. If, on the other hand, the projection is collocation, then the differentiation matrix is equivalent to that found using the Galerkin projection but for a lower-order spline. That is, the differentiation matrix derived for an n-order B-spline basis using a Galerkin projection is the same as the differentiation matrix found for a (2n + 1)-order B-spline basis using a collocation projection, see [11].

6. Equations for finite dimensions

The equations that will be derived in this section will hold for any spline-based wavelet basis defined on a finite-dimensional domain. This includes the equations for both periodic and nonperiodic boundary conditions. Note that in this section all summations are from 0 to d - 1, and consequently the summation limits will often not be shown.

6.1. Quadrature Let g(x) denote the projection of f(x) ~

L2(~) in V0:

d-1

g(x)=Pvof(X)= Y'~ SlC~t(x),

(47)

l=0

where s, = ( f , @).

(48)

Note that in this section the notation bm(x) will be used to denote the B-spline basis functions. For periodic boundary conditions one could have simply used / 3 ( x - m), but for nonperiodic boundary conditions the basis functions near the boundary are truncated B-splines. The notation used here is, therefore, general enough for both periodic and nonperiodic problems. That is, for a nonperiodic domain bo(x) will be only the right-hand portion with support width 1 of the B-spline being used. Now, using the B-spline expansion of the scaling function in (47) we get g ( x ) = Y'.SlY'~Pl,mbm(X). l

(49)

m

Sample g(x) at the integers to get

gk = ~.,SlPt,mbm,k"

(50)

l,m

Utilize matrix and vector notation to get ~'= (pb)Tg '.

(51)

This expression is a relationship between the samples of g ( x ) ~ Vo and the exact scaling

L. Jameson/Applied Numerical Mathematics 17 (1995) 33-45

41

function coefficients. In general we need a relationship between the samples of f ( x ) ~ L2(~) and an approximate set of scaling function coefficients tT, f = (pb)Td.

(52)

6.2. The scaling function and its dual Recall once more the relationship between the scaling function and its dual,

ai,j =

4;;>.

(53)

In terms of the B-spline expansions we get,

6id= ~_, ~.,Pi,mri,, fbm(X)bl(X) dx. m

(54)

I

Let B denote the matrix form of the above integral. That is, B contains the projections of the B-splines projected onto other B-splines. The above equation becomes

Now, rename the expression in the brackets 3, to get, a i d = Eri,,yi,,.

(56)

l

Transpose r to get, T

t~i,j ~ £ Ti,lrl,j •

(57)

l

The above equation has the following matrix form:

I = PBR T.

(58)

This is the matrix form of the requirement that the scaling function be orthogonal to its dual under translation.

6.3. The derivative projection matrix Recall from an earlier section that the elements of the derivative projection matrix come from the following inner product: dj,i = ( ~ i , ~)j).

(59)

Replace the scaling functions with their B-spline expansions to get

d,,i= £pJ,,t f m(x)b,(x) dx. m,l

(60)

42

L. Jameson/Applied Numerical Mathematics 17 (1995) 33-45

Let the matrix A contain the results of the above integration, dj, i = E P i , m r j , l A l , m,l

m.

(61)

The matrix form of this equation is, D = R A P T.

(62)

6.4. The differentiation matrix The differentiation matrix is ..@ = C D C - ' ,

(63)

where

d=(eb) T.

(64)

From above we found that

D = RAP T,

(65)

I = PBRT.

(66)

and

From (66) we get / = R B T p T"

(67)

Solve for R and use the fact the B is symmetric to get, R = (BpT) -1.

(68)

Use this expression in (65) to get

D = (PT)-IB-1ApT.

(69)

Now, using the expression for the quadrature we get the differentiation matrix,

=bTpT(pT)-~B-1ApT(pT)-I(bT)-I,

(70)

or

~ = b T B - 1 A ( b T ) -1.

(71)

Let me restate two notational definitions. First, the elements of the matrix "b" are bm,k which denotes sampling basis "m" at position "k". Second, the elements of the matrix " B " are fbm(X)bk(X) dx. Note that in the case of periodic boundary conditions all three of the matrices b, A, and B are circulant. They, therefore, commute yielding B - ~ I as the differentiation matrix.

L. Jameson/Applied Numerical Mathematics 17 (1995) 33-45

43

7. Example with periodic boundary conditions

The first-order spline is the convolution of two box functions. It is defined a s /~I(x) =X "1- 1 for x ~ [ - 1 , 0 ] , / 3 1 ( x ) = l - x for x ~ ( 0 , 1 ] and 0 otherwise, the matrices A and B for first-order B-spline are the following:

A=

0

1/2

0

0

0

0

0

0

0

- 1/2-

- 1/2 0

0 - 1/2

1/2 0

0 1/2

0 0

0 0

0 0

0 0

0 0

0

0 0

0 0

- 1/2 0

0 - 1/2

1/2 0

0 1/2

0 0

0 0

0 0

0 0 0

0 0 0

0 0 0

0 0 0

- 1/2 0 0

0 - 1/2 0

1/2 0 - 1/2

0 1/2 0

0 0 1/2

0 1/2

0 0

0 0

0 0

0 0

0 0

0 0

-1/2 0

0 - 1/2

1/6 2/3 1/6 0 0 0 0 0 0 0

0 1/6 2/3 1/6 0 0 0 0 0 0

0 0 1/6 2/3 1/6 0 0 0 0 0

0 0 0 1/6 2/3 1/6 0 0 0 0

-2/3 1/6 0 0 0 B = 0 0 0 0 1/6

0 0 0 0 1/6 2/3 1/6 0 0 0

0 0 0 0 0 1/6 2/3 1/6 0 0

0 0 0 0 0 0 1/6 2/3 1/6 0

0

0 0 0 0 0 0 0 1/6 2/3 1/6

0 0 0

(72)

0 0 1/2 0

1/6 0 0 0 0 0 0 0 1/6 2/3

(73)

The derivative projection matrix D, which recall is the same as the differentiation matrix for periodic boundary conditions, is D = B - ~1.

(74)

Note that the matrix D is exactly the fourth-order compact differentiation matrix.

8. Nonperiodic boundary conditions

Here the differentiation matrix will be given for one particular boundary construction. The construction considered here is one of truncated B-splines. That is, away from the boundaries one simply shifts the B-splines by one to generate a spline basis. At the boundaries we simply continue this shifting. If part of the B-spline goes beyond the domain, then simply truncate. This construction maintains the approximation properties of the spline subspace. That is, any nth-order spline can be generated across the domain using nth-order B-splines truncated in this way.

44

L. Jameson/Applied Numerical Mathematics 17 (1995) 33-45

Using the boundary construction outlined in the previous paragraph and the first-order B-spline one gets the following matrices: -

-

A=

1 / 2

1/2 0 0 0 0 0 0 0 0

-

1/2 0 1/2 0 0 0 0 0 0 0

0 1/2 0 -1/2 0 0 0 0 0 0

0 0 1/2 0 -1/2 0 0 0 0 0

0 0 0 1/2 0 -1/2 0 0 0 0

0 0 0 0 1/2 0 -1/2 0 0 0

0 0 0 0 0 1/2 0 -1/2 0 0

0 0 0 0 0 0 1/2 0 -1/2 0

0 0 0 0 0 0 0 l/2 0 -1/2

0 0 0 0 0 0 0 0 1/2 1/2

(75)

and the matrix B is

B =

1/3

1/6

1/6 0 0 0 0 0 0 0 0

2/3 1/6 0 0 0 0 0 0 0

0 1/6 2/3 1/6 0 0 0

0 0 1/6 2/3 1/6 0 0

0 0 0 1/6 2/3 1/6 0

0 0 0 0 1/6 2/3 1/6

0 0 0 0 0 1/6 2/3

0

0

0

0

1/6

0 0

0 0

0 0

0 0

0 0

0 0 0 0 0 0 1/6 2/3 1/6 0

0 0 0 0 0 0 0

0 0 0 0 0 0 0

1/6

0

2/3 1/6

1/6 1/3

(76)

Computationally it can be shown that the differentiation matrix constructed from these two matrices is first-order accurate at the boundaries while maintaining fourth-order accuracy away from the boundaries. This leads to a degradation of accuracy such that the global accuracy is second-order. Similar matrices can be found for third-order splines, and the third-order spline differentiation matrix constructed from truncated B-splines is only third-order accurate at the boundaries. This low accuracy limits the global accuracy to fourth-order.

9. Conclusion

This paper was written in order to illucidate the numerics behind a spline-based wavelet basis. In summary, without data compression the differentiation matrix for a spline-based wavelet basis is a compact finite difference operator with the accuracy of 2n + 2 for an nth-order spline. That is, this accuracy holds regardless of whether the scaling functions are chosen to be orthogonal under translation or not. That is, the scaling functions can be required to be B-splines, C-splines, or orthogonal. The parameters which dictate the properties of the scaling function "divide out" when the differentiation matrix is calculated. Because of this

L. Jameson/Applied Numerical Mathematics 17 (1995) 33-45

45

feature, only one theory is necessary to illustrate the accuracy for the spline-based differentiation matrix when periodic boundary conditions are assumed. When the boundary conditions are, on the other hand, not periodic it has been shown that for boundary functions constructed from truncated B-splines that the superconvergence is lost at the boundary. I speculate that this loss of superconvergence will continue to be a characteristic of wavelet bases defined on an interval. When one begins to eliminate bases functions which have small coefficients, i.e., performs data compression, then one obtains an approximation to the noncompressed solution. In other words, a spline-based wavelet numerical method is equivalent to applying compact finite differencing on a wavelet defined grid. In addition, it is essential to know the accuracy of this numerical method in order to apply appropriate boundary conditions, notably boundary conditions with sufficient accuracy.

Acknowledgment This research was supported by the National Aeronautics and Space Administration under N A S A Contract No. NAS1-19480 while the author was in residence at the Institute for Computer Applications in Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23681. Research was also supported by A F O S R grant 93-1-0090, by D A R P A grant N00014-91-4016, and by NSF grant DMS-9211820.

References [1] [2] [3] [4]

R. Bauer and L. Jameson, Wavelet-based adaptive cell refinement for conservation laws, Preprint. I. Daubechies, Orthonormal basis of compactly supported wavelets, Comm. Pure Appl. Math. 41 (1988) 909-996. P. Davis, Circulant Matrices (Wiley-Interscience, New York, 1979). D. Gottlieb, B. Gustafsson, P. Olsson and B. Strand, On the superconvergence of Galerkin methods for hyperbolic IBVP, SIAM J. Sci. Comput. (submitted). [5] L. Jameson, On the wavelet based differentiation matrix, J. Sci. Comput. (1993). [6] L. Jameson, On the Daubechies-based wavelet interval, ICASE Report No. 93-94, NASA CR-191582 (1993). [7] L. Jameson, On the wavelet-optimized finite difference method, ICASE Report No. 94-9, NASA CR-191601 (1994). [8] I.J. Schoenberg, Contribution to the problem of approximation of equidistant data by analytic functions, Quart. Appl. Math. 4 (1946) 45-99, 112-141. [9] I.J. Schoenberg, Cardinal Spline Interpolation, CBMS-NSF Series in Applied Mathematics 12 (SIAM, Philadelphia, PA, 1973). [10] B. Swartz and B. Wendroff, The relative efficiency of finite difference and finite element methods I: hyperbolic problems and splines, SIAMJ. Numer. Anal 11 (1974) 979-993. [11] B. Swartz and B. Wendroff, The relation between the Galerkin and collocation methods using smooth splines, SIAM J. Numer. Anal 11 (1974) 994-996. [12] M. Unser and A. Aldroubi, Polynomial splines and wavelets--a signal processing perspective, in: C.K. Chui, ed., Wavelets--A Tutorial in Theory and Applications (Academic Press, New York, 1992) 91-112. [13] M. Unser, A. Aldroubi and M. Eden, The L 2 polynomial spline pyramid: a discrete multiresolution representation of continuous signals, IEEE Pattern Anal Mach. Intell. (to appear). [14] Wei Cai and Jian Zhong Wang, Adaptive wavelet collocation methods for initial value boundary problems of nonlinear PDE's, ICASE Report No. 93-48 NASA CR-191507 (1993); also: S/AM J. Numer. Anal (submitted).