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.
%
n = size(p,2);
r = rand;
s = rand;
for i = 1:N
d = [1 -r -s]; % divisor
[q1,r1] = deconv(p,d);
b0 = r1(n) + r*r1(n-1);
b1 = r1(n-1);
p2 = [q1 b1 b0];
[q2,r2] = deconv(p2,d);
c1 = r2(n-1); % define linear system
c2 = q2(n-2);
c3 = q2(n-3);
a = [c2 c3; c1 c2];
b = [-b1; -b0];
drs = a\b;
r = r + drs(1);
s = s + drs(2);
if (norm(drs) < eps)
break;
end;
end;
q = q1;
f = d;
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.
%
n = size(p,2);
q = p;
r = zeros(1,n-1);
for i = 1:floor(n/2)
if (size(q,2) <= 3)
f = q;
q = 0;
else
[q,f] = bairstow1(q,eps,N);
end
if size(f,2) == 2
r(2*i-1) = -f(2)/f(1);
else
disc = f(2)^2 - 4*f(1)*f(3);
r(2*i-1) = (-f(2) - sqrt(disc))/(2*f(1));
r(2*i) = (-f(2) + sqrt(disc))/(2*f(1));
end;
end;
r = r';
We generate a random polynomial of degree five and compute the roots with several different values of eps. Each time we list the error:
+------------------------------------------------------------+ | eps absolute errors on the roots | +============================================================+ | 1.0e-4 -1.591901586106559e-07 + 1.304010855962545e-07i | | -1.591901586106559e-07 - 1.304010855962545e-07i | | 2.924182164698585e-05 | | -1.446172066482676e-05 + 4.286444749590501e-06i | | -1.446172066482676e-05 - 4.286444749590501e-06i | +------------------------------------------------------------+ | 1.0e-8 -5.760947274779937e-13 + 5.523359547510154e-13i | | -5.760947274779937e-13 - 5.523359547510154e-13i | | 5.702327499079729e-12 | | -2.274957999759408e-12 - 8.339551271774326e-12i | | -2.274957999759408e-12 + 8.339551271774326e-12i | +------------------------------------------------------------+ | 1.0e-12 -1.110223024625157e-16 + 1.110223024625157e-16i | | -1.110223024625157e-16 - 1.110223024625157e-16i | | 1.665334536937735e-15 | | -6.106226635438361e-16 + 3.330669073875470e-16i | | -6.106226635438361e-16 - 3.330669073875470e-16i | +------------------------------------------------------------+ | 1.0e-16 -1.110223024625157e-16 + 2.220446049250313e-16i | | -1.110223024625157e-16 - 2.220446049250313e-16i | | 5.551115123125783e-16 | | -8.326672684688674e-17 + 1.110223024625157e-16i | | -8.326672684688674e-17 - 1.110223024625157e-16i | +------------------------------------------------------------+
The propagation of the errors is noticeable with eps = 1.0e-4 and less pronounced with eps = 1.0e-8. For the details we refer to the output of the MATLAB session:
>> format long e
>> p = randn(1,6)
p =
Columns 1 through 3
1.189164201652103e+00 -3.763327659331765e-02 3.272923614086541e-01
Columns 4 through 6
1.746391428209245e-01 -1.867085776814394e-01 7.257905482933027e-01
>> r = roots(p)
r =
-9.156617707168502e-01
-1.877222373959926e-01 + 9.256377302559300e-01i
-1.877222373959926e-01 - 9.256377302559300e-01i
6.613765374439257e-01 + 5.565965671504062e-01i
6.613765374439257e-01 - 5.565965671504062e-01i
>> rb1 = bairstow(p,1.0e-4,100)
rb1 =
6.613766966340843e-01 + 5.565966975514918e-01i
6.613766966340843e-01 - 5.565966975514918e-01i
-1.877077756753278e-01 + 9.256420167006796e-01i
-1.877077756753278e-01 - 9.256420167006796e-01i
-9.156910125384972e-01
>> err = sort(r) - sort(rb1)
err =
-1.591901586106559e-07 + 1.304010855962545e-07i
-1.591901586106559e-07 - 1.304010855962545e-07i
2.924182164698585e-05
-1.446172066482676e-05 + 4.286444749590501e-06i
-1.446172066482676e-05 - 4.286444749590501e-06i
>> sort(rb1)
ans =
6.613766966340843e-01 - 5.565966975514918e-01i
6.613766966340843e-01 + 5.565966975514918e-01i
-9.156910125384972e-01
-1.877077756753278e-01 - 9.256420167006796e-01i
-1.877077756753278e-01 + 9.256420167006796e-01i
>> % we see the largest errors at the roots that were computed later
>> rb2 = bairstow(p,1.0e-8,100)
rb2 =
6.613765374445018e-01 + 5.565965671509585e-01i
6.613765374445018e-01 - 5.565965671509585e-01i
-1.877222373937176e-01 + 9.256377302475904e-01i
-1.877222373937176e-01 - 9.256377302475904e-01i
-9.156617707225525e-01
>> err = sort(r) - sort(rb2)
err =
-5.760947274779937e-13 + 5.523359547510154e-13i
-5.760947274779937e-13 - 5.523359547510154e-13i
5.702327499079729e-12
-2.274957999759408e-12 - 8.339551271774326e-12i
-2.274957999759408e-12 + 8.339551271774326e-12i
>> sort(rb2)
ans =
6.613765374445018e-01 - 5.565965671509585e-01i
6.613765374445018e-01 + 5.565965671509585e-01i
-9.156617707225525e-01
-1.877222373937176e-01 - 9.256377302475904e-01i
-1.877222373937176e-01 + 9.256377302475904e-01i
>> % there is still a slight difference between the accuracy of the roots
>> rb3 = bairstow(p,1.0e-12,100)
rb3 =
6.613765374439258e-01 + 5.565965671504063e-01i
6.613765374439258e-01 - 5.565965671504063e-01i
-1.877222373959920e-01 + 9.256377302559303e-01i
-1.877222373959920e-01 - 9.256377302559303e-01i
-9.156617707168518e-01
>> err = sort(r) - sort(rb3)
err =
-1.110223024625157e-16 + 1.110223024625157e-16i
-1.110223024625157e-16 - 1.110223024625157e-16i
1.665334536937735e-15
-6.106226635438361e-16 + 3.330669073875470e-16i
-6.106226635438361e-16 - 3.330669073875470e-16i
>> % here we see no difference anymore in the errors
>> % due to the order of computation
>> rb4 = bairstow(p,1.0e-16,100)
rb4 =
6.613765374439258e-01 + 5.565965671504064e-01i
6.613765374439258e-01 - 5.565965671504064e-01i
-1.877222373959925e-01 + 9.256377302559301e-01i
-1.877222373959925e-01 - 9.256377302559301e-01i
-9.156617707168507e-01
>> err = sort(r) - sort(rb4)
err =
-1.110223024625157e-16 + 2.220446049250313e-16i
-1.110223024625157e-16 - 2.220446049250313e-16i
5.551115123125783e-16
-8.326672684688674e-17 + 1.110223024625157e-16i
-8.326672684688674e-17 - 1.110223024625157e-16i
>> % same conclusion as before
The code for Newton's was modified so that it computes the ratio of two consecutive errors which reveals (m-1)/m = 0.8, which is then used to estimate the multiplicity. First we show the output of the MATLAB session:
>> p = poly(ones(1,5))
p =
1 -5 10 -10 5 -1
>> [q,r] = newton1(p,0.9,eps,20)
|dx| |f(x)| |dx/previous_dx| m
2.00e-02 3.28e-06 2.00e-02 1.02e+00
1.60e-02 1.07e-06 8.00e-01 5.00e+00
1.28e-02 3.52e-07 8.00e-01 5.00e+00
1.02e-02 1.15e-07 8.00e-01 5.00e+00
8.19e-03 3.78e-08 8.00e-01 5.00e+00
6.55e-03 1.24e-08 8.00e-01 5.00e+00
5.24e-03 4.06e-09 8.00e-01 5.00e+00
4.19e-03 1.33e-09 8.00e-01 5.00e+00
3.36e-03 4.36e-10 8.00e-01 5.00e+00
2.68e-03 1.43e-10 8.00e-01 5.00e+00
2.15e-03 4.68e-11 8.00e-01 5.00e+00
1.72e-03 1.53e-11 8.00e-01 5.00e+00
1.37e-03 5.02e-12 8.00e-01 5.00e+00
1.10e-03 1.65e-12 8.00e-01 5.00e+00
8.79e-04 5.40e-13 8.00e-01 5.00e+00
7.05e-04 1.77e-13 8.01e-01 5.03e+00
5.64e-04 5.82e-14 8.00e-01 5.00e+00
4.54e-04 1.85e-14 8.05e-01 5.13e+00
3.56e-04 6.00e-15 7.85e-01 4.65e+00
2.79e-04 1.22e-15 7.83e-01 4.61e+00
failed requirements after 20 steps
q =
1.0000 -4.0012 6.0035 -4.0035 1.0012
r =
0.9988
The file newton1.m :
function [q,x] = newton(p,x,eps,N)
%
% Applies Newton's method using synthetic division
% on the polynomial p, starting at x.
% On return is the quotient and a new approximation for the root.
% The iteration stops until the accuracy eps is reached.
% No more than N steps are performed.
% This method displays estimates for the multiplicity m.
%
% Examples :
% p = [1 0 -1]; x = 0.9;
% [q,x] = newton(p,x,eps,4)
% p = [1 0 1]; x = 0.1 + 0.9*i;
% [q,x] = newton(p,x,eps,4)
%
n = size(p,2);
d = [1 -x];
[q,r] = deconv(p,d);
prev_dx = 1;
fprintf(' |dx| |f(x)| |dx/previous_dx| m \n');
for i = 1:N
[q1,r1] = deconv(q,d);
dx = r(n)/r1(n-1);
x = x - dx;
d(2) = -x;
[q,r] = deconv(p,d);
dx = abs(dx);
m = prev_dx/(prev_dx - dx);
fprintf('%8.2e %8.2e %12.2e %12.2e\n', dx, abs(r(n)), dx/prev_dx, m);
if (abs(dx) < eps) | (abs(r(n)) < eps)
fprintf('succeeded after %d steps\n', i);
return;
end
prev_dx = dx;
end
fprintf('failed requirements after %d steps\n', N);
If we use the multiplicity information in Newton's method, there is only one step:
>> p = poly(ones(1,5))
p =
1 -5 10 -10 5 -1
>> [q,r] = newton2(p,0.9,5,eps,4)
real(x) imag(x) dx f(x)
9.999999999993783e-01 0.000000000000000e+00 -2.00e-02 0.00e+00
succeeded after 1 steps
q =
1.0000 -4.0000 6.0000 -4.0000 1.0000
r =
1.0000
The file newton2.m :
function [q,x] = newton(p,x,m,eps,N)
%
% Applies Newton's method using synthetic division
% on the polynomial p, starting at x, with m as multiplicity.
% On return is the quotient and a new approximation for the root.
% The iteration stops until the accuracy eps is reached.
% No more than N steps are performed.
%
n = size(p,2);
d = [1 -x];
[q,r] = deconv(p,d);
fprintf(' real(x) imag(x) dx f(x) \n');
for i = 1:N
[q1,r1] = deconv(q,d);
dx = r(n)/r1(n-1);
x = x - m*dx;
d(2) = -x;
[q,r] = deconv(p,d);
fprintf('%.15e %.15e %.2e %.2e\n', real(x), imag(x), dx, r(n));
if (abs(dx) < eps) | (abs(r(n)) < eps)
fprintf('succeeded after %d steps\n', i);
return;
end
end
fprintf('failed requirements after %d steps\n', N);
The polynomial (x-1)*(x-2)* ... *(x-20) was used by Wilkinson to test the routines he was developing to find the roots. The bad behaviour of his routines on this polynomial are due to the nature of the problem.
We see that Bairstow's method fails terribly on the Wilkinson polynomial. This polynomial is very sensitive to changes in the input coefficients. The deflation process add errors to the coefficients in the quotient. What makes the implementation fail completely though is the bad scaling of the coefficients in the linear system we solve. The coefficients of the Wilkinson polynomial vary enormously in magnitude which leads to badly conditioned linear systems.
>> w = diag(wilkinson(41));
>> r = w(1:20,:)'
r =
Columns 1 through 13
20 19 18 17 16 15 14 13 12 11 10 9 8
Columns 14 through 20
7 6 5 4 3 2 1
>> c = poly(r)
c =
1.0e+19 *
Columns 1 through 7
0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000
Columns 8 through 14
-0.0000 0.0000 -0.0000 0.0001 -0.0010 0.0063 -0.0311
Columns 15 through 21
0.1207 -0.3600 0.8038 -1.2871 1.3804 -0.8753 0.2433
>> roots(c)
ans =
20.0003
18.9970
18.0118
16.9695
16.0509
14.9319
14.0684
12.9472
12.0345
10.9836
10.0063
8.9983
8.0003
7.0000
6.0000
5.0000
4.0000
3.0000
2.0000
1.0000
>> rb = bairstow(c,eps,100)
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.253588e-16.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 3.139568e-17.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 7.862779e-18.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.969186e-18.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 4.931684e-19.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.235108e-19.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 3.093241e-20.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 7.746814e-21.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.940136e-21.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 4.858939e-22.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.216888e-22.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 3.047613e-23.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 7.632540e-24.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.911518e-24.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 4.787265e-25.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.198938e-25.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 3.002658e-26.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 7.519953e-27.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.883321e-27.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 4.716648e-28.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.181252e-28.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 2.958366e-29.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 7.409026e-30.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.855540e-30.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 4.647073e-31.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.163828e-31.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 2.914727e-32.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 7.299735e-33.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.828169e-33.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 4.578524e-34.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.146660e-34.
> In /home/6/jan/Courses/MCS471/Project2/bairstow1.m at line 26
In /home/6/jan/Courses/MCS471/Project2/bairstow.m at line 18
rb =
NaN
NaN
NaN
NaN
NaN
NaN
NaN
NaN
NaN
NaN
NaN
NaN
NaN
NaN
NaN
NaN
NaN
NaN
NaN
NaN
Here we show the sensitivity of the roots of the Wilkinson polynomial:
>> w = diag(wilkinson(41)); >> r = w(1:20,:)'; >> c = poly(r); >> c(2) = c(2) + 1.0e-8; >> roots(c) ans = 19.8694 + 0.4832i 19.8694 - 0.4832i 17.8768 + 1.5260i 17.8768 - 1.5260i 15.4027 + 1.7511i 15.4027 - 1.7511i 13.1282 + 1.2783i 13.1282 - 1.2783i 11.2499 + 0.4671i 11.2499 - 0.4671i 9.9403 9.0060 8.0000 6.9999 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000
A change of 10^(-8) to one of the coefficients gives rise to complex roots.