Labs Septic

Author

Eric W. Weisstein
June 6, 2008

This notebook downloaded from http://mathworld.wolfram.com/notebooks/Surfaces/LabsSeptic.nb.

For more information, see Eric's MathWorld entry http://mathworld.wolfram.com/LabsSeptic.html.

©2008 Wolfram Research, Inc. except for portions noted otherwise

Equation

alpha=Root[7*#^3 + 7*# + 1&,1];
a1=-12/7*alpha^2-384/49*alpha-8/7;
a2=-32/7*alpha^2+24/49*alpha-4;
a3=-4*alpha^2+24/49*alpha-4;
a4=-8/7*alpha^2+8/49*alpha-8/7;
a5=49*alpha^2-7*alpha+50;
P=x*(x^6-3*7*x^4*y^2+5*7*x^2*y^4-7*y^6)+7*z*((x^2+y^2)^3-2^3*z^2*(x^2+y^2)^2+2^4*z^4*(x^2+y^2))-2^6*z^7;
U=(z+a5*w)*(a1*z^3+a2*z^2*w+a3*z*w^2+a4*w^3+(z+w)*(x^2+y^2))^2;
LabsW=P-U;

Labs=FullSimplify[LabsW/.w->1]

Labs=(x^7)-21 (x^5) (y^2)+35 (x^3) (y^4)-7 x (y^6)-64 (z^7)+7 ((x^2)+(y^2)) z ((((x^2)+(y^2)-4 (z^2)))^2)-(1/2401)(z+Root[-1+3 #1-52 ((#1)^2)+((#1)^3)&,1]) (((7 (1+z) (-7 (x^2)-7 (y^2)+4 (2+z) (1+2 z))+8 (-1+3 z (-1+z (-1+16 z))) Root[1+7 #1+7 ((#1)^3)&,1]+28 (((1+z))^2) (2+3 z) ((Root[1+7 #1+7 ((#1)^3)&,1])^2)))^2);

AlgebraicDegree[eqn ,{x,y,z}]

Graphic

With[{r=4},
ContourPlot3D[
Evaluate[N[eqn]==0],
{x,-r,r},{y,-r,r},{z,-r,r},Boxed->False,BoxRatios->{1,1,1},Axes->False,Mesh->None,RegionFunction->(#1^2+#2^2+#3^2<r^2&)]
]

g=With[{r=4},
ContourPlot3D[
Evaluate[N[eqn]==0],
{x,-r,r},{y,-r,r},{z,-r,r},Boxed->False,BoxRatios->{1,1,1},Axes->False,Mesh->None,RegionFunction->(#1^2+#2^2+#3^2<r^2&),PlotPoints->25]
]

With[{r=4},
ContourPlot3D[
Evaluate[N[eqn]==0],
{x,-r,r},{y,-r,r},{z,-r,r},Boxed->False,BoxRatios->{1,1,1},Axes->False,Mesh->None,PlotPoints->30,MaxRecursion->4,RegionFunction->(#1^2+#2^2+#3^2<r^2&)]
]

g=With[{r=20},
ContourPlot3D[
Evaluate[N[eqn]==0],
{x,-r,r},{y,-r,r},{z,-r,r},Boxed->False,BoxRatios->{1,1,1},Axes->False,Mesh->None,PlotPoints->30,MaxRecursion->1,RegionFunction->(#1^2+#2^2+#3^2<r^2&)]
]

Show[{g,Graphics3D[{Red,PointSize[.04],Point/@nodes}]},BoxRatios->1]

Properties

Gaussian Curvature

ByteCount/@(k=GaussianCurvature[eqn,{x,y,z}])//Timing

Singular Points

Equations

(gb=GroebnerBasis[{Labs,D[Labs,x],D[Labs,y],D[Labs,z]},x,{y,z},MonomialOrder->EliminationOrder,Method->{"GroebnerWalk","EarlyEliminate"->True}])//Timing

(gb=GroebnerBasis[{Labs,D[Labs,x],D[Labs,y],D[Labs,z]},z,{x,y},MonomialOrder->EliminationOrder,Method->{"GroebnerWalk","EarlyEliminate"->True}])//Timing

GroebnerBasis[{Labs,D[Labs,x],D[Labs,y],D[Labs,z]},{x,y},z,MonomialOrder->EliminationOrder]//Timing

Approximate Attempt

Length[nodes=SingularPoints[Labs,{x,y,z},WorkingPrecision->200]]//Timing

Length[nodes=SingularPoints[Labs,{x,y,z},WorkingPrecision->500]]//Timing

Length[nodes=SingularPoints[Labs,{x,y,z},WorkingPrecision->1000]]//Timing

{Max[Norm/@nodes],Max/@Abs[Transpose[nodes]]}//N

soln=RootApproximant[nodes[[1]]]

N[Labs/.Thread[{x,y,z}->soln],2000]

nodes//N

Length[nodes=SingularPoints[Labs,{x,y,z},WorkingPrecision->1000]]//Timing

(soln=RootApproximant[nodes[[1]],32])//Timing

N[Labs/.Thread[{x,y,z}->soln],2000]

Resultant

zs=Union[RootApproximant/@Last/@nodes]

eqns={Labs,D[Labs,x],D[Labs,y],D[Labs,z]}/.z->zs[[2]];

Length[eqns]

eqnxy=RootReduce[{Labs,D[Labs,x],D[Labs,y],D[Labs,z]}/.z->zs[[2]]]

ByteCount/@eqnxy

GroebnerBasis[Take[eqnxy,2],y,x,MonomialOrder->EliminationOrder,Method->{"GroebnerWalk","EarlyEliminate"->True}]//Timing

GroebnerBasis[eqnxy,x,y,MonomialOrder->EliminationOrder,Method->{"GroebnerWalk","EarlyEliminate"->True}]

labs[z0_]:=Module[{eqns={Labs,D[Labs,x],D[Labs,y],D[Labs,z]}/.z->z0,xs,ys,poss,opts={MonomialOrder->EliminationOrder,Method->{"GroebnerWalk","EarlyEliminate"->True}}},
xs=Select[Union[x/.{ToRules[Roots[First[GroebnerBasis[Take[eqns,2],x,y,Sequence@@opts]]==0,x]]}],Im[#]==0.&];
Print["done gb x"];
ys=Select[Union[y/.{ToRules[Roots[First[GroebnerBasis[Take[eqns,2],y,x,Sequence@@opts]]==0,y]]}],Im[#]==0.&];
Print["done gb y"];
poss=Flatten[Outer[List,xs,ys,{z0}],2];
Print[Length[poss]," possible solutions..."];
Select[poss,(eqns/.{x->#[[1]],y->#[[2]]})==0.{1,1,1,1}&]
(*
(eqns/.{x->#[[1]],y->#[[2]]})&/@Flatten[Outer[List,{z0},xs,ys],2]
*)
]

labs[zs[[1]]]//Timing

Labs/.{x->0,y->0,z->-1}//RootReduce

labs[zs[[2]]]//Timing

(gb=GroebnerBasis[Take[eqns,2],x,y,MonomialOrder->EliminationOrder])//Timing

Roots[First[gb]==0,y]

(gb=GroebnerBasis[Take[eqns,-2],x,y,MonomialOrder->EliminationOrder])//Timing

(gb=GroebnerBasis[Take[eqns,2],y,x,MonomialOrder->EliminationOrder])//Timing

One

FindInstance[Labs==D[Labs,x]==D[Labs,y]==D[Labs,z]==0,{x,y,z},Reals]//Timing

Recognize

<<Utilities`Simplify`

?ExactRoot

?Recognize

?ToExact

soln=ExactRoot[#,32]&/@nodes[[1]]

ExactRoot

N[Labs/.Thread[{x,y,z}->soln],2000]//N

soln=ExactRoot[#,16,10^5]&/@nodes[[1]]

N[Labs/.Thread[{x,y,z}->soln],2000]//N

Exact

Length[nodes=SingularPoints[Labs,{x,y,z},WorkingPrecision->Infinity]]//Timing

...silently died...