- Email: [email protected]
Polynomial factorization and the Q-D algorithm
LINEAR
ALGEBRA
Polynomial L. FOX Oxford
AND
AND
ITS APPLICATIONS
Factorization LINDA
University
and the
445
Q-D Algorithm*
HAYES
Computing
Laboratory
Oxford, England Communicated
by Alan
J. Hoffman
INTRODUCTION
1. The purposes of this note are to give some details and numerical illustrations of a matrix factorization method for finding the zeros of real polynomials, to indicate the connections of this process with the Q-D algorithm, and to outline its advantages and disadvantages compared with other known methods for this problem. The basis of the method for quadratic and quartic polynomials appeared in Fox [l], and Bauer [2] gave a more rigorous general treatment. These accounts do not appear to be very well known, perhaps because there exist more powerful alternative methods for factorizing polynomials, but there may be problems for which the matrix factorization could be valuable, either as a method in its own right or as a preliminary to other well-known methods. As variants of the Q-D algorithm, moreover, we believe that these processes should find their place in the literature of numerical methods. 2. The account of Fox [l] had its motivation in the numerical solution of two-point boundary value problems in ordinary differential equations, and we retain this motivation in the present discussion. Its basis is the fact that for differential equations with constant coefficients the approximating algebraic problem is reduced to the solution of difference equations with constant coefficients. The general equation is given by @NY,+ %-1Yr+1 + * * * + soy,,,
= b
Y = 1,2,3,
. . .)
(1)
and with the support of initial and boundary conditions the problem reduces to the solution of the algebraic equations * Dedicated
to Professor
A. M. Ostrowski
Linear Copyright
0
Algebra
on his 76th birthday.
and Its Applications
1968 by American
Elsevier
1, 445-463
Publishing
Company,
(1908) Inc.
I.. FOX
44ri
Hert) .-I is a “I~tnd”
matrix
of band width
ISD
I,INl).\
H.\YES
N 1- 1, whose general
ro\\’
is 0, 0, . and J’ is the vector One attractive
.)
0, u,y,
N.y.
,’
with components
.
,
a,, 0, 0, .
0)
,
yi, y,, ~‘a, .
method of solving equations
(2) first performs the matrix
decomposition
where
L is a unit lower triangular
matrix.
If the decomposition
is, if every. leading matrices
submatrix
will have the general
In each general
k + 1 adjacent
and
banded
then
are k + 1 adjacent
and the corresponding
nonzero elements
triangular
“interchanges,” the
that
I. ant1 r.
forms
nonzero
elements
on
row of U will hnvc .V
on and above the diagonal,
k depending on the number of initial conditions equation
Ci an upper
without
of A is nonsingular,
row of L there
and below the diagonal,
matrix
is possible
associated
the number
with the general
(1).
The importance provided
of the method
the second
decomposition one coming
reaches
boundary a “steady
from the I, matrix
is that,
under certain
is far enough state”
away
conditions
from
and
the first,
the
which serves to give two factors,
and the other from the U matrix,
of the
polynomial /J\.(_X)= us 4 aR.__ix $- . . . J- aoN?
In \vhat follows we first give some detailed .Y -- 4, matrix
the
quartic
method
polynomial,
and
with the Q-D algorithm.
illustrate
((9
attention the
Subsequently
to the
connection
case
of the
we treat with less
POLYNOMIAL
ZEROS
AND
Q-D
447
ALGORITHM
detail the general polynomial, and give some numerical illustrations the application THE
QUARTIC,
of
of matrix factorization. GENERAL
FACTORIZATION
3. In the case of a quartic polynomial, the general relation (1) has the form U4Yr +
a3yr+1+
UzYr+2 +
%Yr+3
+
UoYr-l-4 =
b
(7)
The number of associated initial conditions may be 1, 2, 3, or 4, of which the last represents the pure initial value problem and has no particular interest in our present context. The relevant types of condition, roughly equivalent to the respective specification of first, second, and third derivatives for a fourth-order differential equation, can be represented by the initial equations "1Yl
+
PlYI+
YlYl
+
=
6 1’
=
62,
Y4Y4 =
63.
E2Y2 P2Y2 +
P3Y3
Y2Y2 +
Y3Y3 +
(8)
With all of (8) given, our matrix and its decomposition look like A \
a3
a2
a1
a0
a4
a3
a2
a1
a0
L
u
(9) %
Linear
Algebra
and
Its Applications
.
%
1, 445-463
(1968)
.-1
-4
-
1
Y3
7’1 Y2
Ys
?'A
a4
a2
a1
a3
6111 a3
a0
fl2 ...
and with the specification
only of the last of (8) vve ha\~
.I Y2
Y3
Y4
u3
a2
al
(14
a3
%
(1, ...
a,
..*
3 iP3
;;4
7'4
1. 1.
It is now useful to consider
A4~t==-0 of the relevant is at infinitv
so that the general
I. is not singular
and this means homogeneous I.lnrar
Algrhra
algebraic
that
successive
md
the
Its Appl~cnlims
of the homogeneous
form
in which the second boundar!,
row of A is repeated
in (4), it follows
form of (l),
the solution
equations,
indefinitely.
Since
in addition
to the
that
values
of J’? satisfy,
“reduced” 1, 445-46s
recurrence (196X)
relations
POLYNOMIAL
ZEROS
““Y” +
Q-D
AND
V,Y,+1
01
=
449
ALGORITHM
“rY, + VrYr+l+
W*Yr+2 = 0, (13)
%Y” + vryr+1+
relevant
respectively
wry,+2
+
GY,+3
to the decompositions
=
01
(9), (lo),
Now if the zeros x1, x2, x3, and xq of the relevant distinct,
the general
solution
and the initial the
arbitrary
conditions arbitrary
constants
almost certainly then
of largest
reach a “steady
(i)
(8) merely in (14).
(14)
give one or more relations Provided
that
distinct
as a result of inevitable
modulus
dominate
in (la),
large r the decompositions
state.”
Provided
moduli,
Decomposition
real “dominant”
A,%,'+ Qql,
are not zero (and if they are exactly
for some sufficiently pairs have
of type
constants
be “introduced”
the roots
(6) are
of (7), with b, = 0, is
y, = ApI'+ A&'+ between
and (11). polynomial
that
If lxil > lxzl >
zero, then
rounding errors), and we infer that
(9), (lo),
and (11) will
any real roots and any complex
we find the following
(9).
the relevant zero they will
for sufficiently
results.
1x3( 3
1~~1, so that
xi is a
large Y we have
U” + v,xi = 0, and the elements (ii)
u, and v, converge
Decomposition
are two “dominant” Y we have
obtained
If
(10). roots,
(iii)
to constant
lx11 3
the quadratic
Decomposition
(11).
If 1x1( >
large
to fixed
lx2[ >
(16)
numbers.
[~a/> 1x11, so that there is the cubic factor
with zeros xi, x2, and x3,
and u,, v,, w,, and Z, all converge
Linear
there
then we have obtained
u, + v,x + wp2 + z$,
In all cases the other factor
1~~1, so that
for sufficiently
with zeros xi and .x2.
modulus,
but is lurking in the L matrix.
then
factor
uI, v~, and w, all converge
one real root of smallest
u and v, respectively.
Ixsl > lxsj >
real or complex,
24, + v,x + wrx2> The elements
(15)
to fixed
of the quartic
numbers. polynomial
The rules for matrix .4lgebrn
and
(17)
is not “lost,”
multiplication
Its Applications
1, 445-
463
easily (1968)
show
that
corresponding
to (15),
(Iti),
and
(IS)
the
other
factors
are
respecti\+-
and the respective
.5 .
I,, nz,, and 12, con\.erge to fixed uumbvrs.
The most interesting
quadratic
factors,
qualifications
already
even including
and useful case is that of decomposition
represented
by the
mentioned
the situation
decomposition
we can “solve”
not present
that is, the modulus of a real or comples
(10).
the quartic
in (le),
into
\Vith the completely,
the case in which
pair svparatcs
the rnoduli of two
real roots. In this case the factors so that numbers
we cannot
(w,, of course,
does dominate
everything,
can show that,
to a,,).
U matrix,
of the “middle”
real root in the
hi1 -L
.y4 is lurking
“steady m,x
-t
state”
in I,~,~L and m,.
pair
On the other
that
for sufficiently
tllat
in fact
to fixrtl hand
xl
large Y
.y2 and .I:~.
in the L matrix,
and in fact \v(
situation,
x2 = (x4 - x)
where x4(“) will also help to determine suffixes
contributions,
large Y. Here x1 I”)depends on Y, but its I-alue is important
for the computation The smaller
equal
and this implies in the
c~lual
II, and 11, in (10) will con\crgc
is always
the zero ,cI is contained
for all sufficiently
.y2” and s3” in (14) makr,
e.xpect that
(.r4(r’-
the other pair.
X))
(23)
We note the different
POLYNOMIAL
6.
ZEROS
451
AND Q-D ALGORITHM
The result (23) is perhaps easier to verify rather than to prove,
but in the process of the verification we can relate the xi(‘) and the XL’) to the other roots x2 and xs. a4 + agx + a,#
We have
+ alx3 + aox = a&x, -
so that A = x2 + x3, B = x2x3. a4 = a,x,x,B,
a3 = - a,{Ax,x,
The other relevant
equations,
x)(B - Ax + x2),
(24)
Then a, = -
a2 = a,{x,x,
x)(x4 -
+
+ B t
a,(xl + x4 + A), w,
+
X4,>>
Ah
+
~4) 1.
coming from the matrix
(25)
decomposition
(lo), are
@r-2 =
lIv,_2+
a4,
&a, + m,v,_, + u, = a2, Equations
m,uu,_l
=
m,a, + v, =
a3,
al,
w, = ao.
(26)
(22) and (23) give Qr,
= x1x1(I),
1r+l =
x4x4
(I),
- v,/a, = x1 + xl(“), (27) -
m, =
x4 +
x4@).
From (26), (27), and the first two of (25) we easily find B = x4(‘)x,(‘-i)
A = x4(y)+ xi(‘)
(28)
and substitution in the second pair of (25) with further use of (26) verifies that everything is satisfied. Equation (28) gives the required formulas for determining the middle pair from the “steady state” situation. EXAMPLES
7. In each of the following three examples we take the first two rows of the matrix A in (10) to be respectively (a,, a,, a,) and (a3, a2, a,, ao). For the first example [3] we consider the quartic p,(x) = 1 - 32x + 160x2 - 256x3 + 128x4, which has four real zeros of distinct moduli. we obtain in 24 steps the quadratic factors Linear
Algebra
(29)
To six significant figures
and Its Applications
1, 445-463
(1968)
0.0117476 and there correct
0346719x
is no further
after
For
-
85.1237
-
211.62OOm $- ldW,
change to this precision.
and in 1X steps
with two complex
we obtain
0.489555
-
In our third
+ AA,
are correct
the second factor,
Some eleven figures are
=
[I]
six-figure 2.04265
after
about
we take
the cluartic
decompositiorr
-- 2.6853i.r
35 steps.
In both casts,
the modulus
which
From
5.9375x" - 4.r” i- x4,
I --- 4.r -t-
three
From
(33)
successive
For
appearance
from
finding of Table
I.
of two real \ve get _ra =
.-\t about
this stage
zeros of the U and IX quadratics
0.536864
0.925692
~~.8%4308
0.669727
1 .ok-N273 of (24) the three
the correct
zeros of a quartic
are
(34)
estimates
value being obvious,
(28) the “correct”
the
implying
_~r= 1.640388,
and from the L matrix
121313x
for the .l
and 1.750000,
the moduli
to the figures gi\ren.
pairs of the other
(28) we deduce
of (24) we have
pair separates
we find from the L’ matrix
to this precision,
which is also correct
1.750000,
8.
of a complex
Y = 27 onward
is correct
0.609612,
of
the quartic
we find that the I,, ML,,u,, and U, do not converge to fixed numbers, roots.
(32)
+- .\.2.
with the larger zeros, belongs to the I: matrix.
example, b4(x)
pairs
the correct
1.31463x
.\gain eleven figures
that
(SO)
35 steps.
an example
course,
_t ~2,
two estimates
I.750002,
and for the
14
of unit\,.
the Q-D array
~3] has
the
POLYNOMIAL
TABLE
AND
Q-D
ALGORITHM
453
I
q,(l)
-
ZEROS
#)
q*w
al
e,(V
0
q,(3)
e,(3)
qr(4)
0
0
a0 0
%
-a3
a4
al
%
a3
q&V
qP
q!f)l
e,!l)
0
0
e!?;
0
qow
q!!;
El(%)
t?,(l) qP
q’“‘,
e,@J ql(“)
qP
0
qp
qP
0
&?,(3)
qllC4)
... The numbers
in the array are generated W’ =
%+1
according
to the “rhombus
rules”
(k’ _ eSf+;C + q,@‘,
e 7
(35)
ejyl = (q,(k+l)/q;$‘l)ep. If the zeros of the quadratic 14
> I%1 > 14
in modulus,
with jxr/ >
then
g(i) I -*Xi)
q(Z) I +x2,
and the e columns converge with
are real and distinct
to zero.
lx11 = lxs\ > [x3\ = 1x41, then B, -
qCS) I +xs,
lj ICi) + * 4p
(36)
If the zeros form two complex th ese belong
A,% + X2,
B, -
to quadratic
pairs,
factors
A,x f x2,
(37)
and
The
case of a complex
pair with modulus
of two real roots has obvious of Eq.
behavior,
(21), to which we have given
LU method.
For this
p (1) --f Xi’ I
greater
than
the moduli
and the other possibility some detailed
attention
is that with the
case we have
@) + qj:S’ ---f‘4, Q7tl Linear
Algebra
ql’d’q;“) + B, md
Ifs
.-lpplicatims
qj”’ + x4. 1, 445-463
(39) (1968)
1,. ITOX .\SI)
454 !I.
Sow
giv~ii
if the
in the
first
that
\vith exact
tion
(lo),
Lli
we should
equations
and
from
(39) \\.v cltduc~,
l‘he
other
roots
of thaw,
10. the
Q-D
unique,
the
coefficients
\vork,
and
perhaps
of rounding
error
avoids
bitrarily,
the
segment It Table
I must
For
w
qy)::
denotetl
in the special
satisf!.
a11~1
ant1
results
if an\-
li.
.I~'"'
(2X)
numcricnl
precision of the
ro\vs
will
usually
to Iiialz
an
polynomial,
LISIWI
“c~tlrck”
of
I i:,
of thcl coefficients ;I new
new
rxxtra
;I minimum The
coefficients. but
ul, that
~)ol\~nomial nirans
to ensure’
suitabl!-,
of tlie matriws
form<9
of Table
0iwurst’,
This,
hi- t.hoosing
tliv
descriptions
it is rc~commentli~cl
to produ(,(,
defect.
that
al-I-ang-ement
(6) is LW).
computation
quartic the
clsist
rewals
published
thr
he shifted
csstra
matrix
possible
that
not
hal.cb this
difficult\.
the
5
111. .I,“’
of (40) , anti &r
tlw
First,
does
should
with
unit
is in fact
and
pol\nomial
do not
first
\\hich
qLfl
13 j) implv
method
this
of the
algorithm.
to tllcs tl(ac,omy)osi-
(39).
(e.g.,
(42, or ug of the quartic in this cast‘ the origin whose
the
advantages.
algorithm
method
I-OI~~ 2s to slio\\~
large: Y, tlkal
of the l,I- ant1 (2-D mr~thocls
significant
so that
for 5ufficirntl~~
from
(‘omparison
some
t\vo
difficult
find, corresponding
cluadratics,
5, are clearly
immecliatel!-
has
tlici tirst
it is not
give) an ~a>\. proof of tilt> rc5iilts oi Sc-rtions i? :; from (40) \vk~ fintl qy)? and y,’
if we eliminate
follow
with
7, then
~l.\S155
identities
l~or
in Section
is p~~rformed
of Section
arithmetic
the
‘l‘ticw
method
sentence
l.lSl)\
I,(.
othtwvise
.i in (!1)--(I I ),
al-=
(Ilearl!-
2
be satisfactory. “arhitrq-”
start
for c,sample, quation
\vith the
the
iiumhcrs
Q-U in
POLYNOMIAL
ZEROS
AND
Q-D
455
ALGORITHM
and we can also establish the identities
Together with the rhombus rules, it turns out that we can start the Q-D table, and build it up to the last two rows given explicitly in Table I, with arbitrary values of q,,(l), e,,(l), and %(‘). The remaining elements in the relevant triangle come successively from the equations ql(l) = qow + eo(l!, qp
= -
elm =
qlN + e,(l),
q2(l) =
[((u&p + %)/q,(1)+ %Nq2(1) + %I/~,~
Ql(2) = q2(l)e2(l)/elU), (&(P) =
e,(‘)q,(2)/q,(‘),
e,(2) =
(@ az...1__ + a-3 -~~ Go92(l)
a, ql(l)q2(1)ql(2)
ql’“) = - ~,/~, - qp
q3C1) - q2C1),
q1t2)+ elil) - q,c’r),
+ q1t2)) +a4
el@) = q2(2)+ ez(l) _ qlt’)’
e2(li =
(45) 40(2)41(2) + qo(“)q2w + q1U)q2(l) qo~l~q~~l~qo~~~q~~~~q~~~~~ ’ a, i
1
40(3)= e,(2)q,P)/eo(2),
- q2r4 - qow,
eoW)=
q1
q.
(3) +
(4) =
$p&$3J
el(2) -
poP),
and we can then continue by the standard method. The corresponding results for the general case, of course, are not very attractive as the basis for a general starting method. 11. Second, both methods can subsequently fail through the vanishing of a denominator, for example, an early &$i in the Q-D equations (35), or an early U, in the LU decomposition. (Clearly, if convergence takes place, U, tends to a nonzero value). Whereas the Q-D failure necessitates a new start with a shifted polynomial or with equations like (45), we can with the LU method replace the offending SC,by an arbitrary nonzero value, say unity, and proceed with the iterative process. The effect can be regarded as that due to some different set of initial conditions, that is, of different initial rows of A, or the result of an isolated “blunder,” neither of which can do more than delay the subsequent convergence. Consider, for example, the polynomial p,(x) = 13 + x + 2x2 - X3 + x4. Linear
Algebra
and
Its Af$lications
(46) 1, 445-463
(1968)
W’ith the “Q-D start” this Ita bJ* unity
of the LI’ method we find ~1~=m0, but if \VC’replacc~
in the LU method
quadratic
factors,
50 steps.
The Q-D algorithm
a complete 12.
restart
Third,
uccww~lntion LU
a six-figure
of rounding
method.
the L and
error,
For suppose 7’ matrices
is reported
to
to obtain in about 0,
qrf’i from
sulfu
this is \irtuallv
we ha\e reached
reached
\Ve can then consider restarting
the corresponding
whereas
that
have
happily
ant1
necessar!‘.
algorithm
Q-D
quite
state result being reached
produces
is apparently
the
we proceed
steady
a steady
some
state
the whole operation
;t slo\\~
abscw t in tllcs stage
at which
to some prrcision.
with the factorization
I,
-1
!. 0
;_
. .
where if, ~1,and ZI’are the stead!, state elements ‘This is simply equivalent representing
the fact
in (14) are already efficients
to starting
that
quite
with
in the LU method
polynomial,
whereas
only in the first Moreover,
those
start
the
One might
are “refreshed”
relevant
say that b!- frequent
of the Q-D method
make
conditions,”
coefficients
Li,
the successive
co-
contact
with the
explicit
contact
few steps.
if we kno\v something
we can use this information start,”
this
small.
of tlir, previous operations.
with very good “initial
about
the distribution
quite deliberately
for which the unwanted
coefficients
It is only fair to add, however, tried this device has not achieved
of thti zeros
to make the “best possible A, in (14) are already
that in many of the examples a spectacular
improvement
small.
WC have in con\-cr-
gence ! Finall!.,
13. properties factor
the LO7 method
belonging
to the I,’ matrix
the rate of convergence of the IJ factor method, of similar
has some slightly
from those of the Q-D algorithm.
modulus.
con\ergencc
has zeros of \-ei-h.similar modulus,
depends 0111~7on the ratio of the smallest
to that of the largest modulus
on the other
different
It does not matter
hand,
has some
trouble
of the L factor. in srpar-sting
if the and
modulus The
Q-D
an>. roots
POLYNOMIAL
ZEROS
AND
Q-II
457
ALGORITHM
For example, with 25 - 56x + 38x2 - 8x3 + x4 = (x -
1)2(x2 - 6x + 25),
(48)
the LU method separates the two relevant quadratic pairs very rapidly. We find I,, = 1.00000058,
ml2 = - 2.00000024, 752 =
-
uu12= 2.49999993, (49)
5.99999976,
which are very accurate, and at Y = 20 these coefficients are all correct to 13 figures. With the Q-D method the individual q/“) and q,‘“), representing the two equal real roots, converge very slowly, being correct to only 1% at about Y = 70. In accordance with Eq. (40), of course, relevant combinations of the qr(k) “settle” at the same time as the L and U coefficients, and we might conclude that the corresponding sum and product checks could be used with advantage even in the case of real roots of nearly equal modulus. Similar results are obtained for the quartic with zeros 5, - 5, 1, and 1. The LU process gives quadratic factors with coefficients correct to five figures at Y = 9, and to 12 figures at Y = 19. In the Q-D process the two equal smaller roots are correct to 1% only at Y = 100, and the columns for the larger roots do not appear to be converging. But at r = 9 we find qJi’ = 13.000007,
qjy
qi4), = 0.899160,
= -
13.000008,
qpil = 1.923077,
qyll = 1.091602,
qPJ2 = 1.100841, q,@) = -
1.923076,
(50)
q;y2 = 0.908397,
and Eqs. (38) give very nearly the correct quadratic pairs x2 - 25 and X2 -2x+ 1. The case of equal real roots belonging to the U matrix is determined quite easily by considering, instead of (14), the relevant general solution y* =
(A, +rA,)x,*+
A,%,'+
of our basic recurrence relation.
A*X4",
I%/ >
(51)
IX3b4~
Here we have, for sufficiently
large
7, Linear
Algebra
and
Ifs
Applicatims
1, 445-
463
(1968)
and we get convergence case of equal 14. have
example,
each
factor
is ultimately
zero in the
roots!
The difficult
equal
because
problems
modulus.
When
arise when at least three of the four roots all the roots
ha\-e the same
modulus,
for
in an!. of the cases
(.53)
they cannot, there
so to speak,
decide which matrix
is no easy way of deducing
(53) gives some hope, because converges,
but
only very
\Vhen only three other
real root.
of which the last factor no such problem state”
factors.
The
and
third
of
of the term rA, in (Til), and this in fact
slowl!-.
roots ha\-e equal modulus
For
with this modulus
they want to inhabit,
the quadratic
example,
out thr
for the quartic
has complex
have the dilemma
and attaches
we can separate
roots of modulus mentioned,
6, the three roots
but the small root has
itself firmly to the L mxtris.
The “stead\.
is gi\ven bTV
where .Y~is the lone small root. we can also show that the steady
Hy methods
similar to those of Section
state remaining
cubic equation
6
is given
b!,
In this results
particular
example
we find
after
a few steps
the se\-en-figure
POLYNOMIAL
ZEROS
AND
Q-D ALGORITHM
1,, + mrex + x2 = 4.187036
-
5.187036x
us = 29.854055,
xq(i’) = 4.187036,
u i. = 17.996859, and our computed x3.
vr,, = -
+ x2,
x4=
zig= -
8.837551,
of course,
separate
125 -
out the roots
55x + 11x2 -
given in the first part of this section.
of equal
modulus
Eq.
and the “reduced
xl(q
+2
POLYNOMIALS
(55) is replaced
cubic”
corresponding
x(.+?z,
+
OF HIGHER
t1
by
by
x2(x1(*) - nz,) - x3 = 0.
(59)
DEGREE
Many of these results can doubtless be extended
ing treatment
(58)
to (56) is given
Z,+,) +
-
with three smaller
x)(Q) - 4,
w,,x2 = (x1 -
v,x +
21, +
15.
(57)
of this cubic for the
In the case when there is a single real root associated roots
1,
6.812964,
cubic is very close to the correct
We cannot,
reasons
459
of polynomials
of other
to the correspond-
and higher degrees.
In any case
we have seen that any splitting is possible provided that all zeros belonging to the U matrix
have moduli greater
than those of the L matrix.
moduli occur only in complex roots, therefore, which will converge of L and
U settle
The best in problems
in the absolute
sense, that is, in which the elements
down to constant
application
values.
of the suggested
splitting
(or smallest) so that
here have
modulus
the certainty,
the “reduced”
lacking
16.
To
test
more the
by Wilkinson
in Newton
method
on fairly
if we seek the k zeros
the
appropriate
the relevant or Bairstow
initial
k zeros.
We
methods,
that
of the distribution
(fastest)
the required
from Olver
is probably of the zeros
and we can then concentrate
the best
or less than
worked on three examples studied
includes
Some knowledge
can help to achieve
if this includes
however, the nature
For example,
our objective,
polynomial.
about
we can make
U (or L) certainly
we have in fact achieved moreover,
method,
in which we know something
and in which we need specific results. of largest
If equal
there is always some splitting
difficult
initial
on
of zeros,
splitting,
even
k zeros. polynomials
we have
[a], two of which have also been
[5]. Linear
Algebra
a?zd Its Applications
1, 445-463
(1968)
The
polynomial
1 + X.Y + 32.8”
t 89.6x” ~,-190.68.x4 -I- 304.08x5
+ 443.,576YG +- 468.88.T;
$mFi24.327,~~ + 378.908.~~ -+-34Fi.072FiB~~~-+ 166.44768~~’ -1. 37,6lilO96x’3
+ 26.1783048.~‘~
11as approximate 0.293.51 f
+ 3.4356048~~‘~ + 2.032.53121 x1’
0.4FiO93i, ~ 0.14762 + O.i7176i,
0.09004 + l.O6119i,
0.0608~ i + 1.2969li,
0.01049 + I .5963Oi, -
(~.00249 -1.: 1.66712i.
that
the
zeros)
order of modulus.
fastest
in which
to the L matrix.
(into
This is (partly)
with
numbers
verified
(the smallest
places in 100 iterations,
L the corresponding
of thts moduli would suggest
factors
an e\~n number
successivelv~ 2, 4, (i, 8, 10, 12, and
\Yith six zeros in the L matrix decimal
~~ 0.02,567 + 1.474383’, (61)
Inspection
decompositions
are those
belong
(60)
zeros
O.l435Oi, ~~~0.22447 f
in ascending
-t 128.218748~~~
of
14 zeros
by our computations.
we tried) we converge
to IO
and with 8, IO, 12, and 14 zeros in the
of correct
decimal places after 100 iterations
are 7, 4, 3, and 2 respectivel\-. To test polynomial
the
speed
and accuracy
we performed
each of the two factors, iteration
was terminated
eight decimal places. results
produced
the (8, 8) splitting, when all relevant
The computed
hv. \I’ilkinson
0.293504?53 & 0.1434993oi3,
-
0.14762378
l 0,7717fi71Ri,
(-
0.14762378
4 0.7717fi72Oi),
-- O.OFiO86440+ 1,2969113i, (-
O.OFiO86444+ 1,2969113i),
--~ 0.01049254
+ 1,5962951i,
(-- O.O10493,5Fi:-i 1,,59629EiAi), I.iwrtrr,Algchvcr
ctnd
method
T/s .?~~licc~linm
coefficients
all the zeros.
aw
& 0,45092796i,
0.22445006 0.09003999
+ 0.4fiO92796i), + 1.0611!~21 i,
0.09003999
% l.O61192li),
~~~0.02t56687.5 3: 1.4743771i, (-
0.02566871
~-- 0.00249022 (1,
of
Each
had convrerged to
in parentheses,
-- 0.22447006
(-
particular
values of the zeros, and the accurate
[fi], given
(-
for this
then the (4, 4) splitting
and so on until we obtained
0.293,504,53 * 0.1434992&‘, (-
of the
+ 1,66712OXi,
0.00248920 44S-463
& 1.474377 li).
+ 1.6671204~‘)
(1968)
(62)
POItYNOMIAL
ZEROS
AND Q-D ALGORITHM
461
The complete computation took 28 seconds on the KDFO machine with an ALGOL compiler. Some of the zeros are fairly close together in modulus, and the results obtained are really quite satisfactory. 17.
The polynomial
420 + 1882x + 4090x2 + 5276x3 + 5644x4 + 1019x5 + 3651x6 - 4900x’ + 3912x8 - 4618x9 + 3086~~~ -
1916+
+ 1056~‘~ - 409x13 + 139x14
_ 36x15 + 4%” has real zeros 34, 3, 24, 2 and twelve complex zeros with moduli < 1.94. This is a rather more ill-conditioned problem. The splitting into two polynomials of degree eight is quite fast, but one of the next splittings (of the U matrix) is relatively slow, and the whole operation took 53 seconds. To seven decimal places we find the zeros 3.5000600,
3.OOOOOO0,2.5000600,
2.0000000,
- 0.0112583 & 1.9332891i,
- 0.0449754 & 1.7381625i,
- 0.1012792 f
- 0.1812607 f
1.429775&,
- 0.2844601 f 0.5987694i,
l.O348308i,
(64)
- 0.3767663 + 0.1888406i.
The reason for the slow convergence of the second U splitting is quite obvious, the first complex root having modulus very near to that of the nearest real root. The maximum error in (64) nowhere exceeds a few units in the seventh decimal. 18.
Our final example, the determination of the zeros of the poly-
nomial 2 + 16x + 224x2 + 1288x3 + 10416x4 + 44184x5 + 267232x6 + 837860x’ + 41782600~s+ 9490840x9 + 41018752~~~ + 64249356~~~ + 247926664xr2 + 240775148~~~ + 845947696x14 + 385455882~~~ + 1250162561+,
(65)
represents an extremely ill-conditioned problem. All the zeros are complex, six pairs have very similar modulus, and the last four pairs have almost equal modulus. Linear
Algebra and Its Afiplications
1, 445- 463 (1968)
For
the first
To test
time
convergence
is here affected
this we ran the program
after 500, 1000, and 2600 steps, of which consumed to five decimal
about
places
results
of our method.
‘I‘.4BLE
11
Correct results
three
times,
respectively,
I10 seconds
the correct
errors.
each splitting
in the three runs, the first
of machine
results
500 iterations
by rounding
truncating time.
Table
[Fi] and the three
1000 itcratwns
II shows computed
2.WO iterations
-0.13245~_0.13601i -0.13246~0.13601i
-0.13245~0.136001
-0.13245 :,O.l359Hi
-0.01869~0.25305i
-0.01870~0.25304i
-0.01870~0.25305i
--0.01869&0.25304i
-O.O0232&0,29258i
-0.00232&0.29256i
-0.00233~0.29:'38~ -0.00233 I 0.29258i
-0.00049+0.30418i
-0.000301_0.30414i -0.00029&0.30418i
-0.00014~0.308612' -0.00146&0.30785i -0.00005~0.31066i -0.00000~0.31220i
+0.00182~0.3099Oi
-0.00161 4 0.31201; -0.00168- 0.31222i
+0.00060~0.312623' +0.00071~0.31319i
The results are interesting,
t~O.OO0811:~0.313ORi
in that there is no significant
with more than 500 iterations. mined
-0.00136_t0.30800i
+0.00182 t 0.31011i +0.00179&0.30998i
--0.00001~0.31170i -0.00134~0.31264i
rounding
0.00026~0.30418i
-0.00130f0.30799i
improvement
This is not due, to any accumulation
of
error, but to the fact that the later roots ~I-C \.er>r poorly deter-
(ill-conditioned)
with respect
The effect of individual can be interpreted, to factorize .1t the next
to small changes
step
these
situation
in the coefficients.
errors in the computed
in the spirit of backward
the given
ill-conditioned
rounding polynomial
with
perturbations
matrix
error analysis,
slightly
perturbed
are slightly
this has a serious
effect
elements
as an attempt coefficients.
different,
and in
on the zeros,
an
and also
on the rate and even the existence
of con\rergence to more than a definite
number
can br alleviated
of figures.
of floating-point scalar
products
this matrix because
method
the device
is apparently
that Linew
is particularly is not included
rather
far more
however,
iteration,
which permits
prior to single-precision
It is perhaps result,
This situation arithmetic
surprising
with double-precision
the last real part .4lgebvn
and Ifs
rounding,
well suited. that
Applicutiors
results.
and treble-precision
had only the figures 1, 445-463
ALGOI,
compiler.)
zeros the real part
the imaginary
by Wilkinson’s
of
(It \vas not used here
in the later than
accumulation
and for this purpose
in the standard
ill-conditioned
is also reflected
by using the variant
the accurate
arithmetic,
-- 0.00000305 (196x1
part.
This
Using Hairstow he found
in common,
POLkCNOMIAL
ZEROS
AND
463
Q-D ALGORITHM
while the imaginary part was the same at 0.312196968, more decimal place and to five more significant
that is, to one
figures.
CONCLUSION
19. We have demonstrated that the method of factorizing a polynomial in terms of a certain matrix decomposition leads to a quite practicable technique for determining some or all of the zeros of the polynomial. In particular we can factorize the polynomial directly into two factors which contain specified numbers of the roots of smallest or largest modulus, and can then concentrate on the part which may be relevant to our particular requirements. Moreover much information is gained at all stages, by observation of the rate of convergence, of the existence of roots of similar modulus, and of their precise positions in the ordering of the zeros. This information is not easily obtained in the early stages of other methods of iteration, such as those of Bairstow or Muller. On the other hand, the Bairstow iteration is probably quite a bit faster, in ill-conditioned cases, than the methods of this paper carried One might, therefore, decide that the to the limits of their accuracy. advantages mentioned in the last paragraph would make this method a useful preliminary to a Bairstow-type “cleaning-up” operation. Various questions require further consideration, such as a satisfactory error analysis and a comprehensive program which would do all that anybody might want to do, and we hope to make some subsequent progress on these and other allied matters. REFERENCES 1 L. Fox, Numerical Solution of Boundary-Value
Problems, Clarendon Press, Oxford,
1957. 2 F. L. Bauer, wiss.
3 P. Henrici, 4 F. W.
Faktorisierung
Elements of Numerical
J. Olver,
244(1952),
eines Polynoms,
Sitz. Ber.
Bayer.
Akad.
Evaluation
Analysis.
Wiley,
New York,
of zeros of high degree polynomials,
1964. Phil.
Trans.
385-415.
5 J. Wilkinson, Num.
Direkte
1966, 163.
The evaluation
Math. 1(1959),
of the zeros of ill-conditioned
polynomials,
Part II,
167-180.
Received Mavch 27, 1968
Linear Algebra and Its Applications
1, 445-463
(1968)
ALGEBRA
Polynomial L. FOX Oxford
AND
AND
ITS APPLICATIONS
Factorization LINDA
University
and the
445
Q-D Algorithm*
HAYES
Computing
Laboratory
Oxford, England Communicated
by Alan
J. Hoffman
INTRODUCTION
1. The purposes of this note are to give some details and numerical illustrations of a matrix factorization method for finding the zeros of real polynomials, to indicate the connections of this process with the Q-D algorithm, and to outline its advantages and disadvantages compared with other known methods for this problem. The basis of the method for quadratic and quartic polynomials appeared in Fox [l], and Bauer [2] gave a more rigorous general treatment. These accounts do not appear to be very well known, perhaps because there exist more powerful alternative methods for factorizing polynomials, but there may be problems for which the matrix factorization could be valuable, either as a method in its own right or as a preliminary to other well-known methods. As variants of the Q-D algorithm, moreover, we believe that these processes should find their place in the literature of numerical methods. 2. The account of Fox [l] had its motivation in the numerical solution of two-point boundary value problems in ordinary differential equations, and we retain this motivation in the present discussion. Its basis is the fact that for differential equations with constant coefficients the approximating algebraic problem is reduced to the solution of difference equations with constant coefficients. The general equation is given by @NY,+ %-1Yr+1 + * * * + soy,,,
= b
Y = 1,2,3,
. . .)
(1)
and with the support of initial and boundary conditions the problem reduces to the solution of the algebraic equations * Dedicated
to Professor
A. M. Ostrowski
Linear Copyright
0
Algebra
on his 76th birthday.
and Its Applications
1968 by American
Elsevier
1, 445-463
Publishing
Company,
(1908) Inc.
I.. FOX
44ri
Hert) .-I is a “I~tnd”
matrix
of band width
ISD
I,INl).\
H.\YES
N 1- 1, whose general
ro\\’
is 0, 0, . and J’ is the vector One attractive
.)
0, u,y,
N.y.
,’
with components
.
,
a,, 0, 0, .
0)
,
yi, y,, ~‘a, .
method of solving equations
(2) first performs the matrix
decomposition
where
L is a unit lower triangular
matrix.
If the decomposition
is, if every. leading matrices
submatrix
will have the general
In each general
k + 1 adjacent
and
banded
then
are k + 1 adjacent
and the corresponding
nonzero elements
triangular
“interchanges,” the
that
I. ant1 r.
forms
nonzero
elements
on
row of U will hnvc .V
on and above the diagonal,
k depending on the number of initial conditions equation
Ci an upper
without
of A is nonsingular,
row of L there
and below the diagonal,
matrix
is possible
associated
the number
with the general
(1).
The importance provided
of the method
the second
decomposition one coming
reaches
boundary a “steady
from the I, matrix
is that,
under certain
is far enough state”
away
conditions
from
and
the first,
the
which serves to give two factors,
and the other from the U matrix,
of the
polynomial /J\.(_X)= us 4 aR.__ix $- . . . J- aoN?
In \vhat follows we first give some detailed .Y -- 4, matrix
the
quartic
method
polynomial,
and
with the Q-D algorithm.
illustrate
((9
attention the
Subsequently
to the
connection
case
of the
we treat with less
POLYNOMIAL
ZEROS
AND
Q-D
447
ALGORITHM
detail the general polynomial, and give some numerical illustrations the application THE
QUARTIC,
of
of matrix factorization. GENERAL
FACTORIZATION
3. In the case of a quartic polynomial, the general relation (1) has the form U4Yr +
a3yr+1+
UzYr+2 +
%Yr+3
+
UoYr-l-4 =
b
(7)
The number of associated initial conditions may be 1, 2, 3, or 4, of which the last represents the pure initial value problem and has no particular interest in our present context. The relevant types of condition, roughly equivalent to the respective specification of first, second, and third derivatives for a fourth-order differential equation, can be represented by the initial equations "1Yl
+
PlYI+
YlYl
+
=
6 1’
=
62,
Y4Y4 =
63.
E2Y2 P2Y2 +
P3Y3
Y2Y2 +
Y3Y3 +
(8)
With all of (8) given, our matrix and its decomposition look like A \
a3
a2
a1
a0
a4
a3
a2
a1
a0
L
u
(9) %
Linear
Algebra
and
Its Applications
.
%
1, 445-463
(1968)
.-1
-4
-
1
Y3
7’1 Y2
Ys
?'A
a4
a2
a1
a3
6111 a3
a0
fl2 ...
and with the specification
only of the last of (8) vve ha\~
.I Y2
Y3
Y4
u3
a2
al
(14
a3
%
(1, ...
a,
..*
3 iP3
;;4
7'4
1. 1.
It is now useful to consider
A4~t==-0 of the relevant is at infinitv
so that the general
I. is not singular
and this means homogeneous I.lnrar
Algrhra
algebraic
that
successive
md
the
Its Appl~cnlims
of the homogeneous
form
in which the second boundar!,
row of A is repeated
in (4), it follows
form of (l),
the solution
equations,
indefinitely.
Since
in addition
to the
that
values
of J’? satisfy,
“reduced” 1, 445-46s
recurrence (196X)
relations
POLYNOMIAL
ZEROS
““Y” +
Q-D
AND
V,Y,+1
01
=
449
ALGORITHM
“rY, + VrYr+l+
W*Yr+2 = 0, (13)
%Y” + vryr+1+
relevant
respectively
wry,+2
+
GY,+3
to the decompositions
=
01
(9), (lo),
Now if the zeros x1, x2, x3, and xq of the relevant distinct,
the general
solution
and the initial the
arbitrary
conditions arbitrary
constants
almost certainly then
of largest
reach a “steady
(i)
(8) merely in (14).
(14)
give one or more relations Provided
that
distinct
as a result of inevitable
modulus
dominate
in (la),
large r the decompositions
state.”
Provided
moduli,
Decomposition
real “dominant”
A,%,'+ Qql,
are not zero (and if they are exactly
for some sufficiently pairs have
of type
constants
be “introduced”
the roots
(6) are
of (7), with b, = 0, is
y, = ApI'+ A&'+ between
and (11). polynomial
that
If lxil > lxzl >
zero, then
rounding errors), and we infer that
(9), (lo),
and (11) will
any real roots and any complex
we find the following
(9).
the relevant zero they will
for sufficiently
results.
1x3( 3
1~~1, so that
xi is a
large Y we have
U” + v,xi = 0, and the elements (ii)
u, and v, converge
Decomposition
are two “dominant” Y we have
obtained
If
(10). roots,
(iii)
to constant
lx11 3
the quadratic
Decomposition
(11).
If 1x1( >
large
to fixed
lx2[ >
(16)
numbers.
[~a/> 1x11, so that there is the cubic factor
with zeros xi, x2, and x3,
and u,, v,, w,, and Z, all converge
Linear
there
then we have obtained
u, + v,x + wp2 + z$,
In all cases the other factor
1~~1, so that
for sufficiently
with zeros xi and .x2.
modulus,
but is lurking in the L matrix.
then
factor
uI, v~, and w, all converge
one real root of smallest
u and v, respectively.
Ixsl > lxsj >
real or complex,
24, + v,x + wrx2> The elements
(15)
to fixed
of the quartic
numbers. polynomial
The rules for matrix .4lgebrn
and
(17)
is not “lost,”
multiplication
Its Applications
1, 445-
463
easily (1968)
show
that
corresponding
to (15),
(Iti),
and
(IS)
the
other
factors
are
respecti\+-
and the respective
.5 .
I,, nz,, and 12, con\.erge to fixed uumbvrs.
The most interesting
quadratic
factors,
qualifications
already
even including
and useful case is that of decomposition
represented
by the
mentioned
the situation
decomposition
we can “solve”
not present
that is, the modulus of a real or comples
(10).
the quartic
in (le),
into
\Vith the completely,
the case in which
pair svparatcs
the rnoduli of two
real roots. In this case the factors so that numbers
we cannot
(w,, of course,
does dominate
everything,
can show that,
to a,,).
U matrix,
of the “middle”
real root in the
hi1 -L
.y4 is lurking
“steady m,x
-t
state”
in I,~,~L and m,.
pair
On the other
that
for sufficiently
tllat
in fact
to fixrtl hand
xl
large Y
.y2 and .I:~.
in the L matrix,
and in fact \v(
situation,
x2 = (x4 - x)
where x4(“) will also help to determine suffixes
contributions,
large Y. Here x1 I”)depends on Y, but its I-alue is important
for the computation The smaller
equal
and this implies in the
c~lual
II, and 11, in (10) will con\crgc
is always
the zero ,cI is contained
for all sufficiently
.y2” and s3” in (14) makr,
e.xpect that
(.r4(r’-
the other pair.
X))
(23)
We note the different
POLYNOMIAL
6.
ZEROS
451
AND Q-D ALGORITHM
The result (23) is perhaps easier to verify rather than to prove,
but in the process of the verification we can relate the xi(‘) and the XL’) to the other roots x2 and xs. a4 + agx + a,#
We have
+ alx3 + aox = a&x, -
so that A = x2 + x3, B = x2x3. a4 = a,x,x,B,
a3 = - a,{Ax,x,
The other relevant
equations,
x)(B - Ax + x2),
(24)
Then a, = -
a2 = a,{x,x,
x)(x4 -
+
+ B t
a,(xl + x4 + A), w,
+
X4,>>
Ah
+
~4) 1.
coming from the matrix
(25)
decomposition
(lo), are
@r-2 =
lIv,_2+
a4,
&a, + m,v,_, + u, = a2, Equations
m,uu,_l
=
m,a, + v, =
a3,
al,
w, = ao.
(26)
(22) and (23) give Qr,
= x1x1(I),
1r+l =
x4x4
(I),
- v,/a, = x1 + xl(“), (27) -
m, =
x4 +
x4@).
From (26), (27), and the first two of (25) we easily find B = x4(‘)x,(‘-i)
A = x4(y)+ xi(‘)
(28)
and substitution in the second pair of (25) with further use of (26) verifies that everything is satisfied. Equation (28) gives the required formulas for determining the middle pair from the “steady state” situation. EXAMPLES
7. In each of the following three examples we take the first two rows of the matrix A in (10) to be respectively (a,, a,, a,) and (a3, a2, a,, ao). For the first example [3] we consider the quartic p,(x) = 1 - 32x + 160x2 - 256x3 + 128x4, which has four real zeros of distinct moduli. we obtain in 24 steps the quadratic factors Linear
Algebra
(29)
To six significant figures
and Its Applications
1, 445-463
(1968)
0.0117476 and there correct
0346719x
is no further
after
For
-
85.1237
-
211.62OOm $- ldW,
change to this precision.
and in 1X steps
with two complex
we obtain
0.489555
-
In our third
+ AA,
are correct
the second factor,
Some eleven figures are
=
[I]
six-figure 2.04265
after
about
we take
the cluartic
decompositiorr
-- 2.6853i.r
35 steps.
In both casts,
the modulus
which
From
5.9375x" - 4.r” i- x4,
I --- 4.r -t-
three
From
(33)
successive
For
appearance
from
finding of Table
I.
of two real \ve get _ra =
.-\t about
this stage
zeros of the U and IX quadratics
0.536864
0.925692
~~.8%4308
0.669727
1 .ok-N273 of (24) the three
the correct
zeros of a quartic
are
(34)
estimates
value being obvious,
(28) the “correct”
the
implying
_~r= 1.640388,
and from the L matrix
121313x
for the .l
and 1.750000,
the moduli
to the figures gi\ren.
pairs of the other
(28) we deduce
of (24) we have
pair separates
we find from the L’ matrix
to this precision,
which is also correct
1.750000,
8.
of a complex
Y = 27 onward
is correct
0.609612,
of
the quartic
we find that the I,, ML,,u,, and U, do not converge to fixed numbers, roots.
(32)
+- .\.2.
with the larger zeros, belongs to the I: matrix.
example, b4(x)
pairs
the correct
1.31463x
.\gain eleven figures
that
(SO)
35 steps.
an example
course,
_t ~2,
two estimates
I.750002,
and for the
14
of unit\,.
the Q-D array
~3] has
the
POLYNOMIAL
TABLE
AND
Q-D
ALGORITHM
453
I
q,(l)
-
ZEROS
#)
q*w
al
e,(V
0
q,(3)
e,(3)
qr(4)
0
0
a0 0
%
-a3
a4
al
%
a3
q&V
qP
q!f)l
e,!l)
0
0
e!?;
0
qow
q!!;
El(%)
t?,(l) qP
q’“‘,
e,@J ql(“)
qP
0
qp
qP
0
&?,(3)
qllC4)
... The numbers
in the array are generated W’ =
%+1
according
to the “rhombus
rules”
(k’ _ eSf+;C + q,@‘,
e 7
(35)
ejyl = (q,(k+l)/q;$‘l)ep. If the zeros of the quadratic 14
> I%1 > 14
in modulus,
with jxr/ >
then
g(i) I -*Xi)
q(Z) I +x2,
and the e columns converge with
are real and distinct
to zero.
lx11 = lxs\ > [x3\ = 1x41, then B, -
qCS) I +xs,
lj ICi) + * 4p
(36)
If the zeros form two complex th ese belong
A,% + X2,
B, -
to quadratic
pairs,
factors
A,x f x2,
(37)
and
The
case of a complex
pair with modulus
of two real roots has obvious of Eq.
behavior,
(21), to which we have given
LU method.
For this
p (1) --f Xi’ I
greater
than
the moduli
and the other possibility some detailed
attention
is that with the
case we have
@) + qj:S’ ---f‘4, Q7tl Linear
Algebra
ql’d’q;“) + B, md
Ifs
.-lpplicatims
qj”’ + x4. 1, 445-463
(39) (1968)
1,. ITOX .\SI)
454 !I.
Sow
giv~ii
if the
in the
first
that
\vith exact
tion
(lo),
Lli
we should
equations
and
from
(39) \\.v cltduc~,
l‘he
other
roots
of thaw,
10. the
Q-D
unique,
the
coefficients
\vork,
and
perhaps
of rounding
error
avoids
bitrarily,
the
segment It Table
I must
For
w
qy)::
denotetl
in the special
satisf!.
a11~1
ant1
results
if an\-
li.
.I~'"'
(2X)
numcricnl
precision of the
ro\vs
will
usually
to Iiialz
an
polynomial,
LISIWI
“c~tlrck”
of
I i:,
of thcl coefficients ;I new
new
rxxtra
;I minimum The
coefficients. but
ul, that
~)ol\~nomial nirans
to ensure’
suitabl!-,
of tlie matriws
form<9
of Table
0iwurst’,
This,
hi- t.hoosing
tliv
descriptions
it is rc~commentli~cl
to produ(,(,
defect.
that
al-I-ang-ement
(6) is LW).
computation
quartic the
clsist
rewals
published
thr
he shifted
csstra
matrix
possible
that
not
hal.cb this
difficult\.
the
5
111. .I,“’
of (40) , anti &r
tlw
First,
does
should
with
unit
is in fact
and
pol\nomial
do not
first
\\hich
qLfl
13 j) implv
method
this
of the
algorithm.
to tllcs tl(ac,omy)osi-
(39).
(e.g.,
(42, or ug of the quartic in this cast‘ the origin whose
the
advantages.
algorithm
method
I-OI~~ 2s to slio\\~
large: Y, tlkal
of the l,I- ant1 (2-D mr~thocls
significant
so that
for 5ufficirntl~~
from
(‘omparison
some
t\vo
difficult
find, corresponding
cluadratics,
5, are clearly
immecliatel!-
has
tlici tirst
it is not
give) an ~a>\. proof of tilt> rc5iilts oi Sc-rtions i? :; from (40) \vk~ fintl qy)? and y,’
if we eliminate
follow
with
7, then
~l.\S155
identities
l~or
in Section
is p~~rformed
of Section
arithmetic
the
‘l‘ticw
method
sentence
l.lSl)\
I,(.
othtwvise
.i in (!1)--(I I ),
al-=
(Ilearl!-
2
be satisfactory. “arhitrq-”
start
for c,sample, quation
\vith the
the
iiumhcrs
Q-U in
POLYNOMIAL
ZEROS
AND
Q-D
455
ALGORITHM
and we can also establish the identities
Together with the rhombus rules, it turns out that we can start the Q-D table, and build it up to the last two rows given explicitly in Table I, with arbitrary values of q,,(l), e,,(l), and %(‘). The remaining elements in the relevant triangle come successively from the equations ql(l) = qow + eo(l!, qp
= -
elm =
qlN + e,(l),
q2(l) =
[((u&p + %)/q,(1)+ %Nq2(1) + %I/~,~
Ql(2) = q2(l)e2(l)/elU), (&(P) =
e,(‘)q,(2)/q,(‘),
e,(2) =
(@ az...1__ + a-3 -~~ Go92(l)
a, ql(l)q2(1)ql(2)
ql’“) = - ~,/~, - qp
q3C1) - q2C1),
q1t2)+ elil) - q,c’r),
+ q1t2)) +a4
el@) = q2(2)+ ez(l) _ qlt’)’
e2(li =
(45) 40(2)41(2) + qo(“)q2w + q1U)q2(l) qo~l~q~~l~qo~~~q~~~~q~~~~~ ’ a, i
1
40(3)= e,(2)q,P)/eo(2),
- q2r4 - qow,
eoW)=
q1
q.
(3) +
(4) =
$p&$3J
el(2) -
poP),
and we can then continue by the standard method. The corresponding results for the general case, of course, are not very attractive as the basis for a general starting method. 11. Second, both methods can subsequently fail through the vanishing of a denominator, for example, an early &$i in the Q-D equations (35), or an early U, in the LU decomposition. (Clearly, if convergence takes place, U, tends to a nonzero value). Whereas the Q-D failure necessitates a new start with a shifted polynomial or with equations like (45), we can with the LU method replace the offending SC,by an arbitrary nonzero value, say unity, and proceed with the iterative process. The effect can be regarded as that due to some different set of initial conditions, that is, of different initial rows of A, or the result of an isolated “blunder,” neither of which can do more than delay the subsequent convergence. Consider, for example, the polynomial p,(x) = 13 + x + 2x2 - X3 + x4. Linear
Algebra
and
Its Af$lications
(46) 1, 445-463
(1968)
W’ith the “Q-D start” this Ita bJ* unity
of the LI’ method we find ~1~=m0, but if \VC’replacc~
in the LU method
quadratic
factors,
50 steps.
The Q-D algorithm
a complete 12.
restart
Third,
uccww~lntion LU
a six-figure
of rounding
method.
the L and
error,
For suppose 7’ matrices
is reported
to
to obtain in about 0,
qrf’i from
sulfu
this is \irtuallv
we ha\e reached
reached
\Ve can then consider restarting
the corresponding
whereas
that
have
happily
ant1
necessar!‘.
algorithm
Q-D
quite
state result being reached
produces
is apparently
the
we proceed
steady
a steady
some
state
the whole operation
;t slo\\~
abscw t in tllcs stage
at which
to some prrcision.
with the factorization
I,
-1
!. 0
;_
. .
where if, ~1,and ZI’are the stead!, state elements ‘This is simply equivalent representing
the fact
in (14) are already efficients
to starting
that
quite
with
in the LU method
polynomial,
whereas
only in the first Moreover,
those
start
the
One might
are “refreshed”
relevant
say that b!- frequent
of the Q-D method
make
conditions,”
coefficients
Li,
the successive
co-
contact
with the
explicit
contact
few steps.
if we kno\v something
we can use this information start,”
this
small.
of tlir, previous operations.
with very good “initial
about
the distribution
quite deliberately
for which the unwanted
coefficients
It is only fair to add, however, tried this device has not achieved
of thti zeros
to make the “best possible A, in (14) are already
that in many of the examples a spectacular
improvement
small.
WC have in con\-cr-
gence ! Finall!.,
13. properties factor
the LO7 method
belonging
to the I,’ matrix
the rate of convergence of the IJ factor method, of similar
has some slightly
from those of the Q-D algorithm.
modulus.
con\ergencc
has zeros of \-ei-h.similar modulus,
depends 0111~7on the ratio of the smallest
to that of the largest modulus
on the other
different
It does not matter
hand,
has some
trouble
of the L factor. in srpar-sting
if the and
modulus The
Q-D
an>. roots
POLYNOMIAL
ZEROS
AND
Q-II
457
ALGORITHM
For example, with 25 - 56x + 38x2 - 8x3 + x4 = (x -
1)2(x2 - 6x + 25),
(48)
the LU method separates the two relevant quadratic pairs very rapidly. We find I,, = 1.00000058,
ml2 = - 2.00000024, 752 =
-
uu12= 2.49999993, (49)
5.99999976,
which are very accurate, and at Y = 20 these coefficients are all correct to 13 figures. With the Q-D method the individual q/“) and q,‘“), representing the two equal real roots, converge very slowly, being correct to only 1% at about Y = 70. In accordance with Eq. (40), of course, relevant combinations of the qr(k) “settle” at the same time as the L and U coefficients, and we might conclude that the corresponding sum and product checks could be used with advantage even in the case of real roots of nearly equal modulus. Similar results are obtained for the quartic with zeros 5, - 5, 1, and 1. The LU process gives quadratic factors with coefficients correct to five figures at Y = 9, and to 12 figures at Y = 19. In the Q-D process the two equal smaller roots are correct to 1% only at Y = 100, and the columns for the larger roots do not appear to be converging. But at r = 9 we find qJi’ = 13.000007,
qjy
qi4), = 0.899160,
= -
13.000008,
qpil = 1.923077,
qyll = 1.091602,
qPJ2 = 1.100841, q,@) = -
1.923076,
(50)
q;y2 = 0.908397,
and Eqs. (38) give very nearly the correct quadratic pairs x2 - 25 and X2 -2x+ 1. The case of equal real roots belonging to the U matrix is determined quite easily by considering, instead of (14), the relevant general solution y* =
(A, +rA,)x,*+
A,%,'+
of our basic recurrence relation.
A*X4",
I%/ >
(51)
IX3b4~
Here we have, for sufficiently
large
7, Linear
Algebra
and
Ifs
Applicatims
1, 445-
463
(1968)
and we get convergence case of equal 14. have
example,
each
factor
is ultimately
zero in the
roots!
The difficult
equal
because
problems
modulus.
When
arise when at least three of the four roots all the roots
ha\-e the same
modulus,
for
in an!. of the cases
(.53)
they cannot, there
so to speak,
decide which matrix
is no easy way of deducing
(53) gives some hope, because converges,
but
only very
\Vhen only three other
real root.
of which the last factor no such problem state”
factors.
The
and
third
of
of the term rA, in (Til), and this in fact
slowl!-.
roots ha\-e equal modulus
For
with this modulus
they want to inhabit,
the quadratic
example,
out thr
for the quartic
has complex
have the dilemma
and attaches
we can separate
roots of modulus mentioned,
6, the three roots
but the small root has
itself firmly to the L mxtris.
The “stead\.
is gi\ven bTV
where .Y~is the lone small root. we can also show that the steady
Hy methods
similar to those of Section
state remaining
cubic equation
6
is given
b!,
In this results
particular
example
we find
after
a few steps
the se\-en-figure
POLYNOMIAL
ZEROS
AND
Q-D ALGORITHM
1,, + mrex + x2 = 4.187036
-
5.187036x
us = 29.854055,
xq(i’) = 4.187036,
u i. = 17.996859, and our computed x3.
vr,, = -
+ x2,
x4=
zig= -
8.837551,
of course,
separate
125 -
out the roots
55x + 11x2 -
given in the first part of this section.
of equal
modulus
Eq.
and the “reduced
xl(q
+2
POLYNOMIALS
(55) is replaced
cubic”
corresponding
x(.+?z,
+
OF HIGHER
t1
by
by
x2(x1(*) - nz,) - x3 = 0.
(59)
DEGREE
Many of these results can doubtless be extended
ing treatment
(58)
to (56) is given
Z,+,) +
-
with three smaller
x)(Q) - 4,
w,,x2 = (x1 -
v,x +
21, +
15.
(57)
of this cubic for the
In the case when there is a single real root associated roots
1,
6.812964,
cubic is very close to the correct
We cannot,
reasons
459
of polynomials
of other
to the correspond-
and higher degrees.
In any case
we have seen that any splitting is possible provided that all zeros belonging to the U matrix
have moduli greater
than those of the L matrix.
moduli occur only in complex roots, therefore, which will converge of L and
U settle
The best in problems
in the absolute
sense, that is, in which the elements
down to constant
application
values.
of the suggested
splitting
(or smallest) so that
here have
modulus
the certainty,
the “reduced”
lacking
16.
To
test
more the
by Wilkinson
in Newton
method
on fairly
if we seek the k zeros
the
appropriate
the relevant or Bairstow
initial
k zeros.
We
methods,
that
of the distribution
(fastest)
the required
from Olver
is probably of the zeros
and we can then concentrate
the best
or less than
worked on three examples studied
includes
Some knowledge
can help to achieve
if this includes
however, the nature
For example,
our objective,
polynomial.
about
we can make
U (or L) certainly
we have in fact achieved moreover,
method,
in which we know something
and in which we need specific results. of largest
If equal
there is always some splitting
difficult
initial
on
of zeros,
splitting,
even
k zeros. polynomials
we have
[a], two of which have also been
[5]. Linear
Algebra
a?zd Its Applications
1, 445-463
(1968)
The
polynomial
1 + X.Y + 32.8”
t 89.6x” ~,-190.68.x4 -I- 304.08x5
+ 443.,576YG +- 468.88.T;
$mFi24.327,~~ + 378.908.~~ -+-34Fi.072FiB~~~-+ 166.44768~~’ -1. 37,6lilO96x’3
+ 26.1783048.~‘~
11as approximate 0.293.51 f
+ 3.4356048~~‘~ + 2.032.53121 x1’
0.4FiO93i, ~ 0.14762 + O.i7176i,
0.09004 + l.O6119i,
0.0608~ i + 1.2969li,
0.01049 + I .5963Oi, -
(~.00249 -1.: 1.66712i.
that
the
zeros)
order of modulus.
fastest
in which
to the L matrix.
(into
This is (partly)
with
numbers
verified
(the smallest
places in 100 iterations,
L the corresponding
of thts moduli would suggest
factors
an e\~n number
successivelv~ 2, 4, (i, 8, 10, 12, and
\Yith six zeros in the L matrix decimal
~~ 0.02,567 + 1.474383’, (61)
Inspection
decompositions
are those
belong
(60)
zeros
O.l435Oi, ~~~0.22447 f
in ascending
-t 128.218748~~~
of
14 zeros
by our computations.
we tried) we converge
to IO
and with 8, IO, 12, and 14 zeros in the
of correct
decimal places after 100 iterations
are 7, 4, 3, and 2 respectivel\-. To test polynomial
the
speed
and accuracy
we performed
each of the two factors, iteration
was terminated
eight decimal places. results
produced
the (8, 8) splitting, when all relevant
The computed
hv. \I’ilkinson
0.293504?53 & 0.1434993oi3,
-
0.14762378
l 0,7717fi71Ri,
(-
0.14762378
4 0.7717fi72Oi),
-- O.OFiO86440+ 1,2969113i, (-
O.OFiO86444+ 1,2969113i),
--~ 0.01049254
+ 1,5962951i,
(-- O.O10493,5Fi:-i 1,,59629EiAi), I.iwrtrr,Algchvcr
ctnd
method
T/s .?~~licc~linm
coefficients
all the zeros.
aw
& 0,45092796i,
0.22445006 0.09003999
+ 0.4fiO92796i), + 1.0611!~21 i,
0.09003999
% l.O61192li),
~~~0.02t56687.5 3: 1.4743771i, (-
0.02566871
~-- 0.00249022 (1,
of
Each
had convrerged to
in parentheses,
-- 0.22447006
(-
particular
values of the zeros, and the accurate
[fi], given
(-
for this
then the (4, 4) splitting
and so on until we obtained
0.293,504,53 * 0.1434992&‘, (-
of the
+ 1,66712OXi,
0.00248920 44S-463
& 1.474377 li).
+ 1.6671204~‘)
(1968)
(62)
POItYNOMIAL
ZEROS
AND Q-D ALGORITHM
461
The complete computation took 28 seconds on the KDFO machine with an ALGOL compiler. Some of the zeros are fairly close together in modulus, and the results obtained are really quite satisfactory. 17.
The polynomial
420 + 1882x + 4090x2 + 5276x3 + 5644x4 + 1019x5 + 3651x6 - 4900x’ + 3912x8 - 4618x9 + 3086~~~ -
1916+
+ 1056~‘~ - 409x13 + 139x14
_ 36x15 + 4%” has real zeros 34, 3, 24, 2 and twelve complex zeros with moduli < 1.94. This is a rather more ill-conditioned problem. The splitting into two polynomials of degree eight is quite fast, but one of the next splittings (of the U matrix) is relatively slow, and the whole operation took 53 seconds. To seven decimal places we find the zeros 3.5000600,
3.OOOOOO0,2.5000600,
2.0000000,
- 0.0112583 & 1.9332891i,
- 0.0449754 & 1.7381625i,
- 0.1012792 f
- 0.1812607 f
1.429775&,
- 0.2844601 f 0.5987694i,
l.O348308i,
(64)
- 0.3767663 + 0.1888406i.
The reason for the slow convergence of the second U splitting is quite obvious, the first complex root having modulus very near to that of the nearest real root. The maximum error in (64) nowhere exceeds a few units in the seventh decimal. 18.
Our final example, the determination of the zeros of the poly-
nomial 2 + 16x + 224x2 + 1288x3 + 10416x4 + 44184x5 + 267232x6 + 837860x’ + 41782600~s+ 9490840x9 + 41018752~~~ + 64249356~~~ + 247926664xr2 + 240775148~~~ + 845947696x14 + 385455882~~~ + 1250162561+,
(65)
represents an extremely ill-conditioned problem. All the zeros are complex, six pairs have very similar modulus, and the last four pairs have almost equal modulus. Linear
Algebra and Its Afiplications
1, 445- 463 (1968)
For
the first
To test
time
convergence
is here affected
this we ran the program
after 500, 1000, and 2600 steps, of which consumed to five decimal
about
places
results
of our method.
‘I‘.4BLE
11
Correct results
three
times,
respectively,
I10 seconds
the correct
errors.
each splitting
in the three runs, the first
of machine
results
500 iterations
by rounding
truncating time.
Table
[Fi] and the three
1000 itcratwns
II shows computed
2.WO iterations
-0.13245~_0.13601i -0.13246~0.13601i
-0.13245~0.136001
-0.13245 :,O.l359Hi
-0.01869~0.25305i
-0.01870~0.25304i
-0.01870~0.25305i
--0.01869&0.25304i
-O.O0232&0,29258i
-0.00232&0.29256i
-0.00233~0.29:'38~ -0.00233 I 0.29258i
-0.00049+0.30418i
-0.000301_0.30414i -0.00029&0.30418i
-0.00014~0.308612' -0.00146&0.30785i -0.00005~0.31066i -0.00000~0.31220i
+0.00182~0.3099Oi
-0.00161 4 0.31201; -0.00168- 0.31222i
+0.00060~0.312623' +0.00071~0.31319i
The results are interesting,
t~O.OO0811:~0.313ORi
in that there is no significant
with more than 500 iterations. mined
-0.00136_t0.30800i
+0.00182 t 0.31011i +0.00179&0.30998i
--0.00001~0.31170i -0.00134~0.31264i
rounding
0.00026~0.30418i
-0.00130f0.30799i
improvement
This is not due, to any accumulation
of
error, but to the fact that the later roots ~I-C \.er>r poorly deter-
(ill-conditioned)
with respect
The effect of individual can be interpreted, to factorize .1t the next
to small changes
step
these
situation
in the coefficients.
errors in the computed
in the spirit of backward
the given
ill-conditioned
rounding polynomial
with
perturbations
matrix
error analysis,
slightly
perturbed
are slightly
this has a serious
effect
elements
as an attempt coefficients.
different,
and in
on the zeros,
an
and also
on the rate and even the existence
of con\rergence to more than a definite
number
can br alleviated
of figures.
of floating-point scalar
products
this matrix because
method
the device
is apparently
that Linew
is particularly is not included
rather
far more
however,
iteration,
which permits
prior to single-precision
It is perhaps result,
This situation arithmetic
surprising
with double-precision
the last real part .4lgebvn
and Ifs
rounding,
well suited. that
Applicutiors
results.
and treble-precision
had only the figures 1, 445-463
ALGOI,
compiler.)
zeros the real part
the imaginary
by Wilkinson’s
of
(It \vas not used here
in the later than
accumulation
and for this purpose
in the standard
ill-conditioned
is also reflected
by using the variant
the accurate
arithmetic,
-- 0.00000305 (196x1
part.
This
Using Hairstow he found
in common,
POLkCNOMIAL
ZEROS
AND
463
Q-D ALGORITHM
while the imaginary part was the same at 0.312196968, more decimal place and to five more significant
that is, to one
figures.
CONCLUSION
19. We have demonstrated that the method of factorizing a polynomial in terms of a certain matrix decomposition leads to a quite practicable technique for determining some or all of the zeros of the polynomial. In particular we can factorize the polynomial directly into two factors which contain specified numbers of the roots of smallest or largest modulus, and can then concentrate on the part which may be relevant to our particular requirements. Moreover much information is gained at all stages, by observation of the rate of convergence, of the existence of roots of similar modulus, and of their precise positions in the ordering of the zeros. This information is not easily obtained in the early stages of other methods of iteration, such as those of Bairstow or Muller. On the other hand, the Bairstow iteration is probably quite a bit faster, in ill-conditioned cases, than the methods of this paper carried One might, therefore, decide that the to the limits of their accuracy. advantages mentioned in the last paragraph would make this method a useful preliminary to a Bairstow-type “cleaning-up” operation. Various questions require further consideration, such as a satisfactory error analysis and a comprehensive program which would do all that anybody might want to do, and we hope to make some subsequent progress on these and other allied matters. REFERENCES 1 L. Fox, Numerical Solution of Boundary-Value
Problems, Clarendon Press, Oxford,
1957. 2 F. L. Bauer, wiss.
3 P. Henrici, 4 F. W.
Faktorisierung
Elements of Numerical
J. Olver,
244(1952),
eines Polynoms,
Sitz. Ber.
Bayer.
Akad.
Evaluation
Analysis.
Wiley,
New York,
of zeros of high degree polynomials,
1964. Phil.
Trans.
385-415.
5 J. Wilkinson, Num.
Direkte
1966, 163.
The evaluation
Math. 1(1959),
of the zeros of ill-conditioned
polynomials,
Part II,
167-180.
Received Mavch 27, 1968
Linear Algebra and Its Applications
1, 445-463
(1968)