Lecture Notes 1 10

LECTURE NOTES: ADVANCED NUMERICAL COMPUTING - CS 543

Michael J. Johnson Spring 2008

1. Notation N := {1, 2, 3, . . . } is the set of natural numbers. N0 := {0, 1, 2, 3, . . . } is the set of non-negative integers. Z := {. . . , −2, −1, 0, 1, 2, . . . } is the set of integers. R denotes the set of real numbers.

Π denotes the space of univariate polynomials. Πk denotes the space of univariate polynomials of degree at most k (the degree of the 0 polynomial is −1). Pk,Ξ denotes the space of piecewise polynomials of degree k with knot sequence Ξ. C k [a, b] denotes the space of functions f which are k-times continuously differentiable on [a, b]. Sk,Ξ denotes the space of splines of degree k with knot sequence Ξ. Cc (R) denotes the space of compactly supported functions in C(R). 2. Piecewise Polynomials For k ∈ N0 := {0, 1, 2, . . . }, the space of univariate polynomials of degree at most k is denoted by Πk . Let Ξ = {a = ξ1 , ξ2 , . . . , ξN +1 = b} be a set of increasing nodes (ξi < ξi+1 ) and let k ∈ N0 . The space of piecewise polynomials of degree k with knots Ξ, denoted Pk,Ξ , is the collection of functions s : [a, b] → R which satisfy s|

[ξi ,ξi+1 )

∈ Πk |

[ξi ,ξi+1 )

for i = 1, 2, . . . , N,

where the above is taken with the closed interval [ξN , ξN +1 ] when i = N . Note that we have adopted the somewhat arbitrary convention that piecewise polynomials (pp) are continuous from the right. Typeset by AMS-TEX

1

2

ADVANCED NUMERICAL COMPUTING

Example 2.1. Verify that s belongs to P3,Ξ if Ξ = {0, 1, 3, 4} and  3 if 0 ≤ x < 1,   x +1 2 s(x) := 2(x − 1) + (x − 1) if 1 ≤ x < 3,   4(x − 3)3 + 2 if 3 ≤ x ≤ 4. Is s continuous?

Theorem 2.2. The dimension of Pk,Ξ is N (k + 1). Piecewise polynomials with Octave In Octave, the polynomial p(x) = p1 xk + p2 xk−1 + · · · pk x + pk+1 is represented by the vector [p1 , p2 , . . . , pk+1 ]. For example, the vector [3, 0, 1, −2] represents the polynomial p(x) = 3x3 + x − 2. Note that the vector [0, 0, 3, 0, 1, −2] also represents this same polynomial while the vector [3, 0, 1, −2, 0, 0] represents the polynomial x2 p(x) = 3x5 + x3 − 2x2 . Some Octave commands which relate to polynomials are polyval (evaluation of a polynomial), conv (multiplication of two polynomials), polyderiv (derivative of a polynomial), and polyinteg (integral of a polynomial). Piecewise polynomials are not a built in part of Octave; however, I have written some routines to handle them. To describe how we represent a piecewise polynomial in Octave, let s ∈ Pk,Ξ . Since s is composed of N polynomial pieces, we can write s as s(x) = si (x − ξi ) for x ∈ [ξi , ξi+1 ), 1 ≤ i ≤ N , where each polynomial si is of degree at most k. The coefficients of these polynomials si , 1 ≤ i ≤ N , are stored in an N × (k + 1) matrix S whose i-th row contains the vector representation of the polynomial s i . Thus the piecewise polynomial s is represented by the pair (Ξ, S). Some commands which I’ve written for piecewise polynomials are pp val (evaluate a pp), pp add (add two pp), pp mult (multiply two pp), pp deriv (differentiate a pp), pp integ (anti-derivative of a pp), pp def integ (definite integral of a pp), and pp knot insert (insert knots into a pp). Example 2.3. For example, the Octave representation of  3 if 0 ≤ x < 1,   x +1 2 s(x) := 2(x − 1) + (x − 1) if 1 ≤ x < 3,   4(x − 3)3 + 2 if 3 ≤ x ≤ 4 

1 is the pair Ξ = [0, 1, 3, 4] and S =  0 4

 0 0 1 2 1 0 . 0 0 2

Knot insertion If two pp’s have identical knots, then adding them or multiplying them is fairly straightforward; however, if they have different knots, then one must first insert knots, as needed, until both pp’s have been rendered with identical knots.

CS-543

3

Example (knot insertion). Let s be as in the previous example and let Ξ1 = {0, 1, 2, 3, 4}. Note that Ξ1 has been obtained from Ξ by inserting the knot 2. Find the representation of s with respect to the knot sequence Ξ1 (ie. as an element of P3,Ξ1 ). In the above example, we see that the computation required is that of expanding s 2 (x+1) in powers of x, when s2 (x) = 2x2 + x. This yields s2 (x + 1) = 2x2 + 5x + 3. In general, the work involved in a knot insertion is simply that of finding, for a given polynomial p(x) = p1 xk + · · · + pk x + pk+1 , the coefficients q1 , . . . , qk+1 such that the polynomial q(x) = q1 xk + · · · + qk x + qk+1 satisfies q(x) = p(x + τ ). In other words, we have to translate the polynomial p by a distance −τ . One can easily write an algorithm for this based on the binomial formula, however the number of flops needed for execution is about 2k 2 . A better algorithm, which uses only k 2 + 1 flops, is the following (assuming that p is the Octave representation of the polynomial p(x) above): q=p s=tau*p(1) for i=k+1:-1:2 q(2)=q(2)+s for j=3:i q(j)=q(j)+tau*q(j-1) end end It is easy repeat the above example using this efficient algorithm and verify that the same result is obtained. 3. Splines and B-splines The space Pk,Ξ contains functions with various smoothness properties (assuming N ≥ 2). The following definition allows us to give a rather precise categorization of the smooth functions contained in Pk,Ξ . Definition 3.1. The space of continuous functions on [a, b] is denoted C[a, b] (or C 0 [a, b]). For a positive integer k, we define C k [a, b] to be the space of functions f : [a, b] → R for which f, f 0 , . . . , f (k) are continuous on [a, b]. Theorem 3.2. For ` = 0, 1, . . . , k, the dimension of Pk,Ξ ∩ C ` [a, b] equals N (k − `) + ` + 1. Moreover, Pk,Ξ ∩ C k [a, b] equals Πk | . [a,b]

Of particular interest is the subspace Pk,Ξ ∩ C k−1 [a, b] which has dimension N + k. Definition 3.3. The subspace Pk,Ξ ∩C k−1 [a, b], denoted Sk,Ξ , is called the space of splines of degree k having knot sequence Ξ. Example 3.4. Determine whether s belongs to S2,Ξ if  2   x −2 s(x) = 2(x − 1)2 + 2(x − 1) − 1   −(x − 2)2 + 6(x − 1) + 3

Ξ = {0, 1, 2, 4} and if 0 ≤ x < 1 if 1 ≤ x < 2 if 2 ≤ x ≤ 4.

4

ADVANCED NUMERICAL COMPUTING

Theorem 3.5. Let s ∈ Sk,Ξ . a) If k > 0, then the derivative s0 belongs to Sk−1,Ξ . b) If se is an anti-derivative of s, then se ∈ Sk+1,Ξ .

The B-splines, which we now define, are important to the theory and applications of splines because they can be used as the elements of a basis for a given spline space Sk,Ξ . Let us begin by assuming, as usual, that we are given knots Ξ = {ξ1 , ξ2 , . . . , ξN +1 } satisfying a = ξ1 < ξ2 < · · · < ξN +1 = b. Additionally, we assume that Ξ has been extended to an infinite sequence of knots {. . . , ξ−1 , ξ0 , ξ1 , ξ2 , ξ3 , . . . } satisfying ξi < ξi+1 for all i ∈ Z. The exact values of the additional knots is of no concern, so long as the entire sequence is increasing. The B-splines of degree 0 are defined by  1 if x ∈ [ξi , ξi+1 ), 0 Bi (x) := 0 if x 6∈ [ξi , ξi+1 ). Recursive Definition 3.6. For k = 1, 2, 3, . . . , the B-splines of degree k are defined by     ξi+k+1 − x x − ξi k−1 k−1 k Bi (x) + Bi+1 (x). Bi (x) := ξi+k − ξi ξi+k+1 − ξi+1 Note that the B-spline is defined on all of R. In the important special case when the knots are simply ξi = i, the B-splines are called Cardinal B-splines and the above recursion reduces to
1 k−1 Bik (x) := (x − i)Bik−1 (x) + ((i + k + 1) − x)Bi+1 (x) ; k moreover, we have (in the cardinal case) Bik (x) = B0k (x+i) so that Bik is simply a translate of B0k . Example 3.7. Find a formula for Bi1 . Partition of Unity Theorem 3.8. Let k ∈ N0 . If lim ξi = ±∞, then i→±∞

∞ X

Bik (x) = 1 for all x ∈ R.

i=−∞

Properties of B-splines. For k ∈ N0 and i ∈ Z, the following hold: (i) Bik is a pp of degree k with knots ξi , ξi+1 , . . . , ξi+k+1 . (ii) Bik (x) > 0 if x ∈ (ξi , ξi+k+1 ), and Bik (x) = 0 if x 6∈ [ξi , ξi+k+1 ]. (iii) If k > 0, then     k d k k k−1 k−1 Bi (x) − Bi+1 (x). B (x) = dx i ξi+k − ξi ξi+k+1 − ξi+1 (iv) Bik is (k − 1)-times continuously differentiable (ie. Bik ∈ C k−1 (−∞, ∞)).   ∞ Z x ξi+k+1 − ξi X k+1 k Bi (t) dt = (v) Bj (x). k+1 −∞ j=i

The support of a function f : R → R is defined to be the closure of the set {x ∈ R : f (x) 6= 0}. One consequence of property (ii) is that the support of Bik is the interval [ξi , ξi+k+1 ].

CS-543

5

Question. Which B-splines have some of their support on the interval (a, b)? The following Corollary extends properties (iii) and (v) to a linear combination of Bsplines. ∞ X ci Bik (x). Corollary 3.9. Let f (x) = i=−∞

(i) If k > 0, then

 ∞  X ci − ci−1 d f (x) = k Bik−1 (x). dx ξ − ξ i+k i i=−∞

(ii) If k ≥ 0 and only finitely many {ci } are nonzero, then Z x ∞ X f (t) dt = di Bik+1 (x), −∞

i=−∞

i X 1 where di = cj (ξj+k+1 − ξj ). k + 1 j=−∞

As mentioned above, the knots in Ξ are denoted a = ξ1 < ξ2 < · · · < ξN +1 = b, and we assume that we have chosen additional knots ξi for i < 1 and i > N + 1, maintaining ξi < ξi+1 for all i ∈ Z. Theorem 3.10. A basis for the space Sk,Ξ is formed by the functions Bik |

[a,b]

for i = 1 − k, 2 − k, . . . , N.

Note that the number of functions in the basis above is N + k, which of course is also the dimension of Sk,Ξ . We also note that the B-splines Bik , for i = 1 − k, 2 − k, . . . , N , are precisely those B-splines whose support has some overlap with the interval (a, b). A consequence of Theorem 3.10 is that every function s ∈ Sk,Ξ can be written in the form N X s(x) = cj Bjk (x), x ∈ [a, b], j=1−k

for some scalars (numbers) {cj }. This form is known as the bb-form of s, where bb is meant to connotate the phrase B-spline basis. Since the B-splines in use are determined by the knots {ξ1−k , ξ2−k , . . . , ξN +k+1 }, the bb-form of s is determined simply by these knots along with the coefficients {cj }, j = 1 − k, 2 − k, . . . , N . As illustrated by the following example, given the bb-form of s one can use Corollary 3.9 to obtain the bb-form of the derivative of s or of an anti-derivative of s. Example. Let ξi = i for all i, and define 2 s(x) := 2B−1 (x) − B02 (x) + 3B12 (x) + B22 (x) − 3B32 (x) − B52 (x) + B62 (x), x ∈ [1, 7]. Rx Find the bb-forms for the derivative s0 (x) and the antiderivative se(x) = −1 s(t) dt, where x ∈ [1, 7].

6

ADVANCED NUMERICAL COMPUTING

4. B-Splines in Octave As usual we assume that we are given the knots Ξ = [a = ξ1 , ξ2 , . . . , ξN +1 = b]. In order to use the B-spline basis for Sk,Ξ we need k additional knots preceeding a and k additional knots following b. These additional knots must be chosen to ensure that the extended knot sequence Ξext = [ξ1−k , ξ2−k , . . . , ξN +k+1 ] is increasing. One suggestion, which is motivated by a desire to maintain numerical stability, is to choose the additional knots to be equispaced with spacing h = (b − a)/N . In Octave, the first index of a vector is always 1, so the vector Ξ ext is stored as [Ξext (1), Ξext (2), . . . , Ξext (N + 2k + 1)]. Thus one has to keep in mind the correspondence ξi = Ξext (i + k) for i = 1 − k, 2 − k, . . . , N + k + 1, or equivalently Ξext (i) = ξi−k for i = 1, 2, . . . , N + 2k + 1. By Theorem 3.10, the restriction to [a, b] of the B-splines Bik , i = 1 − k, 2 − k, . . . , N , form a basis for Sk,Ξ . In order to construct the B-spline Bik in Octave, one must first form the vector x which contains the knots of Bik , namely, x = [ξi , ξi+1 , . . . , ξi+k+1 ]. This can be accomplished with the Octave command x=Xi ext(i+k : i+2*k+1); k The B-spline Bi can then be constructed using the supplementary Octave function B spline. The command C=B spline(x); will produce the (k + 1) × (k + 1) matrix C so that the pair (x,C) represents Bik | [ξi ,ξi+k+1 ]

as a pp in Sk,x (keep in mind that

Bik

= 0 outside [ξi , ξi+k+1 ]).

Warning!. Mathematically speaking, a pp s ∈ Pk,Ξ (as defined in Section 2) is a function defined on the interval [a, b]. However, in Octave one is allowed to evaluate (using pp val) a pp s at any point x ∈ R: If x < a, then the first piece is used (ie. s(x) := s1 (x)) while if x > b, then the last piece is used (ie. s(x) := sN (x)). In other words, the first piece of s is extended all the way down to −∞ and the last piece is extended all the way up to ∞. This extension is wisely chosen and usually very convenient, except for the case of B-splines. The problem is that the B-spline Bik should be 0 outside the interval [ξi , ξi+k+1 ], but the first and last pieces produced by the Octave command C=B spline(x); are non-zero polynomials which extend, to the left and right, as non-zero polynomials. We get around this difficulty with the command [X,C1]=B spline(x); which returns the pp (X,C1) which has an extra 0 piece at the beginning and end. Specifically, C1 is obtained from C by adding a row of zeros to the top and bottom, while X is obtained from x by adjoining an extra knot at the beginning and at the end. In practice, one usually first constructs the bb-form of a spline s ∈ Sk,Ξ (we discuss several methods in the next section). In order to use this spline s, it is desirable to have s

CS-543

7

in its pp form. The supplemental Octave command S=bb2pp(Xi ext,d); returns the coefficient matrix S of the pp (Xi ext,S) which corresponds to s(x) =

N X

cj Bjk (x),

j=1−k

where d = [c1−k , c2−k , . . . , cN ]. If one prefers the restriction of this pp to the interval [a, b], then one should use instead the command [Xi,S]=bb2pp(Xi ext,d,k); 5. Spline Interpolation Let us consider the following interpolation problem: Given points (ξi , yi ), i = 1, 2, . . . , N + 1, we seek a nice function s which satisfies (5.1)

s(ξi ) = yi for i = 1, 2, . . . , N + 1.

The function s is said to interpolate the given data. We assume that the nodes Ξ := {a = ξ1 , ξ2 , . . . , ξN +1 = b} are increasing (ie. ξi < ξi+1 ), and we desire to choose s as an element of the spline space Sk,Ξ , where k ≥ 1 is chosen to reflect the desired degree of smoothness in s. The case k = 1, known as linear spline interplation, is quite easy: The i-th piece of −yi s is simply si (x) = yξi+1 (x − ξi ) + yi . i+1 −ξi Example. Find the linear spline s which passes through the points (1, 2), (3, 6), (4, 4), (7, 1), and write s in both pp-form and bb-form. In addition to the nodes in Ξ, we need nodes ξi for i = 1 − k, 2 − k, . . . , 0 and i = N + 2, N + 3, . . . , N + 1 + k chosen so that ξi < ξi+1 for i = 1 − k, 2 − k, . . . , N + k. For the sake of brevity, we will write the B-spline Bjk simply as Bj . By Theorem 3.10, the restriction to [a, b] of the B-splines Bj , j = 1 − k, 2 − k, . . . , N , form a basis for Sk,Ξ . Hence, we can write s in the bb-form (5.2)

s(x) =

N X

cj Bj (x),

x ∈ [a, b],

j=1−k

where the coefficients {ci } are unknown at present. The condition s(ξi ) = yi becomes s(ξi ) =

N X

cj Bj (ξi ) = yi ,

1 ≤ i ≤ N + 1.

j=1−k

These equations can be written as the matrix equation  B B2−k (ξ1 ) ··· BN (ξ1 )   c1−k   y1  1−k (ξ1 ) B2−k (ξ2 ) ··· BN (ξ2 )   c2−k   y2   B1−k (ξ2 )  .  =  . .  (5.3) . . . ..  .   .   .. .. .. . . . yN +1 cN B1−k (ξN +1 ) B2−k (ξN +1 ) · · · BN (ξN +1 )

8

ADVANCED NUMERICAL COMPUTING

It is very important to note that since Bj (x) = 0 for all x 6∈ (ξj , ξj+k+1 ), the above (N + 1) × (N + k)-matrix has exactly k(N + 1) non-zero entries. Indeed, the i-th row reduces to [0, 0, . . . , 0, Bi−k (ξi ), Bi−k (ξi ), . . . , Bi−1 (ξi ), 0, 0, . . . , 0] which has only k non-zero entries. When k is greater than 1, the dimension of Sk,Ξ (or equivalently, the number of unknowns c1−k , c2−k , . . . , cN ) exceeds the number of conditions in (5.1), and it turns out that there are infinitely many splines s ∈ Sk,Ξ which satisfy the interpolation conditions (5.1). In order to arrive at a unique spline s ∈ Sk,Ξ it is necessary to impose k − 1 additional conditions. In the literature, there are four prominent ways of imposing these additional conditions; these lead to the natural spline, the complete spline, the not-a-knot spline, and the periodic spline. The Natural Spline Theorem 5.4. If k > 1 is an odd integer, then there exists a unique spline s ∈ S k,Ξ which satisfies the interpolation conditions (5.1) along with the additional end conditions k+1 k+3 , , . . . , k − 1. s(`) (a) = 0 = s(`) (b) for ` = 2 2 With s written in the form (5.2), these additional end conditions become N X

j=1−k

(`) cj Bj (a)

=0=

N X

(`)

cj Bj (b).

j=1−k

Adjoining these extra equations to those of (5.3) and setting m = (k + 1)/2, we arrive at the linear system   (m) (m) (m)   B1−k (a) B2−k (a) ··· BN (a) 0   .. .. .. ..  ..    . . . .  .    (k−1) (k−1)    B (k−1) (a)  B2−k (a) · · · BN (a)    0   1−k      c B2−k (ξ1 ) ··· BN (ξ1 )  1−k  y1   B1−k (ξ1 ) c    2−k   .  .. .. .. ..   .  =  . . . . . .   .   .  .      yN +1   B1−k (ξN +1 ) B2−k (ξN +1 ) · · · BN (ξN +1 )  c     N (m) (m) (m)  0   B1−k (b) B2−k (b) ··· BN (b)   .       ..  .. .. .. ..   . . . . 0 (k−1) (k−1) (k−1) B1−k (b) B2−k (b) · · · BN (b)

Let us refer to the above (N + k) × (N + k) matrix as A, and note that those additional equations associated with left end-conditions are placed at the top of A and those associated with right end-conditions have been placed at the bottom of A. It follows from Theorems (`) 5.4 and 3.10 that A is non-singular. Since Bj (x) = 0 = Bj (x) for all x 6∈ (ξj , ξj+k+1 ), it turns out that A is a banded matrix, and the above system can be solved very efficiently using Doolittle’s LU decomposition for banded matrices.

CS-543

9

Example 5.5. Find the natural cubic spline s which passes through the points (1, 2), (3, 6), (4, 4), (7, 1). The natural splines are famous for the following property: Theorem 5.6. Let s be the natural spline as described in Theorem 5.4. If f ∈ C k−1 [a, b] interpolates the given data: f (ξi ) = yi , then

Z

b a

Z (m) 2 s (t) dt ≤

b a

i = 1, 2, . . . , N + 1,

(m) 2 f (t) dt,

where m = (k + 1)/2.

The Complete Spline Theorem 5.7. If k > 1 is an odd integer, then there exists a unique spline s ∈ S k,Ξ which satisfies the interpolation conditions (5.1) along with the additional end conditions s(`) (a) = ya,` s(`) (b) = yb,`

` = 1, 2, . . . ,

k−1 . 2

With s written in the form (5.2), these additional end conditions become N X

j=1−k

(`) cj Bj (a)

= ya,` and

N X

(`)

cj Bj (b) = yb,` .

j=1−k

Adjoining these extra equations to those of (5.3) and setting m = (k + 1)/2, we arrive at the linear system     (1) (1) (1) B1−k (a) B2−k (a) ··· BN (a) ya,1   .. .. .. ..  ..     . . . . .      (m−1) (m−1)   B (m−1) (a)  B2−k (a) · · · BN (a)   ya,m−1   1−k     c   y1  B2−k (ξ1 ) ··· BN (ξ1 )  1−k  B1−k (ξ1 )    c2−k    ..  . .. .. .. ..  .  =    .  . . . .  .     .    yN +1   B1−k (ξN +1 ) B2−k (ξN +1 ) · · · BN (ξN +1 )     cN  (1) (1) (1)  yb,1    B1−k (b) B (b) · · · B (b)   N 2−k     .   .. .. .. .. ..     . . . . (m−1) (m−1) (m−1) yb,m−1 (b) B2−k (b) · · · BN B1−k (b) As with the linear system obtained for the natural spline, the above (N + k) × (N + k) matrix is nonsingular and banded, and consequently, the above system can be solved very efficiently using Doolittle’s LU decomposition for banded matrices.

10

ADVANCED NUMERICAL COMPUTING

Example 5.8. Find the complete cubic spline s which passes through the points (1, 2), (3, 6), (4, 4), (7, 1) and additionally satisfies s0 (1) = 0 and s0 (7) = −2. The complete splines are famous for the following property: Theorem 5.9. Let s be the complete spline as described in Theorem 5.7. If f ∈ C k−1 [a, b] interpolates the given data: f (ξi ) = yi ,

i = 1, 2, . . . , N + 1,

and also satisfies the end conditions f (`) (a) = ya,` f (`) (b) = yb,` then

Z

b a

Z (m) 2 s (t) dt ≤

b a

` = 1, 2, . . . ,

(m) 2 f (t) dt

k−1 , 2

where m = (k + 1)/2.

The complete splines have another property which can be used to find the bb-form of a spline given in pp form. Theorem 5.10. Let f ∈ Sk,Ξ and let s be the complete spline (as in Theorem 5.7) determined by the interpolation conditions s(ξi ) = f (ξi )

i = 1, 2, . . . , N + 1,

along with the end conditions s(`) (a) = f (`) (a) s(`) (b) = f (`) (b) Then s = f .

` = 1, 2, . . . ,

k−1 . 2