6
$\begingroup$

I've got a poor quality when plot the implicit equation singu==0, with

singu=1+ 2. y^5 z - 3. z^2 + 3. z^4 + x^5 (2. y + 2. z) + 
 y^3 (-1. z + 4. z^3) + 
 y^2 (-3. + 9. z^2 - 6. z^4 - 3. z Sqrt[1 - x^2 - y^2 - z^2]) + 
 x^4 (3. - 6. y^2 - 6. z^2 + 2. y Sqrt[1 - x^2 - y^2 - z^2] - 
    2. z Sqrt[1 - x^2 - y^2 - z^2]) + 
 y^4 (3. - 6. z^2 + 2. z Sqrt[1 - x^2 - y^2 - z^2]) + 
 y z (-1. - 1. z^2 + 2. z^4 + 3. z Sqrt[1 - x^2 - y^2 - z^2] - 
    2. z^3 Sqrt[1 - x^2 - y^2 - z^2]) + 
 x^3 (4. y^3 - 1. z + 4. y^2 z + 4. z^3 + y (-1. + 4. z^2)) + 
 x^2 (-3. - 6. y^4 + 4. y^3 z + 9. z^2 - 6. z^4 + 
    y^2 (9. - 12. z^2) + 3. z Sqrt[1 - x^2 - y^2 - z^2] + 
    y (2. z + 4. z^3 - 3. Sqrt[1 - x^2 - y^2 - z^2])) + 
 x (2. y^5 - 2. y^4 Sqrt[1 - x^2 - y^2 - z^2] + y (-1. + 2. z^2) + 
    y^3 (-1. + 4. z^2) + 
    y^2 (2. z + 4. z^3 + 3. Sqrt[1 - x^2 - y^2 - z^2]) + 
    z (-1. - 1. z^2 + 2. z^4 - 3. z Sqrt[1 - x^2 - y^2 - z^2] + 
       2. z^3 Sqrt[1 - x^2 - y^2 - z^2]))

or equivalently,

singu= 1 +2. y^5 z-3. z^2+3. z^4+x^5 (2. y+2. z)+y^3 z (-1.+4. z^2)+y^2 (-3.+9. z^2-6. z^4-3. z Sqrt[1-x^2-y^2-z^2])+x^4 (3. -6. y^2-6. z^2+2. y Sqrt[1-x^2-y^2-z^2]-2. z Sqrt[1-x^2-y^2-z^2])+y^4 (3. -6. z^2+2. z Sqrt[1-x^2-y^2-z^2])+x^3 (4. y^3-1. z+4. y^2 z+4. z^3+y (-1.+4. z^2))+x^2 (-3.-6. y^4+4. y^3 z+9. z^2-6. z^4+y^2 (9. -12. z^2)+3. z Sqrt[1-x^2-y^2-z^2]+y (2. z+4. z^3-3. Sqrt[1-x^2-y^2-z^2]))+y z (-1.+z (3. Sqrt[1-x^2-y^2-z^2]+z (-1.+2. z^2-2. z Sqrt[1-x^2-y^2-z^2])))+x (2. y^5-2. y^4 Sqrt[1-x^2-y^2-z^2]+y (-1.+2. z^2)+y^3 (-1.+4. z^2)+y^2 (2. z+4. z^3+3. Sqrt[1-x^2-y^2-z^2])+z (-1.+z (-3. Sqrt[1-x^2-y^2-z^2]+z (-1.+2. z^2+2. z Sqrt[1-x^2-y^2-z^2]))))

I used the following command:

Table[ContourPlot3D[
  Evaluate@singu == 0, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, 
  PerformanceGoal -> "Quality", MaxRecursion -> 1, 
  PlotPoints -> pp], {pp, {15}}] 

which has a very poor quality, as shown in Fig[Poor quality ContourPlot].

Theoretically the origin shouldn't be part of the surface, moreover, the surface in the figure is very broken. So is there anyway to improve the quality? I have tried increasing the number of MaxRecursion to 2, and more plotpoints, it takes much longer, but does not help much. This expression looks not too difficult to me, so I am a bit puzzled. Could someone give me some suggestions to improve the plot quality, and preferably tips to reduce the evaluation time? Thanks a lot!

Update: I tried to increase the number of plotpoints and maximumRecursion, and I found that the quality becomes better but some of broken parts still exists. Then I plotted the unit sphere, which is the constraint to make Sqrt[1-x^2-y^2-z^2] real. Then I notice that seems the desired surface might intersect the unit sphere, which might be the reason why the plot is broken.

So I wonder if it is possible to plot the equation only within the range of the unit sphere? is there anyway to specify this constraint to the ContourPlot3d command? Or is there any other way to realize this? Thanks a lot!

Update 2: I have tried the FegionFunction to limit the figure within the unit sphere, however, the plot does not change, it seems that this does not affect the evaluation; then I tried collecting the square root term in the equation and put it into the form: a=-b Sqrt[1-x^2-y^2-z^2] then square each side to obtain a new equation free of square root, but then I notice that the new plot includes some more regions corresponding to a=b Sqrt[1-x^2-y^2-z^2], and I can not distinguish them. So is there anyway to deal with this?

Or do you suggest working on the original expression? and is there some way to deal with the broken parts, other than the RegionFunction or increasing the PlotPoints ? Thanks a lot!

Update 3: The solution from Michael E2 is quite satisfying, which solves the problem perfectly. I just have one more question: Is there a direct way to plot the surface generated by an equation involving the spherical coordinates?

I notice that in this algorithm, the graph is plotted first by ContourPlot3d in the CARTESIAN frame first, with the coordinates as r,theta,phi, then transformed in some way to make the frame coordinates into x,y,z. So there is not a direct way to plot the surface generated by an equation involving the spherical coordinates? Thanks a lot!

$\endgroup$
  • 2
    $\begingroup$ You are shooting yourself in the foot by simultaneously restricting the MaxRecursion and using a low PlotPoint count. By setting MaxRecursion -> 1, it can only add one extra layer of points. $\endgroup$ – rcollyer Jan 16 '17 at 21:07
  • $\begingroup$ @rcollyer, thanks for your reply, the reason I use these values is because its evaluation takes long when I increase them, and the examples I found use these values. Since I'm not quite familiar with mathematica yet, could you please give me some recommendations for these values? Also is there anyway to make it faster? E.g., not using the ContourPlot3d command? I know matlab can do the same thing with the trick of isosurface, but I'm not sure if mathematica has a similar way, nor do I know its efficiency. Do you know anything about it? Thanks again! $\endgroup$ – larry Jan 16 '17 at 21:34
  • 1
    $\begingroup$ @larry What is "the trick of isosurface"? $\endgroup$ – C. E. Jan 16 '17 at 21:49
  • 1
    $\begingroup$ @larry It's not a trick then, just a MATLAB function. The corresponding function in Mathematica is ListContourPlot3D. Basically, instead of using adaptive sampling it's using a fixed grid. It should actually be more expensive to get a similar quality figure using that method as compared to what ContourPlot3D is doing. But MATLAB may be more efficient, so I'm not saying that you're wrong. $\endgroup$ – C. E. Jan 17 '17 at 1:55
  • 2
    $\begingroup$ @bbgodfrey WRI beta-tests new versions with some users. I mean that one. I believe internally WRI will have several development versions that they use to evaluate various alternative implementations of fixes/enhancements. $\endgroup$ – Michael E2 Jan 17 '17 at 12:11
9
$\begingroup$

Here's a roundabout way: Plot in spherical coordinates and transform back. There are some slight imperfections, esp. where the boundaries join.

sph = TransformedField[ "Cartesian" -> "Spherical", 
   singu, {x, y, z} -> {r, θ, ϕ}] // Simplify

(* coordinate and vector field (VertexNormals) transformations *)
cXF = CoordinateTransformData["Spherical" -> "Cartesian", "Mapping"];
vXF[{x_, y_, z_}, {a_, b_, c_}] = 
  TransformedField["Spherical" -> "Cartesian", {a, b, c}, {r, θ, ϕ} -> {x, y, z}];

(* spherical plot *)
cpSPH = ContourPlot3D[
   sph == 0, {r, 0, 1}, {θ, 0, Pi}, {ϕ, 0, 2 Pi}, 
   MeshFunctions -> (* cartesian mesh *)
    Thread[cXF /. HoldPattern[Slot[1][[n_]]] :> Slot[n]]];

(* transform back to cartesian *)
cpCAR = Show[
   cpSPH /. 
    GraphicsComplex[p_, g_, rest___] :> 
     With[{x = Transpose@cXF@Transpose@p},
      GraphicsComplex[
       x,
       g,
       VertexNormals -> Transpose[vXF[Transpose@x, Transpose[VertexNormals /. {rest}]]],
       rest
       ]],
   PlotRange -> 1]

Mathematica graphics

$\endgroup$
  • $\begingroup$ Great result (+1). Why do you suppose that ContourPlot3D works so much better in spherical than Cartesian coordinates? $\endgroup$ – bbgodfrey Jan 17 '17 at 5:27
  • 1
    $\begingroup$ @bbgodfrey First, compare cpSPH with cpCAR; the behavior near r == 1 is important. (1) In spherical, we can keep r <= 1 and r == 1 is a "natural" sampling boundary; this doesn't happen in cartesian. I think ContourPlot3D gets confused about how to construct the surface near r == 1 in cartesian, because the function is undefined/complex near it. (2) The geometry of the thin "flanges" near r == 1 looks easier to resolve in spherical. But that's also due to (1) & maybe rounding error. $\endgroup$ – Michael E2 Jan 17 '17 at 12:28
  • $\begingroup$ @MichaelE2, hi, sorry for the delay of reply, because I had to spend some time to understand your code. Thanks a lot for the suggestion, which works perfect. Now there is one last question in the update of the problem. Could you please take a look at it? $\endgroup$ – larry Jan 17 '17 at 20:16
  • 1
    $\begingroup$ @larry I don't know of a way to plot implicit equations directly in other coordinate systems. AFAIK, ContourPlot3D and RegionPlot3D are the only ways to plot implicit surfaces, and they use cartesian coordinates. $\endgroup$ – Michael E2 Jan 18 '17 at 3:47
  • $\begingroup$ @MichaelE2 fine, your code works perfect already. Thanks again for your help! $\endgroup$ – larry Jan 18 '17 at 17:57
7
$\begingroup$

SliceContourPlot3D is not as elegant as ContourPlot3D, but it does display the complexity of the surface, and in less than a minute. It is this complexity that gives ContourPlot3D problems unless PlotPoints is large, in which case the computation crashes, at least on my PC.

s = Collect[singu // Rationalize, Sqrt[_], Simplify];
Show[SliceContourPlot3D[s, {"ZStackedPlanes", {#}}, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, 
    PerformanceGoal -> "Quality", MaxRecursion -> 4, PlotPoints -> 50,
    Contours -> {0}, ContourShading -> Opacity[0], 
    ContourStyle -> Directive[Hue[(# + 1)/2], Thick]] & /@ 
    Range[-.75, .75, .15]]
SwatchLegend[Hue[(# + 1)/2] & /@ Range[-.75, .75, .15], Range[-.75, .75, .15], 
    LegendLayout -> "Row"]

enter image description here enter image description here

Note that s was computed from signu using Collect with Simplify to improves speed. Multiple ``SliceContourPlot3D` plots were generated and combined to give contours on each slice a different color.

$\endgroup$
  • $\begingroup$ Hi, I do not quite catch what you are plotting. Is this the intersection of the surface with the unit sphere? Thanks! $\endgroup$ – larry Jan 17 '17 at 2:10
  • $\begingroup$ @larry It is a stack of nine 2D ContourPlots for z at values given by Range[-.8, .8, .2]. The sphere is at x^2 + y^2 + z^2 == 1 and was added by SliceContourPlot3D to show where s begins to evaluate to complex numbers. By the way, I have experimented with several sets of options for ContourPlot3D but none produce a plot much better than what you have in the question. $\endgroup$ – bbgodfrey Jan 17 '17 at 2:19
  • $\begingroup$ oh, I now understand what you are plotting now. Thanks a lot! $\endgroup$ – larry Jan 17 '17 at 20:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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