2
$\begingroup$

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: enter image description here

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.

$\endgroup$
5
  • 1
    $\begingroup$ 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] $\endgroup$ May 30, 2022 at 20:55
  • 2
    $\begingroup$ 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 $\endgroup$
    – Michael E2
    May 30, 2022 at 23:33
  • $\begingroup$ There are at least three ways to do this, please post the Mathematica code. $\endgroup$
    – cvgmt
    May 31, 2022 at 4:01
  • $\begingroup$ @cvgmt what is your idea to plot the curve in the intesection? I posted the code I'm using, I hope it helps. $\endgroup$
    – User0212
    May 31, 2022 at 20:35
  • $\begingroup$ 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.) $\endgroup$
    – Michael E2
    May 31, 2022 at 23:40

1 Answer 1

2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.