You probably find the answer in the book W. Gautschi, Orthogonal polynomials--Computation and Approximation, Oxford University Press, 2014, or in one of the author's papers on the subject.

Dear @Peyman, if the function has a singularity then it's integral is improper. Is the integral bounded or not?

Yes, assuming the integrals involved exist and are finite, if you have n nodes, you expect a Gaussian quadrature rule to be exact for polynomials of degree 2n+1,  You can compute analytic (or HIGHLY precise approximate) integrals for x^i * w(x) (i = 0 : 2n) on your interval and solve a linear system for weights and node locations.  This is essentially what colleagues and I did 30 years ago in our SVALAQ Self-Validating Adaptive Quadrature package to find weights and nodes for weight functions not in the CRC tables of the day.

@Corliss : It will be exact for polynomials of degree 2n -1. It is possible if you are able to construct family of orthogonal polynomials  satisfying  a three-term recurrence relation. You will form corresponding tri-diagonal matrix and find its eigen values which are nothing but zeros of the orthogonal polynomial giving you nodes. Weights can also be calculated from the integrals. More details you may see in 

Introduction to Numerical Analysis by STOER, J. and BULLRISCH, R

Agreed.  Thanks for correcting and supplying more detail.

ORTHPOL (of Walter Gautschi) does this job. This is in

netlib.ornl.gov  in the library toms the number 726 .

The algorithm behind is all but trivial and there are numerical intricacies. Hence

simply using the well known three term recursion plus a zero finder is strongly

discouraged. Of course you may have only boundary singularities in your weight.

Otherwise you must decompose your interval appropriately.

(1) Please read a book of mathematical physics: integration of complex quantities and how to remove the singularities (e.g. residue theorem). There are different analytical tricks. It will help.

(2) You also can consult the papers:  Gauss-Legendre quadratures are used in JPB 33, 4285 (2000);  JPB 31, 4427 (1998);  JPB 35, 3365 (2002);  Pramana 63, 1063 (2004); Physica stat. sol. 6, 2281 (2009); etc.

Peyman,

There is no general rule to integrate functions in which the weight function has a singularity. We can however, in some circumstances, perform a transformation where the usual quadrature rules become applicable.

please show us the integral you want to evaluate. We can definitely help you to apply the mathematical procedure to follow.

Best of luck

Thank you all. The integration is of the form \int^x_0[(x-t)^{-beta} u(t)]dt.

I think Payman, you can also try the book:

Table on Integrals, Series, and Products, written by I. S. Gradshteyn and L. M. Ryzhik.

Define 'beta', how it varies? If 'beta' is integer, it is a pole type singularity, can be integrated easily; find the residue if the pole is of order 'beta',  you could apply residue theorem. If 'beta' is fraction, then there is diferent trick.  I think it will help. Please read a book for details.

You probably find the answer in the book W. Gautschi, Orthogonal polynomials--Computation and Approximation, Oxford University Press, 2014, or in one of the author's papers on the subject.

Hi Peyman Hessari:

1. Numerical Recipies could address your question. Look: 

http://www.aip.de/groups/soe/local/numres/bookcpdf/c18-3.pdf

2. Adicional references:

http://m.njit.edu/CAMS/Technical_Reports/CAMS09_10/report0910-4.pdf

http://arxiv.org/pdf/1207.4461v2.pdf