LAB #6: Polynomial Interpolation

TABLE OF CONTENTS

• Introduction
• MATLAB Tips
• The Polynomial through Given Data
• Vandermonde's Equation
• Lagrange Polynomials
• Computing Divided Differences
• Evaluating Divided Differences
• ASSIGNMENT

Introduction

Polynomials are often used to approximate data. There are several methods for defining, constructing, and evaluating the polynomial of a given degree that passes through data values. We can:

• solve for the coefficients of the polynomial by setting up a linear system called Vandermonde's equation.
• represent the polynomial as the sum of simple Lagrange polynomials.
• represent the polynomial in the form of Newton's divided difference table.

It can be very helpful to use a consistent notation. We will try to use the names xdata and ydata to refer to known data values, while xval will be a arbitrary point at which we want to evaluate the interpolating polynomial, and pval will be the value of the polynomial at that point.

At the end of this lab is an assignment due by the beginning of the next lab session.

MATLAB Tips

At some point in this lab, you will need to determine if a quantity x is not equal to a quantity y. The way to do this is

if ( x ~= y )

If you have a (square) linear system of equations, of the form

A * x = b

x = A \ b

Exercise: Test the backslash command by setting up and solving the linear system

        [ 1  2 ]         [ 17 ]
        [ 3  4 ]  * x =  [ 39 ]
      

We will want to be able to store certain tables of data in a file. We can have MATLAB read and write these tables using the save and load commands. For now, we are going to want to save two vectors of the same length, so MATLAB can treat them as two rows of one array. Supposing we have two vectors called xdata and ydata, here's how we could store them in a file:

save parabola xdata ydata

load parabola

Exercise: create two vectors, called xdata and ydata, with the values:

        XDATA :  0  1  2
        YDATA :  3  6 11
      

We will want to make plots in which the data values show up with a marker of some sort. One way to do this is to specify 'o', which indicates that the points are to be marked with circles.

plot ( xdata, ydata, 'o' )

Exercise: load the data from the parabola file. On one graph, plot the pairs of points from the data file, and the curve:

y = x² + 1

MATLAB has a command

c = polyfit ( xdata, ydata, n )

Exercise: read the data from the parabola file. Use the polyfit command to determine a set of polynomial coefficients that match the data. Now figure out how to plot the data points as circles, and the polynomial determined by polyfit. Does the formula match the data?

The Polynomial through Given Data

Maybe you didn't notice it, but in the secant method and the regula falsi method, we were solving the problem

Determine the polynomial p(x) of degree 1 that passes through the points (a,fa) and (b,fb).

Determine the point x0 at which a polynomial p(x) of degree 1 is equal to zero.

The idea that is common here is that we are trying to construct a "model function", in this case a polynomial, which would "explain" all the data that we have. We may then want to examine the graph of the polynomial, evaluate it at other points, determine its integral or derivative, or do other things with it.

Vandermonde's Equation

Here's one way to see how to organize the computation of the polynomial that passes through a set of data. It's not the best way, but it will give us a good place to start.

Suppose we wanted to determine the linear polynomial

p(x) = c(1) * x + c(2)

        c(1) * xdata1 + c(2) = ydata1
        c(1) * xdata2 + c(2) = ydata2

        c(1) = ( ydata2 - ydata1 ) / ( xdata2 - xdata1 )
        c(2) = ( xdata2 * ydata1 - xdata1 * ydata2 ) / ( xdata2 - xdata1 )

Compare that situation with the case where we want to determine the quadratic polynomial

p(x) = c(1) * x^2 + c(2) * x + c(3)

        c(1) * xdata1^2 + c(2) * xdata1 + c(3) = ydata1
        c(1) * xdata2^2 + c(2) * xdata2 + c(3) = ydata2
        c(1) * xdata3^2 + c(2) * xdata3 + c(3) = ydata3

You should be able to see that, for any collection of abscissas and data points, it is possible to define a linear system that should be satisfied by the (unknown) polynomial coefficients. If we can solve the system, and solve it accurately, then that is one way to determine the interpolating polynomial.

Interesting fact: Vandermonde's coefficient matrix is nonsingular if and only if all the abscissas xdata are distinct. Why would we never want to have two abscissas the same? What are we saying if two abscissas are very close, and the corresponding ydata values are far apart?

If we look at this problem as an abstract linear system, we need to define the individual entries of the matrix A:

      for i = 1 to 3
        for j = 1 to 3
          A(i,j) = xdata(i)^(3-j)
        end
      end 

c = A \ ydata'

Exercise: write a MATLAB M function file, fitpoly.m, of the form

          function c = fitpoly ( xdata, ydata )
          n = ?
          A = ?
          c = ?

          XDATA:  0   1   2
          YDATA:  3   6  11

Lagrange Polynomials

Suppose we fix the set of n distinct abscissas xdata, and think about the problem of constructing a polynomial which has (not-yet-specified) values ydata at these points. Now suppose I had a polynomial l₇(x) which was zero at each abscissa value, except that at xdata₇, where it was 1. Then ydata₇ * l₇(x) would have the value ydata₇ at xdata₇, and if we could just find a similar polynomial for the other abscissas, I'd be done.

In fact, the Lagrange polynomials are easily constructed for any set of abscissas. Each Lagrange polynomial will be of order n (which is degree n-1). There will be n Lagrange polynomials, one per abscissa, and the i-th polynomial will be named l_i(x), and will have the "special relationship" with the abscissa xdata_i, namely, it will be 1 there, and 0 at the other abscissas.

In terms of Lagrange polynomials, then, the interpolating polynomial has the form:

          p(x) = ydata1 * l1(x)
               + ydata2 + l2(x)
               ...
               + ydatan + ln(x)

Exercise: Let's consider the abscissas of our data set:

            I:  1   2   3
        XDATA:  0   1   2
        YDATA:  3   6  11

Here's the formula for l_i(x):

l_i(x) = f₁(x) * ... * f_i-1(x) * f_i+1(x) * ... * f_n(x)

f_j(x) = ( x - xdata_j ) / ( xdata_i - xdata_j)

Exercise: write a MATLAB M function file lag_poly.m of the form

function pval = lag_poly ( xdata, i, xval )

          pval = 1;
          for j = 1 : n
            pval = pval * SOMETHING
          end 

            I:  1   2   3
        XDATA:  0   1   2
        YDATA:  3   6  11
        L3(X):  0   0   1

If you did the exercise correctly, and you used the proper elementwise operations in your function, then you should be able to issue the following nifty commands:

lag_poly ( xdata, 1, xdata )
lag_poly ( xdata, 2, xdata )
lag_poly ( xdata, 3, xdata )

Exercise: now write a MATLAB M function file lagval.m of the form

function pval = lagval ( xdata, ydata, xval )

            I:  1   2   3
        XDATA:  0   1   2
        YDATA:  3   6  11
      

lagval ( xdata, ydata, xdata )

Computing Divided Differences

If you've ever tried to compute the next number in a sequence by looking at the pattern of differences between successive entries, you've started down the path of divided differences. Newton used this method back when there weren't good tables of function values, so that he had to do a lot of interpolation himself. The divided difference table is a third way to define the interpolating polynomial for a set of data.

Here is the algorithm for computing a divided difference table, as taken from our textbook. Since you need to begin to learn to think these things through on your own, I will present the algorithm as the book does, making little attempt to suggest how to write it in MATLAB, except for the function statement:

function d = divdif ( xdata, ydata )

1. Remark: xdata and ydata are vectors with entries xdata_i and f(xdata_i) respectively, for i from 0 to n. On output, d will contain the divided difference table f[xdata₀,...,xdata_n].
2. Copy the entries of ydata into d.
3. Do through step 5 for i from 1 to n.
4. Do through step 5 for j from n down to i.

5. 
               dj := ( dj - dj-1) /( xdataj - xdataj-i)
             

   (Note the final subscript is j - i, not j - 1!)
6. Exit from the algorithm.

Exercise: write a MATLAB M function called divdif.m as described above. If you test your function on our simple data, the resulting d vector should be as follows:

          XDATA:  0   1   2
          YDATA:  3   6  11
          D:      3   3   1

Evaluating Divided Differences

Once we've gone to all the work of setting up our divided differences, we're only halfway home! Now we have to figure out how to evaluate the polynomial represented by the divided difference table at an arbitrary point xval. The book's algorithm is as follows:

function pval = divval ( xval, xdata, d )

1. Remark: xval is a point at which the polynomial is to be evaluated, xdata is the set of abscissas, and d is the divided difference information computed by divdif. On output, pval is the value of the divided difference polynomial at xval.
2. pval := d_n
3. Do through step 4 for i from n-1 down to 0.
4. pval := d_i + ( xval - xdata_i) * pval
5. Exit from the algorithm.

Exercise: write a MATLAB M function called divval.m as described above. Use the divided difference table d from the previous exercise, and evaluate the divided difference polynomial at xval = 3, at which point you should get pval = 18.

Assignment

We never really thought about it, but we assume that an interpolating polynomial p(x) gives us a good approximation to the function f(x) everywhere, no matter what function we choose. If the approximation is not good, we assume it gets better if we increase the number of data points. And if there was an example where the approximation did not get better, we would assume that happens only because the function is "bad". Well, here is a very simple example of how bad things happen to good functions!

Consider interpolating the function:

f(x) = 1 / ( 1 + x * x )

• Write a MATLAB M file called runge.m to evaluate this function.
• Define a set of n = 5 data values xdata which are evenly spaced from -5 to 5, and set ydata to the value of the Runge function at these points.
• Construct the interpolating polynomial to this data.
• Define a set of 101 sampling values x also evenly spaced between -5 and 5.
• Evaluate the runge function and the polynomial at each of these points, and call the vectors of results y and pval.
• Find the maximum absolute value of the difference between the entries of y and pval.

This should take you no more than about 10 lines, assuming you are using vector operations. Repeat the process for n = 10 and 20. Make a table containing three lines:

        N =  5  Max Diff = ???
        N = 10  Max Diff = ???
        N = 20  Max Diff = ???

Isn't there a theorem that says that on a given closed interval, any continuous function can be uniformly well approximated by a polynomial? What is the relationship between your evidence, and this fact?