Lecture Notes 1 5

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     x − ξi ξi+k+1 − x k−1 k−1 k Bi (x) := Bi (x) + Bi+1 (x). ξi+k − ξi ξi+k+1 − ξi+1 Note that the B-spline is defined on all of R. 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 > 1, then     d k k k k−1 k−1 Bi (x) = Bi (x) − Bi+1 (x). dx ξ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 (v) Bi (t) dt = 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 ]. The following Corollary extends properties (iii) and (v) to a linear combination of B-splines.

CS-543

Corollary 3.9. Let f (x) =

∞ X

5

ci Bik (x).

i=−∞

(i) If k > 1, then

 ∞  X d ci − ci−1 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

f (t) dt = −∞

∞ X

di Bik+1 (x),

i=−∞

i X 1 cj (ξj+k+1 − ξj ). where di = 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).

Find the bb-forms for the derivative s0 (x) and the antiderivative se(x) =

Rx

−1

s(t) dt.