# How to get the correct answer when solve this equation?

I want to find the values of the parameter $m$ to the graph of the funtion $$x^4 - 2 m x^2 + m + 1$$ has three extremal points make an equilateral triangle. I tried

Clear[f];
f[x_] := x^4 - 2 m x^2 + m + 1
sol = Solve[D[f[x], x] == 0, x]
t = x /. sol;
a = {t[[1]], f[x /. sol[[1]]]};
b = {t[[2]], f[x /. sol[[2]]]};
c = {t[[3]], f[x /. sol[[3]]]};
ab = EuclideanDistance[a, b];
bc = EuclideanDistance[b, c];
ac = EuclideanDistance[a, c];
Solve[{ab == bc, bc == ac, ac == ab}, m, Reals]


The correct answer is $m =\sqrt[3]{3}$. But, when I tried

Solve[{ab == bc, bc == ac, ac == ab, m > 0}, m, Reals]


Mathematica out put empty set. How to get the correct answer? And how to reduce my code?

-
Interesting... Try comparing Reduce[Sqrt[m^4 + Abs[m]] == 2 Sqrt[Abs[m]] && m > 0, m, Reals] and Reduce[Sqrt[m^4 + Abs[m]] == 2 Sqrt[Abs[m]] && m > 0, m ], the former one gives False while the latter one gives m == 3^(1/3) on my MMA 9.0.1. –  Silvia May 29 '13 at 3:50
@MichaelE2 I had quickly copied and pasted the first part of the text, which indeed gives that result. Does it do so when you include the m>0 constraint? It didn't for me. –  Jonathan Shock May 29 '13 at 4:45
@JonathanShock D'oh! You're exactly right. I'll remove my comment. –  Michael E2 May 29 '13 at 4:54
This works: Solve[{ab == bc, bc == ac, ac == ab, m > 0}, m] –  Michael E2 May 29 '13 at 4:57
@JonathanShock Math is simpler when it's Complex -- beyond that I'm not sure. It also works if EuclideanDistance is replaced by #.# &[#1 - #2] &, the square of distance formula -- polynomials are easier than algebraic functions. –  Michael E2 May 29 '13 at 5:04

Here is what I know about this. I do not know the answer to the question, but I have some guesses at an explanation. I also give some workarounds and a few simple, curious examples that illustrate some issues.

### Preliminaries

First of all, Solve uses different algorithms for different types of problems, some description of which can be found in the implementation notes and the reference page. For instance it uses different algorithms over Reals and Complexes; it uses different algorithms for polynomials, algebraic functions and transcendental functions. With an inequality, it will assume the variable is real but that function values can be complex, unless the domain is specified to be Reals. Solve also gives only generic solutions (from the reference page):

Solve gives generic solutions only. Solutions that are valid only when continuous parameters satisfy equations are removed. Other solutions that are only conditionally valid are expressed as ConditionalExpression objects.

Aside: As a function of a complex variable Abs[m] is not even algebraic, although it (accidentally?) is over the reals: Abs[m]^2 == m^2. It comes into the OP's system through EuclideanDistance, which can be seen by evaluating this:

{ab == bc, bc == ac, ac == ab, m > 0}

{Sqrt[Abs[m] + Abs[m]^4] == 2 Sqrt[Abs[m]],
2 Sqrt[Abs[m]] == Sqrt[Abs[m] + Abs[m]^4], True, m > 0}


The equations are algebraic as equations in Abs[m]. In terms of a complex variable m, they are transcendental (and not even locally holomorphic).

### Explanation?

This could point to an explanation. Here are two possibilities:

1) If the problem is viewed as primarily to solve the inequality m > 0 but subject to the constraints ab == bc etc., then none of the solutions would be generic.

2) If the problem is viewed as primarily to solve the equations ab == bc etc. but subject to the condition m > 0, then I would expect to get the desired solution.

The question is not what is the user's point of view, but how do the Solve algorithms approach the problem. In the first case, the base problem m > 0 has a one-dimensional solution space and the constraints imply there is no one-dimensional solution space, only some possibly non-generic (zero-dimensional) points left over. I don't know if this is the approach -- perhaps someone here does. (Frankly, it seems a little odd, if it is so.) Comparing the output of GenericCylindricalDecomposition and CylindricalDecomposition we see some support for the viewpoint 1).

GenericCylindricalDecomposition returns False (no solutions). See "Possible Issues".

GenericCylindricalDecomposition[And @@ {ab == bc, bc == ac, ac == ab, Abs[m] > 0}, Abs[m]]

{False, -3 Abs[m] + Abs[m]^4 == 0}


CylindricalDecomposition will handle it, but one needs algebraic equations (here in terms of Abs[m]):

CylindricalDecomposition[And @@ {ab == bc, bc == ac, ac == ab, Abs[m] > 0}, Abs[m]]

Abs[m] == 3^(1/3)


It does not work with the OP's setup (a non-algebraic system):

CylindricalDecomposition[And @@ {ab == bc, bc == ac, ac == ab, m > 0}, m]

False


Solve can also solve the algebraic case, albeit with the solution in a different form:

Solve[{ab == bc, bc == ac, ac == ab, Abs[m] > 0}, Abs[m], Reals]

{{Abs[m] -> Root[-3 + #1^6 &, 2]^2}}


Either solution can be combined with the following step to get the answer:

Solve[{Abs[m] == Root[-3 + #1^6 &, 2]^2, m > 0}, m, Reals]

{{m -> 3^(1/3)}}


### Workarounds

If viewpoint 1) above about generic solutions is correct, then a two-step process of finding all solutions to the equations and then selecting the ones fitting the condition seems appropriate.

First, the equations, as requested in a reduced-code form:

eqns = Equal @@@ Partition[MapThread[EuclideanDistance,
{#, RotateLeft[#]}] &[
{#, f[#]} & /@ (x /. Solve[D[f[x], x] == 0, x])],
2, 1, 1]

{Sqrt[Abs[m] + Abs[m]^4] == 2 Sqrt[Abs[m]],
2 Sqrt[Abs[m]] == Sqrt[Abs[m] + Abs[m]^4], True}

Cases[Solve[eqns, m], pat : Rule[m, sol_?Positive] :> pat, 2]

{m -> 3^(1/3)}


A better approach to me is to avoid EuclideanDistance and use an algebraic distance formula or, when feasible, its square. In this case, the observation about generic solutions does not apply, perhaps because CylindricalDecomposition is chosen for algebraic equations.

dist[p_, q_] := Sqrt[#.#] &[p - q];
{#, RotateLeft[#]}] &[{#, f[#]} & /@
(x /. Solve[D[f[x], x] == 0, x])], 2, 1, 1]

{Sqrt[m + m^4] == 2 Sqrt[m], 2 Sqrt[m] == Sqrt[m + m^4], True}

Solve[Append[eqns, m > 0], m]

{{m -> 3^(1/3)}}

distsq[p_, q_] := #.# &[p - q];
{#, RotateLeft[#]}] &[{#, f[#]} & /@
(x /. Solve[D[f[x], x] == 0, x])], 2, 1, 1]

{m + m^4 == 4 m, 4 m == m + m^4, True}

Solve[Append[eqns, m > 0], m]

{{m -> 3^(1/3)}}


### Simple, curious examples

Simplify[Abs[m]^2 == m^2 && m \[Element] Reals]

m^2 == Abs[m]^2 && m \[Element] Reals

FullSimplify[Abs[m]^2 == m^2 && m \[Element] Reals]

m \[Element] Reals

Solve[{Sqrt[Abs[m]] == 1}, m, Reals] (* generically 0-dimensional? *)

{{m -> -1}, {m -> 1}}

Solve[{Sqrt[Abs[m]] == 1, m > 0}, m] (* generically 1-dimensional? *)

{}

Solve[{Sqrt[Abs[m]] == 1, m > 0}, m, Reals] (* non-algebraic *)

{}

Solve[{Sqrt[Sqrt[m^2]] == 1, m > 0}, m] (* algebraic *)

{{m -> 1}}


Mathematica can handle Abs in simple situations:

Solve[{Abs[m] == 1, m > 0}, m, Reals]

{{m -> 1}}


Curiouser -- the first solution is sort of a contradiction:

Solve[{ab == bc, bc == ac, ac == ab, m > 0}, Abs[m], Reals]

{{Abs[m] -> ConditionalExpression[0, m > 0]},
{Abs[m] -> ConditionalExpression[3^(1/3), m > 0]}}

Simplify[Abs[m] == ConditionalExpression[0, m > 0]]

ConditionalExpression[False, m > 0]

-

Easier if you avoid radicals in the setup of the equations. This can be done by working with squares of distances.

f[x_] := x^4 - 2 m x^2 + m + 1
pts = {x, f[x]} /. Solve[f'[x] == 0, x];
pairs = Partition[pts, 2, 1, {1, -2}]

(* Out[75]= {{{0, 1 + m}, {-Sqrt[m], 1 + m - m^2}}, {{-Sqrt[m],
1 + m - m^2}, {Sqrt[m], 1 + m - m^2}}, {{Sqrt[m], 1 + m - m^2}, {0,
1 + m}}} *)

diffs = Apply[Subtract, pairs, {1}];
sqrdists = Simplify[Expand[Map[#.# &, diffs]], m > 0]

(* Out[77]= {m + m^4, 4 m, m + m^4} *)

ToRadicals[m /. Solve[Apply[Equal, sqrdists], m, Reals]]

(* Out[79]= {0, 3^(1/3)} *)

-

I don't think this is substantially different from what @Daniel Lichtblau proposed, but it directly uses your code up to the final Solve expression, where it also uses the trick of squaring in order to avoid square-roots:

Solve[{ab^2 == bc^2, bc^2 == ac^2, ac^2 == ab^2, m > 0}, m, Reals]
(* Out:  {{m -> 3^(1/3)}}  *)

-