Applications of Symbolic Inequality Solving in Geometry

Adam Strzebonski, Wolfram Research.

Initialization

<< Experimental`;
<< Developer`;
<< Utilities`FilterOptions`;
<< Algebra`InequalitySolve`;
<< Graphics`ImplicitPlot`;
<< Graphics`FilledPlot`;
Off[General::spell1];
Off[General::spell];
Off[$MaxExtraPrecision::meprec];

\(\(\(IneqPlot[ineq_, {x_, x0_, x1_}, {y_, y0_, y1_}, 
        opts___?OptionQ] := 
	
      Module[{cyls, \ plots, \ 
          bds}, 
		cyls = \(GenericCylindricalAlgebraicDecomposition[
              ineq && x0 < x < x1 && y0 < y < y1, {x, 
                y}]\)[\([1]\)]; 
		bds = \(Boundaries /. {opts}\) /. {Boundaries -> False}; 
		plots = 
			If[cyls === False, 
            Plot[y1 + 1, {x, x0, x1}, PlotRange -> {{x0, x1}, {y0, y1}}, 
              Evaluate[FilterOptions[\ FilledPlot, \ opts\ ]], 
              DisplayFunction -> Identity], 
					If[
              Head[cyls] === Or, cyls = List @@ cyls, 
              cyls = {cyls}]; 
					\(cylinderplot[#, x, 
                  FilterOptions[FilledPlot, \ opts]] &\) /@ cyls]; 
		If[
          bds\  && \ cyls =!= False, \ 
          plots = {plots, 	\(boundaryplot[#, x, 
                    FilterOptions[FilledPlot, \ opts]] &\) /@ 
                cyls}]; 
	\ \ Show[\ 
          plots, \ 
\ \ \ \ \ \ FilterOptions[\ Graphics, \ 
            opts], \ 
\ \ \ \ \ \ Axes -> True, \ 
          AxesLabel -> {x, 
              y}, \ 
\ \ \ \ \ \ DisplayFunction -> $DisplayFunction, 
          AspectRatio -> Automatic, 
          PlotRange -> {{x0, x1}, {y0, y1}}\ ]]\)\(
\)
    \)\)
\(\(\(cylinderplot[cyl_, x_, opts___?OptionQ] := 
	 (*\ 
        assumes\ form\ a < x < 
            b\  && \ \((\(\(c < y < d\)\(\ \)\(||\)\)\  ... )\)\ *) 
	
      Module[{x0 = cyl[\([1, 1]\)], x1 = cyl[\([1, 5]\)], ys}, 
		ys = 
          If[Head[cyl[\([2]\)]] === Or, 
            List @@ \((cyl[\([2]\)])\), {cyl[\([2]\)]}]; 
		\(FilledPlot[\ {#[\([1]\)], #[\([5]\)]}, \ {x, \ x0, x1}, \ DisplayFunction -> Identity, 
              opts] &\) /@ ys]\)\(		\)\(
\)
    \)\)
\(\(\(boundaryplot[cyl_, x_, opts___?OptionQ] := 
	 (*\ 
        assumes\ form\ a < x < 
            b\  && \ \((\(\(c < y < d\)\(\ \)\(||\)\)\  ... )\)\ *) 
	
      Module[{x0 = cyl[\([1, 1]\)], x1 = cyl[\([1, 5]\)], ys}, 
		ys = 
          If[Head[cyl[\([2]\)]] === Or, 
            List @@ \((cyl[\([2]\)])\), {cyl[\([2]\)]}]; 
		\({Graphics[
                Line[{{x0, #[\([1]\)] /. x -> x0}, {x0, #[\([5]\)] /. 
                        x -> x0}}]], 
			Graphics[
                Line[{{x1, #[\([1]\)] /. x -> x1}, {x1, #[\([5]\)] /. 
                        x -> x1}}]], 
              Plot[\ {#[\([1]\)], #[\([5]\)]}, \ {x, \ x0, x1}, \ 
                DisplayFunction -> Identity, opts]} &\) /@ 
          ys]\)\(		\)\(
\)
    \)\)
xypoints[{xl_, xu_, yl_, yu_}, pp_] := 
      Block[{x, y}, 
        Module[{ε = 10. ^-6, dx, dy, nx}, δ = 
            1. /pp; nx = 
            Max[Floor[(xu - xl - 2  ε)/δ], 
              1]; δx = (xu - xl - 2  ε)/nx; 
          dx = xu - xl; 
          dy = 0; 
Do[
            x = xl + ε + ix δx; 
            dy = Max[dy, yu - yl], {ix, 0, nx}]; 
{Ceiling[
              dx/δ], Ceiling[dy/δ]}]]

    
makePolygons3D[
        Inequality[xl_, Less, x, Less, xu_] && 
          Inequality[yl_, Less, y, Less, yu_] && 
          Inequality[zl_, Less, z, Less, zu_], pp_] := 
      Block[{x, y}, 
        Module[{ε = 10. ^-12, nx, ny}, {nx, ny} = 
            xypoints[{xl, xu, yl, yu}, pp]; 
          If[nx < 1 || ny < 1, 
            Print[{nx, ny}]]; δx = (xu - xl - 2  ε)/
              nx; 
pointsl = 
            Table[x = 
                xl + ε + 
                  ix δx;  δy = (yu - yl - 
                      2  ε)/ny; 
Table[
                y = yl + ε + iy δy;  {x, y, 
                  Re[zl]}, {iy, 0, ny}], {ix, 0, 
                nx}]; 
pointsu = 
            Table[x = 
                xl + ε + 
                  ix δx;  δy = (yu - yl - 
                      2  ε)/ny; 
Table[
                y = yl + ε + iy δy;  {x, y, 
                  Re[zu]}, {iy, 0, ny}], {ix, 0, 
                nx}]; 
{Table[
              Polygon[{pointsl[[i, j]], pointsl[[i + 1, j]], 
                  pointsl[[i + 1, j + 1]], pointsl[[i, j + 1]]}], {i, 
                nx}, {j, ny}], 
            Table[Polygon[{pointsu[[i, j]], pointsu[[i + 1, j]], 
                  pointsu[[i + 1, j + 1]], pointsu[[i, j + 1]]}], {i, 
                nx}, {j, ny}]}]]

    
CADPlot3D[gcad_, opts___] := 
      With[{gc = 
            gcad[[1]] //. 
              a_ && b_Or :> (a && # & /@ 
                    b)}, 
With[{polyg = makePolygons3D[#, 
                    Evaluate[(PlotPoints /. {opts}) /. 
                        PlotPoints -> 15]] & /@ 
                Switch[Head[gc], Or, List @@ gc, And, {gc}, _, {}]}, 
          Print[Length[Flatten[polyg]]]; 
          Graphics3D[{EdgeForm[], polyg}, FilterOptions[Graphics3D, opts], 
            PlotRange -> All, 
            AspectRatio -> Automatic]]]

\(CADPlot2D[cad_, opts___] :=  (*assumes\ nonempty\ and\ bounded*) 
	
      Module[{cyls, \ plots, \ x, y}, 
cyls = 
          cad //. a_ && b_Or :> \((\(a && # &\) /@ b)\); 
        If[Head[cyls] === Or, cyls = List @@ cyls, cyls = {cyls}]; 
        cyls = Select[cyls, FreeQ[#, Equal] &]; x = cyls[\([1, 1, 3]\)]; 
        y = cyls[\([1, 2, 3]\)]; 
        plots = \(cylinderplot[#, x, FilterOptions[FilledPlot, \ opts]] &\) /@ 
            cyls; 
	\ \ Show[\ 
          plots, \ 
\ \ \ \ \ \ FilterOptions[\ Graphics, \ 
            opts], \ 
\ \ \ \ \ \ Axes -> True, \ 
          AxesLabel -> {x, 
              y}, 
\ \ \ \ \ \ DisplayFunction -> $DisplayFunction, 
          AspectRatio -> Automatic, PlotRange -> All\ ]]\)

clgo = 
        GenericCylindricalAlgebraicDecomposition[
          32 - 216*x^2 + 648*x^2*y - 216*y^2 - 216*y^3 - 150*z + 216*x^2*z + 
                216*y^2*z + 231*z^2 - 113*z^3 < 0 && -3/2 < x < 
              3/2 && -3/2 < y < 3/2 && -5/2 < z < 5/2, {x, y, 
            z}];
Save["ClGo.m", clgo];

<< ClGo.m;

cgplot = CADPlot3D[clgo, PlotPoints -> 30]

Chessboard[gr_, a_, b_] := 
      gr /. Polygon[l_] :> {SurfaceColor[Hue[a], Hue[b],  2.7], 
            Polygon[Apply[Plus, Partition[Append[l, First[l]], 2, 1], {1}]/
                2]}

cgsurf = Chessboard[cgplot, 0.21, 0.5];

Save["ClGoPlot.m", cgplot];
Save["ClGoSurf.m", cgsurf];

<< ClGoSurf.m;

What does it mean to "solve" a system of equations and inequalities.

Solutions of systems of real polynomial equations and inequalities form semialgebraic sets.

ineqs = (((x + (7/4)))^2) + (((y - (9/4)))^2) < 
          (1/9) || (((x - (9/4)))^2) + (((y - (9/4)))^2) < (1/9) || 
        (x^2) + (y^2) < 25 && (((x - 2))^2) + (((y - 2))^2) > 
            1 && (((x + 2))^2) + (((y - 2))^2) > 
            1 && ((x^2) + (((y + 1))^2) > 5 || 
              (x^2) + (((y - 1))^2) < 9); 
    IneqPlot[ineqs, {x, -5, 5}, {y, -5, 5}]

Infinite, non-parametrizable.

Every semialgebraic set can be represented as a finite union of disjoint cells.

A cell in ℝ is an open interval, or a point.

A cell in (ℝ^(n + 1)) can have one of the following two forms:
\(TraditionalForm\`{\((x, 
              y)\) ∈ ℝ\^n×\(ℝ : \ x ∈ C\  ∧ \ y\) = r(x)}\)
\(TraditionalForm\`{\((x, 
            y)\) ∈ ℝ\^n×\(ℝ : \ 
              x ∈ C\  ∧ \ \(r\_1\)(x) < y < \(r\_2\)(x)\)}\),
where r, (Subscript[r, 1]), (Subscript[r, 2]) are continuous functions which are roots of polynomials from ℝ[x][t], (Subscript[r, 1])may also be -∞, and (Subscript[r, 2])may be ∞.

Cells can be directly represented in a computer algebra system.

Mathematica has built-in functions representing roots of polynomials. Root[f, k] represents the k-th root of f, where the roots of f are ordered so that the real roots come first, ordered by their values. So for instance

Root[#(1^3) + a #1 + b &, 1]

represents the smallest real root of (t^3) + a t  + b.

Root functions are an integral part of the Mathematica system, in particular they can be evaluated numerically, and plotted.

Root[#(1^3) + a #1 + b &, 1] /. {a -> 1.23, 
        b -> 3.45}

Plot3D[
      Root[#(1^7) + a #(1^3) + (b^2) #1 + (a^2) b &, 1], {a, -10, 
        10}, {b, -10, 10}, PlotPoints -> 50]

Show[
      Table[Plot3D[
          Re[Root[#(1^3) + #(1^2) + (a + b I) #1 - (1/2) &, 
              i]], {a, -3, 3}, {b, -3, 3}, 
          DisplayFunction -> Identity, PlotPoints -> 25], {i, 3}], 
      DisplayFunction -> $DisplayFunction, BoxRatios -> {1, 1, 2}, 
      ViewPoint -> {1.53, -1.9, 1.15}, Boxed -> False, 
      Axes -> False]

The cylindrical algebraic decomposition algoritm allows to find decompositions of solution sets of systems of polynomial equations and inequalities into finite unions of disjoint cells

CAD = Experimental`CylindricalAlgebraicDecomposition;
CAD[1 ≤ (x^2) + (y^2) < 2, {x, y}]

SeedRandom[21]; 
    IneqPlot[1 ≤ (x^2) + (y^2) < 2, {x, -2, 2}, {y, -2, 2}, 
      Fills :> Hue[Random[]]]

GCAD = 
        Experimental`GenericCylindricalAlgebraicDecomposition;
GCAD[1 ≤ (x^2) + (y^2) < 2, {x, y}]

ParametricPlot[{Cos[5 x], Sin[4 x]}, {x, 0, 
          2 π}];

s = 0; 
    IneqPlot[16 (x^2) - 80 (x^4) + 128 (x^6) - 64 (x^8) - 25 (y^2) + 
          200 (y^4) - 560 (y^6) + 640 (y^8) - 256 (y^10) < 
        0, {x, -(3/2), (3/2)}, {y, -(3/2), (3/2)}, 
      Fills :> Hue[s = 1.27 s + 0.04]]

GCAD[
      16 (x^2) - 80 (x^4) + 128 (x^6) - 64 (x^8) - 25 (y^2) + 200 (y^4) - 
            560 (y^6) + 640 (y^8) - 256 (y^10) < 0 && -(3/2) < x < 
          (3/2) && -(3/2) < y < (3/2), {x, y}]

i = 0; 
    IneqPlot[16 (x^2) - 80 (x^4) + 128 (x^6) - 64 (x^8) - 25 (y^2) + 
          200 (y^4) - 560 (y^6) + 640 (y^8) - 
          256 (y^10) < -(1/20), {x, -(3/2), 
        (3/2)}, {y, -(3/2), (3/2)}, 
      Fills :> Hue[((i++)/25.) ]]; i

Here are some pictures generated using the output of (generic) cylindrical algebraic decompostion.

IneqPlot[((((x^2) + (y^2)))^2) - 2 ((x^2) - (y^2)) < 0, {x, -2, 
        2}, {y, -2, 2}]

lem = ((((x^2) + (y^2) + (z^2) - 1))^2) + 4 (z^2) < 1;

lemcad = GCAD[lem, {x, y, z}]

Show[CADPlot3D[lemcad, PlotPoints -> 30], 
      Boxed -> False]

torus[a_, 
          r_] := ((((x^2) + (y^2) + (z^2)))^2) + ((a^2) - (r^2)) (z^2) - 
            2 ((a^2) + (r^2)) ((x^2) + (y^2)) + ((((a^2) - (r^2)))^2) < 
          0;

torcad1 = GCAD[torus[1, (1/3)], {x, y, z}];

Show[CADPlot3D[torcad1, PlotPoints -> 30], 
      Boxed -> False]

lemtor = 
      GCAD[lem && ! (torus[1, (1/3)]), {x, y, 
          z}]

plt = CADPlot3D[lemtor, PlotPoints -> 20];

Show[Chessboard[plt, 0.1, 0.9], Boxed -> False]

lemtor2 = 
      GCAD[lem && ! (torus[(3/2), (1/3)]), {x, 
          y, z}]

Show[CADPlot3D[lemtor2, PlotPoints -> 30], 
      Boxed -> False]

Clebsch surface:
32 - 216 (x^2) + 648 (x^2) y - 216 (y^2) - 216 (y^3) - 150 z + 
          216 (x^2) z + 216 (y^2) z + 231 (z^2) - 113 (z^3) == 0

Show[cgsurf, ViewPoint -> {-2, -1.6, 1.4}, 
      PlotRange -> {0.7 {-2. , 2. }, 
          0.7 {-2. , 2. }, {-2.2, 2}}]

Geometric theorem proving

A geometric inequality from: D. Mitrinovic, J.E. Pecaric, V. Volenec, "Recent Advances in Geometric Inequalities", Kluwer Academic Publishers, 1989.

triangle = 
        a > 0 && b > 0 && c > 0 && a + c > b && b + c > a && 
          a + b > c;

ineq = 
        a b c (((a^2)/(b^2)) + ((b^2)/(c^2)) + ((c^2)/(a^2))) ≥ 
          (a^3) + (b^3) + (c^3) + a b (b - a) + a c (a - c) + 
            b c (c - b);

Timing[ForAllRealQ[Implies[triangle, ineq], {a, b, c}]]

Conjecture by T.Ono (1914). Counterexample G.Quijano (1915). Proof for acute triangles F.Balitrand (1916).

\(\(\(	\)\(OnoIneq = 
      27\ \((b\^2 + c\^2 - a\^2)\)\^2\ \((a\^2 + c\^2 - b\^2)\)\^2\ \((a\^2 + b\^2 - c\^2)\)\^2 ≤ \((4\ F)\)\^6; 
    F = \@\(s\ \((s - a)\)\ \((s - b)\)\ \((s - c)\)\); \ 
    s = \(1\/2\) \((a + b + c)\); \ 
    triangle = 
      a > 0 && b > 0 && c > 0 && a + b > c && b + c > a && a + c > b; 
    acute = a\^2 + b\^2 > c\^2 && b\^2 + c\^2 > a\^2 && 
        a\^2 + c\^2 > b\^2;\)\)\)

Timing[
      InequalityInstance[! (Implies[triangle, 
              OnoIneq]), {a, b, c}]]

PlotTriangle[3/4, 1/2, 1]

Timing[
      ForAllRealQ[Implies[triangle && acute, OnoIneq], {a, b, c}]]

Global optimization

Find the maximal k for which the following inequality is satisfied for all a, b, c sides of a triangle.

ineq = (a + b) (c^2) + (b^2) (a + c) + 
            (a^2) (b + c) ≥ 
          k (-a + b + c) (a - b + c) (a + b - 
                c);

Infimum[k, 
      k > 0  &&  triangle  &&  ! ineq, {a, b, c, k}]

k>0 && triangle && !ineq is an open set, and therefore so it is its projection on the k-axis, so we know that the infimum cannot be attained. But we can check it anyway...

ForAllRealQ[Implies[triangle,  ineq /. k -> 6]]

Show that 0 < 
        cos(α) + (Sqrt[2]) (cos(β) + 
                cos(γ)) ≤ 2. Is (Sqrt[2]) the optimal constant for both inequalities? How will it change if we change the right hand side of the second inequality?

ca = ((-(a^2) + (b^2) + (c^2))/(2 b c));
cb = (((a^2) - (b^2) + (c^2))/(2 a c));
cc = (((a^2) + (b^2) - (c^2))/(2 a b));

-Infimum[-k, 
        k > 0 && triangle  &&  ca + k (cb + cc) <= 0, {a, b, c, 
          k}]

ForAllRealQ[
      Implies[triangle,  ca + 1/2 (cb + cc) > 0], {a, b, c}]

Infimum[k, 
      k > 0 && triangle  &&  ca + k (cb + cc) > 2, {a, b, c, 
        k}]

Infimum[k, 
      k > 0 && triangle  &&  ca + k (cb + cc) > 3, {a, b, c, 
        k}]

Infimum[k, 
      k > 0 && triangle  &&  ca + k (cb + cc) > 4, {a, b, c, 
        k}]

#[[1]] && 
          With[{t = #[[2]]}, 
            Switch[Head[t], Or, t[[1, 1]], And, t[[1]], _, 
              t]] & /@ GCAD[
          k > 0 && triangle && ca + k (cb + cc) > p, {p, k, a, b, 
            c}][[1]]

ForAllRealQ[
      Implies[triangle, 
        ca + Min[((1 + p)/2), (Sqrt[2 p - 2])] (cb + cc) ≤ 
          p], {a, b, c, p}]

Computing projections a.k.a. quantifier elimination

clproj = 
      Resolve[(Subscript[∃ , x, (x | y | z) ∈ Reals])(-(3/2) < y < (3/2) && -(5/2) < z < (5/2) && 
              32 - 216 (x^2) + 648 (x^2) y - 216 (y^2) - 216 (y^3) - 
                  150 z + 216 (x^2) z + 216 (y^2) z + 231 (z^2) - 
                  113 (z^3) == 0), {y, z}]

CADPlot2D[clproj]

Show[cgsurf, ViewPoint -> {2, 0, 0}, 
      PlotRange -> {0.7 {-2. , 2. }, 
          0.7 {-2. , 2. }, {-2.2, 2}}]

ParametricPlot3D[{Cos[t] (2 + Sin[s]^5), Sin[t] (2 + Sin[s]^5), 
        Cos[s]^5}, {s, 0, 2  Pi}, {t, 0, 2  Pi}, 
      PlotPoints -> 30]

TimeConstrained[
      GroebnerBasis[{x - Cos[t] (2 + Sin[s]^5), 
            y - Sin[t] (2 + Sin[s]^5), z - Cos[s]^5, 
            Sin[s]^2 + Cos[s]^2 - 1, Sin[t]^2 + Cos[t]^2 - 1}, {x, y, 
            z}, {Sin[s], Cos[s], Sin[t], Cos[t]}] // Timing, 3600]

(gb = 
            GroebnerBasis[{x - Cos[t] (2 + Sin[s]^5), 
                y - Sin[t] (2 + Sin[s]^5), z - Cos[s]^5, 
                Sin[s]^2 + Cos[s]^2 - 1, Sin[t]^2 + Cos[t]^2 - 1}, {Sin[s], 
                Cos[s], Sin[t], Cos[t], x, y, z}, 
              MonomialOrder -> DegreeReverseLexicographic];) // 
      Timing

ewm = {{1, 1, 1, 1, 0, 0, 0}, {1, 0, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 
            0, 0}, {0, 0, 1, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 1, 1}, {0, 0, 0, 0, 
            1, 0, 0}, {0, 0, 0, 0, 0, 1, 0}};

(gb1 = 
            Internal`GroebnerWalk[
              gb, {Sin[s], Cos[s], Sin[t], Cos[t], x, y, z}, 
              Internal`FirstMonomialOrder -> DegreeReverseLexicographic, 
              Internal`SecondMonomialOrder -> ewm];) // 
      Timing

Select[gb1, FreeQ[#, s] && FreeQ[#, t] &]