# NIntegrate reacts strangely to changes in WorkingPrecision

I'm trying to integrate a bunch of expressions similar to

integrand =
p^2 (-4 + ((3 - I Sqrt) (1 + 2 p^2) Cot[
Sqrt[-1 - I Sqrt - 2 p^2]/(2 Sqrt)])/(
3 Sqrt[-2 - 2 I Sqrt -
4 p^2]) + ((3 + I Sqrt) (1 + 2 p^2) Cot[
Sqrt[-1 + I Sqrt - 2 p^2]/(2 Sqrt)])/(
3 Sqrt[-2 + 2 I Sqrt - 4 p^2]) + Cot[1/2 Sqrt[-1 - p^2]]/
Sqrt[-1 - p^2] -
4 RootSum[
1 + p^2 + p^4 + 4 \[Pi]^2 #1^2 + 8 p^2 \[Pi]^2 #1^2 +
16 \[Pi]^4 #1^4 &, (
PolyGamma[0, -#1] #1)/(-1 + p^2 + 4 \[Pi]^2 #1^2) &]);

NIntegrate[integrand, {p, 0, \[Infinity]}]


The outcome of the evaluation (both if it succeeds at all and what result it returns if it does) depends strongly on the value of WorkingPrecision. I can see why if I plot the integrand with and without increased WorkingPrecision:

Plot[integrand, {p, 0, 1000}] Plot[integrand, {p, 0, 1000}, WorkingPrecision -> 30] This suggests that with high-enough WorkingPrecision everything should be fine. But for values 16 and higher I get the error SystemException["MemoryAllocationFailure"]. Between 7 and 16 I get warnings about Numerical integration converging too slowly and results spanning many orders of magnitude. Below WorkingPrecision -> 6 or not setting it at all throws the error "The integrand has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries".

I'm thoroughly lost at this point. Can I somehow avoid these problems? And if not, how do I know which result to trust?

• Are you sure the integral converges? Your integrand looks like it has a 1/x dependence. – Carl Woll Oct 30 '17 at 17:58
• @CarlWoll If the integral did not converge, wouldn't I be unable to attain a result regardless of the value of WorkingPrecision? But I do get a finite numerical result for WorkingPrecision $\in \{7,8,...,16\}$. – Casimir Oct 30 '17 at 18:10
• Look at NIntegrate[1/x, {x, 1, Infinity}] which produces a result with the same warnings, and is clearly not convergent. – Carl Woll Oct 30 '17 at 19:19

The memory exception issue is caused by Cot. You are basically feeding Cot an extremely large pure imaginary number. For example:

Cot[416*^14 I]


General::nomem: The current computation was aborted because there was insufficient memory available to complete the computation.

Throw::sysexc: Uncaught SystemException returned to top level. Can be caught with Catch[[Ellipsis], _SystemException].

SystemException["MemoryAllocationFailure"]

I think it's worth reporting this as a bug. You can work around this issue by replacing Cot with:

cot[p_?NumericQ] := Cos[p]/Sin[p]


so that:

integrand = p^2 (-4 + ((3 - I Sqrt) (1 + 2 p^2) cot[
Sqrt[-1 - I Sqrt - 2 p^2]/(2 Sqrt)])/(
3 Sqrt[-2 - 2 I Sqrt -
4 p^2]) + ((3 + I Sqrt) (1 + 2 p^2) cot[
Sqrt[-1 + I Sqrt - 2 p^2]/(2 Sqrt)])/(
3 Sqrt[-2 + 2 I Sqrt - 4 p^2]) + cot[1/2 Sqrt[-1 - p^2]]/
Sqrt[-1 - p^2] -
4 RootSum[
1 + p^2 + p^4 + 4 \[Pi]^2 #1^2 + 8 p^2 \[Pi]^2 #1^2 +
16 \[Pi]^4 #1^4 &, (
PolyGamma[0, -#1] #1)/(-1 + p^2 + 4 \[Pi]^2 #1^2) &]);


Then, we can evaluate the integral with various settings for the WorkingPrecision:

NIntegrate[integrand, {p, 0, Infinity}, WorkingPrecision->16]
NIntegrate[integrand, {p, 0, Infinity}, WorkingPrecision->20]
NIntegrate[integrand, {p, 0, Infinity}, WorkingPrecision->30]
NIntegrate[integrand, {p, 0, Infinity}, WorkingPrecision->40]


NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in p near {p} = {42.53581280075412}. NIntegrate obtained 4.67596994132151916. and 0.00989988117763158216. for the integral and error estimates.

4.675969941321519

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in p near {p} = {107.11269876829140479}. NIntegrate obtained 5.830371115131344422920. and 0.004895820216559319638220. for the integral and error estimates.

5.8303711151313444229

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in p near {p} = {1220.52200359274338453350425272}. NIntegrate obtained 8.6596599770482398526461597489530. and 0.039520892417721033775047475701130. for the integral and error estimates.

8.65965997704823985264615974895

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in p near {p} = {17591.42010943322603964222883582026742128}. NIntegrate obtained 11.4249891106069888138379180795912250891540. and 0.0671623528066706730262224874787473598905740. for the integral and error estimates.

11.42498911060698881383791807959122508915

Note that the integral is not converging when using higher WorkingPrecision. Also, note the first error message, which suspects a singularity as a possible issue. I think this is what is happening, the integrand is not convergent.

• That behaviour of Cot really is strange, especially seeing as N@Cot[4*^14 I] and Cot[4.*^14 I] both don't throw an error. – numbermaniac Oct 31 '17 at 5:05
• Thanks for this great analysis! It helped a lot. I went ahead and reported the unexpected behavior of Cot[] for large imaginary arguments to Wolfram Support. – Casimir Oct 31 '17 at 8:31
• Wolfram Support got back to me saying they were able to reproduce the issue in 11.2. They reported this behavior to the developer team and will update me once the issue is resolved. – Casimir Nov 4 '17 at 6:42

This is due to too less numerical precision maybe jof the RootSum.

Apply ToRadicals and FullSimplify to get a simple espression

integrand2 = integrand // ToRadicals // FullSimplify[#, p > 0] &

(*     -(I p^2 (Sqrt[
2] (1 + p^2) ((-I + Sqrt) Sqrt[-1 + I Sqrt - 2 p^2]
Cot[Sqrt[-1 - I Sqrt - 2 p^2]/(
2 Sqrt)] - (I + Sqrt) Sqrt[-1 - I Sqrt - 2 p^2]
Cot[Sqrt[-1 + I Sqrt - 2 p^2]/(2 Sqrt)]) -
2 I Sqrt[-1 - I Sqrt - 2 p^2]
Sqrt[(-1 + I Sqrt - 2 p^2) (1 + p^2)]
Coth[Sqrt[1 + p^2]/2]))/(2 Sqrt[-1 - I Sqrt - 2 p^2]
Sqrt[-1 + I Sqrt - 2 p^2] (1 + p^2))     *)

Plot[integrand2, {p, 0, 1000}] +++Appendix+++

I just see that NIntegrate still has problems. I try to solve this.

• Thanks for your help! The simplification with ToRadicals was already of some help (even though I too still can't perform the integration without errors or warnings). – Casimir Oct 30 '17 at 20:00

For integrand2 with M8, and WorkingPrecision ->100, i get

The integrand has evaluated to Overflow, Indeterminate, or Infinity for all \
sampling points in the region with boundaries (0.*10^-101   \
4.64781594643957803697370913786877210668235121444751683989453370359290\
1605966095761890817530152292030*10^14
`

With M11 I got the MemoryAllocatonFailure, and during evaluation, in Task Manager I watched the kernel fill up my entire 32GB ram.

• Thanks for the confirmation that it's not just my system. – Casimir Oct 30 '17 at 20:08