6
$\begingroup$

Bug introduced after 9.0, persisting through 13.1. Bug report sent on 12/Nov/2022. WRI confirmed the Bug on 18/Nov/2022.


I'm trying to plot the following implicit equation in 3D:

Sqrt[x^2+y^2-z^2]+Sqrt[-x^2+y^2+z^2]+Sqrt[x^2-y^2+z^2]=Sqrt[2]

The code I used is:

ContourPlot3D[Sqrt[x^2 + y^2 - z^2] + Sqrt[x^2 - y^2 + z^2] + Sqrt[-x^2 + y^2 + z^2] == Sqrt[2], {x, 0, 1}, {y, 0, 1}, {z, 0, 1}, Mesh -> None, ContourStyle -> {Red, Opacity[0.5]}, MaxRecursion -> 3]

However, the resulting contour plot is weird and wrong - it has three extra surfaces:

The Plotting result using Mathematica

The correct result should be something like the following (made with python)

The Plotting result using Python

Could someone help me to identify the issue? Thank you!

Edit: I know the python result is correct because I tried to calculate a few explicit numerical solutions myself.

I also tried to plot an x-y intersection when z=0.4 to see how it looks in Mathematica:

k = 0.4;
ContourPlot[
 Sqrt[x^2 + y^2 - k^2] + Sqrt[x^2 - y^2 + k^2] + 
   Sqrt[-x^2 + y^2 + k^2] == Sqrt[2], {x, 0, 1}, {y, 0, 1}, 
 Mesh -> None, Axes -> False]

x-y intersection

The python code I used was (which is from this answer)

from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import numpy as np

def plot_implicit(fn, bbox=(0,1)):
    ''' create a plot of an implicit function
    fn  ...implicit function (plot where fn==0)
    bbox ..the x,y,and z limits of plotted interval'''
    xmin, xmax, ymin, ymax, zmin, zmax = bbox*3
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    A = np.linspace(xmin, xmax, 100) # resolution of the contour
    B = np.linspace(xmin, xmax, 15) # number of slices
    A1,A2 = np.meshgrid(A,A) # grid on which the contour is plotted

    for z in B: # plot contours in the XY plane
        X,Y = A1,A2
        Z = fn(X,Y,z)
        cset = ax.contour(X, Y, Z+z, [z], zdir='z')
        # [z] defines the only level to plot for this contour for this value of z

    for y in B: # plot contours in the XZ plane
        X,Z = A1,A2
        Y = fn(X,y,Z)
        cset = ax.contour(X, Y+y, Z, [y], zdir='y')

    for x in B: # plot contours in the YZ plane
        Y,Z = A1,A2
        X = fn(x,Y,Z)
        cset = ax.contour(X+x, Y, Z, [x], zdir='x')

    # must set plot limits because the contour will likely extend
    # way beyond the displayed level.  Otherwise matplotlib extends the plot limits
    # to encompass all values in the contour.
    ax.set_zlim3d(zmin,zmax)
    ax.set_xlim3d(xmin,xmax)
    ax.set_ylim3d(ymin,ymax)

    plt.show()

def surface(x,y,z):
    return np.sqrt(-x*x + y*y + z*z) + np.sqrt(x*x - y*y + z*z) + np.sqrt(x*x + y*y - z*z)- np.sqrt(2)

plot_implicit(surface)
$\endgroup$
2
  • $\begingroup$ Why do you think the python result is correct? How do you plot it in python? Can you add the corresponding python code? $\endgroup$
    – xzczd
    Nov 12, 2022 at 2:14
  • $\begingroup$ @xzczd, Thank you for your reply! Please see the edits in the OP. $\endgroup$
    – YYing
    Nov 12, 2022 at 2:37

4 Answers 4

9
$\begingroup$

Edit-2

Make the surface smooth.

Clear[sol, expr, plot, reg];
sol = SolveValues[{u^2 == x^2 + y^2 - z^2, v^2 == -x^2 + y^2 + z^2, 
    w^2 == x^2 - y^2 + z^2}, {x, y, z}];
expr = sol // Last
plot = ContourPlot3D[
    u + v + w == Sqrt[2], {u, 0, Sqrt[2]}, {v, 0, Sqrt[2]}, {w, 0, 
     Sqrt[2]}, MaxRecursion -> 4, 
    PlotPoints -> 80] /. {u_Real, v_Real, w_Real} -> expr;
reg = plot // DiscretizeGraphics;
RegionPlot3D[reg, PlotStyle -> {Opacity[.8], Red}, 
 BoundaryStyle -> {Thick, Green}, Axes -> True]

enter image description here

Clear[sol, plot];
sol = SolveValues[{u^2 == x^2 + y^2 - z^2, v^2 == -x^2 + y^2 + z^2, 
    w^2 == x^2 - y^2 + z^2}, {x, y, z}];
plot = ContourPlot3D[
   u + v + w == Sqrt[2], {u, 0, 2}, {v, 0, 2}, {w, 0, 2}, 
   MaxRecursion -> 4, PlotPoints -> 80, Mesh -> None, 
   ContourStyle -> Directive[Opacity[.8], Red]];
Show[(plot /. {u_Real, v_Real, w_Real} -> # & /@ sol), 
 PlotRange -> 1.1]

Edit-1

Use another change of variables.

Clear[sol, expr, reg, meshreg];
sol = SolveValues[{u^2 == x^2 + y^2 - z^2, v^2 == -x^2 + y^2 + z^2, 
    w^2 == x^2 - y^2 + z^2}, {x, y, z}];
expr = sol // Last;
reg = ParametricRegion[{expr, 
    u + v + w == Sqrt[2]}, {{u, 0, 2}, {v, 0, 2}, {w, 0, 2}}];
meshreg = DiscretizeRegion[reg];
Graphics3D[{EdgeForm[], FaceForm[{Red, Opacity[.8]}], meshreg}, 
 ViewPoint -> {1.36, -1.11, 2.89}, Axes -> True]

enter image description here

Clear[sol, regs, meshregs];
sol = SolveValues[{u^2 == x^2 + y^2 - z^2, v^2 == -x^2 + y^2 + z^2, 
    w^2 == x^2 - y^2 + z^2}, {x, y, z}];
regs = ParametricRegion[{#, 
      u + v + w == Sqrt[2]}, {{u, 0, 2}, {v, 0, 2}, {w, 0, 2}}] & /@ 
   sol;
meshregs = DiscretizeRegion /@ regs;
Graphics3D[{EdgeForm[], meshregs}, Axes -> True]

enter image description here

Edit-0

We can change of variables;

Clear[sol];
sol = SolveValues[{u == x^2 + y^2 - z^2, v == -x^2 + y^2 + z^2, 
    w == x^2 - y^2 + z^2}, {x, y, z}];
sol // Last

enter image description here

After that, we at least have two way to do the original plot.

  • ParametricRegion+ DiscretizeRegion.
Clear[reg,meshreg];
reg = ParametricRegion[{{Sqrt[u + w]/Sqrt[2], Sqrt[u + v]/Sqrt[2], 
     Sqrt[v + w]/Sqrt[2]}, 
    Sqrt[u] + Sqrt[v] + Sqrt[w] == Sqrt[2]}, {{u, 0, 2}, {v, 0, 
     2}, {w, 0, 2}}];
meshreg = DiscretizeRegion[reg];
RegionPlot3D[meshreg, BoundaryStyle -> None, 
 PlotStyle -> {Red, Opacity[0.8]}, ViewPoint -> {1, 1, 1}]

enter image description here

  • ContourPlot3D and change variables.
plot = ContourPlot3D[
   Sqrt[u] + Sqrt[v] + Sqrt[w] == Sqrt[2], {u, 0, 2}, {v, 0, 2}, {w, 
    0, 2}, Mesh -> None, ContourStyle -> White, MaxRecursion -> 2, 
   PlotPoints -> 30];
Show[plot /. {u_Real, v_Real, w_Real} -> {Sqrt[u + w]/Sqrt[2], Sqrt[
    u + v]/Sqrt[2], Sqrt[v + w]/Sqrt[2]}, 
 PlotRange -> {{0, 1}, {0, 1}, {0, 1}}, Lighting -> "ThreePoint"]

enter image description here

  • If we remove the restriction 0<=x<=1,0<=y<=1,0<=z<=1, the full surface is
Clear[sol,regs,L];
sol = SolveValues[{u == x^2 + y^2 - z^2, v == -x^2 + y^2 + z^2, 
    w == x^2 - y^2 + z^2}, {x, y, z}];
regs = ParametricRegion[{#, 
      Sqrt[u] + Sqrt[v] + Sqrt[w] == Sqrt[2]}, {{u, 0, 2}, {v, 0, 
       2}, {w, 0, 2}}] & /@ sol;
L = SolveValues[{Sqrt[x^2 + y^2 - z^2] + Sqrt[-x^2 + y^2 + z^2] + 
        Sqrt[x^2 - y^2 + z^2] == Sqrt[2], x == y == z}, {x, y, z}] // 
    First // Norm;
RegionPlot3D[DiscretizeRegion /@ regs // RegionUnion, 
 ColorFunction -> 
  Function[{x, y, z}, 
   Hue[Rescale[Sqrt[x^2 + y^2 + z^2], {L, Sqrt[3]}]]], 
 ColorFunctionScaling -> False, Axes -> True]

enter image description here

Clear[sol,plot];
sol = SolveValues[{u == x^2 + y^2 - z^2, v == -x^2 + y^2 + z^2, 
   w == x^2 - y^2 + z^2}, {x, y, z}]
plot = ContourPlot3D[
   Sqrt[u] + Sqrt[v] + Sqrt[w] == Sqrt[2], {u, 0, 2}, {v, 0, 2}, {w, 
    0, 2}, Mesh -> None, ContourStyle -> White];
Show[plot /. {u_Real, v_Real, w_Real} -> # & /@ sol, PlotRange -> 1.2]

enter image description here

$\endgroup$
6
  • $\begingroup$ Thank you very much for the two ways to do the parameterization! Do you know how to get the surface as PlotStyle -> {Red, Opacity[0.8]} without mesh? Sorry for such a basic question but I got confused about where to put this option to make it work. $\endgroup$
    – YYing
    Nov 12, 2022 at 13:28
  • $\begingroup$ @YYing Use BoundaryStyle -> None in RegionPlot3D or EdgeForm[] in Graphics3D can remove the mesh. see the updated. $\endgroup$
    – cvgmt
    Nov 12, 2022 at 13:38
  • $\begingroup$ Thank you for the editing with the Faceform and the additional info on specifying the style! Although the solution here may not work for a generic function that does not have a nice parameterization, it does solve the particular problem I have at hand. So I'll accept this answer. $\endgroup$
    – YYing
    Nov 12, 2022 at 13:38
  • $\begingroup$ sorry one last question: how to make the surface smoother? I tried to set MeshQualityGoal -> "Maximal", MaxCellMeasure -> 0.001, AccuracyGoal -> Infinity in the DiscretizeRegion function, but it doesn't seem to work. $\endgroup$
    – YYing
    Nov 12, 2022 at 13:57
  • $\begingroup$ @YYing I have also test several methods, but it is curious to me that non of them work. $\endgroup$
    – cvgmt
    Nov 12, 2022 at 14:02
3
$\begingroup$

Seems to be a bug introduced after v9. I guess it's essentially the same as this one. In v9, though the quality of plot isn't great, it's correct:

ContourPlot3D[
  Sqrt[x^2 + y^2 - z^2] + Sqrt[x^2 - y^2 + z^2] + Sqrt[-x^2 + y^2 + z^2] == 
   Sqrt[2], {x, 0, 1}, {y, 0, 1}, {z, 0, 1}, Mesh -> None, 
  ContourStyle -> {Red, Opacity[0.5]}, PlotPoints -> 50] // AbsoluteTiming

$Version

enter image description here

Please report it to WRI.

The following is a possible way to get the correct visualization in newer version:

dat = 
   ParallelTable[{x, y, z} /. 
      Solve[Sqrt[x^2 + y^2 - z^2] + Sqrt[x^2 - y^2 + z^2] + Sqrt[-x^2 + y^2 + z^2] ==
          Sqrt[2] // N, z], {x, 0, 1, 1/50}, {y, 0, 1, 1/50}] // 
    Flatten[#, 2] &; // AbsoluteTiming
(* {4.56397, Null} *)
ListPointPlot3D[dat, PlotRange -> {0, 1}, BoxRatios -> {1, 1, 1}]

Mathematica graphics

The following is a quick implementation for the idea in the python code:

eq = Sqrt[x^2 + y^2 - z^2] + Sqrt[x^2 - y^2 + z^2] + Sqrt[-x^2 + y^2 + z^2] == 
   Sqrt[2];

{Table[Normal@ContourPlot[eq // Evaluate, {x, 0, 1}, {z, 0, 1}] /. 
     Line[a_] :> Line[{#, y, #2} & @@@ a], {y, 0, 1, 1/25}], 
   Table[Normal@ContourPlot[eq // Evaluate, {y, 0, 1}, {z, 0, 1}] /. 
     Line[a_] :> Line[{x, #, #2} & @@@ a], {x, 0, 1, 1/25}], 
   Table[Normal@ContourPlot[eq // Evaluate, {x, 0, 1}, {y, 0, 1}] /. 
     Line[a_] :> Line[{#, #2, z} & @@@ a], {z, 0, 1, 1/25}]} /. 
  Graphics -> Graphics3D // Show[#, PlotRange -> All] &

enter image description here

Sadly I've no idea how to get the desired surface.

$\endgroup$
3
  • $\begingroup$ Thank you so much! I was thinking to use ListSurfacePlot3D instead of ListPointPlot3D in your mathematical code to get the surface, but it doesn't work - the resulting surface is weird. Do you know why or if it's possible to fix it? $\endgroup$
    – YYing
    Nov 12, 2022 at 12:51
  • $\begingroup$ @YYing The performance issue of ListSurfacePlot3D is a long-standing problem and there's no general enough work-around AFAIK, that's why I decide not to dive into it in my answer 囧. See e.g. mathematica.stackexchange.com/q/30608/1871 mathematica.stackexchange.com/a/109987/1871 $\endgroup$
    – xzczd
    Nov 12, 2022 at 13:25
  • $\begingroup$ I see! Thanks for letting me know! $\endgroup$
    – YYing
    Nov 12, 2022 at 13:29
3
$\begingroup$

Inspired by cvgmt

With[{para = Solve[{x^2 + y^2 - z^2, -x^2 + y^2 + z^2, x^2 - y^2 + z^2} == 
  {u, v, Sqrt[2] - u - v}^2, {x, y, z}][[-1, All, 2]]},
 ParametricPlot3D[para, {u, 0, Sqrt[2]}, {v, 0, Sqrt[2] - u}, 
  PlotStyle -> {Opacity[.8], Red}, BoundaryStyle -> Green, Mesh -> None
 ] /. Line -> (Tube[#, 0.006] &)
]

enter image description here

With[{para = Solve[{x^2 + y^2 - z^2, -x^2 + y^2 + z^2, x^2 - y^2 + z^2} == 
  {u, v, Sqrt[2] - u - v}^2, {x, y, z}][[All, All, 2]]},
 ParametricPlot3D[para, {u, 0, Sqrt[2]}, {v, 0, Sqrt[2] - u}, 
   BoundaryStyle -> Automatic, Mesh -> None
   ] /. Line -> (Tube[#, 0.01] &)
 ]

enter image description here

$\endgroup$
-1
$\begingroup$

Consider the equation:

Sqrt[x^2 + y^2 - z^2] + Sqrt[x^2 - y^2 + z^2] + 
  Sqrt[-x^2 + y^2 + z^2] == Sqrt[2]

This is a function of {x^2,y^2,z^2} and contains no linear terms. Therefore, the x-y, x-z, y-z planes are symmetry planes. And the result from Python can not be correct.

The full surface looks like:

ContourPlot3D[
 Sqrt[x^2 + y^2 - z^2] + Sqrt[x^2 - y^2 + z^2] + 
   Sqrt[-x^2 + y^2 + z^2] == Sqrt[2], {x, -1, 1}, {y, -1, 1}, {z, -1, 
  1}, Mesh -> None, ContourStyle -> {Red, Opacity[0.5]}, 
 MaxRecursion -> 3]

enter image description here

$\endgroup$
1
  • 4
    $\begingroup$ (-1)Notice the box ratio of python graphic is not 1:1:1. See also the discussions in the other 2 answers. $\endgroup$
    – xzczd
    Nov 12, 2022 at 9:23

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.