4

I am new to Mathematica and this forum so apologies in advance if I am not doing anything properly (do let me know). I am trying to visualize the imaginary part of an integral of q w.r.t r when plotted against w in the code bellow. t1 and t2 are the limits of integration and are functions of l and w. The plot looks weird for w < 0.5 (blue patch on the left). What does it mean when this kind of behavior is in the plot? Is this due to the nature of my integrand or is there something I can do to fix it? Ideally for w < 0.5 it should be a curve/line instead of this dark patch. Any help is appreciated.

My code:

q[r_, l_, w_] := 
  Sqrt[(-l (1 + l) (-(2/r^3) + 1/r^2) + 1/r^4 + w^2)/(1 - 2/r)^2 - 1/(4(-2 + r)^2)];

t1[l_, w_] := 1/24 (2 Sqrt[3] Sqrt[(192 w^2 + (49 + (-117649 + 525888 w^2 + 24 Sqrt[3]Sqrt[-73530625 w^2 + 159891584 w^4 - 4096 w^6])^(1/3))^2)/(w^2 (-117649 + 525888 w^2 + 24 Sqrt[3] Sqrt[-73530625 w^2 + 159891584 w^4 - 4096 w^6])^(1/3))] + Sqrt[6] [Sqrt](1/w^2 (392 + (-4802-384 w^2)/(-117649 + 525888 w^2 +24 Sqrt[3]Sqrt[-73530625 w^2 + 159891584 w^4 - 4096 w^6])^(1/3) - 2 (-117649 + 525888 w^2 + 24 Sqrt[3]Sqrt[-73530625 w^2 + 159891584 w^4 - 4096 w^6])^(1/3) - (2304 Sqrt[3])/Sqrt[(192 w^2 + (49 + (-117649 + 525888 w^2 + 24 Sqrt[3]Sqrt[-73530625 w^2 + 159891584 w^4 - 4096 w^6])^(1/3))^2)/(w^2 (-117649 + 525888w^2 + 24 Sqrt[3]Sqrt[-73530625 w^2 + 159891584 w^4 - 4096 w^6])^(1/3))])));

t2[l_, w_] := 1/(4 Sqrt[3]) (Sqrt[(192 w^2 + (49 + (-117649 + 525888 w^2 + 24 Sqrt[3] Sqrt[-73530625 w^2 + 159891584 w^4 - 4096 w^6])^(1/3))^2)/(w^2 (-117649 + 525888 w^2 + 24 Sqrt[3] Sqrt[-73530625 w^2 + 159891584 w^4 - 4096 w^6])^(1/3))] - 1/Sqrt[2] ([Sqrt](1/w^2 (392 + (-4802 - 384 w^2)/(-117649 + 525888 w^2 + 24 Sqrt[3]Sqrt[-73530625 w^2 + 159891584 w^4 - 4096 w^6])^(1/3) - 2 (-117649 + 525888 w^2 + 24 Sqrt[3]Sqrt[-73530625 w^2 + 159891584 w^4 - 4096 w^6])^(1/3) - (2304 Sqrt[3])/Sqrt[(192 w^2 + (49 + (-117649 + 525888 w^2 + 24 Sqrt[3]Sqrt[-73530625 w^2 + 159891584 w^4 - 4096 w^6])^(1/3))^2)/(w^2 (-117649 + 525888 w^2 + 24 Sqrt[3]Sqrt[-73530625 w^2 + 159891584 w^4 - 4096 w^6])^(1/3))]))));

l = 2;

Plot[ γ = NIntegrate[q[r, l, w], {r, t2[l, w], t1[l, w]}]; Im[γ], {w, 0.01, 2}]

enter image description here

m_goldberg
  • 107,779
  • 16
  • 103
  • 257
guestuser
  • 71
  • 3
  • 1
    Neither the definition of t1 nor t2showslon its righthand side, so those functions are not functions oflandw, but only ofw. Good practice would be to removel` from the lefthand side of the definitions as well. – m_goldberg Dec 27 '20 at 04:27
  • I get an error from NIntegrate when I run your code, don't you? Related/duplicate: https://mathematica.stackexchange.com/questions/237154/is-there-a-reason-behind-strange-noise-appearing-for-some-coefficients-of-polyno/237162#237162 – Michael E2 Dec 27 '20 at 14:32

1 Answers1

1

Small round-off errors in t2[l, w] and t1[l, w] lead them to have small imaginary parts, randomly changing as w changes. NIntegrate will integrate between them along a line in the complex plane. In turn, the small, random imaginary parts are propagated through the evaluation of q[r, l, w] = Sqrt[...]. For certain values of r the argument of the square-root should be a negative real number, but it actually carries a small, random imaginary part — in particular, the imaginary part has a random sign (as w changes). Consequently, the argument of Sqrt[...] dances back and forth across the branch cut, causing the imaginary part of q[r, l, w] to discontinuously change signs for certain values of r as w changes.

My recommendation is to get rid of the imaginary parts with either of the following (first one seems better):

Plot[Im@NIntegrate[q[r, l, w], {r, Chop@t2[l, w], Chop@t1[l, w]}, 
   PrecisionGoal -> 4], {w, 0.01, 2}]

Plot[Im@NIntegrate[q[r, l, w], {r, Re@t2[l, w], Re@t1[l, w]}, MaxRecursion -> 20], {w, 0.01, 2}]

MaxRecursion -> 20 gets rid of an error message but yields an unimportant warning. Using PrecisionGoal -> 4 instead of MaxRecursion eliminates all messages and warnings*, is a reasonable lower goal for plotting, and speeds things up a little.

*Note: The error/warning message we get is the result of Plot testing its argument near the beginning of the plot domain (i.e. w near 0.01). After that, Plot routinely suppresses many error messages, which hides from the user problems that may exist. This issue came up recently here: Is there a reason behind strange noise appearing for some coefficients of polynomial in NIntegrate?

Another approach to deal with round-off error is to use arbitrary-precision numbers, which @ciao shows.

Michael E2
  • 235,386
  • 17
  • 334
  • 747