3
$\begingroup$

I need to create some surface plots, based on a function that is computationally expensive (a few minutes per point). The function has some branch cuts.

The command Plot3D does not work because it just takes too long -- never finishes at any acceptable resolution. I can get enough data points (over a day or three), put them in a matrix, and use ListPlot3D to produce a plot. But I don't see how to make Mathematica recognize the branch cut and not connect the two sides at the discontinuity.

Improve this question
$\endgroup$
2
  • $\begingroup$ Have you considered using RegionFunction to exclude a narrow strip around the discontinuity? $\endgroup$
    – Simon Woods
    Commented Oct 30, 2013 at 16:03
  • $\begingroup$ Hmm. Interesting idea. In my case, the branch cuts are not in any analytic form, so calculating them would be nontrivial. I did just think of another solution, making use of the fact that in my case, the branch cuts are closed loops. I could separate the data points "by hand" (values are greater or less than such and such) into two (or, in my case, three) separate matrices, and plot them with ListPlot3D on the same figure. Ad hoc, not elegant, and each figure would take some user effort, but I think doable. There should be a better way. $\endgroup$
    – Richard Tasgal
    Commented Oct 31, 2013 at 7:41

1 Answer 1

Reset to default
2
$\begingroup$

One option is simply to create the ListPlot3D and post-process it to remove polygons which are too steep.

data = With[{g := GaussianFilter[RandomReal[{-1, 1}, {50, 50}], 25]}, ArcTan[g, g]];
plot = ListPlot3D[data, Mesh -> None] // Normal;
toosteep[Polygon[pts_, ___]] := Module[{z = pts[[All, 3]]}, Max[z] - Min[z] > 1]
DeleteCases[plot, _Polygon?toosteep, -1]

enter image description here

Another possibility is to use a RegionFunction to exclude regions where the function gradient is high. Here I use a bit of image processing to locate the branch cuts:

exclude = Image[data] ~GradientFilter~ 1 ~Binarize~ 1.5 ~Dilation~ 1

enter image description here

regionfunc = ListInterpolation[ImageData@exclude, InterpolationOrder -> 0];
ListPlot3D[data, Mesh -> None, RegionFunction -> Function[{x, y, z}, regionfunc[y, x] < 1]]

enter image description here

Improve this answer
$\endgroup$
8
  • $\begingroup$ Thanks! I wonder if I can use the techniques in this library.wolfram.com/infocenter/MathSource/758 function to improve the accuracy of the location of the branch cuts (which in my case are smooth curves). (Apologies for not doing this myself and writing it here. I'm not that proficient in Mathematica. It will take me a few days to figure it out.) $\endgroup$
    – Richard Tasgal
    Commented Oct 31, 2013 at 16:34
  • $\begingroup$ @RichardTasgal, you could, but this of course increases the number of evaluations of your expensive function as you iteratively refine the location of the cuts. This is why Plot3D takes so long - it is increasing the sampling density near the discontinuity. $\endgroup$
    – Simon Woods
    Commented Oct 31, 2013 at 17:16
  • $\begingroup$ One other thing, if you are "not that proficient" in Mathematica it may be that your expensive function is coded inefficiently. What sort of computations are involved? $\endgroup$
    – Simon Woods
    Commented Oct 31, 2013 at 17:18
  • $\begingroup$ The function is as follows: I have a 6th-order polynomial equation where the coefficients depend on a few variables, a0(k,a,b)+a1(k,a,b)y+a2(k,a,b)y^2+a3(k,a,b)y^3+a4(k,a,b)y^4+a5(k,a,b)y^5+a6(k,a,b)y^6=0. The functions a_j(k,a,b) take a few tens of lines to define---involving functions of functions and a determinant---but they are algebraic expressions. Then I find the location k of the largest imaginary root of the 6 solutions y_j (Max,Imag,NSolve) over all k (NMaximize), holding a,b constant. I plot this over a,b. Compile doesn't help the speed. Maybe I'm missing a subtlety. $\endgroup$
    – Richard Tasgal
    Commented Nov 3, 2013 at 8:16
  • $\begingroup$ @RichardTasgal, it sounds like its just a really big calculation. I'd be happy to have a look if you can post a link to the notebook, but it may be that several minutes per point is as fast as it gets. $\endgroup$
    – Simon Woods
    Commented Nov 3, 2013 at 11:30

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.