NIntegrate reacts strangely to changes in WorkingPrecision

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

Cot[4`16*^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[3]) (1 + 2 p^2) cot[
   Sqrt[-1 - I Sqrt[3] - 2 p^2]/(2 Sqrt[2])])/(
 3 Sqrt[-2 - 2 I Sqrt[3] - 
   4 p^2]) + ((3 + I Sqrt[3]) (1 + 2 p^2) cot[
   Sqrt[-1 + I Sqrt[3] - 2 p^2]/(2 Sqrt[2])])/(
 3 Sqrt[-2 + 2 I Sqrt[3] - 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.


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[3]) Sqrt[-1 + I Sqrt[3] - 2 p^2]
       Cot[Sqrt[-1 - I Sqrt[3] - 2 p^2]/(
       2 Sqrt[2])] - (I + Sqrt[3]) Sqrt[-1 - I Sqrt[3] - 2 p^2]
       Cot[Sqrt[-1 + I Sqrt[3] - 2 p^2]/(2 Sqrt[2])]) - 
  2 I Sqrt[-1 - I Sqrt[3] - 2 p^2]
    Sqrt[(-1 + I Sqrt[3] - 2 p^2) (1 + p^2)]
    Coth[Sqrt[1 + p^2]/2]))/(2 Sqrt[-1 - I Sqrt[3] - 2 p^2]
 Sqrt[-1 + I Sqrt[3] - 2 p^2] (1 + p^2))     *)

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

enter image description here

+++Appendix+++

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


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.