# Numerical integration of divergent function

I am having trouble with the numerical integration of a divergent function. For example, $$n= \int f(x)\,dx = \displaystyle\int \dfrac{\Theta(x-\varepsilon)\,dx}{\sqrt{x-\varepsilon}}$$ is not behaving good when calculated numerically: $$n_\text{numerical} \approx f(x_j)(x_j-\varepsilon) + \sum_{k=j+1}^N f(x_k)\Delta x$$ with $$x_i = x_\text{min} + i\Delta x$$ and $x_j$ being the smallest point larger than $\varepsilon$. Due to the divergent behavior of $f(x)$ near $\varepsilon$, the numerical result has discontinuities and jumps at every $\varepsilon=x_i$ when I sweep $\varepsilon$. The smoothness of this integration is crucial in the convergence of my simulations. Can anybody suggest a simple trick to solve the problem? All I need is to have a continuous numerical result by continuously changing $\varepsilon$.

I suppose your function is $f(x) = g(x)/\sqrt{x-\epsilon}$ where $g$ is bounded, and you want to integrate this on $(\epsilon, b)$. The change of variables $x = \epsilon + t^2$ gives you $$\int_\epsilon^b \dfrac{g(x)}{\sqrt{x-\epsilon}}\; dx = 2 \int_0^{\sqrt{b-\epsilon}} g(\epsilon + t^2)\; dt$$ which gets rid of the singularity. Similarly to integrate $g(x)/\sqrt{\epsilon - x}$ on $(a,\epsilon)$, use $x = \epsilon - t^2$.

@RobertIsrael provided the way forward here. I thought that it would be instructive to show a simple "trick" that works well in other situations.

To that end, we assume that $g(x)$ satisfies the Holder condition

$$|g(x)-g(y)|\le C|x-y|^{\alpha} \tag 1$$

with Holder exponent $\alpha \ge 1/2$. Then we can write

\begin{align} \int_{\epsilon}^{b}\frac{g(x)}{\sqrt{x-\epsilon}}dx&=\int_{\epsilon}^{b}\frac{g(x)-g(\epsilon)}{\sqrt{x-\epsilon}}dx+g(\epsilon)\int_{\epsilon}^{b}\frac{1}{\sqrt{x-\epsilon}}dx\\\\ &=\int_{\epsilon}^{b}\frac{g(x)-g(\epsilon)}{\sqrt{x-\epsilon}}dx+2g(\epsilon)\sqrt{b-\epsilon}\tag 2 \end{align}

The integrand $\frac{g(x)-g(\epsilon)}{\sqrt{x-\epsilon}}$ of the integral in $(2)$ will be well-behaved for numerically small $\epsilon>0$ because of the Holder condition of $(1)$ with exponent $\alpha \ge 1/2$.

NOTE:

This approach in which we "move the offending term" works well in a fairly general class of problems. Suppose that we have an improper, convergent Riemann integral

$$\int_a^bf(x)dx$$

for which $f$ has a non-removable singularity at $a<\xi<b$ such that $\lim_{x\to \xi}\left(f(x)-K|x-\xi|^{-\beta}\right)$, $0<\beta<1$, exists and is finite. Then, we can write

\begin{align} \int_a^bf(x)dx&=\int_a^b\,(f(x)-K|x-\xi|^{-\beta})\,dx+\int_a^b\,K|x-\xi|^{-\beta}dx\\\\ &=\int_a^b\,(f(x)-K|x-\xi|^{-\beta})\,dx+K\frac{(b-\xi)^{1-\beta}+(\xi-a)^{1-\beta}}{1-\beta}\tag 3 \end{align}

where the integrand $f(x)-K|x-\xi|^{-\beta}$ is well-behaved near $\xi$ and numerical integration poses no challenge.

• @maziarnoei Please let me know how I can improve my answer. I just want to give you the best answer I can. Commented Jul 15, 2015 at 3:31