I wish to plot a function defined by

F[rx_, ry_, tau_, a_, sx_, sy_,  uf_] := 
   NIntegrate[E^(-qx^2 - a qy^2 - 2 qx (sx - I/2 rx) - 2 qy (a sy - I/2 ry) - I uf tau Sqrt[qx^2 + qy^2]), 
      {qx, -Infinity, +Infinity}, {qy, -Infinity, + Infinity}];

Psi[rx_, ry_, tau_, a_, sx_, sy_,  uf_] := 
   Sqrt[a]/(2 \[Pi]^3) E^(-2 (sx^2 + a sy^2)) Abs[F[rx, ry, tau, a, sx, sy, uf]]^2;

DensityPlot[Psi[rx, ry, 1, 2, 0, 3, 7], {rx, -10, 10}, {ry, -10, 10}]

Values and ranges for all variables are mock, but reasonably close in magnitude to the real values I am interested in.

My problem is, plotting takes a really long time. Is there any way I could possibly speed it up?

I am not overly familiar with the parallelizing of functions. I tried a simple

Parallelize[DensityPlot[Psi[rx, ry, 1, 2, 0, 3, 7], {rx, -10, 10}, {ry, -10, 10}]]

which did not have any effect, with Mathematica telling me it cannot be parallelized, and proceding with sequential evaluation.

Is there a smarter way to do it?


I am now going with the suggestion by @Michael E2. It does work in producing plots quickly, but these are very low quality. In fact, even at MaxRecursion -> 10 the plot is so boxy one can hardly read it. I should say I am interested in the actual Psi2, rather than Log(Psi2). In the example below my problem should be evident.

On the plus side, this image took a very reasonable 70 seconds on MaxRecursion -> 10. So my question now is, is there a way I can improve the quality?


Increasing PlotPoints did do the trick. For reference, this is what I am getting for MaxRecursion -> 8, PlotPoints -> 40in just shy of 4 minutes. The upper bit could use some more precision, I'll decide weighing looks against performance.

Thanks to everyone!

enter image description here

  • 1
    $\begingroup$ Have you tried using compiled function, such as cf = Compile[{{qx, _Real}, {qy, _Real}, {rx, _Real}, {ry, _Real}, {tau, _Real}, {a, _Real}, {sx, _Real}, {sy, _Real}, {uf, _Real}}, E^(-qx^2 - a qy^2 - 2 qx (sx - I/2 rx) - 2 qy (a sy - I/2 ry) - I uf tau Sqrt[qx^2 + qy^2])] for the integrand? $\endgroup$ – anderstood Feb 22 '18 at 15:24
  • 4
    $\begingroup$ DensityPlot is the least of your problems. The real pain is the integration. Simply evaluating Psi[-10, -10, 1, 2, 0, 3, 7] takes about 17 seconds here. Assuming your plot samples the region with 200x200 points, then you need about 8 days! Even with 8 cores in parallel you'll only bring it down to 24h. You need to think hard about this NIntegrate. Can you approximate this integral? Can you use a simpler method? Do you have to integrate to infinity for a vastly dropping function such as Exp? $\endgroup$ – halirutan Feb 22 '18 at 15:55
  • $\begingroup$ @halirutan that's disheartening. The whole reason why I was trying a NIntegrate was to find a result for a function I have no good analytical reason / way to approximate further. I will need to look more into that, then. Thank you for your help. $\endgroup$ – Cello Feb 22 '18 at 16:17
  • 2
    $\begingroup$ Have you tried increasing PlotPoints? Start from low values like 10 and increase slowly. $\endgroup$ – anderstood Feb 23 '18 at 13:51

If you change to polar coordinates and simplify the coefficients, Integrate will calculate the radial integral symbolically. This reduces the numerical integration to one dimension, over the angle, which can efficiently be calculated with the trapezoidal method.

coeffsubs = Thread[{c, b, a} -> CoefficientList[
     (-qx^2 - a qy^2 - 2 qx (sx - I/2 rx) - 2 qy (a sy - I/2 ry) - 
         I uf tau Sqrt[qx^2 + qy^2]) /.
       {qx -> r Cos[t], qy -> r Sin[t]} // Simplify[#, r > 0] &, (* to polar *)
tIntegrand[rx_, ry_, tau_, a_, sx_, sy_, uf_] = 
 Integrate[(E^(c + r  b + r^2 a)) r,
   {r, 0, Infinity},
   Assumptions -> {a < 0}] /. coeffsubs

foo = 0;
F2[rx_?NumericQ, ry_, tau_, a_, sx_, sy_, uf_] := 
  NIntegrate[tIntegrand[rx, ry, tau, a, sx, sy, uf], {t, 0, 2 Pi}, 
   PrecisionGoal -> 6, Method -> {"Trapezoidal"}, 
   EvaluationMonitor :> ++foo];

Psi2[rx_, ry_, tau_, a_, sx_, sy_, uf_] := 
  Sqrt[a]/(2 π^3) E^(-2 (sx^2 + a sy^2)) Abs[F2[rx, ry, tau, a, sx, sy, uf]]^2;

For the density plot, I'd suggest doing a log plot because of the range of Psi:

foo = 0;  (* counts integrand evaluations *)
murf = 0; (* counts Psi2 evaluations *)
PrintTemporary@Dynamic@{foo, murf, Clock[{0, Infinity}]}; (* progress indicators *)
DensityPlot[Log10@Psi2[rx, ry, 1, 2, 0, 3, 7],
  {rx, -10, 10}, {ry, -10, 10},
  PlotLegends -> Automatic,
  MaxRecursion -> 1,  (* save time; increase for greater quality *)
  EvaluationMonitor :> ++murf] // AbsoluteTiming
{foo, murf}

Mathematica graphics

(*  {69278, 348}  *)

Update on accuracy:

I compared the OP's Psi with Psi2 at a few points (rx, ry = ±10, 0). The difference between them seems to grow as ry increases toward 10. The relative difference between them at rx = 10, ry = 10 is


Since no error messages were generated, I decide to check whether the difference was simply due to rounding error. I recalculated the OP's Psi at the same point with the options PrecisionGoal -> 6, WorkingPrecision -> 32 added to NIntegrate in F. (It took a couple of hours and almost 3 GB of memory.) The relative difference dropped to


That is much more in line with the expected truncation error from using NIntegrate with either the trapezoidal or the two-dimensional integrations, as well as with the errors for the lower values ry = -10, 0. My conclusion is that the numerical conditioning of Psi gets worse when ry grows large (around 10 and higher, more or less), and that Psi2 is relatively stable and accurate.

  • $\begingroup$ This looks promising! Thank you, I'll take a good look at it. $\endgroup$ – Cello Feb 22 '18 at 18:05

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.