Aaa

 4.10 TENSION SPLINES

Splines are a wonderful tool for approximation, but they can still exhibit some poor behavior. Consider the data set plotted in Fig. 4.25. Obviously, this represents a function with a severe jump near x — 0.5, but there is no sign of oscillatory behavior. However, a B-spline representation of this data (Fig. 4.26) shows small wiggles on either side of a sharp front. This is fundamentally an artifact of the steep gradient in the data, but in other contexts a spline fit can display behavior that does not match the "sense" of the data. One way to avoid the problem is the notion of a taut spline or tension spline, an idea that appears to have been first published by Schweikert [17], but which also owes a lot to the work of A. K. Cline [4]; we relied heavily on a short paper of Marusic and Rogina [12] in our presentation here.

Imagine that the curve in Fig 4.26 is a piece of string that is constrained to pass through small loops at the data points. If we were to pull the string taut, we would smooth out the spurious oscillations in the curve. This amounts to studying the mechanical properties of a cable hanging between two supports. More prosaically, we construct our spline from the new basis set {1, a;, coshpx, sinhpa;}, where p > 0 is the tension parameter: p = 0 means no tension, and it can be shown that this corresponds to the pure spline approximation; p —» oo gives us a piecewise linear approximation. (We will not attempt to justify either of these statements other than by examples and exercises.)

The reader may well be wondering how practical this scheme might be. After all, we have traded a set of polynomial basis functions for a set of transcendental basis functions. Not only is this going to make execution of any program more expensive, but it leaves open the entire question of how to construct the approximation. The basic idea is the same as in §4.8: First, we construct a primary basis function, as in (4.34):

  (4.56)

where   Confirmation of this formula is deferred to the exercises. A plot of    is given in Fig. 4.27, for p = 4; note that it does not look very different from Fig. 4.13.

Construction of a tension spline follows precisely the same recipe as in §4.8: Given data   we look for an approximation in the form

 

Where the   values are coefficients to be determined, and

 

Note that, as in §4.8, we have added additional grid points, which will again require the imposition of additional conditions. We will explicitly cover the natural spline case, leaving the complete spline case to the exercises. For each original grid point we have the equation

 

Because of the local nature of the   functions, this becomes

 

We quickly have that

 

and

 

for 

 

This certainly looks imposing, but it also looks like some simplification ought to be possible, and if we multiply out the numerator, we quickly get

 

The important thing for our purposes is that   The system of linear equations for the taut spline is, initially, as follows:

This, of course, is a system of   equations in   unknowns. Following the derivation for the natural spline in §4.8, we have the boundary conditions. 

 

and 

 

Continuing with  , we get that 

 

We still need to compute ô" at various points. It is an involved, but not arduous computation, to show that 

 

and 

 

We have plotted   in Fig. 4.29. we therefore have, from (4.58)

 

Similarly,

 

Note that this looks a lot like what we got in §4.8 for the ordinary natural B-spline; we can eliminate   and  , so the system (4.57) becomes (after a bit of work)

 

Where 

 

 

And 

 

With 

 

 

 

and 

 

Note that the matrix is tridiagonal and diagonally dominant, so we know how to solve it.

EXAMPLE

As an obvious illustration of this, let's fit a series of tension splines to the data set in Fig. 4.25. Once the coefficients   and  are computed, the spline is evaluated in the same way as for an ordinary B-spline. Fig. 4.30 shows the data and spline curve for   = 0.01; obviously, with such a small tension value we do not expect much difference, and Fig. 4.31, which plots the difference between the pure polynomial spline and the tension spline, confirms this (although the lower "overshoot" does appear to be significantly affected). If we take   = 1, we get Figs. 4.32-4.33 which shows that the lower "overshoot" is indeed beginning to damp out, but the upper one is much the same. Finally, for   = 6, we get Figs. 4.34-4.35; this is perhaps the smallest value of p for which both overshoots are gone. 

4.11 LEAST SQUARES CONCEPTS IN APPROXIMATION

4.11.1 An Introduction to Data Fitting

An important area in approximation is the problem of fitting a curve to experimental data. Since the data is experimental, we must assume that it is polluted with some degree of error, most commonly measurement error ("noise"), so we do not necessarily want to construct a curve that goes through every data point. (In fact, the material in §4.12.1 suggests that this would be a disastrous way to proceed.) Rather, we want to construct a function that represents the "sense of the data" and which is, in some sense, a close approximation to the data. 

The most common approach is known as least squares data fitting. Consider Figure 4.36; this shows an example set of data for which the general trend is clearly a straight line. But which straight line do we use to represent the data?

The least squares approach defines the "correct" straight line as the line that minimizes the sum of the squares of the distances between the data points and the line. Let the experimental data be defined as pairs  for some n. Thus, we want to find the coefficients m and b in the equation  such that

 

is minimized. This is a straight forward problem from multivariable calculus: We compute the partial derivatives   and   and find where they both vanish, and this will define a critical point. It can be shown that this critical point defines a global minimum for F. Thus, m and b are defined by the two equations

 

 

which can be simplified to a system of two equations in two unknowns

 

 

or

 

 

The solution here is then 

 

 

The notion of a least squares data fit can be generalized beyond simply fitting a straight line to data. We can look at higher degree polynomials and we can also look at higher dimensional data sets. The exercises include some examples of more involved least squares data fit problems.

EXAMPLE

Consider the data in Table 4.24. If we plot this, we get what appears to be a straight line as the general trend of the data, so we look for the equation of the line   which best fits this data in the least squares sense. Forming the separate sums gives us

 

From which it follows that 

 

And the line is plotted, along with the data, in Figure 4.38

Table 4.24 Data for example

x 0.0 1.0 2.00 3.00 4.00 5.00

y 10.0 25.0 51.0 66.0 97.0 118

FIGURE 4.38

EXAMPLE

We don't have to restrict ourselves to linear or quadratic models to make good use of the idea of least squares data fits. Consider the   data in Table 4.25. When we plot this data, we get the curve shown in Figure 4.39. Generally, this looks like an exponential growth curve. Ordinarily this would require us to do a fit to a curve of the form , which will lead to a nonlinear system for the parameters   and  . However, if the raw data is exponential, then the logarithm of the data is linear, since we have

 

and we can fall back on our existing algorithm to do a fit to the log data. 

To verify this, we look at the logarithm of our example data; this data is plotted in Figure 4.40, and the general trend is indeed linear. So we do a least squares fit to the log data, getting the straight line

 

This is plotted, along with the raw (logarithm) data, in Figure 4.41. It follows, then, that our curve fit to the original data is

 

which is plotted in Figure 4.42, along with the original data. Note that, except for the last two points, this is a pretty good fit to the data set.

4.11.2 Least Squares Approximation and Orthogonal Polynomials

The notion of least squares approximation can be extended beyond the data-fitting problem. Consider the problem of finding an approximation to a given function f in terms of a set of basic functions  . How do we find the coefficients in the expansion.

 

One way to do this is to require that the coefficients   produce an approximation that minimizes the error. 

 

In the least square sense, i.e., in sthe sense of the integral 2-norm. Tus we seek   such that

 

is minimized. Before doing this it will be convenient to introduce the notion of an inner product. It is best to do this thoroughly, so we pause in our development of approximations, but only briefly.

Inner Products of Functions The reader should be familiar with the notion of the dot product of two vectors in  :

 

This is an example of a more general operation called an inner product, which can be defined on general vector spaces, including spaces of functions, rather than just Euclidean n-space. The formal definition is as follows.

Definition 4.2 (Inner Product on Real Vector Spaces) Let f and g  be elements of a real vector space V. Let (f,g) denote any operation on f and g that satisfies the following three properties:

1.  

2.  

3.  

Then,   is called an inner product.

If we want to consider our overlying vector space to be   that is, continuous functions on a closed interval, then we can easily establish that the positively weighted integral of a product of two functions will be an inner product.

Theorem 4.8 Let   be integrable on [a,b] and non-negative, i.e.,  is defined and   For given f and g in  , define  as 

 

Then,   defines an inner product on  .

Proof: see the problem 6

Note that it therefore follows that if   is an inner product, then   defines a norm. in the common case when   we will simply write  ;

i.e., we will drop the subscript, and the norm is the ordinary 2-norm for functions that we defined earlier in this chapter. Problem 7 offers some practice with a norm defined by a weighted inner product.

This definition of inner product will allow us to apply a number of ideas from linear algebra to the construction of approximations, as we will soon see. Of more immediate interest is the fact that we can use the inner product notation to write the residual,  , in a very convenient form: 

 

Note that we can regard   as a function of the n variables  , and thus apply ordinary calculus to the problem of minimizing  . After some manipulations (see Problem 8) we find that the   are defined by simultaneously solving the set of equation

 

 

 

 

This system can be organized along matrix-vector lines as 

Solving a system of linear equations is a problem that we do not encounter in the general case until Chapter 7. We can avoid it altogether at this point if our basis functions satisfy the orthogonality condition

  (4.62)

In this case, the matrix in (4.61) is a diagonal matrix and we very easily have

 

So, to summarize what we have done so far, we can construct an approximation to a given function f from a given basis set  , and the construction is very easy, if' the basis satisfies the condition (4.62). So, the question becomes: When can we find a basis that satisfies (4.62), and how good is the resulting approximation?

The answer is that we can always find such a basis if we consider polynomial functions for our basis elements, and the resulting approximations are usually quite good. The special basis functions that satisfy (4.62) are called orthogonal polynomials. To be more specific with this, we have to introduce some new concepts and notation, and recall a major theorem from linear algebra. But first, one more definition.

Definition 4.3 (Vector Space of Polynomials of Degree   For any N ≥ 0 define   as the vector space of polynomials of degree ≤  N. Note that this space has a standard basis consisting of  , and thus is an (N + l)-dimensional space.

And, now, the theorem

Theorem 4.9 Let   be a given non-negative weight function on an interval  , and   the associated inner product, defined as in (4.60). Then there exists a family of orthogonal polynomials  , such that

 

and 

 

In addtition, the   satisfy the following:

1. The set  

2. If    is an arbitrary element of   for k < N, then   (thus, orthogonal polynomials are orthogonal to all polynomials of strictly lower degree);

3. For j ≥ 1, the roots of each   are all in [a, b] and are all distinct.

Proof: The proof is somewhat lengthy, in part because of the length of the theorem, but it is not difficult.

To establish that the family   exists, we will construct it directly, using the Gram- Schmidt process from linear algebra. Take  , and define   for  . Then the subsequent   can be found according to

 

and an inductive argument shows very quickly that the orthogonality holds. In Problem 9 we ask the student to fill in the details of this part of the proof.

Having now proved that the family of orthogonal polynomials exists, we turn our attention to proving each of (l)-(3). 

(1) The space   is finite-dimensional with dimension N + 1. It therefore follows that any set of N + 1 independent elements of  will be a basis. The orthogonality condition (4.62) forces the members of the family   to be independent (see Problem 10); therefore, the set   is a basis for  .

(2) Let   be an arbitrary polynomial of degree  . Then we can write

 

because   is a basis for  . Therefore,

 

since   is orthogonal to each element of  . In fact, we can write (Problem 11)

 

(3) First, suppose that  , has no roots in  . This means that   does not changes signs on the interval, thus 

 

since the integrand does not change sign on   But the orthogonality requires that  ; hence, we have a contradiction, so there must be at least one root in  

Now, let   be any root of   that lies in   and suppose that it is a multiple root. Then it follows that

 

For some polynomial q, thus,

 

Is a polynomial of degree   therefore,  , by another part of this theorem. But

 

so we have another contradiction. Thus, any roots that lie in   must be simple roots.

Suppose now that only some of the roots lie in   Call these roots  , and note that we can write   as

 

where   does not change sign in   and is plynomial of degree  . Therefore,

 

 

is also a polynomial that does not change sign in  . Hence, the integral

 

cannot be zero. However,

 

where   and   is polynomial of degree  . Therefore,  , and we have a contradiction. Thus,   , and since a polynomial of degree   cannot have more than   roots, we must have  .

Families of Orthogonal Polynomials At this point it might be useful to look at some examples of orthogonal polynomial families. Four of the most common ones are discussed below.

Note that the orthogonality condition (4.62) means that an orthogonal polynomial can be multiplied by an arbitrary nonzero constant and still satisfy (4.62). To avoid the problems of nonuniqueness that this can lead to, it is common to impose a specific scaling on the elements of each family. 

1. Legendre polynomials: The Legendre polynomials are the orthogonal polynomials on with no weight function (more correctly, the unit weight function); thus, we have

 

The usual scaling is to take  . The first five Legendre polynomials are:

 

 

 

 

 

2. Chebyshev polynomials: The common notation for the Chebyshev polynomials is , and the interval and weight function are defined in the orthogonality relation

 

The common scaling is to set the leading coefficient equal to . The first five Chebyshev polynomials are:

 

 

 

 

 

It can be shown that the Chebyshev polynomials are related in a very simple way to cosines; see Theorem 4.10.

3. Hermite polynomials: The Hermite polynomials are orthogonal on the entire real line, using the weight function   ; that is,

 

The common scaling is to set the leading coefficient equal to  . The first five Hermite polynomials are:

 

 

 

 

 

4. Laguerre Polynomials: The Laguerre polynomials are orthogonal on the positive real line, using the weight function w  that is,

 

The common scaling is to set the leading coefficient equal to  . The first five Laguerre polynomials are then:

 

 

 

 

 

We can use any of these orthogonal polynomial families to construct approximations to functions defined on the appropriate interval. These approximations are "best possible" in the sense that they minimize the error in the appropriate weighted 2-norm; i.e.,

 

for all   

Consider, as illustrations, the following set of examples.

EXAMPLE 4.14

Let's construct the fourth-degree least squares approximation to the exponential function,  , over the interval   using Legendre polynomials. The approximation is defined by

 

where the   the Legendre polynomials. We thus need to compute the integrals 

 

 

and

 

 

In practice, these integrals would be computed using some type of numerical integration routine, such as the trapezoid rule or the more accurate methods we discuss in Chapter 5. For this simple example, though, it is possible to use direct calculus methods or (much more attractive!) a computer algebra package such as Maple or Mathematica. However it is done, to eight digits the integrals are

 

 

and 

 

 

so the polynomial approximation is

 

 

 which can simplifies to

 

Figure 4.43 shows a plot of the error  . Compare this to the error plots for fourth-degree Taylor approximation and fourth degree Lagrange or Newton interpolation from earlier in this chapter. Note, in particular, that the least squares error oscillates back and forth between its maximum and minimum values (or nearly so), several times. It can be shown that this is a necessary and sufficient condition for the approximating polynomial to be the "best" approximation to the function, and is one reason why least squares approximations are considered valuable: they are close to being the best possible approximations.

Figure 4.43 Error in fourth-degree least squares approximation to the exponential function.

EXAMPLE 4.15

Here we construct a Legendre polynomial approximation to the function   using   , and 32 degree polynomials. Figures 4.44A-C show plots of both / and the least squares approximation; Figure 4.44D shows the error in the 32-degree approximation.

Figure 4.44 Legendre least squares approximation to     error in the   In A and B,   is denoted by the dotted curve. 

EXAMPLE 4.16

This time, we use a Chebyshev polynomial approximation to the same f  as in the previous example, again using  and 32 degree polynomials. Figures 4.45 A – C show plots of both f and the least squares approximation; Figure 4.45D shows the error in the 32-degree approximation. The performance of the Chebyshev and Legendre approximations are very similar, with the Chebyshev being very slightly better. 

Figure 4.45 Chebyshev least squares approximation to  : error in the   = 32 case. In A and B,   is denoted by the dotted curve.

EXAMPLE 4.17

Here we construct a Legendre polynomial approximation to   on the interval   in much the same way as was done for   in Example 2. However, the accuracy here is so much greater that we plot the errors for the  cases in Figure 4.46. 

Figure 4.46 Error in Legendre least square approximation to   A: error for   B: error for  ; C: error for  .

EXAMPLE 4.18

Chebyshev approximation to  ; This is the same as Example 4.17, except we use a Chebyshev expansion instead of a Legendre expansion. 

Figure 4.47 Error in Chebyshev least squares approximation to  

Several comments might be in order here. We first note how much easier it was to obtain a high degree of accuracy for the exponential than it was for  . Second, although it might be difficult to discern in the plots, the results for the Chebyshev approximations were, in each case, very slightly better than for the Legendre approximations. A full discussion of this requires more mathematical machinery than we want to deal with right now, but it is generally true that Chebyshev least squares approximations are superior to those done with any other choice of basis.

Finally, note that to do a least squares approximation, we have to be able to compute the inner products, which are integrals. Thus, we need a tool like the trapezoid rule or, perhaps better, some of the more sophisticated methods to be developed in Chapter 5.