I am trying to plot $J_m(\chi_{mn} r) \cos(m\phi),$ where $J_m(\chi_{mn}) =0,$ as a DensityPlot in polar coordinates (solution to 2D wave equation in polar coordinates). First, some definitions

besselzero = N@Table[BesselJZero[l, n], {l, 0, 25}, {n, 1, 25}];
bzeros[l_, n_] := besselzero[[l + 1, n]];

for speed purposes.

For the purposes of this question consider $m=3$ (but there's essentially the same problem for any $m$):

 Cos[3 ArcTan[y/x] ] BesselJ[3, bzeros[3,1]Sqrt[x^2+y^2]],{x,-1,1},{y,-1,1},

Mathematica graphics

This is nearly right, except that the sign on the left-hand side ($y<0$) of the plot is wrong (it should not be mirror symmetric around $y=0$, rather continue alternating positive and negative after each $\Delta\phi=\pi/3$).

Looking at the Plot3D of ArcTan[y/x] I realized the cause of the problem was because of how ArcTan is defined, so I tried a workaround:

 Cos[3 Piecewise[{{0, x > 0 && y == 0}, {π/2, 
       x == 0 && y > 0}, {(3 π)/2, x == 0 && y < 0}, {(2 π)/2,
        x < 0 && y == 0}, {ArcTan[y/x], 
       x > 0 && y > 0}, {ArcTan[y/x] + π, 
       x < 0 && y > 0}, {ArcTan[y/x] + π, 
       x < 0 && y < 0}, {ArcTan[y/x] + 2 π, 
       x > 0 && y < 0}}] ] BesselJ[3, 
   bzeros[3, 1] Sqrt[x^2 + y^2]], {x, -1, 1}, {y, -1, 1},
 Epilog -> {Thick, Circle[]}, PlotRange -> {-1, 1},
 RegionFunction -> Function[{x, y, z}, Norm[{x, y}] < 1.]]

Mathematica graphics

This is nearly right, except for those white regions. Is there a better way to plot functions of this sort?


Use the two-argument form ArcTan[x, y] instead of ArcTan[y/x]:

ArcTan[x,y] gives the arc tangent of y/x, taking into account which quadrant the point $(x,y)$ is in.

It just works:

enter image description here

  • 1
    $\begingroup$ Yeesh. It's that simple. $\endgroup$
    – Eli Lansey
    Mar 29 '12 at 19:40

Quickly, now, which values are positive and which are negative?

Even when a legend may be present, it can help to have additional visual cues concerning the values in a density plot. ReliefPlot will do this with a "hillshading" algorithm. It also reverses the default colors, because it is more natural to think of light colors as being high or closer to the viewer.

To overcome the technical computing problems, isolate the calculations of the arctangent and the Bessel function:

ClearAll[arctan, besselj];
arctan[0., 0.] = 0;
arctan[x_, y_] := ArcTan[x, y];
besselj[m_, \[Chi]_, x_] /; x > 1 := 0;
besselj[m_, \[Chi]_, x_] := BesselJ[m, \[Chi] x];

This takes care of the special case at the origin and expedites the calculations beyond the unit circle (which are not wanted).

Because ReliefPlot does not appear to have the analog of RegionFunction, a reasonable workaround is to mask the undesired areas:

mask = RegionPlot[x^2 + y^2 > 1, {x, -1, 1}, {y, -1, 1}, 
   PlotStyle -> White, BoundaryStyle -> None];

Finally, we need to remember that ReliefPlot builds the image from the bottom to the top. To generate its array, then, we need to go in reverse order of the first coordinate. Whence the work is done as

Block[{m = 3, n = 1, \[Chi], data},
 \[Chi] = BesselJZero[m, n];
 data = ParallelTable[ 
   besselj[m, \[Chi], Sqrt[x^2 + y^2]] Cos[m arctan[x, y]], 
   {x, 1, -1, -.01}, {y, -1, 1, .01}];
  ReliefPlot[data, DataRange -> {{-1, 1}, {-1, 1}}, 
   Axes -> {True, True}, AxesOrigin -> {-1.05, -1.05}, Frame -> False],


The improvements are subtle but helpful: it is now visually clear where the peaks and valleys are located. If the effect is not strong enough, exaggerate it by lowering the second LightingAngle value (the help page shows how).

  • $\begingroup$ I typically include a legend for the plot, but this is a very nice approach. I also like your approach to dealing with the plot range by defining the Bessel function differently. Other than dealing with the {0,0} coordinate, what is the benefit to using the arctan definition? $\endgroup$
    – Eli Lansey
    Mar 30 '12 at 2:21
  • $\begingroup$ Is there a functional (not array-based) version of ReliefPlot? $\endgroup$
    – Eli Lansey
    Mar 30 '12 at 11:26
  • 1
    $\begingroup$ @Eli: (1) Because ReliefPlot insists on working with an array of machine reals, even a single Indeterminate result throws it off. You are right; the purpose of separately defining arctan was to cope with that problem at $\{0,0\}$. (2) Yes: make a 3D (surface) plot and set the viewpoint from (essentially) infinitely far above the surface. This might actually be a more efficient solution (but it's a little fussier in terms of needing to set many more graphing parameters: get rid of the mesh, set the viewpoint and lighting, etc.). $\endgroup$
    – whuber
    Mar 30 '12 at 13:54

As to the white lines:

Those are caused by small gaps in the plot introduced by the Exclusions option, which is trying to avoid jumps caused by discontinuities in your function. (In this case it's because your function uses Piecewise)

You can turn them off if necessary by using Exclusions -> None:

 Cos[3 Piecewise[{
          {0, x > 0 && y == 0}, 
          {Pi/2, x == 0 && y > 0}, 
          {(3 Pi)/2, x == 0 && y < 0}, 
          {(2 Pi)/2, x < 0 && y == 0},        
          {ArcTan[y/x], x > 0 && y > 0}, 
          {ArcTan[y/x] + Pi, x < 0 && y > 0}, 
          {ArcTan[y/x] + Pi, x < 0 && y < 0}, 
          {ArcTan[y/x] + 2 Pi, x > 0 && y < 0}
         }]] BesselJ[3, bzeros[3, 1] Sqrt[x^2 + y^2]], 
 {x, -1, 1}, {y, -1, 1}, 
 Epilog -> {Thick, Circle[]}, PlotRange -> {-1, 1}, 
 RegionFunction -> Function[{x, y, z}, Norm[{x, y}] < 1.], 
 Exclusions -> None]

enter image description here


The difficulties with either single-argument Arctan[y/x] or two-argument ArcTan[x,y] (the latter being undefined at the origin, as was addressed by @whuber) were also discussed on MathGroup. Heike's suggestion there was to use Arg[x + I y] to get the polar angle. But in fact it turns out to be slightly more efficient (and perhaps more readable) to use the function CoordinatesFromCartesian in the VectorAnalysis package:



 Apply[Cos[3 #2] BesselJ[3, bzeros[3, 1] #1] &, 
  CoordinatesFromCartesian[{x, y, 0}]], {x, -1, 1}, {y, -1, 1}, 
 Epilog -> {Thick, Circle[]}, PlotRange -> {-1, 1}, 
 RegionFunction -> Function[{x, y}, Norm[{x, y}] < 1.]]

(there was an unnecessary z argument in the RegionFunction of your example)

One final remark, related to your setup with BesselJZeros: instead of preparing them in a table, you may want to use the fact that BesselJZeros threads over lists! That makes the initial calculation faster than with a table. Of course you only do that loop once in this application, but I thought it's worth pointing out for the benefit of more complex applications.

  • $\begingroup$ I'm not sure that preparing the zeros in a table is any slower, since I only evaluate numerically at the end. $\endgroup$
    – Eli Lansey
    Mar 30 '12 at 13:06
  • $\begingroup$ @Eli That's probably true. I just wanted to mention it because I've found it to be very useful for my own purposes... Similarly, when doing plots involving large sums of BesselJ I found it faster to use the listability of BesselJ by going to ListDensityPlot instead of DensityPlot. But those things don't matter much in your example. $\endgroup$
    – Jens
    Mar 30 '12 at 15:31

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.