I am having trouble with the numerical integration of a divergent function. For example, \begin{equation} n= \int f(x)\,dx = \displaystyle\int \dfrac{\Theta(x-\varepsilon)\,dx}{\sqrt{x-\varepsilon}} \end{equation} is not behaving good when calculated numerically: \begin{equation} n_\text{numerical} \approx f(x_j)(x_j-\varepsilon) + \sum_{k=j+1}^N f(x_k)\Delta x \end{equation} with \begin{equation} x_i = x_\text{min} + i\Delta x \end{equation} 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$.


2 Answers 2


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$.


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


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.

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

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .