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$.
Total[Table[Abs@integrator[n, 1], {n, 1, 25}]]
$\endgroup$Total
to sum the lists already but it's quite inconvenient for me to have to manually find what cutoff the sum converges at (i.e by continually increasing your 25 until the number stops changing). I would have to find this cutoff for each energy $E$, as I would like to plot a graph against $E$. I guess if there is no other solution, then I will end up doing this for say 10 values of $E$ and doing aListPlot
. $\endgroup$justNSum[integrator[n,1],{n,1,Infinity}]
? $\endgroup$NSum[integrator[n,1],{n,1,5}]
for example fails. Despite the fact I can computeintegrator[1,1],integrator[2,1], integrator[3,1],...
individually (or in aTable
) without issue. $\endgroup$