1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
|
<?xml version="1.0" encoding="UTF-8"?>
<!--
* Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
* Copyright (C) 2006-2008 - INRIA
* Copyright (C) 2011 - DIGITEO - Michael Baudin
*
* This file must be used under the terms of the CeCILL.
* This source file is licensed as described in the file COPYING, which
* you should have received as part of this distribution. The terms
* are also available at
* http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
*
-->
<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns5="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="roots" xml:lang="en">
<refnamediv>
<refname>roots</refname>
<refpurpose>roots of polynomials</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
x=roots(p)
x=roots(p,algo)
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Arguments</title>
<variablelist>
<varlistentry>
<term>p</term>
<listitem>
<para>
a polynomial with real or complex coefficients, or
a m-by-1 or 1-by-m matrix of doubles, the polynomial coefficients in decreasing degree order.
</para>
</listitem>
</varlistentry>
<varlistentry>
<term>algo</term>
<listitem>
<para>
a string, the algorithm to be used (default algo="e").
If algo="e", then the eigenvalues of the companion matrix are returned.
If algo="f", then the Jenkins-Traub method is used (if the polynomial is real and
has degree lower than 100).
If algo="f" and the polynomial is complex, then an error is generated.
If algo="f" and the polynomial has degree greater than 100, then an error is generated.
</para>
</listitem>
</varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
This function returns in the complex vector
<literal>x</literal> the roots of the polynomial <literal>p</literal>.
</para>
<para>
The "e" option corresponds to method based on the eigenvalues of the
companion matrix.
</para>
<para>
The "f" option corresponds to the fast RPOLY algorithm, based on
Jenkins-Traub method.
</para>
<para>
For real polynomials of degree <=100, users may consider the "f" option,
which might be faster in some cases.
On the other hand, some specific polynomials are known to be able to
make this option to fail.
For instance, <literal>p=poly([1.e300,1.e0,1.e-300],'x');</literal>
provokes infinite looping of <literal>roots(p,"f")</literal>
</para>
</refsection>
<refsection>
<title>Examples</title>
<para>
In the following examples, we compute roots of polynomials.
</para>
<programlisting role="example"><![CDATA[
// Roots given a real polynomial
p = poly([1 2 3],"x")
roots(p)
// Roots, given the real coefficients
p = [3 2 1]
roots(p)
// The roots of a complex polynomial
p=poly([0,10,1+%i,1-%i],'x');
roots(p)
// The roots of the polynomial of a matrix
A=rand(3,3);
p = poly(A,'x')
roots(p)
spec(A)
]]></programlisting>
<para>
The polynomial representation can have a significant
impact on the roots.
In the following example, suggested by Wilkinson in the 60s and presented by Moler,
we consider a diagonal matrix with diagonal entries equal to 1, 2, ..., 20.
The eigenvalues are obviously equal to 1, 2, ..., 20.
If we compute the associated characteristic polynomial and compute its roots,
we can see that the eigenvalues are significantly different from the expected
ones.
This implies that just representing the coefficients as IEEE doubles changes the
roots.
</para>
<programlisting role="example"><![CDATA[
A = diag(1:20);
p = poly(A,'x')
roots(p)
]]></programlisting>
<para>
The "f" option produces an error if the polynomial is complex or
if the degree is greater than 100.
</para>
<programlisting role="example"><![CDATA[
// The following case produces an error.
p = %i+%s;
roots(p,"f")
// The following case produces an error.
p = ones(101,1);
roots(p,"f")
]]></programlisting>
<para>
The following script is a simple way of checking that the companion matrix gives the same result as the "e" option.
It explicitly uses the companion matrix to compute the roots.
There is a small step to reverse the coefficients of the polynomial ;
indeed, "roots" expects the coefficients in decreasing degree order,
while "poly" expects the coefficients in increasing degree order.
</para>
<programlisting role="example"><![CDATA[
v= [1.12119799 0 3.512D+13 32 3.275D+27 0 1.117D+41 4.952D+27 1.722D+54 0 1.224D+67 0 3.262D+79 ];
r1 = roots(v,"e"); // With "e" option
dv = size(v,"*");
p = poly(v(dv:-1:1),"x","coeff"); // Reversing v's coefficients
A = companion(p);
r2 = spec(A); // With the companion matrix
max(abs(r1-r2))
]]></programlisting>
</refsection>
<refsection role="see also">
<title>See Also</title>
<simplelist type="inline">
<member>
<link linkend="poly">poly</link>
</member>
<member>
<link linkend="spec">spec</link>
</member>
<member>
<link linkend="companion">companion</link>
</member>
</simplelist>
</refsection>
<refsection>
<title>References</title>
<para>
The RPOLY algorithm is described in "Algorithm 493: Zeros of a Real
Polynomial", ACM TOMS Volume 1, Issue 2 (June 1975), pp. 178-189
</para>
<para>Jenkins, M. A. and Traub, J. F. (1970), A Three-Stage Algorithm for
Real Polynomials Using Quadratic Iteration, SIAM J. Numer. Anal., 7(1970), 545-566.
</para>
<para>Jenkins, M. A. and Traub, J. F. (1970), Principles for Testing Polynomial Zerofinding Programs.
ACM TOMS 1, 1 (March 1975), pp. 26-34
</para>
</refsection>
</refentry>
|
