2
$\begingroup$

I have just got a very good answer to my other question, related to use of NIntegrate and tried to apply it to another situation:

t = 0.2; w = 1; a = 4/w^2; K = k*Sin[t]; k = 2 Pi;

XX[lf_, mf_, l_, mg_, p_, r_] := Abs[YY[lf, mf, l, mg, p, r]]^2;
YY[lf_, mf_, l_, mg_, p_, r_] := WignerD[{lf, mf, l}, t] BesselJ[mf - mg, K r] Sqrt[K/(2 Pi)];

MLG[lf_, l_, mg_, p_, r_] := Sum[XX[lf, ii, l, mg, p, r], {ii, -lf, lf, 1}];
Res1[r_?NumericQ] := (MLG[2, 1, 3, 0, r] - 
 MLG[2, -1, 1, 0, r])/(MLG[2, 1, 3, 0, r] + MLG[2, -1, 1, 0, r]);

NIntegrate[r Res1[r], {r, 0, 1000}, MaxPoints -> 500000, PrecisionGoal -> 10, MaxRecursion -> 50, WorkingPrecision -> 10]

and what I get is

enter image description here

A little bit about function Res1:

  1. the integral over $r^2$ (cylindrical coord. axial vector) on the domain $\{0, \infty \}$ must be zero

This is my ultimate goal. So ideally I do not integrate from 0 to 1000 as it is in the code, I need the integral from 0 to $\infty$ and obtain something close to zero.

  1. the "Bad boy" - Bessel obeys the relation:

$\lim_{R\rightarrow \infty} \int_0^R rdr J^2(ar) = \frac{R}{\pi a}$

so there might be a way to tame the oscillatory behavior within Res1.

The function looks like this

enter image description here

I have tried to use a brute force - Gaussian quadrature, but it seems like depending on the value of the upper integration limit I get very different answers within the range $\pm 0.5$. Is there a trick to make this evaluation work correctly?

$\endgroup$
  • $\begingroup$ @Ms Tais: Executing your code, I obtain NIntegrate::inumr: The integrand r Res1[r] has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,1000}}. $\endgroup$ – user64494 Dec 2 '16 at 18:39
  • $\begingroup$ @user64494 Yeah, my bad, I fixed it. Sorry. Now it does not work in a proper way. $\endgroup$ – MsTais Dec 2 '16 at 18:43
  • $\begingroup$ @Anton Antonov Well, it is possible, but I have just played around with those using Oscilatory method for $r\in \{0, 1000\}$, and the best I get is $~0.07...$ and the message that the integral failed to converge and the best it can do for me is this answer+error estimates. $\endgroup$ – MsTais Dec 2 '16 at 20:34
3
$\begingroup$

Insight can be obtained by plotting r Res1[r] (with t = 1/5 used for increased precision).

Plot[r Res1[r], {r, 0, 20}, AspectRatio -> .5/GoldenRatio, ImageSize -> Large]

enter image description here

After an initial transient, the curve appears to be almost periodic. To investigate this further, consider oscillations at r == 1000 and r == 1000000. The periods there are obtained by

r1 = r /. FindRoot[r Res1[r], {r, 996}, WorkingPrecision -> 30]
(* 996.000670616390157849961696210 *)
r2 = r /. FindRoot[r Res1[r], {r, 998.5}, WorkingPrecision -> 30]
(* 998.517418096292514726985805578 *)
r3 = r /. FindRoot[r Res1[r], {r, 999995.2}, WorkingPrecision -> 30]
(* 999995.189014480314603208960429 *)
r4 = r /. FindRoot[r Res1[r], {r, 999997.7}, WorkingPrecision -> 30]
(* 999997.705759254153466573546910 *)
r2 - r1
(* 2.51674747990235687702410937 *)
r4 - r3
(* 2.51674477383886336458648 *)

So the periods are nearly identical. Moreover, so are the curves themselves, as can be seen from.

Plot[{r  Res1[r] /. r -> s + r1, r  Res1[r] /. r -> s + r3}, {s, 0, 10}, 
    PlotStyle -> {Blue, Red}]

enter image description here

However, they are not precisely identical.

NIntegrate[r Res1[r], {r, r1, r2}, WorkingPrecision -> 30]
(* 0.0000819173063738409068656216604082 *)
NIntegrate[r Res1[r], {r, r3, r4}, WorkingPrecision -> 30]
(* 8.16930968934193286777117962325*10^-8 *)

The integral over a single oscillation is decreasing at just slightly faster than 1/r, based on these two values of r. If this relationship holds as r goes to Infinity, the total integral converges. But, if the integral over a single oscillation eventually decreases as 1/r or slower, the total integral does not converge. It is evident from this simple analysis that integration from 0 to some very large number in order to show that the integral to infinity is zero is impractical.

Addendum

Res1 can be simplified substantially by

lb = BesselJ[n, 2 \[Pi] r Sin[t]] /. n -> Range[0, 5];
Simplify[(MLG[2, 1, 3, 0, r] - MLG[2, -1, 1, 0, r])/
    (MLG[2, 1, 3, 0, r] + MLG[2, -1, 1, 0, r]), r > 0];
resimp = Collect[Numerator[%], lb, FullSimplify]/Collect[Denominator[%], lb, FullSimplify]
(* (4 BesselJ[0, 2 π r Sin[1/5]]^2 (-1 + Cos[3/5]) Sin[1/10]^2 - 
    4 BesselJ[4, 2 π r Sin[1/5]]^2 (-1 + Cos[3/5]) Sin[1/10]^2 + 
      BesselJ[3, 2 π r Sin[1/5]]^2 (2 + 3 Cos[1/5] + 6 Cos[2/5] + 5 Cos[3/5]) Sin[1/10]^2 - 
    2 BesselJ[1, 2 π r Sin[1/5]]^2 (3 - 4 Cos[1/5] + 3 Cos[2/5]) Sin[1/5]^2 + 
    8 BesselJ[5, 2 π r Sin[1/5]]^2 Sin[1/10]^4 Sin[1/5]^2)/
  (16 BesselJ[2, 2 π r Sin[1/5]]^2 Cos[1/10]^4 (1 - 2 Cos[1/5])^2 - 
    2 BesselJ[1, 2 π r Sin[1/5]]^2 (-2 + Cos[2/5] + Cos[4/5]) + 
      BesselJ[3, 2 π r Sin[1/5]]^2 (9 + 4 Cos[1/5] + 7 Cos[2/5]) Sin[1/5]^2 + 
      8 BesselJ[5, 2 π r Sin[1/5]]^2 Sin[1/10]^4 Sin[1/5]^2 + 
      8 BesselJ[0, 2 π r Sin[1/5]]^2 Sin[1/10]^2 Sin[3/10]^2 + 
      8 BesselJ[4, 2 π r Sin[1/5]]^2 Sin[1/10]^2 Sin[3/10]^2) *)

I believe that further simplification is possible by combining Bessel functions, but FullSimplify apparently is unable to do so.

$\endgroup$
  • $\begingroup$ You made a serious progress. Have you tried to find the asymptotics of r*resimp as r approches infinity? $\endgroup$ – user64494 Dec 3 '16 at 21:19
  • $\begingroup$ I did that. The result (see here dropbox.com/s/ivci16zs828d93a/Asymptotic.pdf?dl=0) suggests the integral of r*Res1[r] over the interval $[0,\infty)$ diverges. $\endgroup$ – user64494 Dec 4 '16 at 11:05
  • $\begingroup$ Sorry, I did not succeed to insert @bbgodfrey in the beginnings of my comments. $\endgroup$ – user64494 Dec 4 '16 at 18:14
  • $\begingroup$ @bbgodfrey as far as I got you, the convergence is questionable, correct? Probably, this is a good question for Mathematics forum then. Thank you! $\endgroup$ – MsTais Dec 4 '16 at 21:29
  • $\begingroup$ @MsTais I am reasonably confident that the integral does not converge. And, as you suggest, this now is becoming more a mathematics than a Mathematica problem, which could be posed at Mathematics. I do have a bit more material to add to my answer, which may be helpful. $\endgroup$ – bbgodfrey Dec 4 '16 at 23:17
3
$\begingroup$

If we inspect the function Res1(r), we can see that is the form $(x-y)/(x+y)$, so we can write it as $2x/(x+y)-1$, and in the line of the above post:

AbsoluteTiming[
Total@(NIntegrate[(2 r MLG[2, 1, 3, 0, 
       r]/(MLG[2, 1, 3, 0, r] + MLG[2, -1, 1, 0, r])) - 
   r, {r, #, # + 5}, 
  Method -> {"OscillatorySelection", 
    "TermwiseOscillatory" -> True, 
    Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 10000}}, 
  MaxRecursion -> 100, PrecisionGoal -> 8, AccuracyGoal -> 5] & /@
 Table[j, {j, 0, 1000, 5}])]

 {3.19624, 0.0756571}

From 0 to 10000:

{31.1375, 0.152549}

and from 0 to 100000:

{306.902, 0.266373}

I have realized that Res1 has a numerator that is a sum of oscillatory terms, so I say NIntegrate that separate each term with Method -> {"OscillatorySelection","TermwiseOscillatory" -> True. I have used a PC Desktop Intel i7-3770 QuadCore with Windows 7 (64 bits). I cannot assure that my code is faster, but it seems to be.

On the other hand, your integral seems to increase and not be close to zero, as $r\rightarrow\infty$. That is probably because your integrand is not so symmetrical in its oscillations:

enter image description here

Please check.

$\endgroup$
  • $\begingroup$ This is strange, that you don't get the error message. In my case I still get "The integral failed to converge after 9 recursive bisections...", which, I suppose, in this case means that the result is just garbage. $\endgroup$ – MsTais Dec 2 '16 at 20:51
  • $\begingroup$ @MsTais, please, take a look now to see if useful... $\endgroup$ – José Antonio Díaz Navas Dec 2 '16 at 23:21
  • $\begingroup$ yeah, it seems like all of us get basically the same result: it does not converge. Thank you! $\endgroup$ – MsTais Dec 5 '16 at 16:13
0
$\begingroup$

Somewhat different approach is

Sum[NIntegrate[r *Res1[r], {r, j, j + 50}, {Method -> "GlobalAdaptive"},   WorkingPrecision -> 12, AccuracyGoal -> 5], {j, 0, 950, 50}] // Timing

{150.906, 0.076008561771}

To verify it

Sum[NIntegrate[r *Res1[r] + 0.01, {r, j, j + 50}, {Method -> "GlobalAdaptive"}, 
WorkingPrecision -> 12, AccuracyGoal -> 5], {j, 0, 950, 50}] -   10 // Timing

{144.813, 0.0760084559}

$\endgroup$

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.