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.