Mathematics 171 – Numerical Analysis 1 Lab 3 – Polynomial Interpolation in Maple October 22, 1998

Goals

Today, we want to use Maple to

• compute some interpolating polynomials, both directly and with the interp command,

• continue our study of the error term in polynomial interpolation. Recall that if our polynomial p(x) interpolates points on the graph y = g(x), at x = x_{0}, …, x_{n}, then

g(x) − p(x) = \frac{g^{(n+1)}(ξ(x))}{(n + 1)!}(x − x_{0}) ⋯ (x − x_{n})

for some ξ(x) in the interval containing x and the x_{i}.

• examine how round-off error in floating point arithmetic affects polynomial interpolation,

• study some perhaps unexpected phenomena – a case where polynomial interpolation behaves very badly.

A New Maple Command

Maple has a built-in procedure called interp that computes the unique polynomial of degree less than or equal to n in a given variable which interpolates the points in a given set

S = {(x_{0}, y_{0}), (x_{1}, y_{1}), …, (x_{n}, y_{n})}

You call it like this:

                         interp(xlist,ylist,var);

where xlist and ylist are the lists of x- and y- coordinates of the interpolation points respectively and var is the name of the variable for the interpolating polynomial. The lists are entered in square brackets, with the individual list items separated by commas. For instance, you can put the lists of x- and y-coordinates directly into the call to interp, like this:

             interp([1.0,2.0,3.0,4.0],[-1.3,2.4,3.4,-0.9],x);

Or you can define the lists by name, for example with xlst:=[1.0,2.0,3.0,4.0]; and ylst:=[-1.3,2.4,3.4,-0.9];, then call interp supplying the names as inputs to the procedure:

                           interp(xlst,ylst,x);

Lab Exercises

A) Let g(x) = cos(x^{2}). In the first two parts of this question, use Maple’s default 10 decimal digit rounding floating point number system.

1) Compute the unique polynomial p_{4}(x) of degree ≤ 4 interpolating the values of this function at x = 1.0, 1.05, 1.1, 1.15, 1.2 using interp. Plot the interpolating polynomial together with y = g(x) on the same axes for x in the interval 1.0 ≤ x ≤ 1.2. What do you notice about your plot(s)?

2) What is the maximum of the |g^{(5)}(ξ)| for ξ ∈ [1, 1.2]? (Note: If g is an expression containing the variable x, you can compute g^{(5)} in Maple by entering diff(g,x$5). There is a different form for differentiating functions, though.) Using the error term for interpolation, and the bound on |g^{(4)}(ξ)|, what is the maximum possible error

|g(x) − p_{4}(x)|

for x ∈ [1, 1.2]? Discuss how this relates to your plot in part 1.

3) Now change the value of the Digits system variable to 5, so that computations will be done in a five decimal digit floating point number system. Construct the Lagrange interpolating polynomial “by hand” (that is, by constructing the L_{4,k}(x), k = 0, 1, 2, 3, 4, directly and making the appropriate linear combination to get the interpolating polynomial, but using Maple for all computations). Again plot the interpolating polynomial together with y = g(x) on the same axes for x in the interval 1.0 ≤ x ≤ 1.2. What do you notice about your plots now?

4) Still with Digits set to 4, repeat the computation from part 3, but using interp instead of the direct computation, and plot g(x) together with the polynomial.

The results you see in parts 3 and 4 should be quite different. To get some additional information about why the difference occurs, let’s look at another way to compute the interpolating polynomial (closer to what Maple actually did – Maple did not use the Lagrange interpolation formula as you did in part 3). Writing the points as (x_{0}, y_{0}), …, (x_{4}, y_{4}), and noting that x_{i} − x_{i−1} = .05 for all i, let’s write h = .05. The following polynomial

\begin{matrix}y_{0}+ & \frac{(y_{1} − y_{0})}{h}(x − x_{0}) \\ + & \frac{(y_{2} − 2y_{1} + y_{0})}{2h^{2}}(x − x_{0})(x − x_{1}) \\ + & \frac{(y_{3} − 3y_{2} + 3y_{1} − y_{0})}{6h^{3}}(x − x_{0})(x − x_{1})(x − x_{2}) \\ + & \frac{(y_{4} − 4y_{3} + 6y_{2} − 4y_{1} + y_{0})}{24h^{4}}(x − x_{0})(x − x_{1})(x − x_{2})(x − x_{3}) \end{matrix}

is called the divided difference form of the interpolating polynomial.

5) (For after the lab) Verify by hand that in exact arithmetic, this polynomial also interpolates the given points (x_{i}, y_{i}).

6) Also verify that this expression gives results closer to those from Maple’s interp command with Digits set equal to 5. (Suggestion: There is a Maple command called expand that can be used to multiply out products and collect powers of x.)

7) Discuss your results in parts 1 - 6 of this problem. In particular, how do your results in part 4 and in part 3 compare? How do your results in parts 3 and 4 compare with the results from part 1? What conclusions can you draw here about interpolation in floating point arithmetic? Is the Lagrange formula or the divided difference form preferable? Can you see why one might be better than the other?

B) Before starting in on this problem, be sure Digits is set back to the default value 10. All parts of this problem refer to

f(x) = \frac{1}{1 + x^{2}}

on the interval [−5, 5].

1) Compute the interpolating polynomial of degree ≤ 10 matching the values of f at the 11 distinct x-values

−5.0, −4.0, −3.0, −2.0, −1.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0

and plot its graph together with y = f(x) on this interval. The interpolating polynomial should be a relatively good match on [−1, 1] or so, but much worse toward the ends of the interval.

2) A natural idea to improve the situation you see in part 1 is to try interpolating at more points, for example the points with 21 distinct x-values −5 + i/2, for i = 0, 1, …, 20. Repeat part 1. Did this help?

3) One possible explanation for this might be the error term for polynomial interpolation for this function. Find an upper bound B for \frac{f^{(21)}(ξ)}{21!} valid for all ξ ∈ [−5, 5], and use this to plot the “error bound function”

eb(x) = B(x − (−5))(x − (−4.5)) ⋯ (x − 4.5)(x − 5)

How does the plot here relate to the behavior you were seeing in part 2? (Be careful! Does knowing |eb(x)| is big guarantee the error is large?)

Assignment

Individual write-ups due Thursday, October 29.