M382 – PROJECT PROJECT 1

M382 – PROJECT
PROJECT 1

 I want you submit
1. Diary file
2. Summary file (there is an example in the attached file. It can be txt file or doc
files)
3. Submit m-files and jpeg files .
When you submit your project, I want you to write all of group’s members names who
contribute it.
Notation
This lab is focussed on finding functions that pass through certain specified points.
In customary notation, you are given a set of points and
you want to find a function that passes through each point
( for ). That is, for each of the abscissæ, , the function
value agrees with the ordinate . The given points are variously called the
``given'' points, the ``specified'' points, the ``data'' points, etc. In this lab, we will
use the term ``data points.''
Written in the customary notation, it is easy to see that the quantities and are
essentially different. In Matlab, without kerning and font differentiation, it can be
difficult to keep the various quantities straight. In this lab, we will use the
name xval to denote the values of and the names xdata and ydata to denote the
data and . The variable names fval or yval are used for the value to
emphasize that it is an interpolated value.
Generally, one thinks of the values xval as not equal to any of the data
values xdata although, in this lab, we will often set xval to one or more of
the xdata values in order to test that our interpolation is correct.
Matlab Tips
At some point in this lab, you will need to determine if an integer m is not equal to an integer
n. The Matlab syntax for this is
 if m ~= n
Sometimes it is convenient to enclose the logical expression in parentheses, but this is not
required. Numbers with decimal parts should never be tested for equality! Instead, the
absolute value of the difference should be tested to see if it is small. The reason is that
numbers are rarely known to full 16-digit precision and, in addition, small errors are usually
made in representing them in binary instead of decimal form. To check if a number x is equal
to y, you should use
 if abs(x-y) <= TOLERANCE
where TOLERANCE represents some reasonable value chosen based on the problem.
If you have a (square) linear system of equations, of the form
 A * x = b
where A is an N by N matrix, and b and x are column vectors, then Matlab can solve the linear
system either by using the expression:
 x = inv(A)*b
or, better, via the ``backslash'' command:
 x = A \ b
that you used in the previous labs. The backslash command is strange looking. You might
remember it by noting that the matrix A appears underneath the backslash. The backslash is
used because matrix multiplication is not commutative, so a square matrix should appear to
the left of a column vector. The difference between the two commands is merely that the
backslash command is about three times faster because it solves the equation without
constructing the inverse matrix. You will not notice a difference in this lab, but you might see
one later if you use Matlab for more complicated projects.
The backslash command is actually more general than just multiplication by the inverse. It
can find solutions when the matrix A is singular or when it is not a square matrix! We will
discuss how it does these things later in the course, but for now, be very careful when you use
the backslash operator because it can find ``solutions'' when you do not expect a solution.
Matlab represents a polynomial as the vector of its coefficients. Matlab has several commands
for dealing with polynomials:
c=poly(r)
Finds the coefficients of the polynomial whose roots are given by the vector r.
r=roots(c)
Finds the roots of the polynomial whose coefficients are given by the vector c.
c=polyfit( xdata, ydata, n)
Finds the coefficients of the polynomial of degree n passing through, or
approaching as closely as possible, the points xdata(k),ydata(k),
for , where need not be the same as n.
yval=polyval( c, xval)
Evaluates a polynomial with coefficients given by the vector c at the
values xval(k) for for .
When there are more data values than the minimum, the polyfit function returns the
coefficients of a polynomial that ``best fits'' the given values in the sense of least squares. This
polynomial approximates, but does not necessarily interpolate, the data. In this lab, you will
be writing m-files with functions similar to polyfit but that generate polynomials of the
precise degree determined by the number of data points. (N=numel(xdata)-1).
The coefficients of a polynomial in Matlab are, by convention, defined as
(1)
and in Matlab are represented as a vector c. Beware: With this definition of the vector c of
coefficients of a polynomial:
1. N=numel(c) is one higher than the degree of p(x).
2. The subscripts run backwards! c(1) is the coefficient of the term with degree
N-1, and the constant term is c(N).
In this and some later labs, you will be writing m-files with functions analogous to polyfit
and polyval, using several different methods. Rather than following the matlab naming
convention, functions with the prefix coef_ will generate a vector of coefficients, as by
polyfit, and functions with the prefix eval_ will evaluate a polynomial (or other function) at
values of xval, as with polyval.
In the this lab and often in later labs we will be using Matlab functions to construct a known
polynomial and using it to generate ``data'' values. Then we use our interpolation functions to
recover the original, known, polynomial. This strategy is a powerful tool for illustrative and
debugging purposes, but practical use of interpolation starts from arbitrary data, not contrived
data.
The Polynomial through Given Data
we discussed Newton's method, which can be derived by determining a linear
polynomial (degree 1) that passes through the point with derivative . That
is and . One way of looking at this is that we are
constructing an interpolating function, in this case a linear polynomial, that
explains 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.
Suppose we wanted to determine the linear polynomial that passes
through the data points and . We simply have to solve a set of
linear equations for and constructed by plugging in the two data points into
the general linear polynomial. These equations are


or, equivalently,
which (usually) has the solution
Compare that situation with the case where we want to determine the quadratic
polynomial that passes through three sets of data values.
Then we have to solve the following set of three linear equations for the
polynomial coefficients c:
These are examples of second and third order Vandermonde Equations. It is
characterized by the fact that for each row (sometimes column) of the coefficient
matrix, the succesive entries are generated by a decreasing (sometimes increasing)
set of powers of a set of variables.
You should be able to see that, for any collection of abscissæ and ordinates, 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.
Now, let's see how to construct and solve the Vandermonde equation using Matlab.
This involves setting up the coefficient matrix A. We use the Matlab
variables xdata and ydata to represent the quantities and , and we will
assume them to be row vectors of length (numel) N.
for j = 1:N
 for k = 1:N
 A(j,k) = xdata(j)^(N-k) ;
 end
end
Then we have to set the right hand side to the ordinates ydata, that is assumed to be
a row vector. If we can get all of that set up, then actually solving the linear system
is easy. We just write:
c = A \ ydata';
Recall that the backslash symbol means to solve the system with matrix A and right
side ydata'. Notice that ydata' is the transpose of the row vector ydata in this
equation. (By default, Matlab constructs row vectors unless told to do otherwise.)
Exercise 1: The Matlab built-in function polyfit finds the coefficients of a polynomial
through a set of points. We will write our own using the Vandermonde matrix. (This is the
way that the Matlab function polyfit works.)
1. Write a Matlab function m-file, coef_vander.m with signature
2. function c = coef_vander ( xdata, ydata )
3. % c = coef_vander ( xdata, ydata )
4. % xdata= ???
5. % ydata= ???
6. % c= ???
7. % other comments
8.
9. % your name and the date
that accepts a pair of row vectors xdata and ydata of arbitrary but equal
length, and returns the coefficient vector c of the polynomial that passes
through that data. Be sure to complete the comments with question marks in
them.
Warning: Think carefully about what to use for N.
10. Test your function by computing the coefficients of the polynomial through the
following data points. (This polynomial i is , so you can check your
coefficient vector ``by inspection.'')
11. xdata= [ 0 1 2 ]
12. ydata= [ 0 1 4 ]
13. Test your function by computing the coefficients of the polynomial that passes
through the following points
14. xdata= [ -3 -2 -1 0 1 2 3]
15. ydata= [1636 247 28 7 4 31 412]
16. Confirm using polyval that your polynomial passes through these data points.
17. Double-check your work by comparing with results from the Matlab polyfit
function. Please include both the full polyfit command you used and the
coefficient vector it returned in your summary.
In the following exercise you will construct a polynomial using coef_vander to interpolate
data points and then you will see what happens between the interpolation points.
Exercise 2:
1. Consider the polynomial whose roots are r=[-2 -1 1 2 3]Use the Matlab poly
function to find its coefficients. Call these coefficients cTrue.
2. This polynomial obviously passes through zero at each of these five ``data''
points.We want to see if our coef_vander function can reproduce it. To use
our coef_vander function, we need a sixth point
3. . You can ``read off'' the value of the polynomial at x=0 from its coefficients cTrue.
4. What is this value?
5.
6. Use the coef_vander function to find the coefficients of the polynomial passing through the
``data'' points
7. xdata=[ -2 -1 0 1 2 3];
8. ydata=[ 0 0 ?? 0 0 0];
9. Call these coefficients cVander.
10. Use coef_vander to find the coefficients using only the five roots as xdata. Name these
coefficients something other than cVander. Explain your results.
11. Use the following code to compute and plot the values of the true and interpolant polynomials
on the interval [-3,2]. If you look at the last line of the code, you will see an estimate of the
difference between the two curves. How big is this difference? (We will be using essentially this
same code in several following exercises. You should be sure you understand what it does. You
might want to copy it to a script m-file.)
12. xval=linspace(-3,2,4001); % test abscissae for plotting
13. yvalTrue=polyval(cTrue,xval); % true ordinates
14. yvalVander=polyval(cVander,xval); % our ordinates
15. plot(xval,yvalTrue,'g','linewidth',4); % true curve: thick green
16. hold on
17. plot(xval,yvalVander,'k'); % interpolant curve: thin black
18. hold off
19. max(abs((yvalTrue-yvalVander)))/max(abs(yvalTrue))
Please include a copy of this plot with your summary. (You should observe the two curves are the
same. The second curve appears as a thin black curve overlaying the thick green curve.) Note: The
relative error is used here so that the case where yvalTrue are very large or very small numbers
does not cause difficulty..
Lagrange Polynomials
Suppose we fix the set of distinct abscissæ , and think about the
problem of constructing a polynomial that has (not yet specified) values at these
points. Now suppose I have a polynomial whose value is zero at each
, , and is 1 at . Then the intermediate polynomial would have the
value at , and be 0 at all the other . Doing the same for each abscissa and
adding the intermediate polynomials together results in the polynomial that
interpolates the data without solving any equations!
In fact, the Lagrange polynomials are easily constructed for any set of
abscissae. Each Lagrange polynomial will be of degree . There will
be Lagrange polynomials, one per abscissa, and the polynomial will
have a special relationship with the abscissa , namely, it will be 1 there, and 0 at
the other abscissæ.
In terms of Lagrange polynomials, then, the interpolating polynomial has the form:
(2)
Assuming we can determine these magical polynomials, this is a second way to
define the interpolating polynomial to a set of data.
, Remark: The strategy of finding a function that equals 1 at a distinguished point
and zero at all other points in a set is very powerful. If in (2), then ,
so the are an example of a ``partition of unity.'' One of the places you will
meet it again is when you study the construction of finite elements for solving
partial differential equations.
In the next two exercises, you will be constructing polynomials through the same
points as in the previous exercise. Since there is only one nontrivial polynomial of
degree (N-1) through N data points, the resulting interpolating polynomials are the
same in these and the previous exercise.
In the next exercise, you will construct the Lagrange polynomials associated with
the data points, and in the following exercise you will use these Lagrange
polynomials to construct the interpolating polynomial.
Exercise 3: In this exercise you will construct Lagrange polynomials based
on given data points. Recall the data set for :
 k : 1 2 3
xdata= [ 0 1 2 ]
ydata= [ 0 1 4 ]
(Actually, ydata is immaterial for construction of .) In general, the formula
for can be written as:
(3)
(skipping the factor), where each factor has the form
 1. Explain in a sentence or two why the formula (3) using the factors (4)
yield a
(4
)
2. Write a Matlab function m-file called lagrangep.m that computes the
Lagrange polynomials (3) for any k. (One of the Matlab toolboxes has a
function named ``lagrange'', so this one is named ``lagrangep''.) The
signature should be
function pval = lagrangep( k , xdata, xval )
% pval = lagrangep( k , xdata, xval )
% comments
% k= ???
% xdata= ???
% xval= ???
% pval= ???
% your name and the date
and the function should evaluate the k-th Lagrange polynomial for the abscissæ
xdata at the point xval. Hint, you can implement the general formula using
code like the following.
pval = 1;
for j = 1 : ???
 if j ~= k
 pval = pval .* ??? % elementwise multiplication
 end
end
Note: If xval is a vector of values, then pval will be the vector of
corresponding values, so that an elementwise multiplication (.*) is being
performed.
3. Using (5), determine the values of lagrangep( 1, xdata, xval) for
xval=xdata(1), xval=xdata(2) and xval=xdata(3).
4. Does lagrangep give the correct values for lagrangep( 1, xdata, xdata)?
For lagrangep( 2, xdata, xdata)? For lagrangep( 3, xdata, xdata)?
Exercise 4: The hard part is done. Now we want to use our lagrangep routine as a
helper for a second replacement for the polyfit-polyval pair, called eval_lag, that
implements Equation (2). Unlike coef_vander, the coefficient vector of the
polynomial does not need to be generated separately because it is so easy, and that is
why eval_lag both fits and evaluates the Lagrange interpolating polynomial.
1. Write a Matlab function m-file called eval_lag.m with the signature
function yval = eval_lag ( xdata, ydata, xval )
% yval = eval_lag ( xdata, ydata, xval )
% comments
% your name and the date
This function should take the data values xdata and ydata, and compute the
value of the interpolating polynomial at xval according to (2), using your
lagrangep function for the Lagrange polynomials. Be sure to include
comments to that effect.
2. Test eval_lag on the simplest data set we have been using.
 k : 1 2 3
 xdata= [ 0 1 2 ]
 ydata= [ 0 1 4 ]
by evaluating it at xval=xdata. Of course, you should get ydata back.
3. Test your function by interpolating the polynomial that passes through the
following points, again by evaluating it at xval=xdata.
4. xdata= [ -3 -2 -1 0 1 2 3]
5. ydata= [1636 247 28 7 4 31 412]
6. Repeat Exercise 2 using Lagrange interpolation.
o Return to the polynomial constructed in Exercise 2 with roots r=[-2 -
1 1 2 3]. and coefficients cTrue.
o Reconstruct (using polyval) or recall the ydata values associated with
xdata=[-2 -1 0 1 2 3].
o Compute the values of the Lagrange interpolating polynomial at the
same 4001 test points between -2 and 3. Call these values yvalLag.
o Plot yvalLag and yvalTrue and compute the error between them using
the test plotting approach used in Exercise 2. Include a copy of this plot
with your summary.
Interpolating a function that is not a
polynomial
Interpolating functions that are polynomials and using polynomials to do it is cheating a little
bit. You have seen that interpolating polynomials can result in interpolants that are essentially
identical to the original polynomial. Results can be much less satisfying when polynomials
are used to interpolate functions that are not themselves polynomials. At the interpolation
points, the function and its interpolant agree exactly, so we want to examine the behavior
between the interpolation points. In the following exercise, you will see that some nonpolynomial functions can be interpolated quite well, and in the subsequent exercise you will
see one that cannot be interpolated well. The example used here is due in part to C. Runge,
Über empirische Funktionen und die Interpolation zwischen äquidistanten Ordinaten, Z.
Math. Phys., 46(1901), pp. 224-243.
Exercise 5: In this exercise you will construct interpolants for the hyperbolic
sine function and see that it and its polynomial interpolant are quite
close.
1. We would like to interpolate the function on the
interval , so use the following outline to examine the behavior
of the polynomial interpolant to the exponential function for five
evenly-spaced points. It would be best if you put these commands into
a script m-file named exer5.m.
2. % construct N=5 data points
3. N=5;
4. xdata=linspace(-pi,pi,N);
5. ydata=sinh(xdata);
6.
7. % construct many test points
8. xval=linspace(???,???,4001);
9. % construct the true test point values, for reference
10. yvalTrue=sinh(???);
11.
12. % use Lagrange polynomial interpolation to evaluate
13. % the interpolant at the test points
14. yval=eval_lag(???,???,xval);
15.
16. % plot reference values in thick green
17. plot(xval,yvalTrue,'g','linewidth',4);
18. hold on
19. % plot interpolation data points
20. plot(xdata,ydata,'k+');
21. % plot interpolant in thin black
22. plot(xval,yval,'k');
23. hold off
24.
25. % estimate the approximation error of the interpolant
26. approximationError=max(abs(yvalTrueyval))/max(abs(yvalTrue))
Please send me this plot.
27.By zooming, etc., confirm visually that the exponential and its
interpolant agree at the interpolation points.
28.Using more data points gives higher degree interpolation polynomials.
Fill in the following table using Lagrange interpolation with
increasing numbers of data points.
29. N = 5 Approximation Error = ________
30. N = 11 Approximation Error = ________
31. N = 21 Approximation Error = ________
You should have observed in Exercise 5 that the approximation error becomes
quite small. The exponential function is entire, as are polynomials, so they share
one essential feature. In the following exercise, you will see that attempts to
interpolate functions that are not entire can give poor results.
Exercise 6:
1. Construct a function m-file for the Runge example function .
Name the file runge.m and give it the signature
2. function y=runge(x)
3. % y=runge(x)
4. % comments
5.
6. % your name and the date
Use componentwise (vector) division and exponentiation (./ and .^).
7. Copy exer5.m to exer6.m and modify it to use the Runge example
function. Please send me the plot it generates.
8. Confirm visually that the Runge example function and its interpolant
agree at the interpolation points, but not necessarily between them.
9. Using more data points gives higher degree interpolation polynomials.
Fill in the following table using Lagrange interpolation with
increasing numbers of data points.
10. N = 5 Approximation Error = ________
11. N = 11 Approximation Error = ________
12. N = 21 Approximation Error = ________
13.Are you surprised to see that the errors do not decrease?
Many people expect that an interpolating polynomial gives a good
approximation to the function everywhere, no matter what function we
choose. If the approximation is not good, we expect it to get better if we increase
the number of data points. These expectations will be fulfilled only when the
function does not exhibit some ``essentially non-polynomial'' behavior. You will
see why the Runge example function cannot be approximated well by polynomials
in the following exercise.
Exercise 7: The Runge example function has Taylor series
(6)
as you can easily prove. This series has a radius of convergence of 1 in the
complex plane. Polynomials, on the other hand, are entire
functions, i.e., their Taylor series converge everywhere in the complex plane.
No finite sum of polynomials can be anything but entire, but no entire
function can interpolate the Runge example function on a disk with radius
larger than one about the origin in the complex plane. If there were one, it
would have to agree with the series (6) inside the unit disk, but the series
diverges at and an entire function cannot have an infinite value (a
pole).
1. Make a copy of exer6.m called exer7.m that
uses coef_vander and polyval to evaluate the interpolating polynomial
rather than eval_lag.
2. Confirm that you get the same results as in Exercise 6 when you use
Vandermonde interpolation for the Runge example function.
3. N = 5 Approximation Error = ________
4. N = 11 Approximation Error = ________
5. N = 21 Approximation Error = ________
6. Look at the nontrivial coefficients ( ) of the interpolating
polynomials by filling in the following table.
7. N=5 c( 5)= _____ c( 3)= _____ c( 1)= _____
8. N=11 c(11)= _____ c( 9)= _____ c( 7)= _____ c( 5)= _____
9. N=21 c(21)= _____ c(19)= _____ c(17)= _____ c(15)= _____
10.Look at the trivial coefficients ( ) of the interpolating polynomials
by filling in the following table. (Look carefully at what the colon
notation does.)
11. N=5 max(abs(c(2:2:end)))= _____
12. N=11 max(abs(c(2:2:end)))= _____
13. N=21 max(abs(c(2:2:end)))= _____
You should see that the interpolating polynomials are ``trying'' to reproduce
the Taylor series (6). These polynomials cannot agree with the Taylor series
at all points, though, because the Taylor series does not converge at all
points.
Chebyshev Points
If we have no idea what function is generating the data, but we are allowed to pick
the locations of the data points beforehand, the Chebyshev points are the smartest
choice.
We show that the interpolation error between a function and its polynomial
interpolant at any point x is given by an expression of the form:
(7)
where is an unknown nearby point. This is something like an error term for a
Taylor's series.
We can't do a thing about the expression , because is an arbitrary (smooth
enough) function, but we can affect the magnitude of the bracketed expression. For
instance, if we are only using a single data point ( ) in the interval [10,20],
then the very best choice for the data point would be , because then the
maximum absolute value of the expression
would be 5. The worst value would be to choose at one of the endpoints, in
which case we'd double the maximum value. This doesn't guarantee better results,
but it improves the chance.
The Chebyshev polynomial of degree is given by
(8)
This formula doesn't look like it generates a polynomial, but the trigonometric
formulæ for sums of angles and the relation that can be used to
show that it does generate a polynomial.
These polynomials satisfy an orthogonality relationship and a three-term recurrence
relationship
with and (9)
The facts that make Chebyshev polynomials important for interpolation are
 The peaks and valleys of are smallest of all polynomials of
degree over [-1,1] with
and
 On the interval [-1,1], each polynomial oscillates about zero with peaks and
valleys all of equal magnitude
 Thus, when in (7) are chosen to be the roots of , then
the bracketed expression in (7) is proportional to , and the bracketed
expression is minimized over all polynomials.
Before going on to use the Chebyshev points, it is a good idea to confirm that the
two expressions (8) and (9) do, in fact, refer to the same polynomials.
Mathematically speaking, the best way to approach the problem is to use the
symbolic toolbox, which can produce the identity that can become the central
portion of a proof. Instead, we will approach the problem numerically. This
approach will yield a convincing example, but no proof.
Exercise 8:
1. Write a Matlab function with the signature and first few lines,
2. function tval=cheby_trig(xval,degree)
3. % tval=cheby_trig(xval,degree)
4.
5. % your name and the date
6.
7. if nargin==1
8. degree=7;
9. end
to evaluate as defined using trigonometric functions in (7)
above. If the argument degree is omitted, it defaults to 7, a fact that
will be used below.
10.Write a recursive Matlab function with the signature and first few
lines,
11. function tval=cheby_recurs(xval,degree)
12. % tval=cheby_recurs(xval,degree)
13.
14. % your name and the date
15.
16. if nargin==1
17. degree=7;
18. end
to evaluate as defined recursively in (8) above. If the
argument degree is omitted, it defaults to 7, a fact that will be used
below.
19.Show that cheby_trig and cheby_recurs agree for degree=4 ( ) at
the points x=[0,1,2,3,4] by computing the largest absolute value of
the differences at those points. What is this value?
Remark: If two polynomials of degree 4 agree at 5 points, they agree
everywhere. Hence if (7) defines a polynomial, then (7) and (8)
produce identical values for . This is why only five test points are
required.
20. Plot values of cheby_trig and cheby_recurs on the interval [-1.1, 1.1]
for degree=7 ( ). Use at least 100 points for this plot in order to see
the details. The two lines should lie on top of one another. You should
also observe visually that the peaks and valleys of are of equal
absolute value. Please include this plot with your summary.
Constructing the Chebyshev points
The Chebyshev points are the zeros of the Chebyshev polynomials. They can be
found from (7):




For a given number of data points, then, the Chebyshev points on the
interval can be constructed in the following way.
1. Pick equally-spaced angles , for .
2. The Chebyshev points on the interval are given as
for .
Exercise 9:
1. Write a Matlab function m-file called cheby_points.m with the
signature:
2. function xdata = cheby_points ( a, b, ndata )
3. % xdata = cheby_points ( a, b, ndata )
4. % more comments
5.
6. % your name and the date
that returns in the vector xdata the values of the ndata Chebyshev
points in the interval [a,b]. You can do this in 3 lines using vector
notation:
k = (1:ndata); %vector, NOT A LOOP
theta = (vector expression involving k)
xdata = (vector expression involving theta)
or you can use a for loop. (The vector notation is more compact and
runs faster, but is a little harder to understand.)
7. To check your work, use cheby_points to find the Chebyshev points
on the interval [-1,1] for ndata=7. Do the largest and second largest
roots agree with your roots of computed above?
8. What are the five Chebyshev points for ndata=5 and [a,b]=[-5,5]?
9. You should observe that the Chebyshev points are not uniformly
distributed on the interval but they are symmetrical about its center. It
is easy to see from (7) that this fact is true in general.
Exercise 10: Repeat Exercise 4 but with Chebyshev points on [-5,5]. Fill in
the table. (Note the extra row!) You should see smaller errors than before,
especially for larger ndata.
 Runge function, Chebyshev points
 ndata = 5 Max Error = ________
 ndata = 11 Max Error = ________
 ndata = 21 Max Error = ________
 ndata = 41 Max Error = ________