LAB #8: Piecewise Polynomial Interpolation


TABLE OF CONTENTS

Introduction

Exercise 1

MATLAB Tips

Exercise 2

Piecewise Parabolic Interpolation

Exercise 3

Cubic Hermite Interpolation

Exercise 4

Cubic Spline Interpolation

Exercise 5

Parametric Interpolation

Exercise 6

Introduction

Piecewise linear interpolation has many good properties. In particular, if there was a continuously differentiable function f(x) generating the data, the data points were suitably spread throughout the closed interval, then the interpolant converged to the function.

Now we'd like to consider some common attempts to improve the interpolant, increasing the rate of convergence, or the smoothness. Unfortunately, the formulas for these interpolants get complicated, so I'll give you most of the information you need.

We will also take a quick look at one way that interpolation can move into higher dimensions.

MATLAB Tips

We're using the Runge function to examine the performance of our interpolation, which is fine. It means it's easy to increase the number of data points, get the derivative, and so on. However, the Runge function is pretty smooth, and easy to interpolate. Other functions can be used to give the interpolation routines a few headaches. The following functions are useful for stressing interpolation routines and also to illustrate problems that might arise. We will not be using them explicitly, but you can try the interpolation routines we generate on these functions just for fun.

A very interesting function is 1 at a single place and 0 elsewhere. We'll call this the tack function because its graph looks like a thumb tack. Here's how you could define it in Matlab:

function pval = tack ( x )
n = length ( x );
pval = ( [1:n] == floor(n/2) );

I need to tell you that floor(x) returns the integer part of a real value, and that in Matlab, a logical expression has the value 1 if true, and 0 if false.

Another function is 0 for a while and then 1. We'll call this the step function because it looks like a step and because most other people call it that. Here's how you could define it in Matlab:

function pval = step ( x )
n = length ( x );
pval = floor ( 2 * [0:n-1] / n );

The comb function is 0 and 1 and 0 and 1 and ... . You will meet this function often if you solve differential equations.

function pval = comb ( x )
n = length ( x );
pval = mod ( [1:n], 2 );

If we increase the divisor 2 to a larger value, we might call the result a rasp function.

These functions are artificial; they don't really correspond to a formula that applies for all x. But they do tend to show up the differences and defects of an interpolation method.

In the exercises below, I've asked you to look at the Runge function over and over again. I imagine that must get boring after a while. When you want to do something more interesting, try using one of these artificial functions instead. The tack function is interesting, because it shows how changing one data value affects the entire code.

Piecewise Parabolic Interpolation

The piecewise linear interpolant converges, but it doesn't look very nice. Can we try a similar interpolation using pieces of parabolas, a quadratic interpolant? If we do the right things, we'd expect the convergence rate to increase, and the interpolant to be a little more graceful over the interpolation intervals. We'll still have those unpleasant "kinks" in the derivatives at the interior abscissas, but we'll see if we can deal with them later.

To define a parabolic function, we need three data points. So each piece of the piecewise function will actually be defined over two consecutive data intervals. Keep this in mind.

The following recipe explains how to evaluate the piecewise quadratic interpolant at a point xval. You will use this recipe to help write a M-file in the exercise below.

The first step is to determine the pair of intervals that contain the position xval. We can use the function bracket that we wrote before, but it returns a single interval containing xval. In order to pick out the consectutive pair of intervals containing xval, simply note that each consectutive pair of intervals must start with an odd-numbered point and have an even-numbered point at its center. Thus, if bracket returns an even number, subtract one to get the odd-numbered point that identifies the pair of intervals. You can use the Matlab functions floor (see above), mod or rem to help you establish when a number is even. Think carefully about using these as part of a vector statement.

The next piece of the puzzle is, if we have three data points (x0,y0), (x1,y1), and (x2,y2), what is the value at x of the parabola that goes through these points? The following formula uses divided differences to find and evaluate the parabola.

        c = y0
        b = ( y1 - y0 ) / ( x1 - x0 )
        a = ( ( y2 - y1 ) / ( x2 - x1 ) - ( y1 - y0 ) / ( x1 - x0 ) ) / ( x2 - x0 )
        p = c + b * ( x - x0 ) + a ( x - x0 ) * ( x - x1 )

You'll need to convert these short names to the ones we're using in the program, and you will want to use componentwise multiplications.

Exercise 1: write a Matlab M file called paraval.m which has a declaration of the form:

function pval = paraval ( xdata, ydata, xval )

and which evaluates the piecewise parabolic interpolant at the vector of values xval. Here is a sketch of the routine:

        left = bracket ( xdata, xval );
        left = ???

        c = ???
        b = ???
        a = ???

        pval = ???

This routine should work for a vector of values if you are careful to include the proper operators. If you cannot figure out how to make it work using vector, then use a loop. Here is a sketch of the loop version:

        for i=(1:length(xval))
          left = bracket ( xdata, xval(i) );
          left = ???

          c = ???
          b = ???
          a = ???

          pval(i) = ???
        end

Test your work by choosing 5 points for xdata and using a quadratic polynomial of your choice to generate 5 ydata values. Because you are using quadratic interpolation of a quadratic polynomial, your interpolated results (from paraval) should agree exactly with values from the polynomial you chose. Pick at least three new points for xval and check the agreement. Why is three the right number of points? Be sure to include paraval.m when you send your lab report to me.

Exercise 2: Use 5 equally spaced data points from -5 to 5 for xdata and the Runge function to compute ydata. Estimate the maximum interpolation error by computing the maximum absolute difference at 201 evenly spaced points. Repeat for 11, 21, 41 and 81 data points. As the interval size decreases by a factor of 2, how is the maximum error behaving? I suggest you modify the poly_interpolate.m file you wrote for the last lab to do this exercise.

Cubic Hermite Interpolation

When we look at the results of the piecewise parabolic interpolation, there are a couple of things that might concern you. First, there's the tricky business that the number of data points must be odd, so that the number of intervals is even, so that we can take pairs of intervals to build the parabolas. Second, although the function is nice and smooth within each pair of intervals, there's an obvious "kink" in the derivative at the odd numbered data points, where the form of the parabola changes.

If we were trying to design, say, the shape of the sheet metal pattern for a car door, this kind of kink would not be acceptable. If that's a problem, then perhaps we could try to make sure that our interpolant is smoother and simpler by matching both the function value and the derivative at each point. (This assumes you can get the derivative value. That's not always the case!)

Suppose we have ndata points xdata, at each of which we have both the function value ydata and ypdata. If there really is a function f(x) lurking behind the data, then we may assume, for i from 1 to ndata, that:

xdatai = f ( xdatai )
ypdatai = f' ( xdatai )

Now consider the situation in the i-th interval, [ xdatai, xdatai+1 ]. There is just one cubic function that can be defined on this interval, and which has the appropriate function and derivative values at the two endpoints. As before, we need to determine the interval, (but no extra funny business like with paraval), and then evaluate the function defined by the data. The algebra is a little tricky, so I'll supply that for you.

Exercise 3: create a Matlab M file rungep.m which computes the derivative of the Runge function:

function fp = rungep ( x )

I leave it to you to correctly differentiate the Runge function! Then copy from my web page the file herval.m. We are going to use this file to approximate the Runge function.

The following plan for evaluating the error in Hermite interpolation is very similar to that in the poly_interpolate utility function you wrote for the previous lab, and you will probably want to modify that utility function to perform

Of course, we really want to know if the approximation gets better, or goes crazy the way that the continuous polynomial interpolant did. So repeat this exercise for ndata = 11, 21, 41 and 81.

Cubic Spline Interpolation

One disadvantage of the Hermite interpolation scheme is that you need to know the derivatives of your function. But it's very possible that you don't have any formula for your data, just the values at the data points. A second problem is the Hermite interpolant is smooth, but not smooth enough. The discontinuities in the second derivative at the data points can be noticeable.

Cubic splines are piecewise cubic interpolants which are "very smooth". They are continuous, with continuous first and second derivatives. Only the third derivative is allowed to jump at the "join" points.

The continuity and interpolation conditions are almost enough to specify the cubic. If we add one extra condition at each endpoint, such as the value of the derivative, then the spline is determined.

To use a cubic spline requires several steps. For a given set of data, you compute the spline coefficients. Then, to evaluate the spline at a point, you have to bracket the point with data abscissas (just as we did for the other piecewise functions), and then evaluate the appropriate polynomial.

Exercise 4: From my web page, copy the files cubset.m and cubval.m. These two routines are used in pairs to first compute the coefficients and then to use them. Modify poly_interpolate.m to use these two functions in the following way:

          xtest = linspace ( -5, 5, 201 )
          ypp = cubset ( xdata, ydata )
          ytest = cubval ( xdata, ydata, ypp, xtest )
        

Start with 5 equally spaced data points from -5 to 5 on the Runge function. Estimate the maximum interpolation error in a manner similar to the previous exercise. Repeat for 11, 21, 41 and 81 data points. As the interval size decreases by a factor of 2, how is the maximum error behaving?

Parametric Interpolation

There are many aspects of interpolation we can't cover, especially what happens in multiple dimensions. But let's at least look at a fun aspect, for a change. We're familiar with a parameterized curve, in which a special variable, typically called t, is used to draw objects for which the relationship between x and y is not functional. We are comfortable with the idea that a good definition of a circle is:

          x = cos ( 2 * pi * t )
          y = sin ( 2 * pi * t )
        

Well, suppose we wanted to do interpolation of a complicated curve? If we have a drawing of the curve, we could simply mark points along the curve that are roughly evenly spaced, and make tables of tdata, xdata and ydata, where the values of t are just 1, 2, 3,...

Now if we want to compute intermediate points on the curve, that's really saying, for a given value of t, what would the corresponding values of x and y be? This takes a bit of thought, but just realize that t is the independent variable that we're used to calling x, and both x and y become dependent variables. If you can keep track of that in your mind, you can see how to modify our interpolation methods to apply them to the parameterized case.

Exercise 5. Try the following method for plotting a cycloid.

          ndata = 21;
          tdata = linspace ( 0, 2*pi, ndata );
          r = 0.5 + cos ( tdata );
          xdata = r .* cos ( tdata );
          ydata = r .* sin ( tdata );

          t = linspace ( 0, 2*pi, 5*ndata );
          x = lintval ( tdata, xdata, t ); 
          y = lintval ( tdata, ydata, t );

          plot ( x, y )
         

Try it again with paraval to see a smoother plot.

Exercise 6: on a sheet of graph paper, choose a box of about 10 by 10 cells, and write a large script capital letter that roughly fills the box. You'll want to pick a letter that can be drawn in a single curve. You might use one of your initials, or fancy script letters "E", "M" or "P". Mark 21 (roughly) equally spaced points along the curve, and record the xdata and ydata locations of these points. Now try the following commands:

          ndata = length ( xdata );
          tdata = [ 1 : ndata ];

          t = linspace ( 1, ndata, 3*(ndata-1) + 1 );
          x = lintval ( tdata, xdata, t ); 
          y = lintval ( tdata, ydata, t );

          plot ( x, y )
        

Use a "better" interpolation routine to produce a smoother picture.


Back to the MATH2070 page.

Last revised on October 25, 2000.