At what values does the coefficient computation break down

The question asks to use polyfit to find the interpolating polynomial of y=e^x at N equally spaced points and determine the value of N at which the computation of coefficients break down. What does this mean?

Further it asks to graph it using polyval and find N at which the instability of equally spaced interpolation makes the interpolant useless. I notice a break in the graph, is it what this part of the question is asking for or is it something else?

[John D'Errico]

It means what it says. TRY it for different values of N. LOOK at the coefficients. THINK about what you see. What happens between the points when you plot the interpolant? Does something strange start happening as N grows larger? For N large, and x large or small, what happens to those large powers of x? How many significant digits of precision does a double precision number carry? What does eps tell you in MATLAB?

In MATLAB, what is the result of

1 + eps/10 - 1

Might this have anything to do with the above questions?

Another idea to consider. In general, what is different about the shape of higher order polynomials, compared to that of low order polynomials? (Think about the word flexible here.)

The assignment asks you to do something, and then to think about what you see.

Answer 1

[John D'Errico]

Oh, I suppose it seems as if you made some effort. So let me write a small dissertation on the topic.

There are several issues at work here.

The first one is in terms of finite precision arithmetic. The second is the flexibility of a polynomial as the order of that polynomial grows. The third reason I can think of is related to the linear algebra (or whatever computations are used) to compute the coefficients of the polynomial. I imagine there are other issues, if I thought a bit.

They will all have a serious impact on a polynomial interpolant. As well, they are good reasons to avoid high order polynomials for interpolation, in favor of different tools, like splines.

So what does eps tell you? eps is essentially the smallest number that you can add to 1, and have the result be different from 1. For example, Try these tests in MATLAB:

(1 + eps) == 1  
ans =
   0

(1 + eps/2) == 1
ans =
   1

Think of eps as the least significant bit. Suppose I told you to work in 4 digit floating point arithmetic. You can then represent the numbers 1.111 and 0.001111 equally well, since both have the same mantissa, just differing exponents. However, what happens when you add the two numbers? In mathematics, we know that

1.111 + 0.00111 = 1.112111

However, in 4 digits of precision, we would have the only approximate result of

1.111 + 0.00111 = 1.112

Essentially we have lost those low order bits. This is why we can add eps to 1, and get a different number from 1, but when we add eps/2 to 1, we get 1.

(1 + eps/2) == 1
ans =
   1

EXACTLY 1.

Now, think about these ideas in context of a polynomial. Suppose we will choose to evaluate the polynomial

1 + x + x^2 + x^3 + x^4 + ... + x^19 + x^20

(Yeah, I know, I'm not feeling very creative today.) Do so over a range of values for x that vary from x=0 to x=10.

When x is near 1, things are very well behaved. However, what happens when x is 10? Remember, that just like the example I gave about adding eps to 1, what do you get when you do the computation

1 + 10^20

As you would expect, since a double in MATLAB carries only 52 binary bits in the mantissa, we will lose that stuff that falls below eps(10^20).

So high order polynomials are suspect when we are forced to add and subtract the numbers from each term.

Next, lets look at what happens in terms of flexibilty. In fact, the exponential function should be moderately well behaved for these purposes. I'll talk about this later. What do I mean about flexibility? A zeroth order polynomial (a constant) is rather boring, a straight line a tiny bit less so. A quadratic adds some shape, with an extremum someplace. A cubic will add more. As you go up in order, you can see a wide variety of bends and wiggles though.

When you build an interpolating polynomial, all that matters is what it does at the data points. In between however, it can do as it will. Nothing stops that interior behavior from being quite nasty.

Ugh. Sadly, i've got to run for a bit. I'll return to finish this later today.

[Sachi]

Thanks a lot.

So the computation breaks at N=13?

[John D'Errico]

Oh, I'll admit that MATLAB essentially concedes defeat at that point, AND that is probably the "breaking" point your teacher wanted you to find. But did it truly break at that point? It was probably producing semi-garbage a wee bit before then, and it just did not realize the problem.

This is not truly a question of it working perfectly before then, and producing crapola afterwards. Is near garbage also garbage? Or is it just a bit trashy? Lets just say that I would not trust the last couple of polynomials before that point.

To assert more than that, I would need to know many things. For example, what is the domain over which you have chosen to approximate that exponential function? What is your tolerance for error in the result over that domain?

Yes, in theory, higher order polynomials will approximate an exponential better than lower order ones, since they will reduce the "lack of fit". (Interpolation error if you prefer.) But at some point, numerical trash in the solution will come to dominate that lack of fit.

So not without doing a bit of analysis could I truly know when to quit, and not without a bit more information from you. However, safety (and logic) says you should never push anything all the way to a breaking point. In fact, it is very rarely ever my recommendation to use polynomials anywhere near that high of an order. Splines of some ilk are almost always a more rational choice.