I succeeded proving this now. All comments about contour integration should be forgotten. The result is true for all continuous functions, even for those which are not differentiable. So one should not get into complex analysis. The Wikipedia page about Sokhatsky-Weierstrass theorem outlines a correct way, which begins with:
<br />
\frac{1}{x \pm i\epsilon} = \mp \frac{i\epsilon}{x^2 + \epsilon^2} + \frac{x}{x^2 + \epsilon^2}<br />
We can expect that the reader knows how to prove
<br />
\frac{\epsilon}{x^2 + \epsilon^2} \underset{\epsilon\to 0^+}{\to} \pi \delta(x).<br />
So the real task we are left with is to prove this:
jostpuur said:
<br />
\lim_{\epsilon\to 0^+} \int\limits_a^b \frac{x^2}{x^2+\epsilon^2} \frac{f(x)}{x}dx = \lim_{\epsilon\to 0^+} \int\limits_{[a,-\epsilon]\cup [\epsilon,b]} \frac{f(x)}{x}dx<br />
If we assume it known that the principal value exists, then this equation is equivalent with the equation
<br />
\lim_{\epsilon\to 0^+}\int\limits_a^b\Big( \frac{x^2}{x^2 + \epsilon^2} \;-\; \big(1 - \chi_{[-\epsilon,\epsilon]}(x)\big)\Big)\frac{f(x)}{x} dx =: \lim_{\epsilon\to 0^+} I(\epsilon) = 0<br />
For symmetry reasons, it is sufficient to show that this integral vanishes when it is integrated over the values x>0. So we can set a=0. I split the integral in two pieces, over [0,\epsilon] and over [\epsilon,b].
<br />
\int\limits_0^{\epsilon} \frac{x}{x^2+\epsilon^2} f(x) dx<br />
and
<br />
\int\limits_{\epsilon}^b \Big(\frac{x^2}{x^2 + \epsilon^2} - 1\Big)\frac{f(x)}{x} dx = -\int\limits_{\epsilon}^b \frac{\epsilon^2}{x^2 + \epsilon^2} \frac{f(x)}{x} dx<br />
The integral over the [0,\epsilon] can be calculated like this:
<br />
(\inf_{0<x<\epsilon}f(x)) \int\limits_0^{\epsilon} \frac{x}{x^2 + \epsilon^2} dx \leq \int\limits_0^{\epsilon}\frac{xf(x)}{x^2 + \epsilon^2} dx \leq (\sup_{0<x<\epsilon}f(x)) \int\limits_0^{\epsilon} \frac{x}{x^2 + \epsilon^2} dx<br />
The infimum and supremum are going to approach f(0).
<br />
\int\limits_0^{\epsilon} \frac{x}{x^2 + \epsilon^2} dx = \frac{1}{2}\log(2)<br />
So the integral over [0,\epsilon] approaches the number \frac{1}{2}\log(2)f(0). The integral over [\epsilon, b] can be calculated with a variable change x=\epsilon y. It becomes
<br />
-\int\limits_1^{b/\epsilon} \frac{f(\epsilon y)}{(1+y^2)y} dy<br />
If f is bounded on [a,b], which it is since it is continuous, then we can take the limit f(\epsilon y)\to f(0) outside the integral. We get
<br />
-f(0) \int\limits_1^{\infty} \frac{1}{(1+y^2)y} dy = -f(0)\int\limits_1^{\infty}\Big(\frac{1}{y} \;-\; \frac{y}{1+y^2}\Big)dy = -\frac{1}{2}\log(2)f(0) \;-\; \lim_{R\to\infty}\Big(\log(R) - \log(\sqrt{1 + R^2})\Big) = -\frac{1}{2}\log(2)f(0)<br />
So the integrals over [0,\epsilon] and [\epsilon,b] cancel on the limit \epsilon\to 0^+, and we get the desired result I(\epsilon)\to 0.