# 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 ;-)

• Are you aware of anything clever for Fourier series (which suffer from the same oscillation problem)? Sep 22, 2014 at 11:39
• @IgorRivin unfortunately not, sorry.
– user54300
Sep 22, 2014 at 11:46
• See if you can borrow some ideas in chebfun which is centered around doing numerics with Chebyshev polynomials... Sep 22, 2014 at 22:36
• I am quite unfamiliar with approximations and the like. But, let me add my "one rupee worth": the integral may be converted to one of the form $\int _0 ^1 f(x) \frac{d^n}{dx^n}(\frac{x^n(1-x)^n}{n!})$. Up to $\pm 1$ this is (by integration by parts), the same as $\int _0 ^1 f^{(n)}(x)\frac{x^n(1-x)^n}{n!}$, which may be more tractable. Sep 23, 2014 at 4:51
• @Venkataramana I upvoted, since depending on what $f$ is, this might be helpful to somebody, although this is by no mean a 'general method'.
– user54300
Sep 23, 2014 at 9:27

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.

See:

@article {MR1078802,
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},
}

• so as far as I see this, this paper enables me to evaluate the Legendre polynomials efficiently at the nodes of Chebyshev polynomials. Thus, one could some Gauss quadrature method for example. I will try this.
– user54300
Sep 22, 2014 at 13:09
• @TobiasHurth I find it hard to believe that this algorithm is not already implemented somewhere (in particular, in MATLAB). However, the only thing I could find was this: mathworks.com/matlabcentral/fileexchange/… Sep 22, 2014 at 13:18
• yes, I already looked at this code and as somebody said in the comments, it is highly(!) unaccurate for higher Legendre polynomials near $\pm 1$.
– user54300
Sep 22, 2014 at 13:39

I will make the following assumptions

• 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.

The LegendreSeries command of the OrthogonalExpansions package of Maple does the job.