3
$\begingroup$

I'm trying to plot 2d numeric integral of a complex function which is actually real. First problem - small non zero complex part, I was forced to plot RE of it. Second - I'm using only exact numbers but got General::precw warning and can't tune WorkingPrecision! The output interpolation is very rough and I have no clue how to deal with it. My code:

c[x_, y_] := (x*u + (1 - u)*y + 4*u*(1 - u)*k^2)^(1/2)
a = 1
W[r_, k_] := (2 Exp[-2*I*k*r])/(Pi^3 a^3)*
   NIntegrate[
   Exp[4*I*u*r*k]*
    u*(1 - u)*(3/4 Exp[-2*r*c[1, 1]]/c[1, 1]^5 + 
      3/2 (Exp[-2*r*c[1, 1]]*r)/c[1, 1]^4 + (Exp[-2*r*c[1, 1]]*r^2)/
      c[1, 1]^3), {u, 0, 1}, WorkingPrecision -> 20]
Plot3D[r^2*k^2*W[r, k], {r, 0, 5}, {k, 0, 5}, PlotRange -> All, 
 AxesLabel -> Automatic]
$\endgroup$
2
  • $\begingroup$ You were too quick to accept my answer, recommend that you accept Michael E2's answer as it is much better. $\endgroup$
    – Bob Hanlon
    Apr 9, 2017 at 17:24
  • $\begingroup$ @BobHanlon yeah, I'm newbie in SE and thought changing accepted answer will be not fair $\endgroup$
    – satoru
    Apr 9, 2017 at 18:20

2 Answers 2

5
$\begingroup$

Some recommendations:

  • Since the integral is provably real, use Re on the integral instead of Chop, since Chop will chop small real parts, too.
  • To speed things up slightly, use Re on the integrand, since NIntegrate will compute with real numbers instead complex ones. Simplify it beforehand so that no complex results are generated in computing the integral.
  • To speed things up further, use a higher setting of Gauss points in the "GaussKronrodRule". The integrand is fairly smooth and not very oscillatory, and higher setting leads to a quickly converging integral.
  • Since it is a well-behaved integral, a high WorkingPrecision is unnecessary. Using MachinePrecision will speed things up even more.
  • The plot is ragged because of scaling: The variation in height ($\approx 10^{-3}$) is so small that the plot is virtually flat (in the standard unscaled Euclidean geometry) and Plot3D does no recursive subdivision, no matter what the setting of MaxRecursion.
  • To fixed the raggedness, one can increase the PlotPoints, which is somewhat wasteful, or scale the function, so that the automatic adaptive algorithm works effectively; then rescale the result.
  • For another speed improvement when plotting, lower the PrecisionGoal in NIntegrate. Since the integrand converges quickly, you might risk being bold and assume the actual error is much less than the error estimate NIntegrate calculates. In any case, three or four digits of precision is usually enough for plotting.
  • There's an error in the coding of W[] and c[]. In W[r, k], the value of k will not be passed to the body of c[1, 1,]. It has the appearance of working because Plot3D[] uses Block to temporarily set the global value of k. Thus, while Plot3D is being computed, the parameter k in W[r, k] and the symbol k in c[1, 1] happen to have the same value. Try evaluating W[2, 3] with the OP's code and you get an error.

OP's refactored code (thanks to @xzczd for suggesting "SymbolicProcessing" -> 0):

ClearAll[c, W];
c[x_, y_, k_] := (x*u + (1 - u)*y + 4*u*(1 - u)*k^2)^(1/2);   (* N.B. argument  k *)
a = 1;
$pg = Automatic;  (* symbol  $pg  for  PrecisionGoal  setting in  NIntegrate[] *)
Block[{u, k, r},  (* protect symbols while definition is constructed *)
 With[{integrand = 
    Simplify[(2 Exp[-2*I*k*r])/(Pi^3 a^3) Exp[4*I*u*r*k]*
        u*(1 - u)*(3/4 Exp[-2*r*c[1, 1, k]]/c[1, 1, k]^5 +     (* N.B. argument  k *)
          3/2 (Exp[-2*r*c[1, 1, k]]*r)/
            c[1, 1, k]^4 + (Exp[-2*r*c[1, 1, k]]*r^2)/c[1, 1, k]^3) // 
       Re // ComplexExpand,
     0 < k < 5 && 0 < r < 5 && 0 < u < 1]},
  W[r_?NumericQ, k_?NumericQ] := NIntegrate[
    integrand,
    {u, 0, 1},
    Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0, "Points" -> 11},
    PrecisionGoal -> $pg, WorkingPrecision -> MachinePrecision]
  ]
 ]

Timing comparison for different precision goals:

Block[{$pg = 3},   (* low setting for  PrecisionGoal  in  NIntegrate[] *)
   Table[W[r, k], {r, 0., 5., 0.1}, {k, 0., 5., 0.1}]]; // AbsoluteTiming
(*  {4.4409, Null}  *)

Block[{$pg = 8},    (* ~default setting for  PrecisionGoal  in  NIntegrate[] *)
   Table[W[r, k], {r, 0., 5., 0.1}, {k, 0., 5., 0.1}]]; // AbsoluteTiming
(*  {10.9349, Null}  *)

The unscaled plot. Raising MaxRecursion has no effect:

 Block[{$pg = 3},
 Plot3D[r^2*k^2*W[r, k] // Re, {r, 0, 5},
  {k, 0, 5},
  MaxRecursion -> 15,
  PlotRange -> All, AxesLabel -> Automatic]
 ]

The scaled plot. The default MaxRecursion is nearly perfect when the plot is scaled so that scale of each axis are approximately the same. Update: I tried ScalingFunctions earlier, but I must have done something wrong. It does work, but in this case, the ticks are not as nice as the manual scaling (see below).

Block[{$pg = 3},
 Plot3D[r^2*k^2*W[r, k] // Re, {r, 0, 5},
  {k, 0, 5},
  ScalingFunctions -> {None, None, {10^4 # &, 10^-4 # &}},
  PlotRange -> All, AxesLabel -> Automatic]
 ]

The previous scaled plot code produces nicer ticks:

 Block[{$pg = 3, scale = 1.*^4},
 Plot3D[scale*r^2*k^2*W[r, k] // Re, {r, 0, 5},
   {k, 0, 5},
   PlotRange -> All, AxesLabel -> Automatic] /. 
  gc_GraphicsComplex :> Scale[gc, {1., 1., 1/scale}, {0., 0., 0.}]
 ]

$\endgroup$
4
  • $\begingroup$ Very useful! Why the answer is only 1 post... $\endgroup$
    – satoru
    Apr 8, 2017 at 19:00
  • $\begingroup$ @satoru Thanks! (If I understand your why question, all I can say is that a good answer ought to address all the issues.) $\endgroup$
    – Michael E2
    Apr 8, 2017 at 22:28
  • 3
    $\begingroup$ "SymbolicProcessing" -> 0 will speed up the code even further :) $\endgroup$
    – xzczd
    Apr 9, 2017 at 16:06
  • $\begingroup$ @xzczd Sometimes it helps only a little and is not worth worrying about, but it's significant in this case. I should've checked it earlier. Thanks! :) $\endgroup$
    – Michael E2
    Apr 9, 2017 at 16:14
3
$\begingroup$

EDIT: based on comment by Michael E2 definition of c changed to make dependence on u and k explicit

Clear[c];
c[x_, y_, u_, k_] := (x*u + (1 - u)*y + 4*u*(1 - u)*k^2)^(1/2)    
    a = 1;

Since W uses numerical techniques, its arguments should be restricted with NumericQ. Use Chop to eliminate imaginary artifacts.

W[r_?NumericQ, 
  k_?NumericQ] :=
 (2 Exp[-2*I*k*r])/(Pi^3 a^3)*
   NIntegrate[Exp[4*I*u*r*k]*u*(1 - u)*
     (3/4 Exp[-2*r*c[1, 1, u, k]]/c[1, 1, u, k]^5 +
       3/2 (Exp[-2*r*c[1, 1, u, k]]*r)/c[1, 1, u, k]^4 +
       (Exp[-2*r*c[1, 1, u, k]]*r^2)/c[1, 1, u, k]^3), {u, 0, 1}, 
    WorkingPrecision -> 20] // Chop

You can improve the smoothness of the plot with PlotPoints and/or MaxRecursion although these will slow down the evaluation.

Plot3D[r^2*k^2*W[r, k], {r, 0, 5}, {k, 0, 5},
  PlotRange -> All, AxesLabel -> Automatic,
  WorkingPrecision -> 20, PlotPoints -> 30, MaxRecursion -> 5] // Quiet

enter image description here

$\endgroup$
3
  • $\begingroup$ Thank you! Helps a lot, now i will use this nice interpolation for my further computations. $\endgroup$
    – satoru
    Apr 8, 2017 at 17:30
  • $\begingroup$ You might want to try evaluating W on its own, say W[1, 1]. There's a coding error in the OP that you might want to fix in your answer. $\endgroup$
    – Michael E2
    Apr 9, 2017 at 16:29
  • $\begingroup$ @MichaelE2 - thanks, correction made. $\endgroup$
    – Bob Hanlon
    Apr 9, 2017 at 17:25

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.