Can anyone explain this weird Plot3D error?

Question

This is more of a fun question because I see the source of the error, but cannot explain what Plot3D is doing. Can you figure it out?

L = 111;
X = Table[{0, 0}, {n, 1, L}];
X3D[n_, m_] := X[[Round[Sqrt[n^2 + m^2]]]][[2]];
Plot3D[X3D[n, m], {n, -50, 50}, {m, -50, 50}, PlotRange -> All, 
 ViewPoint -> Top]

Basically, I am taking a line of (0,0)s and swirling it into a plane. The error here is that I at some point, X3D[0,0] gets evaluated which calls for the 0th element in the X table -- and table elements start counting from 1, not 0.

But, here's the strange part. In fact, I ask for the 2nd part of the 0th element which leads to Plot3D generating the following pattern.

You can check for yourself that asking for the 1st part of the element list does NOT produce this cooky plot. It produces a flat plane with a missing point at (0,0), as it should:

This is a peculiarly strange result.

Update: A few people have suggested to use Exclusions->None to get rid of the holes. This smoothes over jumps so it's not a surprise that it works. But it doesn't explain why the interesting pattern appears in Plot3D.

Update again: Also a few people suggested that Round is the source of the errors. But if that was true then taking the first part of the elements would also produce the pattern, but it doesn't, as already mentioned. I also looked at adding a 3rd element and taking those, but this also produces the expected result.

Further update: Michael E2 showed below that changing the sampling for the function produced a different plot. There is still the question -- why does this not appear for the [[1]]-list, only the [[2]]-list. It's like a Mathematica Easter Egg.

You can even add more terms {0,0,0} and take the [[3]]-list (or higher) and this pattern does not appear.

Last update: See Michael E2 and Simon Woods answer below. It seems like Mathematica is using a different order of operations for the [[2]] vs every other possible list ([[1]] or [[3]], [[4]], etc.). Why it does this only Mathematica can answer.

Thanks guys, you and the Mathematica StackExchange site are awesome!

Answer 1

The behavior of Plot3D is due to the discontinuity processing associated to Round. (Note the pattern here has many of the same holes, but not exactly the same pattern.) The discontinuities occur whenever the argument to Round is a half-integer. Given the complexity of the argument to Round in this example, perhaps not all discontinuities are detected.

Plot3D[Round[Sqrt[n^2 + m^2]], {n, -50, 50}, {m, -50, 50}]

One can turn it off with Exclusions -> None.

Plot3D[Round[Sqrt[n^2 + m^2]], {n, -50, 50}, {m, -50, 50}, Exclusions -> None]

Here is what we get when we turn off discontinuity processing in the OP's example:

Plot3D[X3D[n, m], {n, -50, 50}, {m, -50, 50},
 PlotRange -> All, ViewPoint -> Top, Exclusions -> None]

Update. As is usual with such plotting behavior, sampling is the key. While it is often pointed out on the site, sometimes the misbehavior is so interesting that sampling is overlooked. With undersampling of so many discontinuities, it should be expected that gaps and bits of surface are mistakenly connected in a seemingly random, if symmetric, pattern.

With a high-enough setting for PlotPoints, we obtain a resolution that accurately detects all the gaps and connected pieces of the surface.

Plot3D[X3D[n, m], {n, -50, 50}, {m, -50, 50},
 PlotPoints -> 301, PlotRange -> All, ViewPoint -> Top]

(Of course, I'm simply ignoring the mistake of rounding a part index to zero.)

Update 2 -- A glimpse at the meshing of discontinuities

Simon Woods and I investigated some aspects of the meshing algorithms of DensityPlot and Plot3D respectively in our answers to Specific initial sample points for 3D plots. The interaction of exclusions with the meshing appears to be quite complicated at first glance, and I'm not sure I would be able to figure out the details if I tried. The basic idea is fairly simple, but it won't completely explain what we see. From our answers to the linked question, we see that meshing starts with a grid of rectangles subdivided into triangles by their SW-NE diagonals. This introduces a bias into the symmetry that is sometimes seen in the final result. The triangles are then subdivided, if the algorithm decides such a subdivision would improve the graph. This is done MaxRecursion number of times. As the triangles are subdivided, the bias tends to disappear. There appear to be pools of active triangle (being subdivided) and inactive ones (not being subdivided), but the subdivision of a triangle can activate its inactive neighbors.

Now introducing exclusions at first glance seems rather simple. Interpolating along the edges of a triangle, if it is found that two edges cross an exclusion, a narrow strip is cut from the triangle along the line connecting the crossing points. (This might introduce into the mesh some polygons with four or more sides.) If and edge crosses two exclusions or crosses one twice, the triangle is left alone. What is unclear is when this occurs, whether it happens only at the end of meshing the graph or recursively at intermediate points, and what further meshing occurs after it is done. Since parts of the surface are observed to come and go as MaxRecursion increases, the relationship does not seem simple to me. Further, the crossing rule I gave above does seem always to be observed, but that could be because an element became inactive at some stage when it did not appear to contain an exclusion.

Clear[showMesh];
showMesh[opts___] := Module[{a = 4, plotX, ptsX, excl},
   {plotX, {ptsX}} = Reap@Plot3D[X3D[n, m], {n, -a, a}, {m, -a, a},
      EvaluationMonitor :> Sow[{n, m, 0.01}],
      opts, PlotPoints -> 5,
      PlotRange -> All, Mesh -> All, ViewPoint -> Top];
   excl = 
    ParametricPlot3D[
     Evaluate[Table[r {Cos[t], Sin[t], 0.01}, {r, 0.5, 5.5}]],
     {t, 0, 2 Pi}, PlotStyle -> Thin];
   Show[plotX, excl,
    Graphics3D[{Red, PointSize[0.015], Point@ptsX}],
    PlotRange -> PlotRange[plotX]
    ]
   ];

GraphicsGrid@ Partition[Table[showMesh[MaxRecursion -> mr], {mr, 0, 3}], 2]

1. Further recursive subdivision does not change the appearance of the exclusions.

2. I cannot explain why some of the mesh elements appear to recombined in moving from MaxRecursion -> 1 to MaxRecursion -> 2. It is this that makes me think there is a subdivision that occurs after the exclusions are computed.

3. The excluded portion gets wider as MaxRecursion steps from 1 to 3. There are other, bigger changes, but the erosion is noticeable. I don't know if that is connected to the disappearance of big chunks of the surface after each step. As I said, further recursion does not increase the gaps. It seems mainly concentrated on refining the edge of the hole in the center.

4. Sometimes a triangle has darkened shading. This is (usually) because two vertices of the polygon are nearly equal and it confuses the rendering algorithm (due to round-off error perhaps).

5. This all informed my strategy for my initial answer. I tried to get an initial mesh so that any mesh-edge path that crossed two excluded circles would have a segment with vertices lying between the two circles. The distance between the circles is 1 and the length of the diagonal of the domain is less than 142. Doubling that and rounding up gave me 300 (or PlotPoints -> 301), so that the edges of the mesh will be less than 1/2.

6. I view the OP's problem as one of insufficient input. All numerical solvers occasionally need extra help with step size, precision, method selection, etc. This is essentially a case of initial step size insufficient for resolving the details of the graph.

Yet another update -- [[1]] vs [[2]]

The difference between using part [[1]] vs. part [[2]] has to do with the symbolic processing of Plot3D. It evaluates the function symbolically to determine exclusions, for instance. Here is what you get in each case:

X3D1[n_, m_] := X[[Round[Sqrt[n^2 + m^2]]]][[1]];
X3D2[n_, m_] := X[[Round[Sqrt[n^2 + m^2]]]][[2]];

X3D1[m, n]

Part::pkspec1: The expression Round[Sqrt[m^2+n^2]] cannot be used as a part specification. >>

(*  {{0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, ..., {0, 0}}  *)

X3D2[m, n]

Part::pkspec1: The expression Round[Sqrt[m^2+n^2]] cannot be used as a part specification. >>

(*  Round[Sqrt[m^2 + n^2]]  *)

So in the first case Plot3D does not see Round, and in the second case it sees Round. So there are no exclusion in the first case but there are in the second.

Answer 2

If you play with the PlotPoints option you get more examples of this:

Answer 3

Exclusions detection evaluates X3D[n,m] symbolically to analyse it for discontinuities. In your second example this returns the list X:

X[[Round[Sqrt[n^2 + m^2]]]][[1]]

(* {{0, 0}, {0, 0}, {0, 0}, .... *)

That's a meaningless result for discontinuity analysis so it is ignored and the plot continues with no exclusions. But in your first example you get this:

X[[Round[Sqrt[n^2 + m^2]]]][[2]]

(* Round[Sqrt[m^2 + n^2]] *)

This is expected of course - Part returns unevaluated and then you extract the second element of Part[X, Round[Sqrt[m^2 + n^2]]]

So in this case the evaluation of X3D[n,m] with symbolic arguments returns a perfectly good expression which can be analysed for discontinuities, and the plot is created with exclusions at those locations.

You can check the logic by modifying X3D to return a completely different function when evaluated with symbolic arguments:

Clear[X3D]

X3D[n_?NumericQ, m_?NumericQ] := X[[Round[Sqrt[n^2 + m^2]]]][[2]];

X3D[n_Symbol, m_Symbol] := Round[n]

Plot3D[X3D[n, m], {n, -50, 50}, {m, -50, 50}, PlotRange -> All, 
 ViewPoint -> Top]

Here the numerical function is plotted but with exclusions determined from the symbolic function.