20
$\begingroup$

I really can't understand this.

I wrote two functions, one with Compile called funcom, one with librarylink to link a fortran subroutine called funcfortran. And they do exactly the same thing! So I Plot the result first

plotfor=Plot3D[funcfortran[w, 0.06, -1, 1, 1, \[Pi]/2., 2., 3., ky], {w, 0, 
  10}, {ky, -0.1, 0.1}, PlotRange -> All, MaxRecursion -> 7, 
 Mesh -> All]

This is output from topview

enter image description here

Well, this may not seem that wrong, But look at the plot result of funcom with exactly the same parameter

plotcom=Plot3D[funcom[w, 0.06, -1, 1, 1, \[Pi]/2., 2., 3., ky], {w, 0, 
      10}, {ky, -0.1, 0.1}, PlotRange -> All, MaxRecursion -> 7, 
     Mesh -> All]

enter image description here

different !!!

They contains almost the same point data, this can be verified after we extract the data and compared the difference

Sort[DeleteDuplicates@
   Abs@Flatten[
     Cases[%58, x_GraphicsComplex :> x[[1]]][[1]] - 
      Cases[%66, x_GraphicsComplex :> x[[1]]][[1]]], Greater][[1 ;; 3]
 ]

The top three biggest differen

{8.00339*10^-8, 7.99922*10^-8, 7.99703*10^-8}

If I plot it on larger region, the plot of funcfortran is unacceptable, see following and corresponding Mesh->All version

Plot3D[funcfortran[w, 0.06, -1, 1, 1, \[Pi]/2, 2, 3., ky], {w, 0, 
  10}, {ky, -\[Pi], \[Pi]}, PlotRange -> All, MaxRecursion -> 7]

enter image description hereenter image description here

and this is what funcom got

enter image description here

Again, the point data contained in these two plot is almost the same, funcfortran contains 8260 points and funcom contains 8264 points

On the other hand, I can "stupidly" calculation funcfortran on a regular rectangular mesh, and ListPlot3D get pretty nice result.

Block[{datatmp, klist, n1, n2},
  n1 = n2 = 100;
  klist = 
   Tuples[{Subdivide[-N@\[Pi], N@\[Pi], n1], 
     Subdivide[-N@\[Pi], N@\[Pi], n2]}];
  datatmp = 
   Flatten[Outer[
     funcfortran[#1, 0.005, -1., 1., 1., \[Pi]/2., 2., 3., #2] &, 
     Subdivide[0., 10., n1], Subdivide[-N@\[Pi], N@\[Pi], n2]], 1];
  data = Join[klist, Partition[datatmp, 1], 2]];

ListPlot3D[data, PlotRange -> All]

enter image description here

So, strange things, it seems that funcfortran has no problem, because ListPlot3D gives good result. But why Plot3D fails?


update

m_goldberg suggested that this maybe due to loss of precision of my fortran routine. But I want to demonstrate, How Plot3D is defective.

I choose a defective parameter region, and plot it

test = Plot3D[
  iterateG[w, 0.06, -1, 1, 1, \[Pi]/2, 2, 3., ky], {w, 2, 3}, {ky, 0, 
   0.5}, PlotRange -> All, MaxRecursion -> 7, Mesh -> None, 
  PlotPoints -> 50]

it outputs

enter image description here

To notice the two strange dark stripes.

Now we extract the data contained in test plot, and use ListPlot3D. As m_goldberg had pointed out, there is default interpolation, we can turn it off using InterpolationOrder -> 1

ListPlot3D[Cases[test, x_GraphicsComplex :> x[[1]]][[1]], 
 InterpolationOrder -> 1, Mesh -> None]

outputs

enter image description here

It is smooth!! And this time, same data, different behaviour between Plot3D and ListPlot3D!!


update

I attached a zip (download here onedrive). Since librarylink is quite difficult to work with. So I have packed everything: .dll for win, .so for linux, source .f90, .nb etc. Hope everyone extract and open the .nb file should have no problem running it. Thank you for testing.


Summary

The bug hunting is over. The defective plot is due to my fortran code and off course my bad fortran coding. I wrote cmplx instead of dcmplx which cause the rounding of w parameter, loss of precision and finally the weird plot3d.

Many thanks for Jason B's kind help, I learned a lot of skill and insight for tracing such kind of bug. Also I recall that m_goldberg is the first one correctly pointed out that it must be due to loss of precision. I reget not pay enough attention to this. Finally thank all people who have concerned with this post and tried to help me.

$\endgroup$
7
  • 1
    $\begingroup$ It looks to me that your Fortran function is having loss of precision problems. I think you should analyze the numeric behavior of the two functions; it may enlighten you as to the source of your plotting difficulties. I don't think anyone reading your question can help much further without having access to your Fortran code and its running environment. $\endgroup$
    – m_goldberg
    Jan 18, 2016 at 3:10
  • $\begingroup$ Further, I don't think the list plot casts any light on your issue. I suspect it interpolates your data and thus smooths it out nicely. Of course, if you like the look of the list plot, then that's the way for you to plot the Fortran function. $\endgroup$
    – m_goldberg
    Jan 18, 2016 at 3:36
  • $\begingroup$ @m_goldberg Hi, m_goldberg, thank you for comment. But how to "analyze the numeric behavior of the two functions", I have no idea $\endgroup$
    – matheorem
    Jan 18, 2016 at 4:31
  • $\begingroup$ @m_goldberg I updated my post, to show that interpolation is not the reason. $\endgroup$
    – matheorem
    Jan 18, 2016 at 5:27
  • 1
    $\begingroup$ @matheorem - is there any way you can post the code for both funcom and funcfortran? It looks like this has a couple of close votes, but I think this is a legitimate issue you are bringing up. And for the record, rather than being "stupid", I think it is almost always better to generate a rectangular grid and plot based on some form of ListPlot than it is to use the corresponding Plot. Basically, there are some decision making algorithms that decide which points to sample when making a Plot and you can circumvent those by sampling yourself. $\endgroup$
    – Jason B.
    Jan 18, 2016 at 8:24

1 Answer 1

26
+200
$\begingroup$

Source of the problem (possibly)

Here is a clear indication that your Fortran library and the Mathematica function are behaving in fundamentally different ways. I noticed the apparent high frequency oscillations in the difference functional, so I decided to see exactly how quickly they oscillate,

Plot[funcfortran[w, 0.06, -1.0, 1.0, 1.0, N[π/2], 2.0, 3., 0.2] - 
    funcom[w + 0.06*I, π/2., 3., .2], {w, 2.75 - #, 2.75 + #}, 
   PlotPoints -> 800, 
   ImageSize -> 400] & /@ {.001, .0001, .00001, .000001}

enter image description here

So now we zoom in and plot them together on this scale,

Plot[{funcfortran[2.75 + dw, 0.06, -1.0, 1.0, 1.0, N[π/2], 2.0, 
   3., 0.2],
  funcom[2.75 + dw + 0.06*I, π/2., 3., .2]},
 {dw, -.000001, .000001},
 PlotPoints -> 200,
 ImageSize -> 400,
 PlotLegends -> {"Fortran", "Mathematica"}]

enter image description here

So somewhere, Fortran is doing some kind of rounding. Perhaps you have some number defined with single precision? Probably not that simple or you would have caught it, but basically as you vary w in increments of the order $10^{-7}$ then the Fortran function does not vary smoothly. This is not the case for the ky parameter.

I would next check whether you get this behavior from Fortran directly, without using Mathematica. If so, the problem is in your code there. If not, it must have to do with the library linking function.

As I try to show below, I think it is this nonlinear behavior in the Fortran function that leads Mathematica to plot it incorrectly when using the adaptive grid created by ListPlot3D. I assume that Plot3D is trying to come up with derivatives to better plot, but at some points along the w axis the derivative is infinite.

You point out that Plot adaptively samples the plot region, so let's just extract the points that you get for both functions,

fortranlist = 
  Reap[Plot3D[
     Sow[{w, ky, funcfortran[w, 0.06, -1, 1, 1, π/2, 2, 3., ky]}];
      funcfortran[w, 0.06, -1, 1, 1, π/2, 2, 3., ky], {w, 0, 
      10}, {ky, -π, π}, PlotRange -> All, MaxRecursion -> 7, 
     PlotPoints -> 100]][[2, 1]];
funcomlist = 
  Reap[Plot3D[Sow[{w, ky, funcom[w + 0.06*I, π/2., 3., ky]}]; 
     funcom[w + 0.06*I, π/2., 3., ky], {w, 0, 
      10}, {ky, -π, π}, PlotRange -> All, MaxRecursion -> 7, 
     PlotPoints -> 100]][[2, 1]];
fortranlist[[2]]
fortranlist = Delete[evalcoords, 2];

That last line is necessary because the second element of fortranlist looks like this,

enter image description here

For those who want to weigh in, but don't want to evaluate the Fortran library function, you could just download the lists above via:

funcomlist=Import["https://gist.githubusercontent.com/anonymous/34581877899db27e1cf4/raw/409c85ce8d34cd8f90344cbd0291c08f29c68d9e/funcom_pts.dat","Table"];
fortranlist=Import["https://gist.githubusercontent.com/anonymous/1d377c32f675847839d6/raw/00f79e303e5ec787107082516bffe3148af8380e/fortran_pts.dat","Table"];

But we can now check that the adaptive grid for the two functions is identical,

funcomlist[[All, ;; 2]] == fortranlist[[All, ;; 2]]
(* True *)

And we can plot the two lists and their differences

ListPlot3D[#, ImageSize -> 450, PlotRange -> All] & /@ {funcomlist, 
  fortranlist, 
  Transpose[{funcomlist[[All, 1]], funcomlist[[All, 2]], 
    funcomlist[[All, 3]] - fortranlist[[All, 3]]}]}

enter image description here

And here is the density of grid points,

ListPlot[funcomlist[[All, ;; 2]]]

enter image description here

What I notice is that the regions where the differences between the function is largest is also the region with the highest density of grid points. So in that region, Mathematica is trying to get many points so it can get a handle on the curvature in that region. The Fortran function is not behaving well there, but this misbehavior is at a relatively small scale. But the scale isn't that small. The average value for the function is on the order of .1, and the differences are a factor of $10^{-6}$ smaller than that, but $10^{-6}$ is a good deal larger than machine precision. So Plot3D is trying to get a handle on the higher order derivatives in that region, and isn't doing such a great job.

Why are there such large differences between the functions? I'm not going to check through your Fortran code, but somewhere it has a difference.

So in this case, the adaptive sampling is working against you. It is preferentially sampling a region where the precision of the Fortran function is less than optimal. No amount of increasing MaxRecursion or PlotPoints will fix this, because they will just increase the sampling in that region.

Another point is that, at least in version 10, the interpolation algorithms for plotting functions like ListContourPlot, ListPlot3D, etc, do a much better job when given a rectangular grid than they do when given a list of tuples like {x, y, f[x, y]}. See my question here. And all Plot3D seems to do with the data it generates is to give it to ListPlot3D. So you already know that you get a much better plot simply from creating a list and plotting it, but I would suggest not structuring that list in the form of tuples, but instead as a grid.

ListPlot3D[
 Table[funcfortran[w, 0.06, -1, 1, 1, π/2, 2, 3., 
   ky], {ky, -π, π, .05}, {w, 0, 10, .1}], 
 DataRange -> {{0, 10}, {-π, π}}]

It produces a better plot than Plot3D, and does it faster. I'd be interested in hearing what @user21 thinks of this issue, I believe he works directly with the underlying functions here.

Edit

To your final point, I don't think that Cases is giving you every point that is plotted by Plot3D. Consider this:

fortranlist2 = Reap[

    fortranplot = 
      Plot3D[Sow[{w, ky, 
         funcfortran[w, 0.06, -1, 1, 1, π/2, 2, 3., ky]}]; 
       funcfortran[w, 0.06, -1, 1, 1, π/2, 2, 3., ky], {w, 2, 
        3}, {ky, 0, 0.5}, PlotRange -> All, MaxRecursion -> 7, 
       Mesh -> None, PlotPoints -> 50];
    ][[2, 1]];
fortranlist2 = Delete[fortranlist2, 2];
fortranlist2b = Cases[fortranplot, x_GraphicsComplex :> x[[1]]][[1]];

fortranlist2 was generated while plotting, using the Reap and Sow method, while fortranlist2b was generated using Cases off of the graphics object. We can see that there is a lot more information in the first one,

Length /@ {fortranlist2, fortranlist2b}
(* {10337, 2584} *)

They both seem to conform to an almost rectangular grid, but clearly Plot3D corresponds to the data in the Reap and Sow method,

{fortranplot, ListPlot3D[fortranlist2, Mesh -> None], 
 ListPlot3D[fortranlist2b, Mesh -> None]}

enter image description here

If you evaluate the differences between the Mathematica and Fortran functions at the grid points defined above, then you see where your stripes come from

diff = {#1, #2, 
     funcfortran[#1, 0.06, -1, 1, 1, π/2, 2, 3., #2] - 
      funcom[#1 + 0.06*I, π/2., 3., #2]} & @@@ 
   fortranlist2[[All, ;; 2]];

ListPlot3D[diff]

enter image description here

$\endgroup$
20
  • 2
    $\begingroup$ Wow, finally someone gives an answer. Thank you Jason B for spending your time helping me. I am looking into your answer right now. Reap and Sow is showing something, but I think you maybe overlooked those dark stripes in the 3rd plot, they exactly corresponding those defective part of 2nd plot. But on these dark strips, there is no value difference, they are flat. I think this should be the key problem instead of those 10-6 differences, what do you think? $\endgroup$
    – matheorem
    Feb 1, 2016 at 13:38
  • $\begingroup$ @matheorem I notice that the parameters you pass to funcfortran are not all reals (however they should) ... Just to be sure, does it make any difference ? $\endgroup$
    – SquareOne
    Feb 1, 2016 at 14:22
  • $\begingroup$ @SquareOne Sorry for the confusion. The parameter of funcfortran is all real, while the w in funcom is complex. But it is just a sight difference in definition. I should had made it real also. But anyway, I just tried a modified real version of funcom, the result is the same $\endgroup$
    – matheorem
    Feb 1, 2016 at 14:40
  • $\begingroup$ @matheorem I actually meant : is there a difference between calling funcfortran[w, 0.06, -1, 1, 1, π/2, 2, 3., ky] and funcfortran[w, 0.06, -1.0, 1.0, 1.0, π/2., 2., 3., ky] where all the param. are reals ? $\endgroup$
    – SquareOne
    Feb 1, 2016 at 14:51
  • 2
    $\begingroup$ What happens when you have the fortran code natively spit out a list over the values above? is it in the stairstep pattern, or is it continuous? $\endgroup$
    – Jason B.
    Feb 1, 2016 at 16:31

Your Answer

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

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