We saw last time, with the Runge function, that even though we tried harder, the interpolating polynomial could get worse as the degree increased. This means that our strategy of using equally spaced data for high degree polynomial interpolation is a bad idea. We need a sensible way to do interpolation with a guarantee that when we try harder, we can expect good behavior.
In order to be precise about "good behavior", we will define and measure the interpolation error, as some integral of the difference between the interpolant and the mysterious hidden function that is generating the data.
Our first try at a new interpolation algorithm will be to examine how the location of the data points XDATA affects the result.
Then we will try a new idea, which is a desperate attempt to approximate very wiggly functions with a straight line - well, actually, many itty bitty straight lines.
At the end of this lab is an assignment due by the beginning of the next lab session.
So far, when we've done an interpolation, we've started with some data, constructed the polynomial, and graphed it. Because it's an interpolant, it always goes through the data. So it seems like this is a problem for which there is normally no error to be expected.
Suppose we assume, however, that there really is a function f(x) which is generating this data. For some reason, we need to create an interpolant to this function. (There are good reasons for doing this. One example: the Zeta function). We want to compare the "real" function with the interpolating polynomial p(x). Now we have two functions that we can compare. We ask the question:
How close are f(x) and p(x)?
The two functions are identical at the data points, and probably extremely far apart far away from the data points. A sensible question, then, is, how close are the two functions over some closed interval containing the data points? This is the same as asking for the size or "norm" of the difference function:
g(x) = f(x) - p(x) .There are three common measures of the norm of a function g(x) over an interval [A,B]:
We were just trying to approximate a few data points, and now we're going to have to do integrals? We definitely do not want to try to do an exact integration here. There are many ways to estimate integrals; there's a whole chapter of the book coming about this. Let's try a simple approach. The trapezoidal estimate for the integral of f(x) from A to B using N points and stepsize h=(B-A)/(N-1) is the sum:
EST = h * ( 0.5 * f(A) + f(A+h) + f(A+2h) + ... + f(B-h) + 0.5 * f(B) )For now, we will assume that using the value N = 101 is enough for our computations.
Exercise: create a MATLAB M file l1norm.m which computes the L1 norm of a function, that is, the integral of its absolute value. The function should be declared as:
function value = l1norm ( f, a, b )Define a variable called N, which is the number of points used. We'll set N to 101 for now, but might want to change it later. Remember that when we pass an argument like f, which is the name of some other function routine:
fvalue = feval ( f, x )
result = l1norm ( 'sin', 0, 2*pi )The name that's quoted can be the name of a built in function, or the name of an M file.
Exercise: create MATLAB M files l2norm.m and maxnorm.m which compute the L1 and max norms in a similar way. Note that the max norm calculation is easy. To take the maximum of a and b, you write max(a,b). To take the maximum of a vector x, just write max(x).
When you get done, use your routine to verify the following values:
F(X) A B L1 L2 MAX
---- -- --- ---- ----- -------
1 0 10 10 3.16 1
sin 0 2*pi 4.0 1.77 1
cosy = cos(x) - x 0 5 14.26 7.42 4.7163
runge = 1 / ( 1 + x^2 ) -5 5 2.74 1.25 1
You need an M file called one.m to set up the first function.
This is actually tricky to do, if you want to do a vectorized computation.
You can copy the version of one.m from my home page to see
one way to do it. For f(x)=1, you should get exactly the answers
shown. For some of the others, you should be "close".
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.
In class it was shown that the interpolation error at any point x is given by an expression of the form:
f(x) - p(x) = Omega ( xdata1, xdata2, ..., xdatandata ) * fn(xi)where xi is an unknown "nearby" point, and
Omega ( xdata1, xdata2, ..., xdatandata ) = ( x - xdata1 ) * ... * ( x - xdatandata )This is something like an error term for a Taylor's series.
We can't do a thing about the expression fn(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: 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)To check, here's what you should get for NDATA=5, and [A,B]=[10,20]:
theta = ...;
x = ...;
1 2 3 4 5
10.2447 12.0611 15.0000 17.9389 19.7553
Exercise: Let's recompute the interpolating polynomial to the Runge function
runge(x) = 1 / ( 1 + x^2).When we used equally spaced data in [-5,5], our estimates for the maximum absolute difference were:
L-Infinity Interpolation Error, Evenly spaced points
NDATA max ( | runge(x) - p(x) | )
----- ---------------------------
5 0.43
10 0.29
20 8.16
Now let's repeat this exercise, but use the Chebyshev points to
define the abscissas XDATA.
L-Infinity Interpolation Error, Chebyshev points
NDATA max ( | runge(x) - p(x) | )
----- ---------------------------
5 ..........
10 ..........
20 ..........
We need to write a utility routine to help us with our next interesting task. We want to get the utility routine right first, because it's just a little tricky, and 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:
Exercise: write a MATLAB M file called bracket.m which has a declaration of the form:
function left = bracket ( xdata, 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:
ndata = ???
nval = ???
for ival = 1 : nval
Do checks on xval(ival)...
left(ival) = ???
end
Figuring out the correct interval to assign is slightly trickier
to code than you might imagine!
Once you get your code to where you think it's working, here's a nice check:
xdata = linspace ( -5, 5, 11 );
xval = linspace ( -6, 6, 101 );
left = bracket ( xdata, xval );
plot ( xval, left )
What should you see?
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.
The uniform convergence is convergence in the L-infinity or MAX norm, and is a much stronger result than the 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:
Warmup Exercise: 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: 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 = ???
pval = ???
Yup, just two lines!
Exercise: Consider the "bumpy" function which has the form:
f(x) = 1 / ( ( x - 0.3 )^2 + 0.01 ) + 1 / ( ( x - 0.9 )^2 + 0.04 )(You can copy bumpy.m from my web page if you want.) Over the interval [-2,2], use NDATA equally spaced data points, compute the piecewise linear interpolant, and evaluate the maximum interpolation error. This is a pretty nasty function, so we'll actually be happy if the interpolation error just goes down. However, it will not only go down, but after some initial slowness, it should occasionally drop by a factor of 4 when you halve the interval size. To see this, try using some stretch of the values NDATA = 2, 3, 5, 9, 17, 33, 65 129, 255, 513, 1025, comparing interpolation error at successive values. My initial results were:
NDATA Max Error
2 98.01
3 89.23
5 84.30
which is not a promising beginning!
Assignment: Let's do a similar analysis of the Runge function over the interval [-5,5]. Fill out the following table:
NDATA MAXnorm(InterpError)
----- --------------------
5 0.1800
9 ?
17 ?
33 ?
65 ?
129 ?
255 ?
This is a smoother function, so your results should behave.
Mail the results to me.
Back to the MATH2070 page.
Last revised on 26 October 1999.