Advanced Search

Featured Post
dharr
9207

Characterizing polynomial roots with RegularChains and related routines

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 (see cases (a) and (b)) (bug?). Assuming it had worked, RealComprehensiveTriangularize could cover all the cases here, though that will not be true for higher-degree polynomials with 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 displayed 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

   

(f) Two distinct real roots each of multiplicity 2

   

(i) Four distinct real roots of multiplicity 1

   

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

   

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

   

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

   

(b) Two distinct pairs of complex-conjugate roots

   

(a) A duplicate pair of complex-conjugate roots

 

Here we want multiplicity 2 but no real roots. The discriminant is expected to be zero, but in most cases is not.

> 

rct := RealComprehensiveTriangularize({p = 0, p2 = 0, p3 <> 0, p4 <> 0}, np, R, 0, output = cadcell); nops(rct[2]); Display(rct[2][1 .. 4], R)

40

[[PIECEWISE([e[0] < 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, index = real[1]), ``], [e[1] < -(2/9)*(-6*e[2]^3)^(1/2), ``], [e[2] < 0, ``]), []], [PIECEWISE([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, index = real[1]) < e[0], ``], [e[1] < -(2/9)*(-6*e[2]^3)^(1/2), ``], [e[2] < 0, ``]), []], [PIECEWISE([e[0] < -(1/12)*e[2]^2, ``], [e[1] = -(2/9)*(-6*e[2]^3)^(1/2), ``], [e[2] < 0, ``]), []], [PIECEWISE([e[0] = -(1/12)*e[2]^2, ``], [e[1] = -(2/9)*(-6*e[2]^3)^(1/2), ``], [e[2] < 0, ``]), []]]

Consider one of the cases with an equality for e[0], suggesting a zero discriminant.
Find a sample point satisfying the conditions, and see what the roots are like

> 

j := 37; cell := rct[2][j][1]; conds := Info(cell, R); pts := Info(SamplePoints(cell, R), R)[1]; `~`[is](eval(conds, pts))

37

cad_cell

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

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

[true, true, true]

We find two complex roots of multiplicity 2, which are complex conjugates, as expected

> 

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

x^4+(1/2)*x^2+1/16

(1/2)*I, -(1/2)*I, (1/2)*I, -(1/2)*I

Check that discriminant is zero.

> 

eval(Delta, pts)

0

But, for example, the first cell does not have the discriminant zero, and does not give a correct result.

> 

j := 1; cell := rct[2][j][1]; conds := Info(cell, R); pts := Info(SamplePoints(cell, R), R)[1]; `~`[is](eval(conds, pts))

1

cad_cell

[e[2] < 0, e[1] < -(2/9)*(-6*e[2]^3)^(1/2), e[0] < 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, index = real[1])]

[e[0] = 1/2, e[1] = -3/2, e[2] = -1/2]

[true, true, true]

We find a complex conjugate pair and two distinct real roots, which is not expected for this case.

> 

eval(p, pts); fsolve(%, x, complex)

x^4-(1/2)*x^2-(3/2)*x+1/2

-.7473459056-.9001675303*I, -.7473459056+.9001675303*I, .3077440660, 1.186947745

The discriminant is indeed nonzero.

> 

eval(Delta, pts)

-3073/16

We did not yet find the conditions for this case. We can ask for the conditions for different numbers of complex roots (complex in this context includes real).

> 

cmplx := ComplexRootClassification([p], np, R)

[[constructible_set, 1], [constructible_set, 2], [constructible_set, 3], [constructible_set, 4]]

We are interested in the conditions for exactly two distinct roots, which is found from the second constructible_set. There are two subcases.

> 

Display(cmplx[2][1], R)

PIECEWISE([12*e[0]+e[2]^2 = 0, ``], [27*e[1]^2+8*e[2]^3 = 0, ``], [e[2] <> 0, ``]), PIECEWISE([4*e[0]-e[2]^2 = 0, ``], [e[1] = 0, ``], [e[2] <> 0, ``])

Consider the second case

> 

case2 := Info(cmplx[2][1], R)[2]

[[4*e[0]-e[2]^2, e[1]], [e[2]]]

> 

solve(case2[1], {e[0], e[1]}); q1 := eval(p, %); rts2 := [solve(%, x)]

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

x^4+x^2*e[2]+(1/4)*e[2]^2

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

We saw this before when e[2]<0 as case (f) with real roots of multiplicity 2. Now, for e[2]>0 we indeed have two duplicate pairs of complex conjugate roots

> 

`assuming`([simplify(rts2)], [e[2] > 0])

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

Consider the first case

> 

case1 := Info(cmplx[2][1], R)[1]

[[12*e[0]+e[2]^2, 27*e[1]^2+8*e[2]^3], [e[2]]]

> 

ans1 := {solve(case1[1], {e[0], e[1]}, explicit)}; q11 := eval(p, ans1[1]); rts11 := [solve(%, x, explicit)]

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

x^4+x^2*e[2]-(2/9)*x*(-6*e[2])^(1/2)*e[2]-(1/12)*e[2]^2

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

The rts11 subcase polynomial q11 must have real coefficients and therefore only applies for e[2]<0. The roots are real with multiplicity 3 and multiplicity 1 and this case is just case (g) above. The rts12 subcase polynomial q12 (below) also requires e[2]<0 and corresponds to case (g), but with signs reversed.

> 

q12 := eval(p, ans1[2]); rts12 := [solve(%, x, explicit)]

x^4+x^2*e[2]+(2/9)*x*(-6*e[2])^(1/2)*e[2]-(1/12)*e[2]^2

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

Therefore the rts2 case is the solution for case (a); the cell37 example was a special case of this. In fact if we have a duplicate pair of complex-conjugate roots, the polynomial must be a pefect square, as we see it is

> 

factor(q1)

(1/4)*(2*x^2+e[2])^2

> 

NULL

Download RootAnalysis5.mw

 

June 08
5  0
polynomial rootfinding parametric
Read Full Post

Featured Post
pchin
1461

Animating a Polyhedron

A little while ago, I created a video, Engaging and Enlightening Students with Maple Visualizations, that showed a sample of Maple visualizations that would be helpful in teaching math. Doing that allowed me to get reacquainted with some of Maple's plotting features that I hadn't used for a while. As a result, I made a second instructional video for my Maple tips series, Animating a Polyhedron in Maple. 

I chose this topic because I thought it would show several features in Maple that might not be known to all users. I list them below and encourage you to try them out.

  • The plots:-polyhedraplot command allows you to create a 3-D plot of a polyhedron, including one of 138 polyhedra that Maple knows about.

  • The list of named polyhedra available can be obtained by calling the plots:-polyhedra_supported command.

  • The viewpoint option, which allows you to create an animation by varying the viewpoint through a 3D plot, can be used to rotate the polyhedron.

  • Finally, the Export feature allows you to save the plot animation as an animated GIF.

 

June 09
4  0
Read Full Post


MaplePrimes Questions
  • Recent
  • Unanswered
  • Maple
  • MapleSim
    •

    How can I get the correct inverse metric with...

    Maple 2025 asked by Hullzie16 448 June 23
    1  2
    linearize perturbation physics

    What if my permutation has a fixed point?

    Maple 2026 asked by wkehowski 30 June 23
    1  2
    permutation group-theory