Ovals of plane curves (and esoteric surfaces in space)

Consider the curve $y^2 - x^3 + x=0$ with $x, y$ real. This curve has two connected components (one bounded, one not), and the question is: how do you get Mathematica to tell you that (if you do:

Reduce[y^2 - x^3 + x > 0, {x, y}]


It gives something messy with more than two cases.

Of course, the real question is for an arbitrary algebraic curve, where you can't just look at the picture and figure it out...

An edit

The suggestion in the comment worked very well for this case, but here is something weird:

Reduce[x^2 + y^2 + z^2 == 1, {x, y, z}, Reals]

(* (x == -1 && y == 0 && z == 0) || (Inequality[-1, Less, x, Less, 1] &&
Inequality[-Sqrt[1 - x^2], LessEqual, y, LessEqual, Sqrt[1 - x^2]] &&
(z == -Sqrt[1 - x^2 - y^2] || z == Sqrt[1 - x^2 - y^2])) ||
(x == 1 && y == 0 && z == 0) *)


Since last time I looked the unit sphere was connected, there is something goofy going on.

• Reduce[y^2 - x^3 + x == 0, {x, y}, Reals] gives something tidier with only two cases - one with a bounded range of x, the other unbounded. Oct 22 '14 at 20:58
• @wxffles Cool! I was trying assumptions, but they did not work... Oct 22 '14 at 21:07
• LogicalExpand also helps one see what is going on. I don't have any ideas how to automatically interpret the conditions though. Oct 22 '14 at 22:50
• The 3D example in no way conflicts with connectivity, it's the outcome of Reduce using cylindrical decomposition. Oct 22 '14 at 23:26
• (1) Have a look at what CAD does. As far as I am aware, Reduce makes no attempt to post-process to join components. I don't actually think it should, either, since that could be more time consuming than the main computation (which, if it involves CAD, is already of considerable complexity). (2) Thus far I've seen zero indication of why this should be regarded as a misfeature. (3) It absolutely is not a bug (it's a correct result). Oct 23 '14 at 15:28

Not general enough, of course. But:

FullSimplify[ Or @@ Join @@ ((s = Solve[y^2 - x^3 + x == 0, {y}, Reals]) /.
(y -> ConditionalExpression[__, b__]) :> b)]

(* -1 <= x <= 0 || x >= 1 *)

Plot[y /. s, {x, -2, 2}]


The same for your sphere:

s = FullSimplify [ Or @@ Join @@ (Solve[x^2 + y^2 + z^2 == 1, {z},  Reals] /.
(z -> ConditionalExpression[__, b__]) :> b)]

(* x + Sqrt[1 - y^2] > 0 && Sqrt[1 - y^2] > x *)

Quiet@RegionPlot[s, {x, -2, 2}, {y, -2, 2}]


So I wrote a thing. It's a bit ugly, but it seems to work:

Clear[bounds];
bounds[eq_Equal, var_] := bounds[#, var] & /@
List @@ List @@@ LogicalExpand@Reduce[eq, var, Reals];
bounds[cond_List, var_] :=
Module[{collect, ranges = ConstantArray[0, Length@var], n = Length@var},
collect = Table[Select[cond, ! FreeQ[#, var[[i]]] &&
And @@ (Function[v, FreeQ[#, v]] /@ var[[(i + 1) ;;]]) &], {i, n}];
Table[ranges[[i]] = var[[i]] /. {
Last@Quiet@Maximize[{var[[i]], Flatten@collect[[;; i]]}, var[[;; i]], Reals],
Last@Quiet@Minimize[{var[[i]], Flatten@collect[[;; i]]}, var[[;; i]], Reals]}, {i, n - 1}];
ranges[[n]] = {First@Quiet@Maximize[{Last@Last@Last@collect, Flatten@Most@collect},
var, Reals],
First@Quiet@Minimize[{Last@Last@Last@collect, Flatten@Most@collect},
var, Reals]};
]


What it should give you is a bounding box for each section that Reduce finds for you. Let's try it out:

equation = y^2 == x (x + 1) (x - 3) (x + 2) (x - 2);
N@bounds[equation, {x, y}]


{{{∞, 0.}, {3., -∞}}, {{∞, ∞}, {3., 0.}}, {{-1., 0.}, {-2., -2.52703}}, {{2., 0.}, {0., -3.50095}}, {{-1., 2.52703}, {-2., 0.}}, {{2., 3.50095}, {0., 0.}}}

ContourPlot[Evaluate@equation, {x, -3, 4}, {y, -15, 15},
Prolog -> {LightGray, Rectangle @@@ Select[N@bounds[equation, {x, y}],
NumberQ[Quiet@Total@Flatten@#] &]}]


Edit: Some more examples:

trott = 12^2 (x^4 + y^4) - 15^2 (x^2 + y^2) + 350 x^2 y^2 + 81 == 0;
ContourPlot[Evaluate@trott, {x, -1.5, 1.5}, {y, -1.5, 1.5},
Prolog -> {LightGray, Rectangle @@@ Chop[N@bounds[trott, {x, y}], 10^-8]}]


(It just needed a little chopping ;)

eq3 = 9 (x - y^2)^2 + 13 (y - x^2)^2 + 7 z^2 == 1;
Show[ContourPlot3D[Evaluate@eq3, {x, -0.5, 1.5}, {y, -0.5, 1.5}, {z, -1, 1}],
Graphics3D[{Opacity[0.25], Cuboid @@@ N@bounds[eq3, {x, y, z}]}]]


Edit 2: Nested ovals.

If we take a polynomial with some known roots, say 1, 4 and 9. Then use the square distance from the origin as the variable, we should get some nested ovals:

nested = 0 == Expand[(t - 1) (t - 4) (t - 9) /. t -> x^2 + y^2]


$$0 = -36 + 49 x^2 - 14 x^4 + x^6 + 49 y^2 - 28 x^2 y^2 + 3 x^4 y^2 - 14 y^4 + 3 x^2 y^4 + y^6$$

ContourPlot[Evaluate@nested, {x, -3, 3}, {y, -3, 3},
Prolog -> {EdgeForm[Red], FaceForm[LightGray], Rectangle @@@ bounds[nested, {x, y}]}]


So it does handle nested ovals, although it's cut up into even more components.

• Cool! But it doesn't like unbounded components? I think another potential problem (I don't have an example in mind right now, but I will think of one tomorrow, hopefully) is of nested ovals... Oct 23 '14 at 3:10
• You might want to check this with the Trott curve: Oct 23 '14 at 3:16
• Which is: $144(x^4+y^4)-225(x^2+y^2)+350x^2y^2+81=0$ Oct 23 '14 at 3:17
• It handles unbounded components. That was the original point - trying to distinguish those. You just can't make rectangles out of them. As for nested ovals, no idea! Oct 23 '14 at 3:33

Maybe these two versions of Plot3D can give some insight:

Plot3D[{0, y^2 - x^3 + x}, {x, -2, 2}, {y, -2, 2}]

Plot3D[Boole[y^2 - x^3 + x > 0], {x, -2, 2}, {y, -2, 2}]

(* Pictures not given *)


Regards, Wolfgang

• Cool pictures! But I already knew what the curve looked like, I am just trying to get Mathematica to give the right number of components (or, as they are known in the trade, "ovals"). Alas, not trivial. Oct 22 '14 at 22:34
• @ Igor Rivin: Could you provide a more complex example, please? Oct 22 '14 at 23:27
• Here is the "Trott curve": $144(x^4+y^4)-225(x^2+y^2)+350x^2y^2+81=0$ Oct 23 '14 at 3:12
• The Trott curve has four ovals (according to Wikipedia), but the Reduce disjunction has nine cases. Oct 23 '14 at 3:14
• @ Igor Rivin: Why not simplify the problem? In polar coordinates your new example is gr = 81 - 225*r^2 + (607*r^4)/4 - (31/4)*r^4*Cos[4*f], and Reduce[gr>0&&r>=0, Reals] gives r == 0 || (r > 0 && Cos[4*f] < (324 - 900*r^2 + 607*r^4)/(31*r^4)). And gr==0 gives us 4 nontrivial solution r(f), as it should. Oct 23 '14 at 8:11