How to plot Clebsch surface with those 27 lines?

I can get its Algebraic Equation by this code:

Entity["Surface", "ClebschDiagonalCubic"][EntityProperty["Surface", "AlgebraicEquation"]]


Then get its graphics:

ContourPlot3D[81 (x^3 + y^3 + z^3) - 9 (x^2 + y^2 + z^2) -
189 (x^2 y + x^2 z + x y^2 + x z^2 + y^2 z + y z^2) + 54 x y z -
9 (x + y + z) + 126 (x y + x z + y z) - 1 == 0, {x, -1, 1}, {y, -1,
1}, {z, -1, 1}, Boxed -> False, BoxRatios -> {1, 1, 1},
AxesOrigin -> {0, 0, 0}, AxesStyle -> {Red, Green, Blue}, Mesh -> None]


As we know, there are 27 lines in its surface:

But I don't know how to draw these lines onto the surface

• AMS blog on the topic for reference.
– Syed
Feb 1, 2023 at 10:59
• Is this a question about plotting the lines of finding these 27 lines? Feb 1, 2023 at 12:23
• @kirma Yeah, wasn't that clear in my post? Could you edit the question for me if I'm not making myself clear? My English is not very good.
– yode
Feb 1, 2023 at 12:28
• @yode The added link on 27 lines makes this clearer. Feb 1, 2023 at 13:00
• @yode could you perhaps consider adding this link to the OP? I checked that the author has included all the equations for the lines and there's a mathematical explanation on how to derive said equations. I think this should make the task at hand much easier
– bmf
Feb 1, 2023 at 13:16

EDIT:

With help from @cvgmt, both fixing a wrong term in the equation(!) and adding constraints for finding lines:

With[{expr =
1 + 54 x y z - 9 (x + y + z) + 126 (x y + x z + y z) -
9 (x^2 + y^2 + z^2) -
189 (x^2 y + x y^2 + x^2 z + y^2 z + x z^2 + y z^2) +
81 (x^3 + y^3 + z^3) == 0},
(* Expression over a parametric infinite line. *)
expr /. Thread[{x, y, z} -> {a, b, c} + t {u, v, w}] //
(* Find lines which lie on the surface. *)
Solve[
(* Expression must be true for all t,
that is, it must lie on the surface everywhere. *)
Resolve[ForAll[t, #]] &&
(* Constrain solutions:
line offset and line direction vectors must be at right angle,
and line direction vector must be unit length. *)
{a, b, c} . {u, v, w} == 0 && u^2 + v^2 + w^2 == 1 &&
(* Prevent finding lines with "direction" negated. *)
((w > 0) || (w == 0 && v > 0) || (w == 0 && v == 0 && u > 0)),
{a, b, c, u, v, w}, Reals] & //
(* Line segments inside the plotted area. *)
RegionIntersection[
InfiniteLine[{a, b, c}, {u, v, w}], Ball[]] /. N[#] & //
(* Create animated graphics. *)
Animate[
Show[
(* Plot the surface. *)
ContourPlot3D[expr, Element[{x, y, z}, Ball[]],
PlotPoints -> 100, ContourStyle -> Opacity[9/10], Boxed -> False,
BoxRatios -> Automatic, Axes -> None, Mesh -> None,
RegionBoundaryStyle -> None,
SphericalRegion -> True, ViewAngle -> 10 Degree,
ViewPoint -> Dynamic[{5 Sin[a], 5 Cos[a], 5/2}]],
(* Plot line segments. *)
Graphics3D[{Thick, CapForm["Butt"], #}]],
{a, 0, 2 Pi}] &]


This is more of an extended comment - but anyway, it would seem that Mathematica has at least some success finding those lines on the surface:

With[{eqn = -1 + 54 x y z - 9 (x + y + z) + 126 (x y + x z + y z) -
9 (x^2 + y^2 + z^2) -
189 (x^2 y + x y^2 + x^2 z + y^2 z + x z^2 + y z^2) +
81 (x^3 + y^3 + z^3)},
CoefficientList[eqn /. Thread[{x, y, z} -> {a, b, c} + t {u, v, w}],
t] & // Map[# == 0 &] //
Solve[
And @@ # && {a, b, c} . {u, v, w} == 0 &&
u^2 + v^2 + w^2 == 1, {a, b, c, u, v, w}, Reals] & // N //
RegionIntersection[InfiniteLine[{a, b, c}, {u, v, w}],
Cuboid[{-1, -1, -1}, {1, 1, 1}]] /. # & //
ListAnimate@
Table[
Show[ContourPlot3D[eqn == 0, {x, -1, 1}, {y, -1, 1}, {z, -1, 1},
PlotPoints -> 50, ContourStyle -> Opacity[9/10],
Boxed -> False, BoxRatios -> Automatic, Axes -> None,
Mesh -> None, SphericalRegion -> True,
ViewPoint -> {5 Sin[a], 5 Cos[a], 5/2}],
Graphics3D[{Thick, #}]],
{a, 0, 2 Pi - Pi/64, Pi/64}] &]


• (+1) If we add the restriction to the direction {u,v,w} to be u^2 + v^2 + w^2 == 1 && (( w > 0) || ( w == 0 && v > 0 ) || (w == 0 && v == 0 && u > 0) ) and replace the orginal as below, we can get exact 27 solutions.( replace -1 to 1. Feb 1, 2023 at 16:20
• sols = With[{eqn = 1 + 54 x y z - 9 (x + y + z) + 126 (x y + x z + y z) - 9 (x^2 + y^2 + z^2) - 189 (x^2 y + x y^2 + x^2 z + y^2 z + x z^2 + y z^2) + 81 (x^3 + y^3 + z^3)}, CoefficientList[eqn /. Thread[{x, y, z} -> {a, b, c} + t {u, v, w}], t] & // Map[# == 0 &] // Solve[And @@ # && {a, b, c} . {u, v, w} == 0 && u^2 + v^2 + w^2 == 1 && ((w > 0) || (w == 0 && v > 0) || (w == 0 && v == 0 && u > 0)), {a, b, c, u, v, w}, Reals] &]; sols // Length Feb 1, 2023 at 16:22
• @cvgmt The unit length approach was the first that came to my mind to parameterise solution search in a sensible way. (Thus the "extended comment.") Thanks for the tip! Feb 1, 2023 at 16:24
• @cvgmt BTW, what you intended to tell about -1 and 1? :) Feb 1, 2023 at 16:28
• The right equation should be eqn = 1 + 54 x y z+... instead of eqn=-1+54 x y z+... It is a typo of the original equations. Feb 1, 2023 at 16:30

Provide a more concise solution:

expr = 1 + 54 x y z - 9 (x + y + z) + 126 (x y + y z + z x) -
9 (x^2 + y^2 + z^2) -
189 (x^2 y + x y^2 + y^2 z + y z^2 + z^2 x + z x^2) +
81 (x^3 + y^3 + z^3);
sols = Solve[{Resolve[
ForAll[t, (expr /.
Thread[{x, y, z} -> {a, b, c} + t {u, v, w}]) == 0]], {a, b,
c} . {u, v, w} == 0, u^2 + v^2 + w^2 == 1,
w > 0 || w == 0  && v > 0 || w == 0 && v == 0 && u > 0}, {a, b, c,
u, v, w}, Reals];
Show[ContourPlot3D[expr == 0, {x, -2, 2}, {y, -2, 2}, {z, -2, 2},
Mesh -> None, PlotPoints -> 60, MaxRecursion -> 4,
ContourStyle -> Cyan],
Graphics3D[{White, InfiniteLine[{a, b, c}, {u, v, w}] /. sols}],
Boxed -> False, Axes -> False, ViewPoint -> Front]


Explanation

1. We can get the lines in the surface by the parametric equation of a line to combine the original equation to get eqns.
2. Since the original equation is $$0$$, We can let the equation of the substituted line be constant $$0$$
3. I assumed that the line was translated through the origin to the surface:

So it has to satisfy:{a,b,c}.{u,v,w}==0

1. Since {u,v,w} is a vector, we can assume that it's distributed along the circumference of the unit circle. So it meet:u^2+v^2+w^2==1

2. Now, one last problem, we solve the equation this way and we get $$54$$ solutions. Because MMA can't tell the difference between $$-vec$$ and $$vec$$ being the same vector. We have to limit it with some conditions:

 ((w > 0)(*At an acute Angle to the w axis*)
|| (w == 0 && v > 0)(*Perpendicular to the w axis but at an acute Angle to the v axis*)
|| (w == 0 && v == 0 && u > 0)(*Perpendicular to both the w and v axes,but at an acute Angle to the u axis*))


But I have an additional puzzle: Why can I only solve three lines on a cubic surface $$x^3 + 3 y^3 + z^3-2 x^2 + 5 x y - x + 7=0$$? Why can't I get $$27$$ too as this statement?

• Yes, this is definitely a more thorough explanation. :) I'd say that the most complicated part of finding this kind of a way to find the lines which lie on the surface is to use a parameterised line whose defining vectors are constrained in a way amenable to processing and getting result from Mathematica. There are many ways to define a line... Feb 3, 2023 at 8:27
• The CoefficientList methods (polynomial is zero everywhere iff all coefficients are zero) can be replaced by Resolve and ForAll, which might be more descriptive from the mathematical viewpoint, although they're the same... Feb 3, 2023 at 8:30