# Why doesn't Roots work on a certain quartic polynomial equation?

Bug introduced in 9.0 or earlier and fixed in 10.0

In everything that follows I am using Roots to get the solutions to the equations.

I start off with this equation: $$\frac{A}{3}x^4 - x^2 + 2ax - b^2 = 0$$ Now Roots gives me back $4$ complicated looking solutions to this.

Now if one stares at the solutions it appears that if one tunes to $a =b = \sqrt{\frac{3}{16 A}}$ then these 4 solutions collapse to $3$ namely, $x = \frac{1}{2}\sqrt{\frac{3}{A}}, -\frac{1}{2}\sqrt{\frac{3}{A}} \pm \sqrt{\frac{6}{A} }$

however

• when I from the beginning itself set $a = b = \sqrt{\frac{3}{16 A}}$ and run Roots on the equation the system just hangs!

• And of these $3$ solutions the "$\pm$" doesn't seem to satisfy the original equation. (with $a = b = \sqrt{\frac{3}{16 A}}$)

Can someone kindly explain what is happening?

• Including the actual code you use in the question eliminates ambiguities. – Yves Klett Jun 26 '14 at 12:39
• It's a bug. Will file a report and investigate. Thanks for the example. – Daniel Lichtblau Jun 26 '14 at 16:09

## Exposition of the problem

Roots works well on the following equation (obviously it yields generic formulae and it doesn't deal with exceptional cases e.g. when A == 0 qurtics expressions are not valid):

eq = A/3 x^4 - x^2 + 2 a x - b^2 == 0;


while in a special case of eq it fails (however in Mathematica 9 the system doesn't hang) returning the following warning:

eq1 = eq /. {a -> Sqrt[3/(16 A)], b -> Sqrt[3/(16 A)]}

-(3/(16 A)) + 1/2 Sqrt[3] Sqrt[1/A] x - x^2 + (A x^4)/3 == 0

Roots[ eq1, x]

\$IterationLimit::itlim: Iteration limit of 4096 exceeded. >>

Hold[ Roots[ -9 Sqrt[1/A] + 24 Sqrt[3] x - (48 x^2)/Sqrt[1/A]
+ (16 x^4)/(1/A)^(3/2) == 0, x]]


## Remedy

We can rewrite eq1 with:

Simplify[ eq1, A > 0 && x != 0]

24 Sqrt[3] Sqrt[A] x + 16 A^2 x^4 == 9 + 48 A x^2


now Roots works well:

TraditionalForm @ Roots[ Simplify[ eq1, A > 0 && x != 0], x]


## Explanation

In the latter case (of eq1) the system automatically rewrites the equation to another form not purely polynomial (the constants involving a square root of A are treated like another variables) involving more solutions than ordinary polynomials.

## Discussion and subtleties

In general solving symbolically equations we encounter various problems related to basic concepts, space of variables, range of interesting solutions, constants, etc. Here we tackle with an issue of constants and when A is not restrected somehow Roots is not able to provide a satisfactory solution, even tough with a more general case it works reasonably (yielding generically correct answers). There are another functions like Reduce which can provide full information on the solution space.

Whenever one doesn't know why an equation cannot be solved with one method it is reasonable to try another ways, and since symbolic capabilities for solving equations are related behind the scenes one should try more reliable aproaches. Roots basically works on polynomial equations, while Reduce and Solve can tackle with transcendental ones, moreover Reduce is more capable to provide symbolic results (an example: How to get intersection values from a parametric graph?). These issues were discussed more extensively in another post What is the difference between Reduce and Solve?.

Let's provide results from more reliable functions. Solve was updated in version 8 and now it can tackle with much larger class of equations, however it fails:

Solve[ eq1, x]

{ ToRules[ Hold[ Roots[-9 Sqrt[1/A] + 24 Sqrt[3] x - (48 x^2)/Sqrt[1/A]
+ (16 x^4)/(1/A)^(3/2) == 0, x]]]}


on the other hand Reduce shows what happens behind:

Reduce[ eq1, x]

 Reduce::useq: The answer found by Reduce contains unsolved equation(s)
{ 0 == (-1-Sqrt[1/A] Sqrt[A])/Sqrt[A], 0 == (-1-Sqrt[1/A] Sqrt[A])/Sqrt[A],
0 == (-1-Sqrt[1/A] Sqrt[A])/Sqrt[A], 0==(-1-Sqrt[1/A] Sqrt[A])/Sqrt[A],
0 == (1-Sqrt[1/A] Sqrt[A])/Sqrt[A], 0 == (1-Sqrt[1/A] Sqrt[A])/Sqrt[A],
0 == (1-Sqrt[1/A] Sqrt[A])/Sqrt[A], 0 == (1-Sqrt[1/A] Sqrt[A])/Sqrt[A]}.

A likely reason for this is that the solution set depends on branch cuts
of Mathematica functions. >>


and then it yields also solutions in terms of

Root[ 9/A^2 - 24 (1/A)^(3/2) ReduceCADAlgVar[1] #1 + (48 #1^2)/A - 16 #1^4 &, k]


This information clarifies what the problem is, namely the equation is transformed to another form which is not a pure polynomial one and such an equation has more solutions than an ordinary polynomial, see e.g.

Simplify[ eq1]

8 x (3 Sqrt[3] Sqrt[1/A] - 6 x + 2 A x^3) == 9/A


More on this issue one can find in another post which I linked above. But for a quick reference one should realize that Solve in similar problems yields this warning:

Solve::fdimc: When parameter values satisfy the condition A ∈ Reals,
the solution set contains a full-dimensional component;
use Reduce for complete solution information. >>


This obviously happens when we deal with symbolic equations, however choosing specific numeric coefficients it works well. Ahother way is to supplement the equation eq1 by a condition e.g. A > 0, however Roots is not capable to deal with boolean formulae, nonetheless it can be done with Solve or Reduce (see e.g. this answer Solve an equation in R+), e.g.:

Solve[ eq1 && A > 0, x]

{{ x -> ConditionalExpression[1/2 Sqrt[3] Sqrt[1/A], A > 0]},
{ x -> ConditionalExpression[ Root[9 - 72 A #1^2 + 16 A^2 #1^4 &, 1], A > 0]},
{ x -> ConditionalExpression[Root[9 - 72 A #1^2 + 16 A^2 #1^4 &, 3], A > 0]}}


In conclusion we can say that the underlying issue is not a desirable feature of the system(it works on more general ones while fails on special cases), however Roots is not intended to tackle with exceptional cases. Moreover we can expect that similar issues will be present in next versions, since there can always appear some problems related to appropriate restricting constants in algebraic equations.

    Clear[A, a, b, x];
eq1 = A/3*x^4 - x^2 + 2*a*x - b^2 == 0;


Now, as you propose:

    eq2 = Simplify[eq1 /. {a -> Sqrt[3/(16 A)], b -> Sqrt[3/(16 A)]},
A > 0]

(*   24 Sqrt[3] Sqrt[A] x + 16 A^2 x^4 == 9 + 48 A x^2  *)


With this one finds:

    Solve[eq2, x]

(*  {{x -> Sqrt[3]/(2 Sqrt[A])}, {x -> Sqrt[3]/(
2 Sqrt[A])}, {x -> (-Sqrt[3] - Sqrt[6])/(
2 Sqrt[A])}, {x -> (-Sqrt[3] + Sqrt[6])/(2 Sqrt[A])}}  *)


Exactly as you expect to obtain. So, what?

• Thanks for your efforts. (1) When I set the values to those special points and evaluated "roots" on it the system hung up. Did it happen for you too? (2) And aren't you also seeing that those last 2 of your roots are not a solution of the equation? – user6818 Jun 26 '14 at 12:38
• @user6818 1) No the solution has been returned immediately. 2) eq2 /. x -> (-Sqrt[3] + Sqrt[6])/(2 Sqrt[A])returns True as well as eq2 /. x -> (-Sqrt[3] - Sqrt[6])/(2 Sqrt[A]). Thus, no, I do not see that. – Alexei Boulbitch Jun 26 '14 at 14:54