Numerical integration of divergent function 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$.
 A: 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$.
A: @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.
