11
$\begingroup$

I'm plotting an ellipse 1/(1 - e Cos[x]) with e=0.9 using PolarPlot. Taking $x\in[0,2\pi)$, i.e. one revolution, everything looks fine. But taking n=500 revolutions produces something very vaguely resembling an ellipse:

e = 0.9;
n = 500;
PolarPlot[1/(1 - e Cos[x]), {x, 0, 2 n π}, Frame -> True, Axes -> False, ImageSize -> 700, FrameTicks -> False]

enter image description here

The solution is quite straightforward, and uses PlotPoints:

PolarPlot[1/(1 - e Cos[x]), {x, 0, 2 n π}, Frame -> True, 
Axes -> False, ImageSize -> 700, FrameTicks -> False, 
PlotPoints -> 10^5]

enter image description here

Without PlotPoints generating the image took 0.05 sec, with PlotPoints -> 10^5 it needed almost 4 sec - rather understandable as it has more to plot.

Finally, consider a different approach: generate 500 images for each revolution separately, without PlotPoints, and Show them together:

Show@Table[PolarPlot[1/(1 - e Cos[x]), {x, 2 (i - 1) π, 2 i π}, 
Frame -> True, Axes -> False, ImageSize -> 700, 
FrameTicks -> False], {i, 1, 500}]

enter image description here

That looks neat, but takes an enourmous 9 sec.

The thing is that I'm going to plot something like this hundreds of times, so PlotPoints will need a few minutes to generate a set of images like this. My question is: can this be sped up while keeping the plot looking nice? The third approach, with Show, is meant to be a demonstration that in fact just PolarPlot works fine for each revolutions, but fails in a longer run. Of course plotting this takes such a long time because it creates, stores and stacks 1000 images.

NOTES: Taking 10^4 instead of 10^5 won't do because (i) it only improves by a factor of two, (ii) in fact I'm going to have a slightly wiggly ellipse (perturbations in celestial mechanics) so that won't give enough detail. This is also the reason why I don't think that sampling the ellipse and using a ListPlot will work as (i) if the sampling will be too sparse, the details will be lost, and (ii) using a large enough number of points will take time to generate, store and plot them.

EDIT: Using MaxRecursion:

enter image description here

and the timings for generation of each case:

enter image description here

It looks like MaxRecursion=11 or 12 is sufficient, and speeds up the thing by a factor of between 4 and 8 - quite neat.

Still, is there a way to speed it up even more? It's not particularly crucial, but now I'm just curious how fast this can be done.

$\endgroup$
  • $\begingroup$ It's generally discouraged to edit your question with posted answers. If you want to encourage more answers that's fine of course. meta.mathematica.stackexchange.com/questions/1918/… That said, you did add a lot yourself, so I don't know. $\endgroup$ – Feyre Sep 4 '16 at 9:41
  • $\begingroup$ @Feyre Your link is about incorporating someone's answer into the question. I made an edit because before it was just incorrect. I don't see any reason why a mathematically incorrect statement couldn't be corrected. I wouldn't mind if you edited your answer too. That would make this thread clear for any potential future reader. $\endgroup$ – corey979 Sep 4 '16 at 10:00
  • $\begingroup$ I meant the thing about MaxRecursion. I at least wouldn't do it if you do get a more definitive answer. $\endgroup$ – Feyre Sep 4 '16 at 10:09
  • $\begingroup$ @Feyre Ah, sorry, I misunderstood you. Yeah, I indeed took the idea from your answer to estimate how can this improve my issue, and at this point I agree with you comment. But I explicitly write "Still, is there a way to speed it up even more?" so indeed I hope that someone maight present another idea here. I apologize if something was unclear. $\endgroup$ – corey979 Sep 4 '16 at 10:17
5
$\begingroup$

Note that with n=1000 you have 500 revolutions. MaxRecursion is the go to here, I tend to use PlotPoints only when there are tiny features like narrow peaks which are missing.

AbsoluteTiming[
 PolarPlot[1/(1 - e Cos[x]), {x, 0, n π}, Frame -> True, 
   Axes -> False, ImageSize -> 700, FrameTicks -> False, 
   MaxRecursion -> 12];]
AbsoluteTiming[
 PolarPlot[1/(1 - e Cos[x]), {x, 0, n π}, Frame -> True, 
   Axes -> False, ImageSize -> 700, FrameTicks -> False, 
   PlotPoints -> 10^5];]

{0.975767, Null} {3.65704, Null}

$\endgroup$
  • $\begingroup$ Yeah, I messed up with the n. I edited my post so that it's consistent now. $\endgroup$ – corey979 Sep 4 '16 at 9:14
9
$\begingroup$

It's a listable expression, so you can plot it with a high precision fixed step size very quickly:

e = 0.9;
n = 500;
Graphics@Line@With[{x = Range[0, 2 n Pi, 0.01]},
    Transpose[{Cos[x]/(1 - e Cos[x]), Sin[x]/(1 - e Cos[x])}]
    ] // AbsoluteTiming

Mathematica graphics

With a fixed step size the number of points required is determined by the part of the curve that requires the most careful sampling. For this curve the curvature is rather uniform and no part of the curve requires a lot more careful sampling than any other part of the curve. So for this example, the number of points required with a fixed step size will be comparable to the number of points that PolarPlot uses. If this is also true for the curve you are actually trying to draw, then this approach may also be preferable in that case.

$\endgroup$
4
$\begingroup$

The following tables show that the major part of the time to plot is proportional to the number of function evaluations, which for the OP's function is about a 10 microseconds per evaluation, or about 100 times as slow as C.E.'s method. Recursive subdivision adds a small overhead, which is not tremendously significant. However, comparing the extremes,

MaxRecursion -> 12, PlotPoints -> 50      (* 50 is the default setting *)
MaxRecursion -> 1,  PlotPoints -> 100000

we see that for roughly the same number of sample points, the extra recursion adds almost 30% to the execution time over the high number for PlotPoints.

Mathematica graphics

Mathematica graphics

The "error" was from an estimate of the maximum distance between the line segments of the plot and the actual ellipse. An error of 0.1% represents one pixel in an image of width 1000, which should give a pretty good plot. (It was also chosen because of the OP's remark that MaxRecursion of 11 or 12 for the default number of PlotPoints seemed sufficient.)

The general rule is that when you need denser sampling everywhere raise PlotPoints. When denser sampling is needed in a limited part of the domain, where the bending of the curve is not being rendered accurately, then raise MaxRecursion. Note that for raising MaxRecursion to be effective, the PlotPoints sampling has to be dense enough to detect the bending that needs further refinement. (Worst case example: Sampling Sin[x] every 2 Pi results in a straight-line plot with no bending and no apparent need for further refinement.) How much bending is tolerated is controlled by Method -> {Refinement -> {ControlValue -> (*radian angle*) }} or by Method -> {"MaxBend" -> 1} (now apparently deprecated). Note that the recursive refinement algorithm is not publicly known, as far as I know. It generally examines the angle change at each segment and subdivides if it is too great, but that is not the whole story. Subdivision occurs, shown in the tables above, each time MaxRecursion is increased. Careful inspection shows that it sometimes occurs in places where it is unneeded, well after having satisfied the bending tolerance.

Code dump

Note the execution times in the tables above. Expect to take the sum of them to reproduce the table, depending on machine speed.

(* Collect data on a plot: {time, npts, error, plot}
   Memoized so that one can play with the results; takes < 0.5 GB *)
Clear[pdata];
pdata[pp_, mr_] := pdata[pp, mr] = 
   Module[{e = 0.9, n = 500, plot, time, pts, npts = 0, radii, angles,
      error},
    {time, plot} = 
     AbsoluteTiming@
      PolarPlot[1/(1 - e Cos[x]), {x, 0, 2*n* π}, 
       PlotStyle -> Directive[Black, Thin], Frame -> True, 
       Axes -> False, FrameTicks -> False, MaxRecursion -> mr, 
       PlotPoints -> pp];
    {pts} = Cases[plot, Line[p_] :> p, Infinity];
    npts = Length@pts;
    radii = Norm /@ pts;
    angles = ArcCos[(radii - 1)/(e radii)];
    error = 
     Max@Abs[1/(1 - e Cos[Mean /@ Partition[angles, 2, 1]]) - 
        Norm /@ Mean /@ Partition[pts, 2, 1]];
    {time, npts, error, plot}];
{$time, $npts, $error, $plot} = Range[4]; (* parts of pdata[] *)

mypts = Round[50*((10^4/5)^(1/5))^Range[0, 5]];
mymr = Range[0, 12]~Append~15;
Labeled[
 TableForm[
  Table[
   Style[
    pdata[pp, mr][[$time]],
    ColorData["RedGreenSplit"][
     UnitStep[0.01 - pdata[pp, mr][[$error]]]]
    ],
   {mr, mymr},
   {pp, mypts}
   ],
  TableHeadings -> {mymr, mypts}],
 {"MaxRecursion", "PlotPoints", 
  Row[{"Evaluation time, by options settings, colored by error ( ", 
    ColorData["RedGreenSplit"][0], " > 0.1%, ", 
    ColorData["RedGreenSplit"][1], " ≤ 0.1%)"}, 
   BaseStyle -> {Medium}]},
 {Left, Top, Bottom}, RotateLabel -> True,
 BaseStyle -> {"Label", 10}, LabelStyle -> "Label"
 ]

Labeled[
 TableForm[
  Table[
   Style[
    pdata[pp, mr][[$npts]],
    ColorData["RedGreenSplit"][
     UnitStep[0.01 - pdata[pp, mr][[$error]]]]
    ],
   {mr, mymr},
   {pp, mypts}
   ],
  TableHeadings -> {mymr, mypts}],
 {"MaxRecursion", "PlotPoints", 
  Row[{"Evaluations, by options settings, colored by error ( ", 
    ColorData["RedGreenSplit"][0], " > 0.1%, ", 
    ColorData["RedGreenSplit"][1], " ≤ 0.1%)"}, 
   BaseStyle -> {Medium}]},
 {Left, Top, Bottom}, RotateLabel -> True,
 BaseStyle -> {"Label", 10}, LabelStyle -> "Label"
 ]
$\endgroup$
  • $\begingroup$ So, the compact message is that taking a careful combination of MaxRecursion and PlotPoints will do the job best as possible? Btw, +1 for the great tables! $\endgroup$ – corey979 Sep 4 '16 at 21:28
  • $\begingroup$ @corey979 Yes, and the best settings depend on the function being plotted. And if, as in your ellipse, uniform, dense sampling does a good job, then C.E.'s method will be the best for speed. $\endgroup$ – Michael E2 Sep 4 '16 at 21:37

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.