# NIntegrate with BesselJ, slow convergence

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

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

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?

• @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}}. – user64494 Dec 2 '16 at 18:39
• @user64494 Yeah, my bad, I fixed it. Sorry. Now it does not work in a proper way. – MsTais Dec 2 '16 at 18:43
• @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. – MsTais Dec 2 '16 at 20:34

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]


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}]


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.

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.

• You made a serious progress. Have you tried to find the asymptotics of r*resimp as r approches infinity? – user64494 Dec 3 '16 at 21:19
• 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. – user64494 Dec 4 '16 at 11:05
• Sorry, I did not succeed to insert @bbgodfrey in the beginnings of my comments. – user64494 Dec 4 '16 at 18:14
• @bbgodfrey as far as I got you, the convergence is questionable, correct? Probably, this is a good question for Mathematics forum then. Thank you! – MsTais Dec 4 '16 at 21:29
• @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. – bbgodfrey Dec 4 '16 at 23:17

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:

• 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. – MsTais Dec 2 '16 at 20:51
• @MsTais, please, take a look now to see if useful... – José Antonio Díaz Navas Dec 2 '16 at 23:21
• yeah, it seems like all of us get basically the same result: it does not converge. Thank you! – MsTais Dec 5 '16 at 16:13

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}