# How to numerically calculate residues?

I need to calculate the following integral: $${1\over 2\pi i} \int_C f(E) \, d E$$ $$f(E) = {\rm Tr}\,\left(({\bf h} + E)\,{\bf G}(E) \right)$$ Where $\bf h$ is a matrix (one particle kinetic and potential energy expressed in a basis), $\bf G$ is a matrix which depends on $E$ (one-particle many-body Green's function) and the contour integral is a left semicircle. The integrand $f(E)$ has poles on the negative real axis and it is expensive to evaluate. What is the most effective way to calculate such integral?

Here is my research so far:

1) I use Gaussian integration, my integration path is a rectangle. I fixed the left and right side (i.e. width) and played with height (above and below the real axis) such that for the given integration order I get the highest accuracy. For example for order 20, if the height is too large, the accuracy goes down (obviously), but if it is too small, it also goes down (my theory is that it needs more and more points around the poles as height goes to 0). I settled with optimal height 0.5 for my function.

2) Then I set the right side of the rectangle at E0, typically E0=0, but it could be E0=-0.2 or something similar.

3) I start moving the left side of the rectangle to the left and for each step I do integration order convergence to make sure my integral is fully converged for each rectangle. By increasing the width, I eventually get a converged value in the limit of the infinite left semicircle.

The computation is really slow and also not very accurate for large widths. One improvement is to simply partition the long width into "elements" and use Gaussian integration on each element (just like in FE).

Another option would be to integrate a small circle around each pole and sum it up. Problems:

a) How to numerically find the poles of the function $f(E)$? It should be robust. The only thing I know is that they are on the negative real axis. For some of them (but not all) I also know a pretty good initial guess. Does there exist a method that works for any analytic function $f(E)$? Or does it depend on the actual form of $f(E)$?

b) Once we know the poles, what numerical scheme is the best for integrating the small circle around it? Should I use Gaussian integration on a circle? Or should I use some uniform distribution of the points?

Another option might be that once I know the poles thanks to a), there might be some semi-analytic way to get the Residues without the need of the complex integration. But for now I'd be happy to just optimize the contour integration.

• Have you checked the book "Numerical Methods for Laplace Transform Inversion" by Cohen (2007)? IIRC, Robert Piessens (of QUADPACK fame) also worked on this topic. – GertVdE Jun 23 '12 at 16:47

I can offer a suggestion for your first question: If you know your poles are somewhere along the real axis, you could localize them quite efficiently using Rational interpolation/approximation. This amounts to finding polynomials $p(x)$ and $q(x)$ such that

$$f(x) \approx \frac{p(x)}{q(x)}$$

for $x$ in some interval. The poles of $f(x)$ should then match the roots of $q(x)$.

Rational interpolation/approximation can be a tricky thing, but I've recently co-authored a paper on a stable algorithm to compute them using the SVD. The paper contains Matlab code implementing the algorithm, and a more extensive version thereof is available as the function ratinterp in the Chebfun project, of which I am one of the developers.

For your second question, this paper may be useful.

• Thanks for all the tips! Here is the code netlib.org/toms/579 of the Bengt Fornberg paper. Unfortunately, there is some numerical bug, as this is the output that I am getting: gist.github.com/2942970#file_output. So I'll have to reimplement or debug it. The Chebfun link gives me 404 (I tried it couple months ago with the same results, so maybe it simply doesn't work from the USA). – Ondřej Čertík Jun 17 '12 at 0:27
• @OndřejČertík: I've never used the TOMS 579 code myself, so I don't know what to tell you about the errors. As for the Chebfun homepage, could you try "googling" it and seeing if it works then? – Pedro Jun 17 '12 at 1:46
• Google finds the Chebfun homepage and shows cached versions. But when I click on the page, this is what I get: pastehtml.com/view/c1ts4h3ct.html – Ondřej Čertík Jun 17 '12 at 20:32
• Try a different browser? Or from a different ISP. The website works fine from here (in USA.) – Costis Jun 23 '12 at 22:05
• I tried Firefox and Chrome. So it must by my ISP. Weird. – Ondřej Čertík Jun 25 '12 at 21:38