Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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$.

share|improve this question
    
Not an answer to your question ("why does it give this message"), but how about Total[Table[Abs@integrator[n, 1], {n, 1, 25}]] –  acl Dec 13 '12 at 16:19
    
I have been using 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 a ListPlot. –  fpghost Dec 13 '12 at 16:36
    
you mean that in your problem you use justNSum[integrator[n,1],{n,1,Infinity}]? –  acl Dec 13 '12 at 17:01
    
@acl yes, that would be the ultimate aim, but even NSum[integrator[n,1],{n,1,5}] for example fails. Despite the fact I can compute integrator[1,1],integrator[2,1], integrator[3,1],... individually (or in a Table ) without issue. –  fpghost Dec 13 '12 at 19:29

1 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).

share|improve this answer
    
thanks. Accumulate is a nice function I had yet to come across. I would still kind of like to know the reason NSum is failing for curiosity's sake. As for replacing with just Sum I had tried this and got the error regarding the derivative which stopped me using it further, how did you get around this? –  fpghost Dec 13 '12 at 19:27
1  
I think the only way you could get the derivative error with Sum is if you had called N[Sum[...]], but with Sum alone I had to do nothing to make it work. NSum is the one that gives trouble because it tries to be too clever. It's designed to avoid cancellation erros and to use extrapolation of the summand so it can do infinite series. But those methods require it to take the derivative of the summand. On the other hand, Sum by itself does not do that at all. –  Jens Dec 13 '12 at 19:44
    
Oh how odd, I just tried Sum alone and it did work as you say. I could have sworn that I tried that and got into bother. Maybe I was still trying the Infinite sum at that stage of my experimentation. OK so if NSum cannot work to give me the infinite sum, then it looks like my best bet is manual cutoff that I spot via eye looking at the Accumulate results and a ListPlot of points in energy. –  fpghost Dec 13 '12 at 19:54
    
"which is part of its evaluation method" - yes, the Euler-Maclaurin method involves taking a few derivatives. A change of method (e.g. use "WynnEpsilon" instead of "EulerMaclaurin") might be useful in this case. –  J. M. Apr 11 '13 at 0:23

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.