Numerical Methods – Numerical Integration of Divergent Function

Question

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

Answer 1

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