Is there a way to deal with this singularity in numerical integration? I would like to compute numerically, e.g., using the standard method of trapezes the following definite integral over the interface $[0,1]$:
$$
I = \int_0^1 \frac{f(x)}{\sqrt{1-x^2}} \, \mathrm{d} x \, , 
$$
where $f(x)$ is continuous in the interval [0,1]. 
It can readily be shown analytically that the integral is well convergent. 
However, when proceeding numerically, difficulties arise since the integrand diverges at the upper limit of integration for $x=1$.
I was wondering whether there exists a procedure that can help to remove the singularity in this integral.
Any help is highly appreciated  
Thank you,
hartmut 
 A: Since your integral is over a subset of the interval $(-1,1)$, your $f$ is continuous on that interval, and your integrand only has singularities on the boundary of that interval, tanh-sinh quadrature is a candidate method for you.
This will change your interval of integration to $[0,\infty]$, but the integrand will decay doubly-exponentially, (e.g., like $\mathrm{e}^{-\mathrm{e}^n}$), so truncation at finite right endpoint need not be significant.
You can continue your plan of using the trapezoidal method with step size $h$ by using the sample points $x_k = \tanh(\frac{\pi}{2} \sinh kh)$ and weights $w_k = \frac{h \pi}{2} \cdot \frac{\cosh kh}{\cosh^2( \frac{\pi}{2} \sinh kh)}$.  For instance, the area of your first trapezoid would be 
$$  A_0 = \frac{x_1 - x_0}{2} \left(\frac{w_0 f(x_0)}{\sqrt{1 - x_0^2}} + \frac{w_1 f(x_1)}{\sqrt{1 - x_1^2}}  \right)  \text{.}  $$
Just using $x_k$ and $w_k$ in a left Riemann sum, (and halving $w_0$ when summing this way) convergence is exponential: halving $h$ doubles the number of correct digits.  Error analysis suggests the trapezoidal rule applied here should improve convergence slightly (although perhaps is not worth the "effort" compared to just computing $\sum_{k=0}^\infty w_k \frac{f(x_k)}{\sqrt{1-x_k^2}}$, stopping when the terms are small enough to meet your precision and accuracy requirements or when either $1-x_k$  or $w_k$ is so small you cannot easily represent it).
Example:  $f(x) = 1 $ : 
The value of the integral is $\frac{\pi}{2} = 1.57079\dots$.  Using left-Riemann summation (and halving $w_0$), with $h = 1/10$, and stopping when $k=20$ (because $1-x_k$ and $w_k$ are about $10^{-5}$), the sum is $1.5659 \dots$.  Continuing on to $k = 27$ obtains $1.57079\dots$.


*

*With $h = 1/20$, stopping at $k = 40$ ($1 - x^k$ around $10^{-5}$ and $w_k$ around $10^{-36}$), obtains $1.56504\dots$.

*With $h = 1/20$, stopping at $k = 54$ ($1 - x^k < 10^{-6}$ and $w_k$ around $10^{-150}$), obtains $1.57078\dots$.  (An arbitrary precision calculation gives an error of about $5 \times 10^{-12}$ when we sum $70$ terms.)


We're only summing tens of terms to get these results...
A: It's easy to convert your integral into one without any singularities via the substitution $x=\sin\theta$:
$$I=\int^{\pi/2}_0f(\sin\theta)\,d\theta.$$
If you're lucky to have a function $f$ that is even (so that $I$ is just a quarter of the integral over the full period, from $-\pi$ to $\pi$), you might be pleasantly surprised by the speed of convergence of the trapezoidal rule in this special situation.
A: I like Professor Vector's method, but here is a way to take advantage of the integrability of $\frac1{\sqrt{1-x^2}}$ on $(0,1)$.
Using the integral
$$
\int_0^1\frac{a+bx}{\sqrt{1-x^2}}\,\mathrm{d}x=\frac{\pi a}2+b
$$
we can write
$$
\int_0^1\frac{f(x)}{\sqrt{1-x^2}}\,\mathrm{d}x=\int_0^1\frac{f(x)-(1-x)f(0)-xf(1)}{\sqrt{1-x^2}}\,\mathrm{d}x+\left(\frac\pi2-1\right)f(0)+f(1)
$$
If $f$ is smooth at $0$ and $1$, the integrand on the right vanishes at $0$ and $1$

A Lesson in Overthinking
I was wondering why $f(0)$ and $f(1)$ had different weights. Then I realized that $\frac1{\sqrt{1-x^2}}$ doesn't have a singularity at $0$. We don't need to subtract $(1-x)f(0)$ at all. While the formula I give above is correct, I was thinking of $\sqrt{x(1-x)}$ in the denominator, which does have a singularity at $0$ and $1$. Doh!
A: A general technique is to remove the singularity by moving it to an analytically integrable function:
$$\int_0^1\frac{f(x)}{\sqrt{1-x^2}}dx=\int_0^1\frac{f(x)-f(1)+f(1)}{\sqrt{1-x^2}}dx=-\int_0^1\frac{f(x)-f(1)}{\sqrt{1-x^2}}dx+f(1)\left.\arcsin(x)\right|_0^1.$$
The limit of the integrand at $x=1$ should be finite.

For instance with $f(x):=x$,
$$\frac{x-1}{\sqrt{1-x^2}}=-\sqrt{\frac{1-x}{1+x}},$$ which is well-behaved.
A: The other answers so far focus on how to transform the integral so that you don't have a singularity. But that's not always viable so it's important to know how to deal with a singularity when one exists.
The trapezoid rule is a poor choice for this case, because it uses the values of the function at the endpoints, which can be infinite. Simposon's rule will have a similar problem.
More appropriate is the more primitive rectangle rule. That is, divide the domain to segments, and multiply each segment's width by the value of the function at the midpoint of the segment. This will work since you never try to sample the function at the infinite endpoints.
Alternatively, if you want faster convergence, you can use the more sophisticated and extremely powerful Gaussian Quadrature. More complicated but well worth knowing. Here, too, you divide the domain to intervals. For each interval, you sample the function at a few carefully chosen points to obtain surprisingly high accuracy.
There are different methods depending on the number of points in each interval, but the simplest is the 2-point method. And again, it will work with singularities, since the endpoints of intervals aren't sampled.
There is even a variant of Gaussian quadrature tailored specifically for a well-behaved function multiplied by $\frac{1}{\sqrt{1-x^2}}$, called Chebyshev–Gauss quadrature. To use it, you don't even have to divide to intervals - you just choose $n$ (the higher $n$, the better the accuracy) and you have
$$\int_{-1}^1\frac{f(x)}{\sqrt{1-x^2}}\approx\frac{\pi}{n}\sum_{i=1}^nf\left(\cos\left(\frac{2i-1}{2n}\pi\right)\right)$$
Which is very similar to transforming the integral as some other answers suggested. It will maximize the accuracy you can have for a given number of samples from the function $f$, and convergence is exponential.
A: The standard way is $x = 1 - t^{2}$:
$$
\int_{0}^{1}{\operatorname{f}\left(x\right) \over
\,\sqrt{\,{1 - x^{2}}\,}\,}\,\mathrm{d}x =
2\int_{0}^{1}{\operatorname{f}\left(1 - t^{2}\right) \over
\,\sqrt{\,{2 - t^{2}}\,}\,}\,\mathrm{d}t
$$
