Isolate all real roots of a real univariate polynomial
MuPAD® notebooks are not recommended. Use MATLAB® live scripts instead.
MATLAB live scripts support most MuPAD functionality, though there are some differences. For more information, see Convert MuPAD Notebooks to MATLAB Live Scripts.
polylib::realroots(p) returns intervals isolating
the real roots of the real univariate polynomial
polylib::realroots(p, eps) returns refined
intervals approximating the real roots of
the relative precision given by
All coefficients of
p must be real and numerical,
i.e., either integers, rationals or floating-point numbers. Numerical
symbolic objects such as
are accepted, if they can be converted to real floating-point numbers
same holds for the precision goal
The isolating intervals are ordered such that their centers are increasing, i.e., ai + bi < ai + 1 + bi + 1.
nops(realroots(p)) of intervals
is the number of real roots of
p. Multiple roots
are counted only once. Cf. Example 3.
Isolating intervals may be quite large. The optional argument
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 ,
i.e., each center approximates
a root with a relative precision
Some care should be taken when trying to obtain highly accurate
approximations of the roots via small values of
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
It may be faster to use moderate values of
obtain first approximations of the roots via
These approximations may then be improved by a fast numerical solver
an appropriately high value of
DIGITS. Cf. Example 6. However, note that
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 .
polylib::realroots returns rigorous bounds
for the real roots of the rationalized polynomial. Despite the fact
that all coefficients are approximated correctly to
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. See Example 4 and Example 5.
The function is sensitive to the environment variable
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.
We use a polynomial expression as input to
p := (x - 1/3)*(x - 1)*(x - 4/3)*(x - 2)*(x - 17):
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
because no floating-point arithmetic is used:
We convert these exact bounds for the real roots to floating
point approximations. Note that with the default value of
ignore 2 of the 12 correct digits the rational bounds could potentially
map(%, map, float)
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
polylib::realroots(orthpoly::legendre(5, x), 10^(-DIGITS)):
map(%, float@op, 1)
We consider a polynomial with a multiple root:
p := poly((x - 1/3)^3*(x - 1), [x])
Note that only one isolating interval
1] is returned for the triple root :
We consider a polynomial with non-rational roots:
p := (x - 3)^2*(x - PI)^2:
Converting the result of
floating-point numbers one sees that the exact roots
PI, PI are approximated only to 3 decimal places:
map(polylib::realroots(p, 10^(-10)), map, float)
This is caused by the internal rationalization of the coefficients
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.
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)
Now we introduce roundoff errors by replacing one entry by a floating-point approximation:
p2 := (x - 1.0/3)^3*(x - 4/3):
In this example rationalization caused the triple root
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
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)
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
DIGITS := 1000:
roots := map(approx, x0 -> numeric::fsolve([p = 0], [x = x0]))
[[x = 1.489177598846870281338916114673844643894...], [x = 1.752191733304413195335101727880090131407...], [x = 3.255184555797733438479691333705558491124...]]
delete approx, DIGITS, roots, x0:
A univariate polynomial: either an expression or a polyomial
of domain type
A (small) positive real number determining the size of the returned intervals.
List of lists [[a1, b1],
…] with rational numbers ai ≤ bi is
returned. Lists with ai = bi represent
exact rational roots. Lists with ai < bi represent
open intervals containing exactly one real root. If the polynomial
has no real roots, then the empty list
[ ] is returned.