Numerical integration of legendre polynomials I hope that numerical questions are also permitted here.
I want to expand a smooth functions $f \in C^{\infty}$in terms of Legendre polynomials. Thus I need to calculate integrals of the form $\int_{-1}^{1} f(x) P_n(x)$, where $n$ becomes sufficiently large (between 40 and 80). In that regime, the Legendre polynomials oscillate heavily, so my question is whether you are aware of a particular method that is good for integrating such things, cause the standard MATLAB method ( probably some low-level Newton Cotes method) cannot do it sufficiently accurate. 
Although this is not a pure math question, I think that this question is of some particular value in applications, so please be not to hard with me ;-)
 A: See:
@article {MR1078802,
    AUTHOR = {Alpert, Bradley K. and Rokhlin, Vladimir},
     TITLE = {A fast algorithm for the evaluation of {L}egendre expansions},
   JOURNAL = {SIAM J. Sci. Statist. Comput.},
  FJOURNAL = {Society for Industrial and Applied Mathematics. Journal on
              Scientific and Statistical Computing},
    VOLUME = {12},
      YEAR = {1991},
    NUMBER = {1},
     PAGES = {158--179},
      ISSN = {0196-5204},
     CODEN = {SIJCD4},
   MRCLASS = {65D20 (41A50)},
  MRNUMBER = {1078802 (91i:65042)},
       DOI = {10.1137/0912009},
       URL = {http://dx.doi.org/10.1137/0912009},
}

A: It seems fundamentally ill-conditioned.
Since $\int_{-1}^{+1} x^rP_n(x),dx=0$ for $r=0,1,\dots,n-1$, your integral is unchanged if you subtract a polynomial of degree $n-1$ from $f(x)$.  I'm guessing that if you can do that very accurately for some polynomial that approximates $f(x)$, the resulting integral won't be so ill-conditioned as the original.
A: I will make the following assumptions


*

*you can evaluate your function wherever your like.

*you are unconcerned with machine precision (this should not be an issue for the number of coefficients you are after).


Then the solution is to use an n-point Gauss-Legendre quadrature in order to optimally extract n Legendre coefficients.
To obtain this result, the key point is that the function space spanned by the first $n$ Legendre polynomials is simply the space of the zero-th to n-th polynomial moments. More specifically:
$$ \mathbb{P}_n\equiv\mathrm{span}(P_0(x),P_1(x),P_2(x),\ldots,P_n(x)) \equiv \mathrm{span}(1,x,x^2,\ldots,x^n).$$
Hence to perfectly integrate the $\{0,\ldots,n-1\}$ Legendre moments, we are really looking to integrate an order-$2n-2$ polynomial exactly. We are looking for an answer to the following.
Question. Given that $f$ is a polynomial of at most order $n$, is there a numerical integration scheme that will integrate the product $f$ exactly?
In this form, the question has a known, provable solution.
Theorem 1 (Gauss-Legendre Quadrature). The $n$-node Gaussian quadrature scheme, whose nodes are defined at the roots of the order $n$ Legendre polynomial, and whose weights are defined via Lagrange interpolants, will peform the integration exactly$$\int_0^1f(x)\,dx=\sum_{i=1}^m w_i f(x_i)$$for all $f\in\mathbb{P}_{2n-1}$, where $\mathbb{P}_{2n-1}$ denotes the space of order $2n-1$ polynomials.
In terms of where to find code for Gauss-Legendre Quadrature, I would look for it on MATLAB central file exchange.
A: The LegendreSeries command of the OrthogonalExpansions package of Maple does the job.
