Finite difference end conditions and semi-cardinal interpolation with quintic B-splines

Finite difference end conditions and semi-cardinal interpolation with quintic B-splines

Applied Mathematics and Computation 237 (2014) 318–329 Contents lists available at ScienceDirect Applied Mathematics and Computation journal homepag...

411KB Sizes 1 Downloads 111 Views

Applied Mathematics and Computation 237 (2014) 318–329

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Finite difference end conditions and semi-cardinal interpolation with quintic B-splines Aurelian Bejancu ⇑, Michael J. Johnson, Hamid Said 1 Department of Mathematics, Kuwait University, PO Box 5969, Safat 13060, Kuwait

a r t i c l e

i n f o

Keywords: Semi-cardinal interpolation Quintic B-spline Approximation order End condition Not-a-knot Superconvergence

a b s t r a c t For quintic spline interpolation on a uniform grid, we formulate a class of end conditions in terms of homogeneous finite differences of B-spline coefficients. The analysis of the resulting schemes, including their approximation order, is based on the model of semi-cardinal interpolation at the set of non-negative integers. Numerical results are also presented to illustrate the implementation on a finite grid. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction An important factor that determines the approximation quality of a univariate spline interpolation scheme over a bounded interval is the choice of end conditions. To study ‘natural’ end conditions for a spline scheme with a single endpoint, Schoenberg [1] introduced the ‘semi-cardinal’ model for interpolation at the set of non-negative integers. In the cubic spline case, Bejancu [2] obtained an alternative construction of semi-cardinal interpolation via B-spline representations. This approach formulated a finite hierarchy of end conditions, including de Boor’s ‘not-a-knot’ end condition [3], in terms of finite differences of B-spline coefficients, and allowed a simplified analysis that established approximation orders of the resulting schemes ranging from 0 to 4. A similar approach previously employed by Bejancu [4], Bejancu and Sabin [5] led to the introduction and analysis of ‘natural’ and ‘not-a-knot’-type boundary conditions for bivariate interpolation with the quartic 3-direction box-spline M 222 and was also applied in [6] to the construction of subdivision surfaces for CAGD. Our paper is a twofold extension of the B-spline approach of [2]: not only do we treat here the case of univariate quintic splines, but we also consider an infinite family of end conditions expressed in terms of homogeneous finite difference (FD) operators of arbitrary integer order s P 0 (see Section 2). In particular, our FD conditions for s ¼ 6, 7; 8; 9, turn out to be equivalent to certain end conditions proposed via a different method by Behforooz and Papamichael [7] for the so-called ‘Eða; b; cÞ quintic splines’ on a finite grid. For each s, our construction and error analysis of semi-cardinal quintic spline interpolation are based on explicit computation of B-spline coefficients satisfying five-term linear recurrence relations (see Section 3). In comparison, the error bounds for interpolation with Eða; b; cÞ quintic splines on a finite uniform grid were established in [7] under the extra assumption that the norm of a certain inverse matrix is bounded independently of the number of interpolation points — a condition which requires separate verification for each relevant triple ða; b; cÞ. Note that a B-spline construction of quintic spline interpolation on a finite grid was proposed by Schoenberg and Sharma [8] for three types of end conditions (natural, complete, Euler–Maclaurin), but it did not include an error analysis. ⇑ Corresponding author. 1

E-mail address: [email protected] (A. Bejancu). Present address: Department of Applied Mathematics and Statistics, Stony Brook University, NY 11794-3600, USA.

http://dx.doi.org/10.1016/j.amc.2014.03.117 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.

A. Bejancu et al. / Applied Mathematics and Computation 237 (2014) 318–329

319

The implementation and numerical approximation order of quintic B-spline interpolation with FD end conditions on a bounded interval are demonstrated in Section 4. The results obey the approximation orders of the semi-cardinal model and also exhibit, for s ¼ 7; 8; 9, the additional property of gradient superconvergence at the grid points, which was first established theoretically by Behforooz and Papamichael [7, (3.11)] for a sub-class of Eða; b; cÞ quintic splines. We also verify this property numerically for s ¼ 10 and conjecture that it holds for any integer s P 7. Our final result, inspired by numerical evidence, proves the pointwise convergence of semi-cardinal Lagrange functions as s ! 1. 2. Semi-cardinal quintic spline interpolation Recall that a cardinal quintic spline is a function S : R ! R such that: (i) S is a quintic polynomial on ½k; k þ 1 for any k 2 Z, and (ii) S 2 C 4 ðRÞ. Denote by S 5 the space of all cardinal quintic splines. Also, let M 2 S 5 be the centered quintic B-spline defined by the binomial formula

MðxÞ ¼

  3 6 1 X 5 maxf0; ðx  kÞ g; ð1Þkþ1 5! k¼3 3þk

x2R

and note that M vanishes outside the interval ½3; 3. It is well known (see Schoenberg’s monograph [9]) that any cardinal spline S 2 S 5 has a B-spline series representation

SðxÞ ¼

X ak Mðx  kÞ;

x2R

ð1Þ

k2Z

for a unique sequence of coefficients fak gk2Z . To construct a quintic spline taking prescribed values at the set Zþ of non-negative integers, we consider the restriction s of a cardinal spline (1) to the positive axis Rþ ¼ ½0; 1Þ. This semi-cardinal restriction admits the B-spline representation

sðxÞ ¼

1 X

ak Mðx  kÞ;

x 2 Rþ

ð2Þ

k¼2

For k 2 Zþ , we will make use of the following evaluation formulae:

1 ðak2 þ 26ak1 þ 66ak þ 26akþ1 þ akþ2 Þ; 120 1 s000 ðkÞ ¼ ðak2 þ 2ak1  2akþ1 þ akþ2 Þ; 2 sð4Þ ðkÞ ¼ ak2  4ak1 þ 6ak  4akþ1 þ akþ2 ;

sðkÞ ¼

ð3Þ

ð5Þ

sþ ðkÞ ¼ sð5Þ  ðk þ 1Þ ¼ ak2 þ 5ak1  10ak þ 10akþ1  5akþ2 þ akþ3 ; ð5Þ

where sþ and sð5Þ  denote the right and left fifth derivatives of s. 2.1. Finite difference end conditions We now turn to the formulation of certain end conditions at 0 for a semi-cardinal quintic spline s in terms of its B-spline coefficients. First, note that the well-known ‘natural’ end conditions for quintic splines, namely s000 ð0Þ ¼ sð4Þ ð0Þ ¼ 0, are equivalent via (3) to the system

a2  3a1 þ 3a0  a1 ¼ 0; a1  3a0 þ 3a1  a2 ¼ 0:

ð4Þ

Alternatively, in analogy with the cubic spline ‘not-a-knot’ condition, we may require that s have a continuous fifth ð5Þ ð5Þ ð5Þ derivative at the interior knots k ¼ 1 and k ¼ 2, i.e. sð5Þ  ð1Þ ¼ sþ ð1Þ and s ð2Þ ¼ sþ ð2Þ. Therefore s becomes a single quintic polynomial on ½0; 3 and each of the knots at 1 and 2 is ‘not-a-knot’. Using the last line of (3), this is equivalent to the system

a2  6a1 þ 15a0  20a1 þ 15a2  6a3 þ a4 ¼ 0; a1  6a0 þ 15a1  20a2 þ 15a3  6a4 þ a5 ¼ 0:

ð5Þ

If D is the standard forward difference operator defined by Dak :¼ akþ1  ak , it follows that (4) and (5) can be written as:

Ds a2 ¼ Ds a1 ¼ 0;

ð6Þ

where s ¼ 3 and 6, respectively. We are thus motivated to consider these automatic endpoint conditions for an arbitrary integer s P 0.

320

A. Bejancu et al. / Applied Mathematics and Computation 237 (2014) 318–329

Definition 1. Let s be a fixed non-negative integer. A function s : Rþ ! R of the form (2) that satisfies (6) is called a semicardinal quintic spline with FD end conditions of order s. The space of all semi-cardinal quintic splines with FD end conditions of order s will be denoted by S þ 5;s . Note that, since sð4Þ is piecewise linear, imposing conditions (6) for s ¼ 4 amounts to sð4Þ  0, i.e. s is a cubic polynomial, on the interval ½0; 1. Also, for s ¼ 5, conditions (6) are equivalent to the property that s is a single quartic polynomial over the interval ½0; 2. On the other hand, as the next result shows, the above ‘not-a-knot’ conditions corresponding to s ¼ 6, along with three more end conditions of our hierarchy, have previously been introduced via different arguments in [7]. Proposition 1. The homogeneous FD end conditions (6) of order s ¼ 6; 7; 8, and 9 are equivalent to the left end conditions  21 1 satisfied by the Eða; b; cÞ quintic splines of [7] for ða; b; cÞ ¼ 33 5 ; 5 ; 5 ; ð9; 9; 1Þ; ð17; 33; 9Þ, and ð25; 61; 21Þ, respectively. Proof. In the context of quintic spline interpolation on a finite uniform grid, each parameter triple ða; b; cÞ generates two left end conditions, initially formulated in [7, (3.1)] by means of gradient values at the first five interpolation points. These conditions are then equivalently expressed [7, (3.18) & (3.32)] in terms of differences of fourth order derivative values at knots, as well as fifth order derivative jumps. In our B-spline setup, the jump discontinuity of sð5Þ at the knot k can be written via (3) as ð5Þ

dk :¼ sþ ðkÞ  sð5Þ ðkÞ ¼ D6 ak3 : Then, by [7, (3.34)], the left end conditions satisfied by the Eða; b; cÞ quintic splines with  21 1 ða; b; cÞ ¼ 33 ; 5 ; 5 ; ð9; 9; 1Þ; ð17; 33; 9Þ, and ð25; 61; 21Þ are equivalent to Dm d1 ¼ Dm d2 ¼ 0, for m ¼ 0; 1; 2, and 3, respectively. 5 In turn, these correspond to our FD end conditions of order s ¼ 6; 7; 8, and 9, respectively. h Remark 1. For semi-cardinal interpolation with cubic splines, a similar hierarchy of FD end conditions was formulated in [2] for s 2 f0; 1; 2; 3; 4g, where s ¼ 2 and s ¼ 4 correspond to ‘natural’ and ‘not-a-knot’ conditions, respectively. If, however, that hierarchy is also extended for arbitrary s, it can be shown that the end condition corresponding to s ¼ 5 is equivalent to condition ‘H3’ of [10] (also satisfied by the ‘Eð3Þ cubic splines’ of [11]), while the end condition for s ¼ 6 is equivalent to condition ‘H4’ of [10]. 2.2. Main results Let s P 0 be a fixed integer. The following result establishes, for power growth data, existence and uniqueness of semicardinal interpolation from the space S þ 5;s introduced in Definition 1. Theorem 1. Let the sequence fyj g1 of real numbers satisfy, for some constants K 1 > 0 and c P 0, the polynomial growth j¼0 condition c

j yj j6 K 1 ð1 þ j Þ;

j 2 Zþ :

ð7Þ

Then there exists a unique quintic spline s 2 S þ 5;s such that

sðjÞ ¼ yj ;

j 2 Zþ

ð8Þ

and, for each m 2 f0; 1; 2; 3; 4g,

j sðmÞ ðxÞ j6 K 2 ð1 þ xc Þ;

x 2 Rþ ;

ð9Þ

where K 2 is a constant that depends on K 1 and c. Let f : Rþ ! R be a continuous function having at most polynomial growth at 1. For each fixed parameter h > 0, define the interpolant sf ;h to f on the scaled grid hZþ by the formula 1

sf ;h ðxÞ :¼ sðh xÞ;

x 2 Rþ ;

ð10Þ

where s is the quintic spline of Theorem 1 corresponding to the data sequence yj :¼ f ðhjÞ; j 2 Zþ . In the statement of the next result we make the following convention: ð5Þ sf ;h ðxÞ

n  o 8   < max sf ;h ð5Þ ðkÞ; sf ;h ð5Þ ðkÞ ; if x ¼ k P 1; k 2 Z;  þ :¼   ð5Þ : sf ;h þ ð0Þ; if x ¼ 0:

Theorem 2. Let

s P 1 and j :¼ min fs; 6g. If f 2 C j ðRþ Þ and f ðjÞ is bounded on Rþ , then there exists cf ;s > 0 such that

    jm ðmÞ supf ðmÞ ðxÞ  sf ;h ðxÞ 6 cf ;s h x2Rþ

ð11Þ

321

A. Bejancu et al. / Applied Mathematics and Computation 237 (2014) 318–329

for 0 6 m 6 j  1 and all sufficiently small h. This establishes the approximation order minfs; 6g for semi-cardinal quintic spline interpolation with FD end conditions of order s. In particular, the quintic spline schemes with ‘natural’ and ‘not-a-knot’ end conditions have approximation orders 3 and 6, respectively. 3. Proofs 3.1. The uniqueness property in Theorem 1 To prove the uniqueness part of Theorem 1, we employ Schoenberg’s space [12, p. 413] of cardinal ‘null-splines’

S 05 :¼ fS 2 S 5 : SðkÞ ¼ 0; Due to (1), we have S 2

S 05

8k 2 Zg:

if and only if the B-spline coefficients of S 2 S 5 satisfy the homogeneous recurrence relations

ak2 þ 26ak1 þ 66ak þ 26akþ1 þ akþ2 ¼ 0;

8k 2 Z:

ð12Þ

To solve the associated characteristic equation

k4 þ 26k3 þ 66k2 þ 26k þ 1 ¼ 0;

ð13Þ

pffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi let k þ k ¼ z, so that z þ 26z þ 64 ¼ 0, with roots z1 ¼ 13  105; z2 ¼ 13 þ 105. Hence, Eq. (13) has the four roots 1 k1 ; k1 1 , k2 ; k2 , with 1 < k2 < k1 < 0, where 1

ka ¼

2

 qffiffiffiffiffiffiffiffiffiffiffiffiffi 1 za þ z2a  4 ; 2

a 2 f1; 2g:

From the theory of linear recurrences, it follows that the linear space S 05 is spanned by the four ‘eigenspline’ functions

Sa ðxÞ :¼

X k ka M ðx  kÞ;

x 2 R; a 2 f1; 2g:

k2Z

Due to the compact support of M, it also follows that, for a 2 f1; 2g; Sa ðxÞ decays exponentially as x ! 1 and has exponentially increasing oscillations as x ! 1, while the reversed behavior at 1 holds for Sa ðxÞ ¼ Sa ðxÞ. þ Now consider two splines in S þ 5;s , each satisfying (8) and (9). Then their difference s0 2 S 5;s inherits the growth behavior c s0 ðxÞ ¼ Oðx Þ, as x ! 1, and interpolates zero-data:

8k 2 Zþ :

s0 ðkÞ ¼ 0;

Let s0 replace s in (2). Clearly, this B-spline representation of s0 can be extended to that of a cardinal null-spline es 0 2 S 05 with es 0 ðxÞ ¼ s0 ðxÞ; 8x P 0, by choosing inductively the coefficients a3 ; a4 ,. . ., such that (12) holds. Hence there exist constants c1 ; c2 , for which

es 0 ðxÞ ¼ c2 S2 ðxÞ þ c1 S1 ðxÞ þ c1 S1 ðxÞ þ c2 S2 ðxÞ;

x 2 R:

Since the growth behavior of s0 at 1 implies c2 ¼ c1 ¼ 0, we obtain

s0 ðxÞ ¼

1  X c1 kk1 þ c2 kk2 M ðx  kÞ;

x P 0:

k¼2

Now the endpoint conditions (6) imply the system s s 2 c1 k2 1 ð1  k1 Þ þ c 2 k2 ð1  k2 Þ ¼ 0; s s 1 c1 k1 1 ð1  k1 Þ þ c 2 k2 ð1  k2 Þ ¼ 0;

ð14Þ

2 whose determinant is ð1  k1 Þs ð1  k2 Þs k2 h 1 k2 ðk2  k1 Þ – 0. Therefore c 1 ¼ c2 ¼ 0 and s0 ðxÞ ¼ 0; 8x P 0.

3.2. Lagrange scheme Let s 2 Zþ . The existence of the quintic spline s 2 S þ 5;s with the properties (8) and (9) of Theorem 1 will be obtained via a Lagrange representation. Specifically, for each j 2 Zþ , we construct in this subsection a bounded function Ls;j 2 S þ 5;s , which satisfies the Lagrange interpolation conditions

Ls;j ðjÞ ¼ 1 and Ls;j ðkÞ ¼ 0 for k 2 Zþ n fjg:

ð15Þ

Then s will be defined as a convergent Lagrange series

sðxÞ ¼

1 X j¼0

yj Ls;j ðxÞ;

x 2 Rþ :

ð16Þ

322

A. Bejancu et al. / Applied Mathematics and Computation 237 (2014) 318–329

Theorem 3. Given j 2 Zþ , there exists a unique sequence of real coefficients 1 X

Ls;j ðxÞ :¼

n

lj;k

o kP2

such that the function

lj;k Mðx  kÞ; x 2 Rþ ;

ð17Þ

k¼2

is bounded on Rþ , satisfies (15), and Ls;j 2 S þ 5;s . In addition, there exist positive constants K 3 and b depending only on s such that, for all integers j 2 Zþ and 0 6 m 6 5,

   ðm Þ  Ls;j ðxÞ 6 K 3 ebjxjj ;

8x 2 R þ ;

ð18Þ

   o    n  ð5Þ ð5Þ ð5Þ       where the left side is replaced by max  Ls;j  ðkÞ;  Ls;j þ ðkÞ if m ¼ 5 and x ¼ k P 1; k 2 Z, and by  Ls;j þ ð0Þ if m ¼ 5 and x ¼ 0. Proof. n o The well-known stability property of a cardinal B-spline series [13] implies that Ls;j is bounded on Rþ if and only if lj;k is a bounded sequence. To determine this sequence, we impose the end conditions (6), namely kP2

    s s X X s s ð1Þk lj;k2 ¼ 0 ¼ ð1Þk lj;k1 ; k k k¼0 k¼0

ð19Þ

as well as the interpolation conditions (15). The latter imply the system of recurrence relations:

lj;k2 þ 26lj;k1 þ 66lj;k þ 26lj;kþ1 þ lj;kþ2 ¼ 120djk ; k 2 Zþ ;

ð20Þ

where djk is the Kronecker delta. For each j P 0, all Eqs. (20) corresponding to k > j are homogeneous recurrences with characteristic Eq. (13), therefore their bounded solution has the form

lj;k ¼ Akk1 þ Bkk2 ; k P j  1

ð21Þ

for some constants A; B. On the other hand, if j P 1, the homogeneous Eqs. (20) corresponding to 0 6 k < j have the general solution k k lj;k ¼ C 1 kk1 þ D1 kk k 2 f2; 1; . . . ; j þ 1g 1 þ C 2 k2 þ D2 k2 ;

ð22Þ

for some constants C 1 ; C 2 ; D1 ; D2 . Hence, due to overlapping, for j P 1 we obtain three compatibility conditions: k k Akk1 þ Bkk2 ¼ C 1 kk1 þ D1 kk 1 þ C 2 k2 þ D2 k2 ;

k 2 fj  1; j; j þ 1g:

ð23Þ

Eliminating A and B from two of these conditions for k ¼ j  1, we get j j A ¼ C 1 þ Ek2j 1 D1 þ Fk1 k2 D2 ;

ð24Þ

j 2j B ¼ C 2 þ Gkj 1 k2 D1  Ek2 D2 ;

where E ¼

1k21 k22 k21 k22

1k4

1k4

1 2 1 ; F ¼ k1 k1 2 k2 k2 ; G ¼ k2 k1 k2 k2 . Then substituting A and B into condition (23) for k ¼ j and simplifying gives 1

2

2

    kj1 1  k21 D1 þ kj1 1  k22 D2 ¼ 0: 1 2

1

Further, consider the remaining interpolation condition (20) corresponding to k ¼ j, in which we use (22) to substitute lj;k for k ¼ j  2; . . . ; j þ 1, while lj;jþ2 is substituted from (21) for k ¼ j þ 2, with A; B from (24). Thus, using the fact that k1 and k2 satisfy the characteristic relation (13), after simplification we obtain

   1  k21 k21 k2 þ k1 þ k2 kjþ2 1

 D1 þ

  1  k22 k1 k22 þ k1 þ k2 kjþ2 2

D2 ¼

120ðk1 þ k2 Þ : k1 k2  1

From the last two displays, a straightforward calculation gives

D1 ¼ a1 kj1 ; where a1 ¼ ð1k

D2 ¼ a2 kj2 ; 120k21 k2

ð1k21 Þ

1 k2 Þðk1 k2 Þ

ð25Þ

; a2 ¼ ð1k

j 2 F Þk1 j 2 EÞk2

120k1 k22

ð1k22 Þ

1 k2 Þðk1 k2 Þ

. Updating (24) and simplifying, we obtain

j 1 k1 ; j 2 k2 :

A ¼ C 1 þ ða1 E þ a

¼ C1 þ a

B ¼ C 2 þ ða1 G  a

¼ C2 þ a

ð26Þ

Hence, for j P 1, by using (25) and (26) in (22) and (21), we can write

lj;k ¼ C 1 kk1 þ a1 kj1kjj þ C 2 kk2 þ a2 kj2kjj ; 8k P 2;

ð27Þ

A. Bejancu et al. / Applied Mathematics and Computation 237 (2014) 318–329

323

where C 1 and C 2 are to be determined from the end conditions (19). Indeed, on substituting (27) into (19), the resulting system for C 1 ; C 2 has the same nonsingular 2  2 matrix as that of (14). Thus, for each j P 1, there is a unique sequence n o of coefficients lj;k : k P 2 satisfying the first statement of the theorem. In the remaining case j ¼ 0, we use formula (21) to write the two boundary conditions (19) as s s 2 2 2 l0;2  Ak2 1  Bk2 þ Ak1 ð1  k1 Þ þ Bk2 ð1  k2 Þ ¼ 0;

ð28Þ

s s 1 Ak1 1 ð1  k1 Þ þ Bk2 ð1  k2 Þ ¼ 0:

Also, substituting (21) into the interpolation condition (20) for k ¼ 0 and recalling that k1 and k2 are roots of (13), we obtain 2 l0;2  Ak2 1  Bk2 ¼ 120;

therefore the first equation of (28) becomes s s 2 Ak2 1 ð1  k1 Þ þ Bk2 ð1  k2 Þ ¼ 120:

Thus A and B in this case satisfy a system of linear equations with the same non-singular matrix as that of (14), which detern o l0;k : k P 2 .

mines uniquely the corresponding sequence of coefficients

The localization property (18) now follows from a straightforward estimate of the series (17) and its term-wise derivatives, using the compact support of M and the next Lemma. h Lemma 1. There exists K 4 ¼ K 4 ðsÞ > 0, such that

    jkjj lj;k  6 K 4 jk2 j ;

8k P 2; 8j 2 Zþ :

ð29Þ

Proof. If j P max f1; s  1g, employing (27), the system (19) becomes:

C1 C1

ð1  k1 Þs k21

þ C2

ð1  k2 Þs k22

¼ a1

ðk1  1Þs k1s2

kj1  a2

ðk2  1Þs k2s2

kj2 ;

ð1  k1 Þs ð1  k2 Þs ðk1  1Þs j ðk2  1Þs j þ C2 ¼ a1 k1  a2 k2 : s 1 k1 k2 k1 k2s1

ð30Þ

This has a solution of the form

C 1 ¼ b1 ðsÞkj1 þ b2 ðsÞkj2 ; C 2 ¼ c1 ðsÞkj1 þ c2 ðsÞkj2 ; with constants b1;2 ðsÞ and c1;2 ðsÞ depending only on s. Substituting C 1 ; C 2 into (27) and using jk1 j < jk2 j, as well as j þ k ¼ jk  jj þ 2 min fk; jg, we obtain (29). Next, assume that 0 6 j 6 s  2. Then C 1 ; C 2 verify a system with the same left-hand side as (30) and with a modified right-hand side. Since j takes only a finite range of values in this case, it follows that, in absolute value, C 1 ¼ C 1 ðs; jÞ and C 2 ¼ C 2 ðs; jÞ admit a common upper bound which depends on s, but not on j. Hence, for k P j, a direct estimate in (27) implies again the decay property (29). On the other hand, there are only a finite number of pairs ðj; kÞ satisfying 2 6 k < j,     for which the corresponding expressions lj;k jk2 jjkjj admit a common upper bound depending only on s. Finally, if j ¼ 0, then (29) follows from a direct estimate in (21), using the coefficients A; B determined above. h Proof of existence in Theorem 1. As stated at the begining of the section, we define s by the Lagrange series (16). Due to the localization estimate (18) and the hypothesis (7), this series and its term-wise derivatives of order  4 are absolutely and uniformly convergent on compact subsets of Rþ . Further, (8) is a consequence of (15), while the growth property (9) holds via a standard argument [12, p. 416]. It remains to show that s 2 S þ 5;s . For an integer k0 P 0 and a fixed value x 2 ½k0 ; k0 þ 1, the representation (17) implies

sðxÞ ¼

1 X j¼0

yj

kX 0 þ3

lj;k Mðx  kÞ ¼

k¼k0 2

kX 0 þ3

ak M ðx  kÞ;

k¼k0 2

P where the coefficient ak :¼ 1 j¼0 yj lj;k is defined by an absolutely convergent series due to (29) and (7). Since, for each j 2 Zþ , the coefficients lj;2 ; . . . ; lj;s1 satisfy (19), we obtain the same end conditions for a2 ; . . . ; as1 , which completes the proof. h 3.3. Polynomial reproduction and approximation order The following result describes the polynomial reproduction properties of the Lagrange scheme (16).

324

A. Bejancu et al. / Applied Mathematics and Computation 237 (2014) 318–329

Proposition 2. If

s P 1, then for any polynomial p of degree 6 min fs  1; 5g,

1 X pðxÞ ¼ pðjÞLs;j ðxÞ;

x 2 Rþ :

ð31Þ

j¼0

This identity fails to hold for

s P 0 and pðxÞ ¼ xmin fs;6g .

Proof. Note that the right-hand side of (31) represents the interpolant s 2 S þ 5;s of (16) corresponding to data values yj ¼ pðjÞ; j 2 Zþ . By the uniqueness property established in Theorem 1, the identity (31) holds if and only if p 2 S þ 5;s . To verify this condition, we recall, for each m 2 f0; 1; 2; 3; 4; 5g, the classical Marsden identities [9, p. 16], written here on the positive axis:

xm ¼

1 X

qm ðkÞM ðx  kÞ;

8x 2 R þ ;

k¼2

where

8 1; > > > > > k; > > > > 2 1 >
qm ðkÞ ¼

if if if

3

k  32 k; if > > > > > 4 2 4 > k  3k þ ; if > 5 > > > : 5 3 2 k  5k þ 5 k; if

m ¼ 0; m ¼ 1; m ¼ 2; m ¼ 3; m ¼ 4; m ¼ 5:

 m Since Ds k  0 if 0 6 m < s, it follows that p 2 S þ 5;s whenever p is a polynomial of degree not greater than min fs  1; 5g.  s 6 Further, Ds k  s! implies xs R S þ R Sþ h 5;s if 0 6 s  5, and it is obvious that x 5;s for any s P 6. Proof of Theorem 2. This proof follows along the lines of the cubic spline case treated in [2]. Recall that s P 1 and j :¼ min fs; 6g. For a fixed x 2 Rþ , denote by px the Taylor polynomial of degree j  1 that matches f and its derivatives of order 6 j  1 at x. Since f ðjÞ is bounded on Rþ , Taylor’s formula with Lagrange remainder implies

jf ðtÞ  px ðt Þj 6 K 5 jt  xjj ;

8t 2 R þ ;



 with K 5 :¼ j1! sup f ðjÞ ðnÞ : n 2 Rþ . Hence jf ðt Þj is of polynomial growth as t ! 1 and, for each scaling parameter h > 0, the interpolant sf ;h to f is obtained via formulae (10) and (16):

sf ;h ðxÞ ¼

1  X 1 f ðhjÞLs;j h x ;

x 2 Rþ :

j¼0

Next, the exponential decay property (18) implies the existence of a constant K 6 ¼ K 6 ðsÞ > 0 such that, for 0 6 m 6 j  1,

   ðm Þ  j2 ; Ls;j ðt Þ 6 K 6 ð1 þ jt  jjÞ

8t 2 Rþ ; 8j 2 Zþ :

Then, using the reproduction formula (31), we obtain:

     X  j2      X       m  ðmÞ  1 1 m j ðmÞ ðmÞ ðm Þ jhj  xj 1 þ h x  j f ðxÞ  sf ;h ðxÞ ¼ pðxmÞ ðxÞ  sf ;h ðxÞ ¼  ½px ðhjÞ  f ðhjÞLs;j h x h 6 K 5 K 6 h  j2Z j2Zþ þ   2 X  1  jm jm 6 K5K6h 1 þ  h x  j 6 c f ;s h ; j2Zþ

where cf ;s :¼ K 5 K 6 K 7 and K 7 :¼ supt2½0;1Þ

P

j2Zþ ð1

þ jt  jjÞ

2

< 1. h

   j Remark 2. By an argument of [2], if pðxÞ ¼ xj , then supx2Rþ pðxÞ  sp;h ðxÞ – o h , as h ! 0, i.e.the approximation order j is optimal. Theorem 2 can also be trivially stated for s ¼ 0. 4. Implementation and conclusions 4.1. FD end conditions on a bounded domain

n We consider interpolation at n þ 1 uniformly spaced knots in the interval ½0; 1. Let h :¼ n1 . Given a data sequence yj j¼0 , we seek a quintic B-spline interpolant of the form

sn ðxÞ ¼

nþ2 X k¼2

 1 ak M h x  k ;

x 2 ½0; 1:

ð32Þ

325

A. Bejancu et al. / Applied Mathematics and Computation 237 (2014) 318–329

The n þ 5 coefficients a2 ; . . . ; anþ2 are determined from n þ 1 interpolation conditions

sn ðjhÞ ¼ yj ;

j 2 f0; 1; . . . ; ng;

augmented by four homogeneous FD end conditions of order

    s s X X s s akþi2 ¼ 0 ¼ anþ2ki ; ð1Þk ð1Þk k k k¼0 k¼0

s 6 n þ 3:

i 2 f0; 1g:

We now demonstrate the implementation for the ‘not-a-knot’ boundary conditions corresponding to the ðn þ 5Þ  ðn þ 5Þ system of boundary and interpolation conditions takes the form:

s ¼ 6. In this case,

a2  6a1 þ 15a0  20a1 þ 15a2  6a3 þ a4 ¼ 0; a1  6a0 þ 15a1  20a2 þ 15a3  6a4 þ a5 ¼ 0; aj2 þ 26aj1 þ 66aj þ 26ajþ1 þ ajþ2 ¼ 120yj ;

j ¼ 0 : n;

ð33Þ

an5  6an4 þ 15an3  20an2 þ 15an1  6an þ anþ1 ¼ 0; an4  6an3 þ 15an2  20an1 þ 15an  6anþ1 þ anþ2 ¼ 0:

Proposition 3. For any n P 6 and any data sequence

n on yj , the linear system (33) admits a unique solution. j¼0

Proof. Denote by LðkÞ the k-th line of system (33). We apply ten elementary operations to the lines of this system, in the following order:

ð1Þ Lð3Þ :¼ Lð3Þ  Lð1Þ

& Lðn þ 3Þ :¼ Lðn þ 3Þ  Lðn þ 5Þ;

ð2Þ Lð3Þ :¼ Lð3Þ  32Lð2Þ

& Lðn þ 3Þ :¼ Lðn þ 3Þ  32Lðn þ 4Þ;

ð3Þ Lð4Þ :¼ Lð4Þ  Lð2Þ

& Lðn þ 2Þ :¼ Lðn þ 2Þ  Lðn þ 4Þ;

ð4Þ Lð4Þ :¼ Lð4Þ 

32 Lð3Þ 243

32 & Lðn þ 2Þ :¼ Lðn þ 2Þ  243 Lðn þ 3Þ;

1 1 Lð3Þ & Lðn þ 1Þ :¼ Lðn þ 1Þ  243 Lðn þ 3Þ: ð5Þ Lð5Þ :¼ Lð5Þ  243

Hence lines Lð3 : n þ 3Þ of the updated system (33) can be written as:

243a0  434a1 þ 626a2  474a3 þ 191a4  32a5 ¼ R0 ; 26281 8854 3922 4654 781 a1  a2 þ a3  a4 þ a5 ¼ R1 ; 243 243 81 243 243 6752 15412 2264 52 32 a1 þ a2 þ a3 þ a4 þ a5 ¼ R2 ; 243 243 81 243 243 aj2 þ 26aj1 þ 66aj þ 26ajþ1 þ ajþ2 ¼ 120yj ;

j ¼ 3 : n  3;

ð34Þ

32 52 2264 15412 6752 an5 þ an4 þ an3 þ an2 þ an1 ¼ Rn2 ; 243 243 81 243 243 781 4654 3922 8854 26281 an5  an4 þ an3  an2 þ an1 ¼ Rn1 ; 243 243 81 243 243  32an5 þ 191an4  474an3 þ 626an2  434an1 þ 243an ¼ Rn ; 40 where R0 :¼ 120y0 ; R1 :¼ 120y1  1280 y0 ; R2 :¼ 120y2  81 y0 ; Rn2 :¼ 120yn2  40 y ; Rn1 :¼ 120yn1  1280 yn ; Rn :¼ 120yn . 81 81 n 81 We now consider the sub-system obtained by ignoring the first and last lines of (34). It can be observed that the ðn  1Þ  ðn  1Þ matrix of this sub-system is strictly diagonally dominant, hence, by the classical Levy–Desplanques theorem (see [14]), this matrix is nonsingular. Therefore there exists a unique solution a1 , . . ., an1 of the sub-system. The remaining values a2 ; a1 , a0 ; an ; anþ1 ; anþ2 are then uniquely determined via back-substitution from the first three and the last three equations of the updated system (33). h

Remark 3. The above method of proof has previously been used in [7] for establishing the nonsingularity of similar matrices of linear relations satisfied by spline derivative values at the knots of a finite grid. Further, for several particular end conditions considered in [7], the error analysis of the resulting Eða; b; cÞ quintic splines is crucially based on the fact that the inverses of such matrices are bounded in the 1-norm independently of the number of interpolation points. However, since the details of this method depend on the specific numerical coefficients for each set of end conditions, a similar error analysis does not appear as possible for our FD end conditions of an arbitrary order s. This serves to underscore the relevance of the semi-cardinal model.

326

A. Bejancu et al. / Applied Mathematics and Computation 237 (2014) 318–329

4.2. Numerical results We work with the data function f ðxÞ ¼ ex and set yj :¼ f ðhjÞ, for j ¼ 0; 1; . . . ; n. The number of subintervals considered is n 2 f32; 64; . . . ; 2048g and the selected end conditions correspond to s 2 f3; 6; 7; 8; 9; 10g. For each s and n, we solve the ðn þ 5Þ  ðn þ 5Þ system of interpolation and end conditions. We then evaluate sn at 31 equi-spaced points between each pair of consecutive interpolation nodes, and compute

  

 hl  En :¼ max ðf  sn Þ hj þ : j ¼ 0; . . . ; n  1; l ¼ 1; . . . ; 31 32 

j

as a numerical approximation to the uniform error kf  sn kL1 ½0;1 . We expect the values of En to be proportional to h , where j :¼ min fs; 6g. Hence

En  2j E2n

)

j  jn :¼

The computed values of En and

E0n

1 En log : log 2 E2n

jn are displayed in the second and third columns of Tables 1–6. In addition, we evaluate



  :¼ max  f 0  s0 ðhjÞ : j ¼ 0; . . . ; n ; n

where, using (32),

s0n ðhjÞ ¼

 1  aj2  10aj1 þ 10ajþ1 þ ajþ2 : 24h r

Expecting that E0n is proportional to h , for some value

r  rn :¼

1 log 2

E0 log 0n E2n

r satisfying

;

we display the computed values of E0n and

rn on the fourth and fifth columns of Tables 1–6.

Remark 4. The above numerical computations are performed in quad precision (104-bit mantissa), using GNU Octave (similar to Matlab) in concert with GNU’s multiple precision library gmp. Note that En closely reflects the approximation order established in Theorem 2: for ‘natural’ end conditions (s ¼ 3) we have jn  3 in Table 1, while the rest of the tables, corresponding to s P 6, exhibit jn  6. The numerical error E0n of first derivatives at the knots also obeys Theorem 2 in Tables 1 and 2, where rn  2 and rn  5, respectively. f denote a sufficiently smooth extension of f on ½0; 1Þ, such that e f Remark 5. Given the smooth data function f on ½0; 1, let e has polynomial growth at 1. Further, let s be the semi-cardinal interpolant of Theorem 1 corresponding to the data sequence yj :¼ e f ðj=nÞ, j 2 Zþ . Then the observed agreement of the above tables with Theorem 2 justifies the view that, for large n, the behavior of the finite-grid interpolant sn ðxÞ is well-modeled by the scaled semi-cardinal interpolant sðnxÞ when x is close to the left-end point 0. Table 1

s = 3: ‘natural’ end conditions. n

jn

En

E0n

rn

32

2:2889  10

6

2.9911

4

4:5856  10

1.9922

64

2:8789  107

2.9956

1:1526  104

1.9961

128

3:6097  108

2.9978

2:8894  105

1.9981

256

4:5191  109

2.9989

7:2331  106

1.9990

512

5:6532  1010

2.9994

1:8095  106

1.9995

1024

7:0692  1011

2.9997

4:5253  107

1.9998

2048

8:8382  1012



1:1315  107



Table 2

s = 6: ‘not-a-knot’ end conditions E n

33 5

 ; 21 ; 1 of [7]. 5 5

jn

En

E0n

rn

32

11

3:2638  10

5.9589

7:7628  10

4.9598

64

5:2468  1013

5.9795

2:4944  1010

4.9799

128

8:3157  1015

5.9897

7:9043  1012

4.9899

256

1:3086  1016

5.9949

2:4874  1013

4.9950

512

2:0520  1018

5.9974

7:8003  1015

4.9975

1024

3:2120  1020

5.9987

2:4418  1016

4.9987

2048

5:0232  1022



7:6374  1018



9

327

A. Bejancu et al. / Applied Mathematics and Computation 237 (2014) 318–329 Table 3

s = 7: end conditions E (9, 9, 1) of [7]. n

jn

En

E0n

rn

32

13

9:5168  10

6.7648

2:0823  10

10

5.9489

64

8:7516  1015

6.6434

3:3709  1012

5.9744

128

8:7546  1017

6.4671

5:3612  1014

5.9872

256

9:8955  1019

6.2945

8:4515  1016

5.9936

512

1:2606  1020

6.1725

1:3264  1017

5.9968

1024

1:7477  1022

6.0905

2:0771  1019

5.9984

2048

2:5648  1024



3:2491  1021



jn

E0n

rn

Table 4

s = 8: end conditions E (17, 33, 9) of [7]. n

En 32

13

1:7951  10

6.0952

5:2190  10

12

7.0713

64

2:6256  1015

6.0235

3:8807  1014

7.2796

128

4:0363  1017

6.0047

2:4977  1016

6.7816

256

6:2861  1019

6.0005

2:2702  1018

6.1164

512

9:8187  1021

5.9998

3:2722  1020

6.0625

1024

1:5344  1022

5.9998

4:8960  1022

6.0322

2048

2:3979  1024



7:4813  1024



Table 5

s = 9: end conditions E (25, 61, 21) of [7]. n

jn

En 32

13

1:6263  10

64

2:5558  1015

128

E0n

rn

5.9916

5:0337  10

13

6.0130

5.9947

7:7948  1015

5.9981

4:0082  1017

5.9972

1:2195  1016

5.9968

256

6:2750  1019

5.9986

1:9097  1018

5.9971

512

9:8143  1021

5.9993

2:9899  1020

5.9985

1024

1:5342  1022

5.9996

4:6765  1022

5.9999

2048

2:3978  1024



7:3077  1024



En

jn

E0n

rn

Table 6

s = 10. n 32

13

5.9883

13

1:6219  10

64

2:5549  1015

4:9769  10

5.9885

5.9942

7:8387  1015

128

5.9984

4:0080  1017

5.9971

1:2261  1016

5.9998

256

6:2750  1019

5.9986

1:9161  1018

6.0000

512

9:8143  1021

5.9993

2:9939  1020

6.0000

1024

1:5342  1022

5.9996

4:6781  1022

6.0003

2048

2:3978  1024



7:3082  1024



Moreover, Tables 3 to 6 display the numerical result rn  6, wich is one order higher than that predicted for the uniform norm of the first derivative error in Theorem 2. For each of the Eða; b; cÞ end conditions considered in Tables 3–5, i.e. s ¼ 7; 8; 9, this ‘superconvergence’ result was established theoretically by Behforooz and Papamichael [7, (3.21)], using the method of analysis indicated in Remark 3. However, the FD end conditions corresponding to s ¼ 10 in Table 6 do not belong to the class Eða; b; cÞ studied in [7]. These observations, as well as similar numerical results for other smooth data functions f, support the following conjecture. Conjecture 1. For each integer s P 7 and each sufficiently smooth function f of polynomial growth at 1, the quintic interpolant sf ;h defined by (10) satisfies the superconvergence property

328

A. Bejancu et al. / Applied Mathematics and Computation 237 (2014) 318–329

     6 supf 0 ðhjÞ  s0f ;h ðhjÞ ¼ O h ;

as h ! 0:

j2Zþ

As proved in [15], gradient superconvergence at knots also implies several other superconvergence properties that occur at special points within a bounded domain, for all derivatives of order up to five. 4.3. Convergence of Lagrange functions as

s!1

For s; j 2 Zþ , recall that Ls;j 2 S þ 5;s is the semi-cardinal Lagrange function satisfying the interpolation conditions (15). The pointwise convergence result stated below has initially been conjectured based on Remark 5 and on numerical approximations of Ls;j via Lagrange functions constructed on a finite grid ½0; n \ Z, for large n. Let l0 ¼ 120; l1 ¼ l2 ¼ l3 ¼ 0, and, for k P 1, define lk recursively via

lk þ 26lk1 þ 66lk2 þ 26lk3 þ lk4 ¼ 0; 8k P 1: Letting, for each j 2 Zþ ,

Lj ðxÞ :¼

j2 X

ljk2 Mðx  kÞ; x 2 Rþ ;

k¼2

it follows that Lj is the unique semi-cardinal quintic spline of the form (2), such that: Lj ðjÞ ¼ 1; Lj ðkÞ ¼ 0 for k 2 f0; . . . ; j  1g, and Lj ðxÞ ¼ 0 for all x P j þ 1. Proposition 4. For each j 2 Zþ and x > 0; lim Ls;j ðxÞ ¼ Lj ðxÞ. s!1

Proof. Fix j 2 Zþ . We employ the form (21) of the B-spline coefficients lj;k of Ls;j for k P j  1. Substituting these coefficients into (20) written successively for k ¼ j; j  1, . . ., 0, and using the recursive definition of lk ; k P 0, and the fact that k1 ; k2 are roots of the characteristic Eq. (13), we obtain

lj;k ¼ Akk1 þ Bkk2 þ ljk2 ; k ¼ j  2; j  1; . . . ; 1; 2:

ð35Þ

Next, substituting (35) and (21) into the boundary conditions (19), it follows that

  j X s s s 2 ð1Þk Ak2 ljk ; 1 ð1  k1 Þ þ Bk2 ð1  k2 Þ ¼  k k¼0   j1 X s s k s 1 ð 1  k Þ þ Bk ð 1  k Þ ¼  ð 1 Þ ljk1 : Ak1 1 2 1 2 k k¼0 In the special case j ¼ 0, this system reduces to (28). Solving the above system for A and B via Cramer’s rule and using the fact that 1  k1 and 1  k2 are greater than 1, we deduce A; B ! 0 as s ! 1. Hence (35) and (21) imply

lim lj;k ¼

s!1



ljk2 ; if  2 6 k 6 j  2; 0;

if j  1 6 k:

The Proposition now follows via the B-spline representations of Ls;j and Lj . h Note that the limiting Lagrange functions Lj ; j 2 Zþ , exhibit large oscillations inside their support and do not generate a    practical Lagrange interpolation scheme. Indeed, the magnitude of lj  grows exponentially with j, causing Lj to loose the localization property (18) of Ls;j that holds for a fixed s. 5. Conclusion For quintic spline interpolation, in this paper we formulate a hierarchy of end conditions in terms of B-spline coefficients, using finite differences of arbitrary order s. The case s ¼ 3 corresponds to ‘natural’ end conditions, while s ¼ 6; 7; 8; 9 generate end conditions that have previously been described by Behforooz and Papamichael [7] in terms of derivative values at knots. We propose the semi-cardinal interpolation model as a framework for the convergence analysis of the resulting schemes. This model allows the treatment of the entire infinite hierarchy of end conditions and may also be adapted to other types of end conditions that can be expressed in terms of B-spline coefficients. Our theoretical results are confirmed by numerical tests performed in quad precision. We conjecture that the superconvergence of first-order derivatives at the knots (proved in [7] for s ¼ 7; 8; 9 and also observed numerically in our paper for s ¼ 10) holds for any s P 7.

A. Bejancu et al. / Applied Mathematics and Computation 237 (2014) 318–329

329

Acknowledgements A.B. is grateful to David Levin (Tel-Aviv University) and Nicolas Papamichael (University of Cyprus) for bringing to his attention references [7,11] back in 2003. The result of subsection 4.3 was prompted by a query of one of the referees regarding a preliminary conjecture. References [1] I.J. Schoenberg, Cardinal interpolation and spline functions VI. Semi-cardinal interpolation and quadrature formulae, J. Anal. Math. 27 (1974) 159–204. [2] A. Bejancu, On semi-cardinal interpolation with cubic splines, Kuwait J. Sci. Eng. 38 (1A) (2011) 25–37. [3] C. de Boor, The Method of Projections as Applied to the Numerical Solution of Two Point Boundary Value Problems Using Cubic Splines (Ph.D. thesis), University of Michigan, 1966. [4] A. Bejancu, Semi-cardinal interpolation and difference equations: from cubic B-splines to a three-direction box-spline construction, J. Comput. Appl. Math. 197 (2006) 62–77. [5] A. Bejancu, M.A. Sabin, Maximal approximation order for a box-spline semi-cardinal interpolation scheme on the three-direction mesh, Adv. Comput. Math. 22 (2005) 275–298. [6] M.A. Sabin, A. Bejancu, Boundary conditions for the 3-direction box-spline, in: M.G. Wilson, R.R. Martin (Eds.), Mathematics of Surfaces, Proc. 10th IMA Intl. Conf., LNCS, vol. 2768, Springer, Berlin, 2003, pp. 244–261. [7] G.H. Behforooz, N. Papamichael, End conditions for interpolatory quintic splines, IMA J. Numer. Anal. 1 (1981) 81–93. [8] I.J. Schoenberg, A. Sharma, The interpolatory background of the Euler–Maclaurin quadrature formula, Bull. Am. Math. Soc. 77 (1971) 1034–1038. [9] I.J. Schoenberg, Cardinal spline interpolation, CBMS-NSF Series in Applied Mathematics, vol. 12, SIAM, Philadelphia, PA, 1973. [10] T.R. Lucas, Error bounds for interpolating cubic splines under various end conditions, SIAM J. Numer. Anal. 11 (1974) 569–584. [11] G.H. Behforooz, N. Papamichael, End conditions for cubic spline interpolation, J. Inst. Maths. Appl. 23 (1979) 355–366. [12] I.J. Schoenberg, Cardinal interpolation and spline functions II. Interpolation of data of power growth, J. Approx. Theory 6 (1972) 404–420. [13] C. de Boor, Splines as linear combinations of B-splines. A survey, in: G.G. Lorentz et al. (Eds.), Approximation Theory II, Academic Press, New York, 1976, pp. 1–47. [14] R.A. Horn, C.R. Johnson, Matrix Analysis, Cambridge University Press, Cambridge, 1985. [15] G.H. Behforooz, N. Papamichael, Overconvergence properties of quintic interpolatory splines, J. Comput. Appl. Math. 24 (1988) 337–347.