Show description

We first form the (n + 1) × (n + 1) Gram matrix M = ( xi , xj )0≤i,j≤n = (mi+j )0≤i,j≤n which is symmetric positive definite, and use the Cholesky decomposition to write M = Rt R where R = (ri,j )0≤i≤n is an upper triangular matrix. R=A$. 1. (1) Note that here the indices vary from 1 to n instead of from 0 to n. (2) The program works only when the matrix is nonsingular, more precisely when the ri,i do not become too small. Then if the (pn ) are orthogonal polynomials associated to the weight w with moments mk , the pn satisfy the following 3-term relation with βn+1 pn+1 (X) = (X − αn )pn (X) − βn−1 pn−1 (X) rn+1,n+2 rn,n+1 rn+1,n+1 and αn = − , rn,n rn+1,n+1 rn,n is interpreted as 0.

Gaussian Integration With Polynomials in 1/x When computing integrals on a compact interval [a, b], using Gaussian integration based on orthogonal polynomials (possibly multiplied by functions such as (x − a)α (b − x)β ) is quite reasonable, in particular since polynomials are dense in continuous functions. However the same is not true on noncompact intervals such as [a, ∞[: when the function to be integrated tends slowly to 0, such as a rational function, changes of variable such as x → 1/x transform the integral into an integral over a compact interval of a reasonable function, which can then sometimes be reasonably Gaussian integrated.

To be able to use this formula, we need an estimate on the growth of bm so as to see where we can truncate the series. This is easily given by a standard theorem of complex analysis: the radius of convergence of the series for F (1/X) around 0 is equal to the smallest modulus of a pole of F (1/X), so equal to 1/ max |α|, where as above α ranges over all the poles of F . Thus bm is of the order of (max |α|)m , and neglecting the coefficient N/(m − 1), to obtain D decimals of accuracy we need (max |α|/N )m < 10−D , in other words N > max |α| and m > D log(10)/ log(N/ max |α|).