I'ld like to calculate the PV of an integral with the form $$ \tilde{G}_l(\omega) = -\frac{2\omega}{\pi} PV\int_0^\infty \frac{\tilde{G}_d(\omega^\prime)}{\omega^2 - {\omega^\prime}^2}d\omega^\prime$$ in MATLAB. It is not obvious to me how to do this. I've tried using factorization to change the power in the denominator to allow for use with the Hilbert function in MATLAB, but that hasn't worked. Any insight would be appreciated.

$\tilde{G}_d(\omega)$ is an even function so I feel like I could do a contour integration and pick up the residue in the right-half plane, but I would rather do it numerically, as $\tilde{G}_d(\omega)$ exists more completely as a numerical function than an analytic one.


1 Answer 1


Away from the singular point the standard quadrature methods should be fine, so the problem is really about computing an integral of the form $$ \mathrm{P.V.}\int_{-h}^{h} \frac{f(x)}{x}\,dx. $$

I would suggest trying the following things:

If you know a closed-form Taylor series for $f(x)=G(\omega+x)$, the you can choose $h$ to small enough for the series approximation to be accurate, and calculate the principal value in closed form for each series term.

Alternatively, you can try to interpolate $f$ numerically, using something like polyfit, which will give you the Taylor series numerically. Depending on how well-behaved your function is, it is difficult to say from the start whether any of this will definitely work.

Still the simplest way (although numerically unstable) is to write the integral as $$ \int_0^h \frac{f(x)-f(-x)}{x}\,dx, $$ and apply the standard quadrature algorithm to this. The issue is that the calculation $f(x)-f(-x)$ will lead to inaccuracies on the order of $\int_0^h (\epsilon_{\mathrm{mach}}/x)\,dx$, which may or may not be tolerable. If you have a way to simplify the expression $f(x)-f(-x)$ to make it numerically stable, this is probably the best option.

There are ways to derive quadrature nodes and weights for computing principal values of such integrals (just search for them, even GSL has them), but I couldn't find this in matlab.

These are numerical methods for principal value integrals on the real line. You can also try to compute the integral as $$ \mathrm{P.V.}\int_{-h}^{h}\frac{f(x)}{x}\,dx = \int_\gamma \frac{f(x)}{x}\,dx + \pi \mathrm{i}f(0), $$ where $\gamma$ is a contour going from $-h$ to $h$ passing above the real line. If you can evaluate $G$ for complex $\omega$, then even standard quadrature routines will be able to compute the integral over $\gamma$.

  • $\begingroup$ I'm going to try some of this and look into the GSL implementation. $\endgroup$
    – user7257
    Oct 17, 2014 at 14:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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