reasonable to use instead. It assumes that any small imaginary part cannot be determined with an
accuracy better than eps relative to the real part.
In the meantime, I had tested polyroot() with about 30 polynomials from textbooks and numerical
methods books, and in all cases it returned both real and complex roots, as appropriate. While
empirical, the criteria in polyroot() seems to work pretty well. Perhaps the best approach is to try the
criteria now used in polyroot(), but if you can't find the roots you want, then try one of the other two
criteria.
Method 3: Eigenvalue methods
We can find the roots of a polynomial by finding the eigenvalues of this m x m companion matrix:
−
a
m
a
m+1
−
a
m−1
a
m+1
... −
a
2
a
m+1
−
a
1
a
m+1
10...00
01...00
... ... ... ... ...
00...10
where the eigenvalues will be the roots of the polynomial
P
(
x
)
=
✟
i=1
m+1
a
i
x
i−1
cybernesto has written a function called proots() which implements this method:
proots(aa)
Func
Local nn,am
dim(aa)→nn
aa[1]→am
right(aa,nn-1)→aa
list▶mat(aa/(⁻am))→aa
augment(aa;augment(identity(nn-2),newMat(nn-2,1)))→aa
eigVl(aa)
EndFunc
The input argument aa is the list of polynomial coefficients. For our example polynomial,
y1
(
x
)
= x
3
+ 4.217E17 $ x
2
− 3.981E20 $ x − 6.494E22
proots() returns the candidate roots in about 1 second, compared to 6.4 seconds required for
Laguerre's method as implemented in polyroot() above. Also, polyroot() is about 1773 bytes, while
proots() is only 156. Unfortunately, the root at about -141.8 returned by proots() fails the sign test.
As another example, consider the polynomial
p
(
x
)
=
(
x + 1
)
20
=
x20 + 20x
19
+ 190x
18
+ 1140x
17
+ 4845x
16
+ 15504x
15
+ 38760x
14
+ 77520x
13
+ 125970x
12
+
167960x
11
+ 184756x
10
+ 167960x
9
+ 125970x
8
+ 77520x
7
+ 38760x
6
+
6 - 123
To Next Page
Loading...
Need help?
General
