|
|
We saw last time, with the Runge function, that the interpolating polynomial could get worse as its degree increased. This means that our strategy of using equally spaced data for high degree polynomial interpolation is a bad idea. In this lab, we will test other possible strategies for spacing the data and then we will look at an alternative way to do interpolation with a guarantee that the error will get smaller when we take more points.
This lab will take three sessions. The only plot you need to send to me is the loglog plot in Exercise 9.
When you encounter an error while using an m-file look at the first error message you receive, not the last. This message will identify the file from which it came, and the line number. Be careful: the file that caused the error may not be the file you think you are using--it could be called from the one you think you are using.
The error message "??? Input argument 'x' is undefined." means that the m-file you are using needs to have the variable x included inside parentheses, and you did not include it when you called it. This is the error message you will get if your m-file starts out "function y=func(x)" but you tried to call it as func instead of func(1.03).
In the previous lab, we examined Newton's method of divided differences for computing polynomial interpolants. Part of this lab will be to generate polynomial interpolants for a few different functions on different sets of points. We will be comparing the accuracy of the interpolating polynomials, just as we did last lab. It will be convenient to automate this process in a single m-file so that we are not continually repeating the same steps.
Exercise 1: Construct a utility function with the calling sequence
function max_error=poly_interpolate(f,x)
This function will take as input the name of a function (m-file) to interpolate and a set of points (x) on which to perform the interpolation. The interpolation error will be computed by picking a much larger number of points than are used for interpolation and computing the maximum error on this larger set of test points. An outline of this function follows:
function max_error=poly_interpolate(f,x)
% Choose the number of the test points and generate them
% Last lab we used 101, and this will be good enough for now.
NTEST=101;
% construct NTEST points evenly spaced so that
% xtest(1)=x(1) and xtest(NTEST)=x(length(x))
xtest= ...
% we need the results of f at the points x to do the interpolation
% WARNING: this is a vector statement
y=feval(f,x);
% use divdif and divval from lab6 to do the interpolation
% WARNING: these use vector statements.
d=divdif( ...
ytest=divval( ...
% we will be comparing ytest with the exact results of f at xtest
yexact= ...
% plot the exact and interpolated results on the same plot
% this gives assurance that everything is reasonable
plot( ...
% compute the error. The easiest way is to use the
% Matlab norm function. The "inf" indicates a maximum norm.
max_error=norm(yexact-ytest,inf);As a test, repeat the table from last lab for the Runge function, with N evenly spaced points between -5 and +5.
Runge function, evenly spaced points
N = 5 Max Error = ________
N = 11 Max Error = ________
N = 21 Max Error = ________You should get the values .4382, 1.9156, and 58.278.
After looking at the plotted comparisons between the Runge and its polynomial fit, one possibility for the poor fit is that the points are evenly spaced. Maybe they should be concentrated near the endpoints (so they don't oscillate wildly there) or near the center (maybe the oscillations are caused by so many points near the endpoints). Let's examine these two hypotheses.
Exercise 2: We have been looking at evenly-spaced points between -5 and 5. One way to concentrate these points near zero would be to recall that, for numbers less than one, |x|^3 is smaller than x. Hence, if x represents a sequence of numbers between -5 and 5, then 5*( (x./5).^3 is a similar sequence of numbers concentrated near 0. Fill in the following table for the Runge function with N points concentrated near 0.
Runge function, points concentrated near 0
N = 5 Max Error = ________
N = 11 Max Error = ________
N = 21 Max Error = ________Exercise 3: Points can be concentrated near the endpoints in a similar manner, by forcing them away from zero using the |x|.^(1/3) function. Taking care to get the signs of negative x's correct, fill in the following table with points concentrated near the endpoints.
Runge function, points concentrated near endpoints
N = 5 Max Error = ________
N = 11 Max Error = ________
N = 21 Max Error = ________You could ask why I chose x^3 and x^(1/3). Well, I chose them because I was working with negative numbers and x^2 would not do the trick. Basically, it was an educated guess. It turns out that there is a systematic way to pick 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. (Note: there is more on Chebyshev points and interpolation on p. 228 of Atkinson.)
Theorem 3.2 in Atkinson shows that the interpolation error at any point x is given by an expression of the form:
f(x) - p(x) = Omega ( xdata1, xdata2, ..., xdatandata ) * fndata(xi)
where xi is an unknown "nearby" point, and
Omega ( xdata1, xdata2, ..., xdatandata ) = ( x - xdata1 ) * ... * ( x - xdatandata ) /(n+1)!
This is something like an error term for a Taylor's series.
We can't do a thing about the expression fndata(xi), but we can affect the magnitude of Omega. For instance, if we are only using a single data point in the interval [10,16], then the very best choice for the data point would be xdata1=13, because then the maximum absolute value of Omega would be 3. The worst value would be to choose it at one of the endpoints, in which case we'd double the maximum value of Omega. We're not saying this guarantess better results, just that we're trying what we can to improve our chances.
For a given number ndata of data points, the Chebyshev points minimize the Omega factor in the error estimate. We start by defining some equally spaced angles:
thetai = (2 * i - 1 ) * pi / ( 2 * ndata )
Then, to pick points in the interval [-1,1], the Chebyshev points are:
xi = cos ( thetai )
while for a general interval [a,b], we use the formula:
xi = 0.5 * ( a + b + (a-b) * cos ( thetai ) )
Exercise 4: Write a MATLAB M file called cheby_points.m with the calling statement:
function x = cheby_points ( a, b, ndata )
which returns in the vector x the values of the ndata Chebyshev points in the interval {A,B]. You can do this in 3 lines using vector notation:
i =...; (a range of values)
theta = ...;
x = ...;
To check, here's what you should get for ndata=5, and [a,b]=[10,20]:
1 2 3 4 5
10.2447 12.0611 15.0000 17.9389 19.7553
Exercise 5: The same as exercises 1, 2, and 3 but with Chebyshev points. Fill in the table. (Note the extra row!)
Runge function, Chebyshev points
N = 5 Max Error = ________
N = 11 Max Error = ________
N = 21 Max Error = ________
N = 41 Max Error = ________ We need to write a utility routine to help us with our next interesting task. We want to get the utility routine right, because we need to be able to rely on it. As usual, it's one of those things that's very easy to describe, and not quite so easy to program.
Suppose we have a set of ndata abscissas xdata. We will assume that these are strictly increasing. Now suppose that we are given a number xval and we are asked:
What interval [ xdata(left), xdata(left+1) ] contains xval?
Here's some clarifications we need to make:
If xval is exactly equal to a data value, it really belongs to two intervals. We don't care which one is chosen as the answer;
If xval is less than xdata(1), we would like to receive the value left=1;
If xval is greater than xdata(ndata), we would like to receive the value left=ndata-1 (read that last value carefully!);
Exercise 6: write a MATLAB M file called bracket.m which has a declaration of the form:
function left = bracket ( xdata, xval )
We really want to be able to handle vectors of input data xval, but there is no way to do this using vector notation, so we need to loop over the components of the vector xval. The easiest way to write this kind of routine is to first write a helper function in which xval is a scalar, and then go back and wrap a loop around it. Call the helper function scalar_bracket because it takes a scalar argument.
function left=scalar_bracket(xdata,xval)
ndata = number of x data points
% first check to see if xval is less than xdata(1)
if (xval < xdata(1) )
left = ???
return
elseif (xval >= xdata(ndata))
left = ???
return
else
% loop down through the xdata values to find the
% largest one that is to the left of xval.
for idata = (ndata-1):-1:1
if (xdata(idata) <= xval)
left = ???
return
end
end
end
error('Scalar_bracket: fell through loop.')Now, make up a set of values for xdata (don't forget to make them ascending!), and then test scalar_bracket with selected values that span the range of possible values.
Next, complete the bracket function.
function left_vec = bracket(xdata,xval)
nval = number of intervals to be determined
for i = 1 : nval
left_vec(i)=scalar_bracket(xdata,xval(i));
end
You can put these two functions in the same file, called bracket.m with bracket coming before scalar_bracket in the file. If you do so, you will not be able to call scalar_bracket directly, but it will keep your directory neater. A more compact and efficient way of doing the same task is available in the file bracket.m on my web page.
Again, make up a set of values for xdata and a set of values for xval and see if your program works correctly.
Now we are ready to consider piecewise linear interpolation. The idea is that our interpolating function is not going to be a "smooth" polynomial defined by a formula. Instead, it will simply be defined by patching together linear interpolants that go through each consecutive pair of data points. To handle values below the first data point, we'll just extend the first linear interpolant all the way to the left. Similarly, we'll extend the linear function in the last interval all the way to the right.
Our graph may not look pretty, but one thing we know is this: it will never go through unexpected oscillations between the data values. And the interpolation error is probably going to involve the size of the intervals, and the norm of the second derivative of the function, both of which are easy to understand. As the number of points gets bigger, the interval size will get smaller, and the norm of the second derivative won't change, so we have good reason to hope that we can assert a CONVERGENCE RESULT:
given an interval [a,b] and a function f(x) with a continuous second derivative over that interval, and any sequence Ln(x) of linear interpolants to f with the property that hmax, the size of the maximum interval, goes to zero as n increases, then the linear interpolants converge to f both pointwise, and uniformly.
If C is a bound for the maximum absolute value of the second derivative of the function over the interval, then in fact the pointwise interpolation error, and the L-infinity norm are bounded by C*hmax^2, while the L1 and L2 norms are bounded by (b-a)*C*hmax^2.
Uniform convergence is convergence in the L-infinity or MAX norm, and is a much stronger result than pointwise convergence.
Thus, if the only thing you know is that your function has a bounded second derivative on [a,b], then you are guaranteed that the maximum interpolation error you make decreases to zero as you make your intervals smaller.
Stop and think: The convergence result for piecewise linear interpolation is so easy to get. Is it really better than polynomial interpolation? Why did we have such problems with polynomial interpolation before? One reason is that the error result for polynomial interpolation couldn't be turned into a convergence result. Using 10 data points, we had an error estimate in terms of the 10-th derivative. Then using 20 points, we had another error estimate in terms of the 20-th derivative. These two quantities can't be compared easily, and for nonpolynomial functions, successively higher derivatives can easily get successively and endlessly bigger in norm. The linear interpolation error bound involved a particular fixed derivative and held that fixed, so then it was easy to drive the error to zero. because the other factors in the error term depended only on interval size.
To evaluate this piecewise linear interpolant function at a point xval, we need to:
determine the interval (x(left),x(left+1)) containing xval;
determine the equation of the line which passes through the points (x(left),y(left)) and (x(left+1),y(left+1));
evaluate that line at xval.
Exercise 7: We need to compute the equation through an abitrary pair of points (x1,y1) and (x2,y2). What is the equation of the line through the points (1,17) and (7,47)? If you compute this correctly, you should be able to see that the value of this line at x=4 is 32. Now, can you write a general one-line formula for the linear function through an abitrary pair of points?
Exercise 8: write a MATLAB M file called lintval.m which has a declaration of the form:
function pval = lintval ( xdata, ydata, xval )
Since we really want to be able to handle vectors of input data xval, we need to do some things carefully. Let me sketch one way to set up the routine:
left = something involving bracket
pval = something like exercise 7
Yup, just two lines! This is the power of vector notation.
Exercise 9: Consider the same Runge function we used above. Over the interval [-5,5], use ndata=5 equally spaced data points, compute the piecewise linear interpolant, and evaluate the maximum interpolation error. Plot both the Runge function and the interpolant. Choosing a sequence of trial values of ndata that halve the point spacing from each previous trial, fill out the following table:
Runge function, Piecewise linear
N = 5 Max Error = ________
N = 11 Max Error = ________
N = 21 Max Error = ________
N = 41 Max Error = ________
N = 81 Max Error = ________
N =161 Max Error = ________
N =321 Max Error = ________Plot the values in the table using a log-log plot (the Matlab function loglog is used just like plot, but results in a log-log plot). Your points
should be roughly a straight line on the plot. Why? Roughly what is the slope of the line?
Back to the MATH2070 page.
Last revised on October 22,2000.