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.

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.

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.

Determine the two consecutive intervals that contain the position

**xval**.Determine the parabola which passes through the three data points that define the intervals; and,

Evaluate that parabola at

**xval**.

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.

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:

xdata_{i}= f ( xdata_{i})

ypdata_{i}= f' ( xdata_{i})

Now consider the situation in the **i**-th interval, **[
xdata _{i}, xdata_{i+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

**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

define 5 data points

**xdata**equally spaced from -5 to 5;define

**ydata**and**ypdata**, the values of the Runge function and its derivative at the data points;define 201 sampling points

**xtest**equally spaced from -5 to 5;compute the value

**ytest**of the Hermite interpolant at 201 points;compute the value

**yexact**of the Runge function at 201 points.compare the plots of the

**(x,yexact)**and**(x,ytest)**;note the maximum norm of

**ytest - yexact**.

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.

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?

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. *