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.
How to numerically calculate residues?
numerics;complex analysis
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.