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

  • 1
    $\begingroup$ Are you aware of anything clever for Fourier series (which suffer from the same oscillation problem)? $\endgroup$
    – Igor Rivin
    Sep 22, 2014 at 11:39
  • $\begingroup$ @IgorRivin unfortunately not, sorry. $\endgroup$
    – user54300
    Sep 22, 2014 at 11:46
  • $\begingroup$ See if you can borrow some ideas in chebfun which is centered around doing numerics with Chebyshev polynomials... $\endgroup$
    – Suvrit
    Sep 22, 2014 at 22:36
  • 1
    $\begingroup$ 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. $\endgroup$ Sep 23, 2014 at 4:51
  • $\begingroup$ @Venkataramana I upvoted, since depending on what $f$ is, this might be helpful to somebody, although this is by no mean a 'general method'. $\endgroup$
    – user54300
    Sep 23, 2014 at 9:27

4 Answers 4


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.



@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},
  • $\begingroup$ 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. $\endgroup$
    – user54300
    Sep 22, 2014 at 13:09
  • $\begingroup$ @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/… $\endgroup$
    – Igor Rivin
    Sep 22, 2014 at 13:18
  • $\begingroup$ 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$. $\endgroup$
    – user54300
    Sep 22, 2014 at 13:39

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.


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy