polylib::realroots
-- isolate
all real roots of a real univariate polynomialpolylib::realroots
(p)
returns intervals
isolating the real roots of the real univariate polynomial
p
.
polylib::realroots
(p, eps)
returns refined
intervals approximating the real roots of p
to the
relative precision given by eps
.
polylib::realroots(p)
polylib::realroots(p, eps)
p |
- | a univariate polynomial: either an expression or a
polyomial of domain type DOM_POLY . |
eps |
- | a (small) positive real number determining the size of the returned intervals. |
A list of lists [[a1,b1], [a2,b2], ...] with rational
numbers a.i <= b.i is returned. Lists with
a.i=b.i represent exact rational roots. Lists with
a.i<b.i represent open intervals containing exactly one
real root. If the polynomial has no real roots, then the empty list
[ ]
is returned.
The function is sensitive to the environment variable DIGITS
, if there are non-integer or
non-rational coefficients in the polynomial. Any such coefficient is
replaced by a rational number approximating the coefficient to DIGITS
significant decimal
places.
numeric::fsolve
,
numeric::polyroots
,
numeric::realroot
,
numeric::realroots
p
must be real and numerical,
i.e., either integers, rationals or floating point numbers. Numerical
symbolic objects such as sqrt(2)
, exp(10*PI)
etc. are accepted, if they can be converted to real floating point
numbers via float
. The
same holds for the precision goal eps
.nops(realroots(p))
of intervals is the
number of real roots of p
. Multiple roots are counted only
once. Cf. Example 3.eps
may be used to refine the intervals such that they
approximate the real roots to a relative precision eps
.
With this argument the returned intervals satisfy b.i-a.i <=
eps
* abs((a.i + b.i)/2), i.e., each center
(a.i+b.i)/2 approximates a root with a relative precision
eps/2
.Some care should be taken when trying to obtain
highly accurate approximations of the roots via small values of
eps
. Internally, bisectioning with exact rational
arithmetic is used to locate the roots to the precision
eps
. This process may take much more time than determining
the isolating intervals without using the second argument
eps
in polylib::realroots
. It may be faster
to use moderate values of eps
to obtain first
approximations of the roots via polylib::realroots
. These
approximations may then be improved by a fast numerical solver such as
numeric::fsolve
with
an appropriately high value of DIGITS
. Cf. Example 6. However, note that polylib::realroots
will
always succeed in locating the roots to the desired precision
eventually. Numerical solvers may fail or return a root not belonging
to the interval which was used for the initial approximation.
Unexpected results may be obtained when the
polynomial contains irrational coefficients. Internally, any such
coefficient c is converted to a floating point number. This
float is then replaced by an approximating rational number r
satisfying abs(r-c) <= 10^(-DIGITS)* abs(c). Finally,
polylib::realroots
returns rigorous bounds for the real
roots of the rationalized polynomial. Despite the fact that all
coefficients are approximated correctly to DIGITS
decimal
places this may change the roots drastically. In particular, multiple
roots or clusters of poorly separated simple roots are very sensitive
to small perturbations in the coefficients of the polynomial. Cf.
Examples 4 and 5.
We use a polynomial expression as input to
polylib::realroots
:
>> p := (x - 1/3)*(x - 1)*(x - 4/3)*(x - 2)*(x - 17):
>> polylib::realroots(p)
[[0, 1], [1, 1], [1, 2], [2, 2], [16, 32]]
The roots 1 and 2 are found exactly: the corresponding
intervals have length 0. The other isolating intervals are quite large.
We refine the intervals such that they approximate the roots to 12
decimal places. Note that this is independent of the current value of
DIGITS
, because no
floating point arithmetic is used:
>> polylib::realroots(p, 10^(-12))
[[1466015503701/4398046511104, 733007751851/2199023255552], [1, 1], [1466015503701/1099511627776, 733007751851/549755813888], [2, 2], [17, 17]]
We convert these exact bounds for the real roots to
floating point approximations. Note that with the default value of
DIGITS=10
we ignore 2 of the 12 correct digits the
rational bounds could potentially give:
>> map(%, map, float)
[[0.3333333333, 0.3333333333], [1.0, 1.0], [1.333333333, 1.333333333], [2.0, 2.0], [17.0, 17.0]]
>> delete p:
Orthogonal polynomials of degree n have
n simple real roots. We consider the Legendre polynomial of
degree 5, available in the library orthpoly
for orthogonal
polynomials:
>> polylib::realroots(orthpoly::legendre(5, x), 10^(-DIGITS)):
>> map(%, float@op, 1)
[-0.9061798459, -0.5384693101, 0.0, 0.5384693101, 0.9061798459]
We consider a polynomial with a multiple root:
>> p := poly((x - 1/3)^3*(x - 1), [x])
4 3 2 poly(x - 2 x + 4/3 x - 10/27 x + 1/27, [x])
Note that only one isolating interval [
0,
1] is returned for the triple root 1/3:
>> polylib::realroots(p)
[[0, 1], [1, 1]]
>> delete p:
We consider a polynomial with non-rational roots:
>> p := (x - 3)^2*(x - PI)^2:
Converting the result of polylib::realroots
to floating point numbers one sees that the exact roots 3, 3, PI,
PI
are approximated only to 3 decimal places:
>> map(polylib::realroots(p, 10^(-10)), map, float)
[[2.998807805, 2.998807805], [3.001213582, 3.001213583], [3.140323518, 3.140323519], [3.142840401, 3.142840401]]
This is caused by the internal rationalization of the
coefficients of p
. Information on the rationalized
polynomial may be optained by a builtin userinfo
command:
>> setuserinfo(polylib::realroots, 1):
>> polylib::realroots(p, 10^(-10))
Info: polynomial rationalized to poly(x^4 - 12283185307/1000..) ...
The intervals returned by
polylib::realroots
(p, 10^(-10))
correctly
locate the 4 exact roots of this rationalized polynomial to a precision
of 10 digits. However, because all 4 roots are close, the small
perturbations of the coefficients introduced by rationalization have a
drastic effect on the location of the roots. In particular,
rationalization splits the two original double roots into 4 simple
roots.
>> setuserinfo(polylib::realroots, 0): delete p:
We consider a further example involving non-exact coefficients. First we approximate the roots of a polynomial with exact coefficients:
>> p1 := (x - 1/3)^3*(x - 4/3):
>> map(polylib::realroots(p1, 10^(-10)), map, float)
[[0.3333333333, 0.3333333333], [1.333333333, 1.333333333]]
Now we introduce roundoff errors by replacing one entry by a floating point approximation:
>> p2 := (x - 1.0/3)^3*(x - 4/3):
>> map(polylib::realroots(p2, 10^(-10)),map,float)
[[0.3332481323, 0.3332481323], [1.333333333, 1.333333333]]
In this example rationalization caused the triple root
1/3
to split into one real root and two complex conjugate
roots.
>> delete p1, p2:
We want to approximate roots to a precision of 1000 digits:
>> p := x^5 - 129/20*x^4 + 69/5*x^3 - 14*x^2 + 12*x - 8:
We recommend not to obtain the result directly by
polylib::realroots
(p,10^(-1000))
, because the
internal bisectioning process for refining crude isolating intervals
converges only linearly. Instead, we compute first approximations of
the roots to a precision of 10 digits:
>> approx := map(polylib::realroots(p, 10^(-10)), float@op, 1)
[1.489177599, 1.752191733, 3.255184556]
These values are used as starting points for a numerical
root finder. The internal Newton search in numeric::fsolve
converges
quadratically and yields the high precision results much faster than
polylib::realroots
:
>> DIGITS := 1000: Roots := []:
>> for x0 in approx do Roots := append(Roots, numeric::fsolve([p = 0], [x = x0])); end_for
[[x = 1.489177598846870281338916114673844643894...], [x = 1.752191733304413195335101727880090131407...], [x = 3.255184555797733438479691333705558491124...]]
>> delete approx, Digits, Roots, x0: