Advanced Search
All Posts : Applications, Examples and Libraries

Recently @salim-barzani asked a question about a paper that involved analysing the different types of roots of polynomials. The appendix in that paper gave the example of the roots of x^4+x^2*e[2]+x*e[1]+e[0]using the analysis in Lu et al, "A complete discrimination system for polynomials", Science in China (Ser. E), 39 (1996) 628-646. The analysis uses the discriminant sequence and extensions. Maple provides this through RegularChains:-ParametricSystemTools:-DiscrminantSequence. For example for this polynomial we find there is a real root of multiplicity 2 and a complex conjugate pair when D__2*D__3 < 0 and D__4 = 0 where the D__i are the ith entries in the discriminant sequence [1, -e[2], -2*e[2]^3+8*e[0]*e[2]-9*e[1]^2, 16*e[0]*e[2]^4-4*e[1]^2*e[2]^3-128*e[0]^2*e[2]^2+144*e[0]*e[1]^2*e[2]-27*e[1]^4+256*e[0]^3].

 

The problem with these conditions is that they are in terms of the D__i and not directly in terms of the e__i parameters. One can derive these conditions and then solve them to find the conditions on the parameters, but Maple has various routines in the RegularChains, RootFinding:-Parametric and SolveTools packages that directly find conditions on parameters to find when there are specified numbers of real or complex roots for polynomial systems. So this post is my attempt to use these tools to find the conditions on the parameters of the above polynomial that give various types of roots. One immediate difficulty is that generally these routines count distinct roots irrespective of multiplicity, and so some indirect analysis is required. There are several different types of commands and analyses that could be used, and my choices here are more to do with my learning experience than an optimum analysis.

 

The first conclusion is that it is possible, although RealComprehensiveTriangularize did not work as I expected when asking for zero real roots (bug?). Assuming it had worked, RealComprehensiveTriangularize could cover all the cases here, though that will not be true when there are more parameters. There doesn't seem to be an obvious systematic way of doing this analysis, which is a downside. Another downside is the large number of subcase conditions, which look as if they could be combined into fewer subcases. CellDecomposition works well for cases without multiplicity.


Main worksheet [not all is displaying below]:

Download RootAnalysis4.mw

> 

restart
> 

with(RegularChains); with(ParametricSystemTools); with(RootFinding:-Parametric)

Consider the following polynomial in x with the three real parameters e[0], e[1], e[2]. We would like to know the conditions on these parameters that lead to different numbers of real and complex-conjugate pairs of roots of different multiplicities.

> 

p := x^4+x^2*e[2]+x*e[1]+e[0]

x^4+x^2*e[2]+x*e[1]+e[0]

Consider first how many cases there are. We can set this up as a combinatorial problem in the combstruct package.

> 

sys := {C = Atom, R = Atom, realrts = Set(multiplereal), rts = Prod(realrts, complexrts), complexpr = Prod(C, C), complexrts = Set(multiplecomplex), multiplecomplex = Sequence(complexpr, card > 0), multiplereal = Sequence(R, card > 0)}
> 

Draw := proc (q) options operator, arrow; eval(q, {Epsilon = NULL, Prod = `[]`, Set = (proc () options operator, arrow; args end proc), Sequence = `*`}) end proc

For a degree 4 polynomial there are 9 different cases to consider. Here [C, C] means a (non-real) complex-conjugate pair of roots and R means a real root; the exponents indicate the multiplicities.

> 

all := combstruct:-allstructs([rts, sys], size = degree(p, x)); nops(%); `~`[Draw](all)

9

[[[C, C]^2], [[C, C], [C, C]], [R^2, [C, C]], [R^4], [R, R, [C, C]], [R^2, R^2], [R^3, R], [R, R, R^2], [R, R, R, R]]

These are (in order)
(a) A duplicate pair of complex-conjugate roots

(b) Two distinct pairs of complex-conjugate roots

(c) A real root of multiplicity 2 and a pair of complex-conjugate roots

(d) A real root of multiplicity 4

(e) Two distinct real roots of multiplicity 1 and a complex-conjugate pair

(f) Two distinct real roots each of multiplicity 2

(g) A real root of multiplicity 3 and a real root of multiplicity 1

(h) Three distinct real roots, of multiplicities 2, 1, and 1

(i) Four distinct real roots of multiplicity 1

Declare the variables first and parameters last. np is the number of parameters. Use the suggested order.

> 

vp := SuggestVariableOrder([p = 0], [x]); R := PolynomialRing(vp); np := nops(`minus`({vp[]}, {x}))

[x, e[0], e[1], e[2]]

polynomial_ring

3

Define derivative polynomials. p2 = 0 when there is a root of multiplicity 2 or more; p3 = 0 when there is a root of multiplicity 3 or more and p4 = 0when there is a root of multiplicity 4.

> 

p2 := diff(p, x); p3 := diff(p2, x); p4 := diff(p3, x)

4*x^3+2*x*e[2]+e[1]

12*x^2+2*e[2]

24*x

Discriminant is zero if and only if there are repeated roots.

> 

Delta := discrim(p, x)

16*e[0]*e[2]^4-4*e[1]^2*e[2]^3-128*e[0]^2*e[2]^2+144*e[0]*e[1]^2*e[2]-27*e[1]^4+256*e[0]^3

(d) A real root of multiplicity 4

 

This is perhaps the simplest case and can be done using RealComprehensiveTriangularize. Specifying np = 3 means use the last 3 variables in the PolynomialRing as parameters. The argument 1 means we want the cases where there is one distinct real root. We specify that all of p, p2, p3, p4 are zero so that the 1 real root is the common root of these polynomials, i.e., is a root on multiplicity 4. (I find the cadcell output a little easier to use, but it is not critical.)

> 

rct := RealComprehensiveTriangularize({p = 0, p2 = 0, p3 = 0, p4 = 0}, np, R, 1, output = cadcell); Display(rct, R)

PiecewiseTools:-Is, "Wrong kind of parameters in piecewise"

This means that the conditions on the parameters to get a single real root of multiplicity 4 are

> 

conds := Info(rct[2][1][1], R)

[e[2] = 0, e[1] = 0, e[0] = 0]

and that the polynomial to solve to find this root under these conditions is

> 

poly := Info(rct[1][][2], R)

[x = 0]

which we can check:

> 

eval(p, conds); solve(%, x)

x^4

0, 0, 0, 0

(g) A real root of multiplicity 3 and a real root of multiplicity 1

 

For multiplicity 3 we need p, p2, p3 all zero ,but p4 nonzero to exclude multiplicity 4. We want there to be one real root meeting these conditions. The other root must be real with multiplicity 1.

> 

rct := RealComprehensiveTriangularize({p = 0, p2 = 0, p3 = 0, p4 <> 0}, np, R, 1, output = cadcell); Display(rct, R)

PiecewiseTools:-Is, "Wrong kind of parameters in piecewise"

There are two conditions, differing only in the sign of e[1]

> 

conds1 := Info(rct[2][1][1], R); eqs1, ineqs1 := selectremove(type, conds1, `=`); conds2 := Info(rct[2][2][1], R); eqs2, ineqs2 := selectremove(type, conds2, `=`)

[e[2] < 0, e[1] = -(2/9)*(-6*e[2]^3)^(1/2), e[0] = -(1/12)*e[2]^2]

[e[1] = -(2/9)*(-6*e[2]^3)^(1/2), e[0] = -(1/12)*e[2]^2], [e[2] < 0]

[e[2] < 0, e[1] = (2/9)*(-6*e[2]^3)^(1/2), e[0] = -(1/12)*e[2]^2]

[e[1] = (2/9)*(-6*e[2]^3)^(1/2), e[0] = -(1/12)*e[2]^2], [e[2] < 0]

The polynomial to be solved to find the root of multiplicity 3 is

> 

poly := Info(rct[1][][2], R)

[4*e[2]*x+3*e[1] = 0]

so the root of multiplicity under the two sets of conditions is (for any negative e[2])

> 

`assuming`([simplify(solve(eval(poly, eqs1), x))], [ineqs1[]]); `assuming`([simplify(solve(eval(poly, eqs2), x))], [ineqs2[]])

{x = -(1/6)*6^(1/2)*(-e[2])^(1/2)}

{x = (1/6)*6^(1/2)*(-e[2])^(1/2)}

which we can check in the original polynomial

> 

`assuming`([simplify([solve(eval(p, eqs1), x, explicit)])], [ineqs1[]]); `assuming`([simplify([solve(eval(p, eqs2), x, explicit)])], [ineqs2[]])

[(1/2)*6^(1/2)*(-e[2])^(1/2), -(1/6)*6^(1/2)*(-e[2])^(1/2), -(1/6)*6^(1/2)*(-e[2])^(1/2), -(1/6)*6^(1/2)*(-e[2])^(1/2)]

[(1/6)*6^(1/2)*(-e[2])^(1/2), (1/6)*6^(1/2)*(-e[2])^(1/2), (1/6)*6^(1/2)*(-e[2])^(1/2), -(1/2)*6^(1/2)*(-e[2])^(1/2)]

 

(f) Two distinct real roots each of multiplicity 2

 

This can be done similarly. For multiplicity 2 we need p, p2 both zero ,but p3, p4 nonzero to exclude multiplicity 3 or 4. There needs to be two such roots.

> 

rct := RealComprehensiveTriangularize({p = 0, p2 = 0, p3 <> 0, p4 <> 0}, np, R, 2, output = cadcell); Display(rct, R)

PiecewiseTools:-Is, "Wrong kind of parameters in piecewise"

The condition, and the polynomial to find the roots of multiplicity 2 are

> 

conds := Info(rct[2][1][1], R); eqs, ineqs := selectremove(type, conds, `=`); poly := Info(rct[1][][2], R)

[e[2] < 0, e[1] = 0, e[0] = (1/4)*e[2]^2]

[e[1] = 0, e[0] = (1/4)*e[2]^2], [e[2] < 0]

[2*x^2+e[2] = 0]

so the roots of multiplicity 2 under the condition are (for any negative e[2])

> 

`assuming`([simplify([solve(eval(poly, eqs), x)])], [ineqs[]])

[{x = (1/2)*2^(1/2)*(-e[2])^(1/2)}, {x = -(1/2)*2^(1/2)*(-e[2])^(1/2)}]

which we can check in the original polynomial

> 

`assuming`([simplify([solve(eval(p, eqs), x, explicit)])], [ineqs[]])

[(1/2)*2^(1/2)*(-e[2])^(1/2), -(1/2)*2^(1/2)*(-e[2])^(1/2), (1/2)*2^(1/2)*(-e[2])^(1/2), -(1/2)*2^(1/2)*(-e[2])^(1/2)]

(i) Four distinct real roots of multiplicity 1

 

This is also a simple case; we just ask for p = 0 to have 4 distinct real roots.

> 

rct := RealComprehensiveTriangularize({p = 0}, np, R, 4, output = cadcell); Display(rct, R)

PiecewiseTools:-Is, "Wrong kind of parameters in piecewise"

Of course the polynomial to solve is just the original one.

> 

poly := Info(rct[1][][2], R)

[x^4+e[2]*x^2+e[1]*x+e[0] = 0]

We recognize the discriminant appearing here in the first and third conditions, where e[0] is chosen so that it is not zero to ensure there are repeated roots; there are additional conditions to distinguish this from cases (b) and (e).

> 

solve(Delta, e[0], explicit = false)

RootOf(256*_Z^3-128*e[2]^2*_Z^2+(16*e[2]^4+144*e[1]^2*e[2])*_Z-4*e[1]^2*e[2]^3-27*e[1]^4)

When e[1] = 0, (as in the second condition), the discriminant is

> 

eval(Delta, e[1] = 0); solve(%, e[0])

16*e[0]*e[2]^4-128*e[0]^2*e[2]^2+256*e[0]^3

0, (1/4)*e[2]^2, (1/4)*e[2]^2

With a little thought, the three conditions may be combined into one rule:
Choose any negative e[2]. Find a value of e[1] between -2*sqrt(-6*e[2]^3)*(1/9)and 2*sqrt(-6*e[2]^3)*(1/9). Then find a value of e[2] between the first and second real solutions of the discriminant.

For cases where the polynomials to solve do not have multiple roots (as here), an alternative to RealComprehensiveTriangularize is CellDecomposition in the RootFinding:-Parametric package. The same information about the conditions is given in a more cryptic form using CellDescription for the cells with 4 real roots (cell numbers 4 and 8) from CellDecomposition. The two cells differ in the sign of e[1].

> 

cells := CellDecomposition([p = 0], [x], [e[0], e[1], e[2]]); NumberOfSolutions(cells); cells4 := CellsWithSolutions(cells, 4); CellDescription(cells, cells4[1]); CellDescription(cells, cells4[2])

[[1, 2], [2, 0], [3, 2], [4, 4], [5, 2], [6, 0], [7, 2], [8, 4], [9, 2], [10, 0], [11, 2], [12, 0], [13, 2], [14, 0], [15, 2], [16, 0]]

[4, 8]

[[-infinity, 0, e[2], e[2], 1], [8*e[2]^3+27*e[1]^2, 1, e[1], e[1], 1], [16*e[0]*e[2]^4-4*e[1]^2*e[2]^3-128*e[0]^2*e[2]^2+144*e[0]*e[1]^2*e[2]-27*e[1]^4+256*e[0]^3, 1, e[0], 16*e[0]*e[2]^4-4*e[1]^2*e[2]^3-128*e[0]^2*e[2]^2+144*e[0]*e[1]^2*e[2]-27*e[1]^4+256*e[0]^3, 2]]

[[-infinity, 0, e[2], e[2], 1], [e[1], 1, e[1], 8*e[2]^3+27*e[1]^2, 2], [16*e[0]*e[2]^4-4*e[1]^2*e[2]^3-128*e[0]^2*e[2]^2+144*e[0]*e[1]^2*e[2]-27*e[1]^4+256*e[0]^3, 1, e[0], 16*e[0]*e[2]^4-4*e[1]^2*e[2]^3-128*e[0]^2*e[2]^2+144*e[0]*e[1]^2*e[2]-27*e[1]^4+256*e[0]^3, 2]]

The SamplePoints and SampleSolutions commands automate finding a valid solution in a cell, which exhibits the 4 distinct real roots.

> 

cells:-SamplePoints[4]; SampleSolutions(cells, %); cells:-SamplePoints[8]; SampleSolutions(cells, %)

[e[0] = 0, e[1] = -2570530724767170816689/9444732965739290427392, e[2] = -1]

[[x = -.8164965809], [x = -.2988584907], [x = 0.], [x = 1.115355072]]

[e[0] = 0, e[1] = 20564245798137366533639/75557863725914323419136, e[2] = -1]

[[x = -1.115355072], [x = 0.], [x = .2988584907], [x = .8164965809]]

 

 

(h) Three distinct real roots, of multiplicities 2, 1, and 1

 

This is the only case with three distinct real roots, so it can also be done with RealComprehensiveTriangularize similarly to the previous case.

> 

rct := RealComprehensiveTriangularize({p = 0}, np, R, 3, output = cadcell); Display(rct, R)

PiecewiseTools:-Is, "Wrong kind of parameters in piecewise"

The conditions on e[2] and e[1] are the same as those of case (i), but now e[0] is chosen to make the discriminant zero, since there are repeated roots.
Choose e[2] < 0, e.g., e[2] = -1, and then -2*sqrt(-6*e[2]^3)*(1/9) < e[1] and e[1] < 2*sqrt(-6*e[2]^3)*(1/9), e.g., e[1] = (1/2)*(2*sqrt(-6*e[2]^3)*(1/9))

> 

params := {e[2] = -1}; params := `union`(params, {e[1] = eval((1/2)*(2*sqrt(-6*e[2]^3)*(1/9)), params)})

{e[2] = -1}

{e[1] = (1/9)*6^(1/2), e[2] = -1}

Now choose e[0] to make the discriminant zero

> 

fsolve(eval(Delta, params), e[0], digits = 20); params := `union`(params, {e[0] = %[1]})

-0.1889013308e-1, 0.6793395576e-1, .4509561773

{e[0] = -0.1889013308e-1, e[1] = (1/9)*6^(1/2), e[2] = -1}

The three distinct roots are

> 

poly := Info(rct[1][][2], R); eval(poly, params); [fsolve(%[], complex)]

[((512*e[2]^5+6912*e[1]^2*e[2]^2)*e[0]^2+(-256*e[2]^7-2304*e[2]^4*e[1]^2-7776*e[2]*e[1]^4)*e[0]+32*e[2]^9+400*e[1]^2*e[2]^6+1728*e[1]^4*e[2]^3+2916*e[1]^6)*x^3+((-256*e[2]^4*e[1]+6912*e[2]*e[1]^3)*e[0]^2+(128*e[2]^6*e[1]+576*e[2]^3*e[1]^3-3888*e[1]^5)*e[0]-16*e[1]*e[2]^8-192*e[1]^3*e[2]^5-648*e[1]^5*e[2]^2)*x^2+((512*e[2]^6+4992*e[1]^2*e[2]^3+5184*e[1]^4)*e[0]^2+(-256*e[2]^8-1856*e[2]^5*e[1]^2-4320*e[2]^2*e[1]^4)*e[0]+32*e[2]^10+408*e[1]^2*e[2]^7+1692*e[1]^4*e[2]^4+2430*e[1]^6*e[2])*x+(256*e[2]^5*e[1]+8640*e[2]^2*e[1]^3)*e[0]^2+(-128*e[2]^7*e[1]-1440*e[2]^4*e[1]^3-7776*e[2]*e[1]^5)*e[0]+16*e[1]*e[2]^9+204*e[1]^3*e[2]^6+972*e[1]^5*e[2]^3+2187*e[1]^7 = 0]

[-13.08460370*x^3-.7573725292*6^(1/2)*x^2+12.82157099*x-.7116973772*6^(1/2) = 0]

[-1.121474463, .1417831433, .8379081760]

And the original polynomial has these 3 roots, with one repeated

> 

eval(p, params); simplify(fnormal([fsolve(%, complex)], 7), zero)

x^4-x^2+(1/9)*x*6^(1/2)-0.1889013308e-1

[-1.121474, .1417831, .1417831, .8379082]

(e) Two distinct real roots of multiplicity 1 and a complex-conjugate pair

 

For this case we have no repeated roots (discriminant nonzero) and two distinct real roots

> 

rct := RealComprehensiveTriangularize({p = 0, Delta <> 0}, np, R, 2, output = cadcell); Display(rct, R)

PiecewiseTools:-Is, "Wrong kind of parameters in piecewise"

There are two types of solutions. (i) Those with e[2] < 0 are similar to those before, but now e[1] can be on the boundaries: -2*sqrt(-6*e[2]^3)*(1/9) <= e[1] and e[1] <= 2*sqrt(-6*e[2]^3)*(1/9) and e[2] is variously related to the real roots of the discriminant. (ii) e[2] >= 0, e[1] is anything and e[0] is less than the first real root of the discriminant.
Similar information is available from CellDecompostion for the cells with two distinct real solutions. CellDecomposition does not deal with values on the boundary. It also does not deal with multiple roots so all the cells with two roots are for this case, case (e).

> 

NumberOfSolutions(cells); cells2 := CellsWithSolutions(cells, 2)

[[1, 2], [2, 0], [3, 2], [4, 4], [5, 2], [6, 0], [7, 2], [8, 4], [9, 2], [10, 0], [11, 2], [12, 0], [13, 2], [14, 0], [15, 2], [16, 0]]

[1, 3, 5, 7, 9, 11, 13, 15]

> 

for i in cells2 do CellDescription(cells, i) end do

[[-infinity, 0, e[2], e[2], 1], [-infinity, 0, e[1], 8*e[2]^3+27*e[1]^2, 1], [-infinity, 0, e[0], 16*e[0]*e[2]^4-4*e[1]^2*e[2]^3-128*e[0]^2*e[2]^2+144*e[0]*e[1]^2*e[2]-27*e[1]^4+256*e[0]^3, 1]]

[[-infinity, 0, e[2], e[2], 1], [8*e[2]^3+27*e[1]^2, 1, e[1], e[1], 1], [-infinity, 0, e[0], 16*e[0]*e[2]^4-4*e[1]^2*e[2]^3-128*e[0]^2*e[2]^2+144*e[0]*e[1]^2*e[2]-27*e[1]^4+256*e[0]^3, 1]]

[[-infinity, 0, e[2], e[2], 1], [8*e[2]^3+27*e[1]^2, 1, e[1], e[1], 1], [16*e[0]*e[2]^4-4*e[1]^2*e[2]^3-128*e[0]^2*e[2]^2+144*e[0]*e[1]^2*e[2]-27*e[1]^4+256*e[0]^3, 2, e[0], 16*e[0]*e[2]^4-4*e[1]^2*e[2]^3-128*e[0]^2*e[2]^2+144*e[0]*e[1]^2*e[2]-27*e[1]^4+256*e[0]^3, 3]]

[[-infinity, 0, e[2], e[2], 1], [e[1], 1, e[1], 8*e[2]^3+27*e[1]^2, 2], [-infinity, 0, e[0], 16*e[0]*e[2]^4-4*e[1]^2*e[2]^3-128*e[0]^2*e[2]^2+144*e[0]*e[1]^2*e[2]-27*e[1]^4+256*e[0]^3, 1]]

[[-infinity, 0, e[2], e[2], 1], [e[1], 1, e[1], 8*e[2]^3+27*e[1]^2, 2], [16*e[0]*e[2]^4-4*e[1]^2*e[2]^3-128*e[0]^2*e[2]^2+144*e[0]*e[1]^2*e[2]-27*e[1]^4+256*e[0]^3, 2, e[0], 16*e[0]*e[2]^4-4*e[1]^2*e[2]^3-128*e[0]^2*e[2]^2+144*e[0]*e[1]^2*e[2]-27*e[1]^4+256*e[0]^3, 3]]

[[-infinity, 0, e[2], e[2], 1], [8*e[2]^3+27*e[1]^2, 2, e[1], infinity, 0], [-infinity, 0, e[0], 16*e[0]*e[2]^4-4*e[1]^2*e[2]^3-128*e[0]^2*e[2]^2+144*e[0]*e[1]^2*e[2]-27*e[1]^4+256*e[0]^3, 1]]

[[e[2], 1, e[2], infinity, 0], [-infinity, 0, e[1], e[1], 1], [-infinity, 0, e[0], 16*e[0]*e[2]^4-4*e[1]^2*e[2]^3-128*e[0]^2*e[2]^2+144*e[0]*e[1]^2*e[2]-27*e[1]^4+256*e[0]^3, 1]]

[[e[2], 1, e[2], infinity, 0], [e[1], 1, e[1], infinity, 0], [-infinity, 0, e[0], 16*e[0]*e[2]^4-4*e[1]^2*e[2]^3-128*e[0]^2*e[2]^2+144*e[0]*e[1]^2*e[2]-27*e[1]^4+256*e[0]^3, 1]]

Verify that all these cells have two distinct real roots and a complex-conjugate pair

> 

for i in cells2 do fsolve(eval(p, cells:-SamplePoints[i]), complex) end do

-.8774388331-.7448617666*I, -.8774388331+.7448617666*I, .7548776662, 1.

-1.208211169, -0.6103890665e-1-.7864141726*I, -0.6103890665e-1+.7864141726*I, 1.330288982

-.7135149791-.3232096177*I, -.7135149791+.3232096177*I, .4198192076, 1.007210751

-1.330288982, 0.6103890665e-1-.7864141726*I, 0.6103890665e-1+.7864141726*I, 1.208211169

-1.007210751, -.4198192076, .7135149791-.3232096177*I, .7135149791+.3232096177*I

-1., -.7548776662, .8774388331-.7448617666*I, .8774388331+.7448617666*I

-.3627365530-1.149253336*I, -.3627365530+1.149253336*I, .1226790229, .6027940830

-.6027940830, -.1226790229, .3627365530-1.149253336*I, .3627365530+1.149253336*I

   

 
polynomial rootfinding parametric celldecomposition regularchains

Please Wait...