Project Two : Bairstow's Method
The goal of this project is to investigate a numerical
algorithm to find all roots of a polynomial equation: Bairstow's method.
Since this method involves synthetic division readily available in
MATLAB's deconv command, we will use MATLAB in our numerical
0. MATLAB and Polynomials
Polynomials in MATLAB are represented by their coefficient vectors.
For example, the cubic polynomial 3.3 x^3 - 2.23 x^2 - 5.1 x + 9.8
is entered in MATLAB typing c = [3.3 -2.23 -5.1 9.8]
at the prompt. With y = polyval(c,x), we evaluate the polynomial
at some point x. After typing r = roots(c), the vector r
will contain approximations for the roots of the polynomial.
We can group commands into .m files, defining functions.
The name of the function should match the file name. For example,
place the definition of the function ``newton'' in the file
``newton.m''. If the path is set correctly (do help path),
then you can call this function just as a regular MATLAB command.
See the .m file ``newton.m'' in
1. Study and implementation of Bairstow's method
Read section 1.8 in the text book. A good starter is to go through
example 1.3 in an interactive MATLAB session, using the command
deconv and some extra manipulations. The book suggests
Cramer's rule to solve Delta r and Delta s, but you may
also use the \ command. For example, to compute the
solution x to the linear system a x = b, we type
x = a \ b.
In a next step, we could write
one intermediate routine for splitting off one quadratic factor.
The prototype for this routine could be
function [q,f] = bairstow1(p,eps,N)
% Computes one quadratic factor of the polynomial p.
% On entry :
% p coefficient vector of polynomial, size(p,2) > 2;
% eps accuracy requirement on the coefficients of result;
% N maximal number of iterations.
% On return :
% q coefficient vector of the quotient: p = q*f
% f coefficient vector of a quadratic factor.
We use bairstow1 repeatedly until the quotient has degree two
or less. On this last quotient and on all factors we apply the
quadratic formula to find the roots. We could use the following prototype
function r = bairstow(p,eps,N)
% Computes all roots of the polynomial p with Bairstow's method.
% On entry :
% p coefficient vector of a polynomial;
% eps accuracy requirement on the quadratic factors;
% N maximal number of iterations to compute one factor.
% On return :
% r approximations for all the roots of p.
2. Numerical experiments on random and special polynomials
In principle, we cannot start with the second part of the project
without a correct implementation of Bairstow's method.
However, if you fail to complete the implementation and/or feel
you cannot trust your code, then you may resort to Program 1.2
in the text book.
In the following assignments we study the finding of the roots for
three classes of polynomials.
1. Random Polynomials.
We generate a coefficient vector c of length 6 typing
c = randn(1,6), defining a polynomial of degree 5.
Compare the output of your implementation with the output of the
roots command. Experiment with different levels of the
accuracy requirement eps, for example 1.E-4, 1.E-8, 1.E-12,
and 1.E-16. Make a table to show the relation between eps
and the errors on the roots.
Do you notice the effect of the propagation of errors caused
by the deflation? In particular, is the first root computed
more accurately than the last root?
2. Multiple Roots.
Type r = ones(1,5); followed by c = poly(r) to
generate the coefficient vector c of the polynomial (x-1)^5.
Apply the implementation of Newton's method used in Lecture 8
to answer the following questions:
In your answer to these questions, show the numerical experiments,
the code for the extended Newton's method you used and the numerical output.
State your conclusions briefly.
- Investigate how to show the geometric convergence rate
to the root 1 of multiplicity 5, starting sufficiently
close to 1 and using sufficiently many iterations.
- Extend the implementation of Newton's method to take
advantage of the multiplicity to restore the quadratic
3. A famous polynomial.
Consider the polynomial (x-1)(x-2) ... (x-20).
In MATLAB, we can generate its coefficient vector c by
w = diag(wilkinson(41)); r = w(1:20,:)'
(observe the prime) and c = poly(r).
Compare the output of the roots command and the output
of your implementation of Bairstow's method with the exact roots.
Explain why it is difficult to find the roots of this polynomial
with a numerical method.
3. Deadline is Wednesday 13 February, at 1PM
Bring your project solution to class. It should contain the following:
In addition, you could also hand in the log files created via the
diary command, eventually after some proper editing.
Typing diary followed by the name of a file in a MATLAB session
creates a new file with the given name which will contain everything
you see on the screen during the session. While this file created by
diary is by no means a substitute for the other items in the
project solution, it could help you to gain some (partial) credit when
some answers are (partially) wrong.
- The listing of the MATLAB code to implement Bairstow's method.
If you resorted to the C code in the text book and extended the
code, then you may get some partial credit for completing the code.
- Numerical output, with numbers formatted in scientific format.
- Answers to the questions in the assignments. Please write complete
grammatically correct sentences and avoid spelling mistakes.
If you have questions, comments, or difficulties, feel free to come
to my office for help.