## `numeric::polysysroots` -- numerical roots of a system of polynomial equations

#### Introduction

`numeric::polysysroots``(eqs, ..)` returns numerical approximations of all real and complex roots of the polynomial system of equations `eqs`.

#### Call(s)

`numeric::polysysroots(eqs)`
`numeric::polysysroots(eqs, vars)`

#### Parameters

`eqs` | - | a polynomial equation or a list or a set of such equations. Also expressions and polynomials of domain type `DOM_POLY` are accepted wherever an equation is expected. They are interpreted as homogeneous equations.
`vars` | - | an unknown or a list or set of unknowns. Unknowns may be identifiers or indexed identifiers.

#### Returns

a set of lists of equations. The set `{[]}` containing an empty list is returned, if no solutions can be computed.

#### Side Effects

The function is sensitive to the environment variable `DIGITS`, which determines the numerical working precision.

#### Related Functions

`linsolve`, `numeric::fsolve`, `numeric::linsolve`, `numeric::polyroots`, `numeric::realroot`, `numeric::realroots`, `polylib::realroots`, `solve`

#### Details

- The coefficients of the polynomials may contain symbolic parameters.
- If no unknowns are specified by `vars`, then `numeric::indets``(eqs)` is used in place of `vars`.
- The solution is a set of lists of solved equations of the form

      [x1=value1,x2=value2,..] ,

  where x1,x2,.. are the unknowns. These simplified equations should be regarded as constraints on the unknowns. E.g., if an unknown x[1], say, does not turn up in the form x[1]=.. in the solution, then there is no constraint on this unknown and it is an arbitrary parameter. This holds true in general for all unknowns that do not turn up on the left hand side of the solved equations. Cf. example 2.
- The ordering of the unknowns in `vars` determines the ordering of the solved equations. If a set `vars` is used, then an internal ordering is used.

- If the solution set of `eqs` is not finite, then `numeric::polysysroots` may return solutions with some of the unknowns remaining as free parameters. In this case the representation of the solution depends on the ordering of the unknowns! Cf. example 3. Further, if higher degree polynomials are involved, then `numeric::polysysroots` may fail to compute roots. Cf. example 5. The same may happen, when `eqs` contains symbolic parameters.

- You may try `numeric::fsolve` to compute a single numerical root, if `numeric::polysysroots` cannot compute all roots of the system. Note, however, that `numeric::fsolve` does not accept symbolic parameters in the equations.
- We recommend to use `numeric::polyroots` to compute all roots of a single univariate polynomial with numerical coefficients.
- `numeric::polysysroots` is a hybrid routine: it calls the symbolic solver

                solve(eqs, vars, BackSubstitution = FALSE)

  and processes its symbolic result numerically via `solvelib::allvalues`. Cf. example 4.

#### Example 1

Equations, expressions as well as `DOM_POLY` objects may be used to specify the polynomials:

    >> numeric::polysysroots(x^2 = PI^2, x)

                      {[x = -3.141592654], [x = 3.141592654]}

    >> numeric::polysysroots({x^2 + y^2 - 1, x^2 - y^2 = 1/2}, [x, y])

          {[x = -0.8660254038, y = -0.5], [x = -0.8660254038, y = 0.5],

             [x = 0.8660254038, y = -0.5], [x = 0.8660254038, y = 0.5]}

    >> numeric::polysysroots({poly(x^2 + y + 1), y^2 + x = 1}, [x, y])

          {[x = -0.4533976515, y = -1.20556943], [x = 0.0, y = -1.0],

             [x = 0.2266988258 - 1.467711509 I,

             y = 1.102784715 + 0.6654569512 I],

             [x = 0.2266988258 + 1.467711509 I,

             y = 1.102784715 - 0.6654569512 I]}

Symbolic parameters are accepted:

    >> numeric::polysysroots(x^2 + y + exp(z), [x, y])

                                      1/2
          {[x = (- 1.0 y - 1.0 exp(z))   ],

                                              1/2
             [x = - 1.0 (- 1.0 y - 1.0 exp(z))   ]}

#### Example 2

The returned solutions may contain some of the unknowns remaining as free parameters:

    >> numeric::polysysroots({x^2 + y^2 = z}, [x, y, z])

                               2 1/2                        2 1/2
               {[x = (z - 1.0 y )   ], [x = - 1.0 (z - 1.0 y )   ]}

#### Example 3

The ordering of the unknowns determines the representation of the solution, if the solution set is not finite. First, the following equation is solved for `x` leaving `y` as a free parameter:

    >> numeric::polysysroots({x^3 = y^2}, [x, y])

                  2 1/3                                   2 1/3
          {[x = (y )   ], [x = - (0.5 + 0.8660254038 I) (y )   ],

                                             2 1/3
             [x = - (0.5 - 0.8660254038 I) (y )   ]}

Reordering the unknowns leads to a representation with `x` as a free parameter:

    >> numeric::polysysroots({x^3 = y^2}, [y, x])

                               3 1/2                3 1/2
                       {[y = (x )   ], [y = - 1.0 (x )   ]}

#### Example 4

The symbolic solver produces a `RootOf` solution of the following system:

    >> eqs := {y^2 - y = x, x^3 = y^3 + x}:

    >> solve(eqs, BackSubstitution = FALSE)

                                 2
          {[x = 0, y = 0], [x = y  - y, y =

                        3       2            4     5
             RootOf(3 X3  - 2 X3  - X3 - 3 X3  + X3  + 1, X3)]}

Internally, `numeric::polysysroots` calls `solve` and processes this result numerically:

    >> numeric::polysysroots(eqs, [x, y])

          {[x = -0.238328984, y = 0.6080324762], [x = 0.0, y = 0.0],

             [x = 0.8911259621, y = -0.5682349751],

             [x = 2.237302267, y = 2.077118343],

             [x = - 1.445049623 + 0.1279930535 I,

             y = 0.441542078 - 1.094745154 I],

             [x = - 1.445049623 - 0.1279930535 I,

             y = 0.441542078 + 1.094745154 I]}

    >> delete eqs:

#### Example 5

The following equation is solved for the first of the specified unknowns:

    >> eqs := y^5 - PI*y = x: solve(eqs, [x, y])

                                        5
                                 {[x = y  - y PI]}

`numeric::polysysroots` processes this output numerically:

    >> numeric::polysysroots(eqs, [x, y])

                                   5
                            {[x = y  - 3.141592654 y]}

The equation is solved for `y` when the unknowns are reordered. However, no simple representation of the solution exists, so a `RootOf` object is returned:

    >> solve(eqs, [y, x])

                                                   5
                        {[y = RootOf(x + X5 PI - X5 , X5)]}

The roots represented by the `RootOf` expression cannot be computed numerically, because the symbolic parameter `x` is involved:

    >> numeric::polysysroots(eqs, [y, x])

                                                         5
                 {[y = RootOf(x + 3.141592654 X6 - 1.0 X6 , X6)]}

    >> delete eqs:

#### Changes

- new function

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