Why does NSum fail here?

Question

I want to evaluate a sum of integrals; each integral has a pole on the real axis and I handle this via the Cauchy Principal Value

$ f(E)=-\frac{a^2}{2q}\sum^{\infty}_{n=-\infty}\int^{\infty}_{-\infty}dr \exp{(-2iEr/q)}\frac{1}{\left[\sinh^2{(r)}-(\frac{rp}{q}+\frac{npP}{2})^2\right]} $

$a,p,q, E$ are constants. I should say first of all the following applies for $n\neq 0$, for zero $n$ I handle differently. For each $n$ and $E$ the denominator of the integrand has two poles on the real axis (amongst others in the complex plane I shall not be concerned with). I will use Mathematica's FindRoot to get the locations of these poles, since the equation whose root we need is transcendental. I will then use NIntegrate with the CauchyPrincipalValue method to evaluate the principle parts of my integral. To get the additional contribution from the $i\pi \times \text{residue}$ at each pole, I first use partial fractions on my integrand to obtain

$ f(E)=-\frac{a^2}{4q}\sum^{\infty}_{n=-\infty}\int^{\infty}_{-\infty}dr \frac{\exp{(-2iEr/q)}}{\left(\frac{rp}{q}+\frac{npP}{2}\right)}\left(\frac{1}{\left[\sinh{(r)}-(\frac{rp}{q}+\frac{npP}{2})\right]}-\frac{1}{\left[\sinh{(r)}+(\frac{rp}{q}+\frac{npP}{2})\right]}\right). $

Now the two poles for each $n,E$ are simple poles and providing we know their location (numerically via FindRoot) the residue is easy to compute for each. For example, considering the first term in the integrand

$-\frac{a^2}{4q}\frac{\exp{(-2iEr/q)}}{\left(\frac{rp}{q}+\frac{npP}{2}\right)}\frac{1}{\left[\sinh{(r)}-(\frac{rp}{q}+\frac{npP}{2})\right]}$

and say we have numerically found the root to be at $r_0$ for this piece, then

$\text{Res}=\lim_{r\to r0}\left[-\frac{a^2}{4q}\frac{\exp{(-2iEr/q)}}{\left(\frac{rp}{q}+\frac{npP}{2}\right)}\frac{(r-r_0)}{\left[\sinh{(r)}-(\frac{rp}{q}+\frac{npP}{2})\right]}\right]$ $\text{Res}=\lim_{r\to r0}\left[-\frac{a^2}{4q}\frac{\exp{(-2iEr/q)}}{\left(\frac{rp}{q}+\frac{npP}{2}\right)}\frac{1}{\left[\cosh{(r)}-\frac{p}{q}\right]}\right]$

the second equality following from L'Hopital's rule. Note the zeroes of $\left(\frac{rp}{q}+\frac{npP}{2}\right)$ are not really poles when the full expression is considered together as further above. There will be a similar residue from the other piece, and we must add $i\pi\times$ both of these residues to our PV.

Mathematica Implementation

(*some constants and definitions*)
M = 1;
R = 10 M;
wp = 40;
ag = wp - 8;
pg = wp/2;
a = 1/(4 M) Sqrt[R/(R - 2 M)];
p = 1/(4 M) Sqrt[(M R)/((R - 2 M) (R - 3 M))];
q = 1/(4 M) Sqrt[R/(R - 3 M)];
P = 2 Pi Sqrt[(R^2 (R - 2 M))/M];

(*find the poles for each n,E*)

findNroots[n_?IntegerQ, x0_?NumericQ] := Block[{rM, rP, roots}, 
rP = r /. FindRoot[Sinh[r] + ((r p)/q + (n p P)/2) == 0, {r, x0}, PrecisionGoal -> pg, WorkingPrecision -> wp][[1]];
rM = r /. FindRoot[Sinh[r] - ((r p)/q + (n p P)/2) == 0, {r, x0}, PrecisionGoal -> pg, WorkingPrecision -> wp][[1]];
roots = {rP, rM};
roots  
]

(*integrate for given n, E, summing PVs and residues*)

 integrator[n_ /; n != 0, En_?NumericQ] := Block[{rtM, rtP, lowRoot, highRoot, mid, tmp1, tmp2, tmp3, tmp4, sum},

 rtP = findNroots[n, 0][[1]];
 rtM = findNroots[n, 0][[2]];
 lowRoot = Sort[{rtM, rtP}][[1]];
 highRoot = Sort[{rtM, rtP}][[2]];
 mid = (highRoot + lowRoot)/2;  (*const somewhere in middle of two poles on real line*)

 tmp1 = -a^2/(2 q) NIntegrate[ Exp[-2 I En r /q] (1/((Sinh[r])^2 - ((r p)/q + (n p P)/2)^2)), {r, -Infinity, lowRoot, mid}, Method -> {"PrincipalValue"}];
 tmp3 = -a^2/(2 q) NIntegrate[Exp[-2 I En r /q] (1/((Sinh[r])^2 - ((r p)/q+ (n p P)/2)^2)), {r, mid, highRoot, Infinity}, Method -> {"PrincipalValue"}];

 (*add residues too:*)
 tmp2 = I Pi ( -a^2/(4 q) Exp[-2 I En rtM /q] 1/((rtM p)/q + (n p P)/2) ) Limit[1/( Cosh[r] -  p/q), r -> rtM];
 tmp4 = I Pi ( -a^2/(4 q) Exp[-2 I En rtP /q] 1/((rtP p)/q + (n p P)/2) ) Limit[-1/( Cosh[r] +  p/q), r -> rtP];

 sum = tmp1 + tmp2 + tmp3 + tmp4;
 {sum}
 ]

This seems to work OK, I can do Table[integrator[n,1],{n,1,100}] and Mathematica seems to compute this list just fine. However if I do justNSum[integrator[n,1],{n,1,2}] then Mathematica will give all manner of complaints such as NIntegrate::deorela. Why is this when it is perfectly happy to compute the values of integrator in a list individually? Is there another way I should be performing this sum?

A separate issue: I also find that at certain values of $n,E$ for example integrator[1, 24] I also get errors like NIntegrate::deodiv: DoubleExponentialOscillatory returns a finite integral estimate, but the integral might be divergent and I can't even compute the result individually let alone in NSum. I suspect this is because things become too oscillatory at higher values of $E$.

Answer 1

If I replace NSum by Sum, your code works fine. With NSum, you'd have to first find a way to suppress forming the derivate with respect to n, which is part of its evaluation method.

But reading your comment, I would suggest doing something slightly different anyway. As you noted, the Table gives no problems, but you would also like to adjust the number of terms in the sum manually by looking at the intermediate results. Then why not do this:

t = Table[integrator[n, 1], {n, 1, 100}];
sums = Accumulate[t]

This outputs a long list, and each entry sums[[m]] is equal to the sum $\sum_{n=1}^{m} \text{integrator}(n,1)$

By the way, you may also want to replace the Limit by Residue (but that's not related to your actual question).