## `polylib::realroots` -- isolate all real roots of a real univariate polynomial

#### Introduction

`polylib::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`.

#### Call(s)

`polylib::realroots(p)`
`polylib::realroots(p, eps)`

#### Parameters

`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.

#### Returns

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.

#### Side Effects

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.

#### Related Functions

`numeric::fsolve`, `numeric::polyroots`, `numeric::realroot`, `numeric::realroots`

#### Details

- All coefficients of `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`.
- The isolating intervals are ordered such that their centers are increasing, i.e., a.i+b.i < a.(i+1)+b.(i+1).
- The number `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 `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.

#### Example 1

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:

#### Example 2

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]

#### Example 3

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:

#### Example 4

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:

#### Example 5

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:

#### Example 6

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:

#### Changes

- new function

Do you have questions or comments?
Copyright © SciFace Software GmbH & Co. KG 2000