not being able to find the a good picture of the intersection of two surfaces

Question

I'm having trouble finding a good picture around the origin of the intersection curve of the following surfaces $z - \dfrac{x^2}{2} - \dfrac{x^3}{12 \sqrt3} + \dfrac{y^2}{2} + \dfrac{x y^2}{4 \sqrt 3}=0$ and $x^2 + y^2 - 2 (1 + x + y) z + 3 z^2$. So far, the best try was using the code:

ContourPlot3D[{x^2 + y^2 - 2 (1 + x + y) z + 3 z^2, z - (x^2/2 - x^3/(12 Sqrt[3]) + y^2/2 + (x y^2)/(4 Sqrt[3]))}, {x, -1, 1}, {y, -1, 1}, {z, -1, 1},BoundaryStyle -> {1 -> None, 2 -> None, {1, 2} -> {Blue}},ContourStyle -> None, Mesh -> None]

which results in the following picture:

I believe the curve has a triple point at the origin, with this plot I can't identify if that is the case. Please, share some ideas on how to improve this.

It's a little annoying that I had to retype your equations because you didn't provide code that can be easily copied, but if you are using the technique from this answer, then you have something like ContourPlot3D[{z - x^2/2 - x^3/(12 Sqrt[3]) + y^2/2 + (x y^2)/(4 Sqrt[3]), x^2 + y^2 - 2 (1 + x + y) z + 3 z^2}, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, Contours -> {0}, ContourStyle -> None, Mesh -> None, BoundaryStyle -> {1 -> None, 2 -> None, {1, 2} -> {Blue}}, MaxRecursion -> 0, PlotPoints -> 45] — J. M.'s missing motivation, May 30 '22 at 20:55
People here generally like users to post code as Mathematica code instead of just images or TeX, so they can copy-paste it. It makes it convenient for them and more likely you will get someone to help you. You may find the meta Q&A, How to copy code from Mathematica so it looks good on this site, helpful — Michael E2, May 30 '22 at 23:33
There are at least three ways to do this, please post the Mathematica code. — cvgmt, May 31 '22 at 04:01
@cvgmt what is your idea to plot the curve in the intesection? I posted the code I'm using, I hope it helps. — User0212, May 31 '22 at 20:35
Adding the options MaxRecursion -> 8, PlotPoints -> 35 helps quite a bit. (Note that the Jacobian is of deficient rank at the origin, which explains the poor result near the origin.) — Michael E2, May 31 '22 at 23:40

score 2 · Accepted Answer · answered Jun 01 '22 at 01:30

The rank of the Jacobian of

eqns = {x^2 + y^2 - 2 (1 + x + y) z + 3 z^2,
   z - (x^2/2 - x^3/(12 Sqrt[3]) + y^2/2 + (x y^2)/(4 Sqrt[3]))};

as a mapping $F:{\bf R}^3 \rightarrow {\bf R}^2$ drops to $1$ at the origin, which is why numerical error in ContourPlot3D makes the intersection, given by the inverse image $F^{-1}(\{(0,0)\})$, appear to flatten out near the origin.

One way to trace the path is to start at some point (found with FindRoot) and use the cross product of the normals to the surface as the tangent velocity to intersection curve. The origin is still a problem, for the velocity is indeterminate there; but we get lucky and the solver steps over the bad point.

eqns = {x^2 + y^2 - 2 (1 + x + y) z + 3 z^2,
   z - (x^2/2 - x^3/(12 Sqrt[3]) + y^2/2 + (x y^2)/(4 Sqrt[3]))};
vars = Variables@eqns;
vel = Cross @@ D[eqns, {{x, y, z}}] // #/Sqrt[# . #] & // Simplify;
vel = vel /. v : x | y | z :> v[t];
boundaryEVT[x_, a_, b_] := {
   WhenEvent[x < a, "StopIntegration"], 
   WhenEvent[x > b, "StopIntegration"]};
curve = NDSolveValue[{
    D[Through[vars[t]], t] == vel,
    Through[vars[0]] == ({0.5, y, z} /. 
       FindRoot[eqns /. x -> 0.5, {y, -0.2}, {z, 0.3}]),
    boundaryEVT[x[t], -1, 1], boundaryEVT[y[t], -1, 1], 
    boundaryEVT[z[t], -1, 1]
    }, Through[vars[t]], {t, -3, 3}];
tdom = Flatten@{t, First[curve] /. t -> "Domain"};
pp = ParametricPlot3D[curve, Evaluate@tdom,
  PlotRange -> {{-1, 1}, {-1, 1}, {-1, 1}}, PlotStyle -> Magenta]

An alternative is to increase sampling in ContourPlot3D and hope it resolves the jaggedness of the curve. One might have to raise WorkingPrecision to handle the numerical problem at the origin, but again we get lucky (?).

cp = ContourPlot3D[
 {x^2 + y^2 - 2 (1 + x + y) z + 3 z^2 == 0,
  z - (x^2/2 - x^3/(12 Sqrt[3]) + y^2/2 + (x y^2)/(4 Sqrt[3])) == 0},
 {x, -1, 1}, {y, -1, 1}, {z, -1, 1},
 BoundaryStyle -> {1 -> None, 2 -> None, {1, 2} -> {Blue}},
 ContourStyle -> None, Mesh -> None,
 MaxRecursion -> 8, PlotPoints -> 35, AxesLabel -> Automatic]

Compare:

Show[cp, pp /. l_Line :> {Dashed, l}]

Note that

Method -> {"Projection", 
  "Invariants" -> (eqns /. v : x | y | z :> v[t])}

might normally be used in NDSolve to make the solution track the intersection curve with greater accuracy. However, if this is tried, one finds that the bad behavior of the Jacobian causes the numerical integration to stop when the solution reaches the origin.

not being able to find the a good picture of the intersection of two surfaces

1 Answers1