This is a long, complicated project. We will allow 3 lab days for it.
There will be no lab on Monday, 11 October.
Polynomials occur so often in mathematical calculations that it is important to have a good idea of:
At the end of this lab is an assignment which must be turned in before Wednesday, 13 October.
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).
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:
cp = polyder ( c )We'll be writing our own versions of these routines today!
pp = polyval ( cp, x )
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 ];MATLAB responds with
c = poly [ v ];
c = [ 1, -6, 11, -6, 0 ]Exercise: Make a polynomial with roots -2, 1 + 2i, 1 - 2i, and 3. Assuming you do this correctly, your coefficient vector will be 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. (Isn't this impossible?) Try the commands:
c = [ 1, -6, 11, -6, 0 ]and see if you get back your original set of roots!
v = roots ( c )
Suppose someone asks you to evaluate the polynomial
p(x) = x^3 - 2 * x^2 - 5 * x + 6Lucky 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 ) - 4Then, 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 ) + constantthen 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. Remember synthetic division? This is the method by which you can simplify a polynomial fraction, in that beloved method of partial fractions from calculus.
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.
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. Remember that, by default, everything in MATLAB is "really" a two dimensional array. Your polynomial coefficient vector is actually a matrix, with 1 row and n columns. To ask MATLAB how many columns a data item c has, we write:
n = size ( c, 2 )(When we need the number of rows of something, of course, that will turn out to be m = size ( c, 1 )).
Class Discussion: 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 convince yourself of this, determine the number of multiplications required to evaluate
p(x) = 2 * x^4 + 7 * x^3 + 6 * x^2 + 8 * x + 12using the "natural" way to evaluate the polynomial. Now determine the number of multiplications involved in Horner's method. If p(x) were a polynomial of degree 20, how many multiplications would be involved?
Exercise: Working from the outline above, write an M file called horner.m of the form
function px = horner ( c, x )which 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 - 4which has the following values:
X: -1 0 1 2 3
P(X): -10 -4 -2 2 14
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:
pp = 0.0
px = c(1)
for i = 2:n
pp = pp * x + px
px = px * x + c(i)
end
Class Discussion: Actually, it may be easier to type in the change that computes the derivative, than it is to understand it. Can we see what's going on here? Let's write out the steps for a simple polynomial. It's worth talking about!
Exercise: Working from the outline above, modify your M file horner.m so it now has the form:
function [ px, pp ] = 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 - 4which has the following values:
X: -1 0 1 2 3
P(X): -10 -4 -2 2 14
P'(X): 10 3 2 7 18
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 such an algorithm would look like:
function root = polynew ( c, x )
FATOL = ?
STEPMAX = ?
for step = 1, STEPMAX
px = ?
pp = ?
if ( | px | <= FATOL )
...we're happy...
end
if ( pp is zero )
'Badness'
end
x = x - px / pp
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!
Programming Problem: your M file for Horner's method now returns two values at once, px and pp. What do you think is the right way to call your Horner function, and "catch" both values it wants to return?
Exercise: 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 - 4There are at least two ways to check your answer!
Exercise: Newton's method is flexible in ways that bisection is not. Start at x = 2+3i and use your polynew routine to find a root of the polynomial
p(x) = x^2 - 6 * x + 10
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 any 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 real 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.
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).
Exercise 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 good 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.
Exercise:
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, x )
Set roots to an "empty" vector;
Set n to the size of the coefficient vector.
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: make a version of this routine, called mulpolynew.m, that works properly. You can check it out on the polynomial
x^2 - x - 6Does your program also properly handle the polynomial
x^3 - 3*x^2 + 4*x - 12
Run your mulpolynew code as follows:
diary ex5and mail the resulting file ex5 to me.
c = [ 1, -5, -47, -39, 90, 0 ]
x = 5
mulpolynew ( c, x )
diary off