LAB #5: Polynomials and Newton's Method

TABLE OF CONTENTS

Introduction

This is a long, complicated project. We will allow 3 lab days for it.

Polynomials occur so often in mathematical calculations that it is important to have a good idea of:

Vocabulary - A polynomial of degree n has a highest term of x^n, which means it has n+1 coefficients. A polynomial of order n has n coefficients, for terms from the constant through x^(n-1). The term degree is used in Atkinson and is the more common one. Unfortunately, MATLAB deals with polynomials as vectors of coefficients, and the length of the vector of coefficients is the order of the polynomial. We will be using both terms.

MATLAB Tips

Let's recall what we already know about MATLAB and polynomials, and discuss some other MATLAB commands that will be useful.

MATLAB will not let us use zero as a vector index. So, unfortunately, if we want to store the coefficients of a polynomial in a vector, the index of a coefficient can't be made to match the exponent of x. Our convention will be that if the coefficients of a polynomial p(x) of order n are stored in the vector c, then

p(x) = c(1) * x^(n-1) + c(2) * x^(n-2) + ... + c(n-1) * x + c(n)

MATLAB has a built in command that, given the coefficients c of a polynomial, will evaluate it at a point x

px = polyval ( c, x )

The value c will be a vector, and x can be a vector too, in which case we evaluate the polynomial at all the given points. We can even enter the values of c or x as part of the command:

polyval ( [ 1, -2, 12 ], [ 0, 1, 2, 3, 4 ] )

The polyder command computes the coefficients of the derivative polynomial. So to evaluate the derivative at a point requires two steps:

cprime = polyder ( c )
pp = polyval ( cprime, x )

We'll be writing our own versions of these routines today!

MATLAB can take a collection of roots ( real or imaginary ) and construct the coefficients of the corresponding polynomial. The command is poly(v). If I set

v = [ 0, 1, 2, 3 ];
c = poly ( v );

MATLAB responds with

c = [ 1, -6, 11, -6, 0 ]

Exercise 1: Make a polynomial with roots -2, 1 + 2i, 1 - 2i, and 3. Name the coefficient vector c and check to be sure it is real. Evaluate your polynomial at the first two roots to verify that it is zero there. MATLAB can also solve for the roots, given the coefficients. Try the command:

v = roots ( c )

to see if you get back your original set of roots!

Horner's Method

Suppose someone asks you to evaluate the polynomial

p(x) = x^3 - 2 * x^2 - 5 * x + 6

Lucky for you, they want to know its value at x = 0. Presto, we're done already, and can simply read off the answer. Now suppose we wanted to know the value at x = 2? That's not so easy. We'd need a pen and paper, or a calculator. But suppose I were able to rewrite the polynomial as:

p(x) = ( x^2 - 5 ) * ( x - 2 ) - 4

Then, again, the value of p(2) is easily seen.

The important idea to see here is that whenever we can rewrite a polynomial as

p(x) = q(x) * ( x - a ) + constant

then for the special argument x=a, we know p(a)=constant. From this idea, Horner determined a method that can be used to evaluate a polynomial, to evaluate its derivative, or, for any desired value a, to rewrite a polynomial as

p(x) = q(x) * ( x - a ) + p(a)

This is essentially a form of synthetic division, in which we are always dividing by a linear factor of the form x-a.

Horner's Method: To evaluate a polynomial with coefficients c at the point x, do the following:

        px = c(1)
        for i = 2:n
          px = px * x + c(i)
        end

Please don't get confused here. I've switched gears, and now I'm thinking of x as a particular value, like 7, rather than as a symbol standing for any value. The variable px stands for "p of x".

Programming tip: to use this method, we need to determine the order n of the polynomial. We could pass this in as an extra argument, but isn't MATLAB smart enough to know how many entries there are in the coefficient vector? Yes, it is, and we have seen the length function for a vector n = length(c) .

In one sense, this is no big deal. Anybody can evaluate a polynomial. But think of this in computational terms. Horner's method uses a significantly smaller number of multiplications. To see this, let's first determine the number of multiplications required to evaluate

p(x) = 2 * x^4 + 7 * x^3 + 6 * x^2 + 8 * x + 12

using the "natural" way to evaluate the polynomial. To make it more obvious, write it out as

p(x) = 2*x*x*x*x + 7*x*x*x + 6*x*x + 8*x + 12
and counting operation signs gives 10 multiplications and 4 additions.

Exercise 2: Determine the number of multiplications and additions involved in Horner's method.

Exercise 3: If p(x) were a polynomial of degree 20, how many multiplications would be involved?

Exercise 4: Working from the outline above, write an M file called horner.m of the form

function px = horner ( c, x )

that evaluates a polynomial at a point x. When you're done, try it out on the polynomial

p(x) = x^3 - 2 * x^2 + 3 * x - 4

and confirm the following values: 

        X:       -1   0   1   2   3
        P(X):   -10  -4  -2   2  14

Horner's Method Plus Derivative

If we want to use Newton's method, we're going to need not just the polynomial, but its derivative as well. We could do this many wrong ways, but it turns out there's a very efficient way that "piggybacks" on top of Horner's method.

Horner's Method Plus Derivative: 

        pprimex = 0.0
        px = c(1)
        for i = 2:n
          pprimex = pprimex * x + px
          px = px * x + c(i)
        end

Exercise 5: Actually, it may be easier to type in the change that computes the derivative than it is to understand it. Explain what is going on in this algorithm in words, if you can, or by illustrating it on a simple polynomial.

This algorithm computes both the value of the polynomial p(x), px, and the value of its derivative p'(x), pprimex. Matlab provides the capability of returning both of these values in the result of a function. The way this is accomplished is to make a vector of the two values and put it on the left of the equal sign. When you call the function, you need to remember that the function returns a vector and not a scalar result. The exercise below uses this syntax. You can get detailed information on how to do these tasks by getting help for the function statement by typing help function at the MATLAB prompt and by getting help for the eval statement by typing help eval. As an alternative, you can go to the MATLAB help site and get help on these statements by typing their names in the box one at a time.

Exercise 6: Working from the outline above, modify your M file horner.m so it now has the form:

function [ px, pprimex ] = horner ( c, x )

and evaluates a polynomial and its derivative at a point x. When you're done, try it out on the polynomial

p(x) = x^3 - 2 * x^2 + 3 * x - 4

which has the following values: 

        X:       -1   0   1   2   3
        P(X):   -10  -4  -2   2  14
        P'(X):   10   3   2   7  18

A Newton-Horner Method

Now if Horner's method makes it easy to compute the value and derivative of a polynomial at any point x, then we are all set to use Newton's method! Instead of writing two functions that evaluate the function and its derivative, we just pass in the coefficients of the polynomial.

Here is a sketch of what an algorithm to compute one root of a polynomial starting from an initial guess x would look like: 

      function root = polynew ( c, x )

      FATOL = ?
      STEPMAX = ?

      for step = 1: STEPMAX

        px = ?
        pprimex = ?

        if ( | px | <= FATOL ) 
          ...we're happy...
        end 

        if ( pprimex is zero )
          'Badness'
        end

        x = x - px / pprimex
        
      end

      'Too many steps'

I'm deliberately being very sloppy in this outline, for two reasons. You need to think, as you fill in the correct steps, and I want to encourage you to develop your programming ideas by starting out with a vague sketch like this, and then filling in the details. Get the overall picture right first, and then fill it in!

Exercise 7: Working from the outline above, write an M file called polynew.m that accepts the coefficients of a polynomial and a starting point, and seeks a root. You will need to invoke your Horner's method file as part of the solution. Start at x = -4 and try to find a root of the polynomial

p(x) = x^3 - 2 * x^2 + 3 * x - 4

Check your answer! Include the file polynew.m when you send me your lab work. You can get a copy of polynew.m into your diary with the Matlab command type polynew.m .

Exercise 8: Newton's method is flexible in ways that bisection is not. Start at the complex value x = 2+3i and use your polynew routine to find a root of the same polynomial as above. This root is different from the one you found above. What is the third root of the polynomial?

Deflation

Now Newton's method gives us a way to find one root of a polynomial. But the thing about a polynomial of degree n is that we know that it could well have as many as n real roots. Supposing we've found one root of a polynomial, is there a systematic way to look for another?

Obviously, there are some very simple ways to try this. We could for instance, just pick a different starting point. If that still converged to the original root, we could keep trying many other starting points. But this is a disorganized way to proceed. Moreover, it could be that the polynomial has only one root, with multiplicity greater than 1. There would be no way to tell that this was so using this haphazard method.

Instead, we can try to be methodical by using an idea known as deflation. The idea of deflation is very simple: if p(x) has n roots, and we know one of them, r1, then to search for the others we might as well divide out the factor x-r1 and consider the simpler polynomial p2(x):

p2(x) = p(x) / ( x - r1 )

The polynomial p2(x) is well defined at x = r1 by invoking continuity, and has all the roots of p(x), except for r1. More precisely, the multiplicity of r1 has decreased by 1. Now we look for a root r2 of p2(x), and if we find it, we factor that out as well, and proceed to an even simpler polynomial p3(x) and so on.

Presumably, if our calculations are reasonably accurate, we could have a good chance of finding every real root of the original polynomial, and we would also account for their multiplicities. That is, if 3.5 is a root with multiplicity 4, we would have found that value 4 separate times.

Horner Factoring

Assuming that r1 is a root of p(x), we can use a variation of Horner's method to determine the polynomial factor p2(x) so that

p(x) = p2(x) * ( x - r1 ) + p(r1)

Horner's Method for Factoring with Remainder: Given the coefficients c of the polynomial p(x) determine the quantities b:

        b(1) = c(1)
        for i = 2:n-1
          b(i) = c(i) + r1 * b(i-1)
        end
        rem = c(n) + r1 * b(n-1)

Then the polynomial p2(x) is defined by the coefficients b and the remainder rem is the value p(r1). This algorithm is presented in Atkinson in Section 2.9 under the topic Nested multiplication. This algorithm is sometimes presented as "synthetic division."

Exercise 9 Working from the outline above, create the M file horner_factor.m, with the form

function [ b, rem ] = horner_factor ( c, r1 )

Before we move on to finding multiple roots, let's make sure this piece of code not only looks good, but works! We'll take a simple problem, and see if the code can pull out the roots as we specify them. Be sure to include the file horner_factor.m when you send your work to me.

Exercise 10:

Multiple Polynomial Roots

So let's see if we can sketch a method of getting some or all of the real roots of a polynomial p(x). Well, the first root r1 is easy, because we can just use Newton's method on p(x). But then we can use Horner's method to rewrite p(x) as

p(x) = p2(x) * ( x - r1 ) + p(r1)

Since p(r1) is presumably 0, or "close enough", we can assume that we can now concentrate on the polynomial p2(x). If we can pass the coefficients of the new polynomial to the Newton method, then we're all set for our next step!

Here is a sketch of the routine you need. 

      function roots = mulpolynew ( c )

      Set roots to an "empty" vector;
      Set n to the size of the coefficient vector.
      Pick a reasonable initial guess.

      Loop n-1 times.

        Call polynew to seek a root;

        If the function value at the new root is too high, then something's
        wrong, so print a warning and return now.

        Add the root to your roots vector;
       
        Call horner_factor to divide out the new root;

      end

Exercise 11: make a version of this routine, called mulpolynew.m, that works properly. You can check it out on the polynomial

x^2 - x - 6

Does your program fail on the polynomial

x^3 - 3*x^2 + 4*x - 12 ?

If it fails, why does it fail?

Exercise 12: Run your mulpolynew code for the two similar cases

c = [ 1, -5, -47, -39, 90, 0 ]
c = [ 1, -11,1,99,-90,0 ]

Do they both work or does one fail and the other work? Explain.

Back to the MATH2070 page.

Last revised on 30 September 2000.