Numerical Integration

One of the most widely used methods of numerical integration is Gauss-Legendre quadrature. It posses very attractive property of to be exact on polynomials of degree up to 2n-1, while using only n integrand evaluations (n-point quadrature).
The algorithm consists in approximation of initial definite integral by the sum of weighted integrand values sampled at special points called abscissas:

    \[ \begin{matrix} \displaystyle \int_a^b f(x)\,\mathrm{d}x\approx\frac{b-a}{2}\sum_{i=1}^{n}{w_i\,f\left(\frac{b-a}{2}\xi_i+\frac{b+a}{2}\right)}&\\ &\\ \displaystyle w_i = \frac{2}{\left(1-\xi_i^2\right)\,\left[P'_n(\xi_i)\right]^2}&(n=1,2,\dots) \end{matrix} \]

where values \xi_i are n zeroes of the n^{th}-degree Legendre polynomial P_n(\xi).

Accuracy of the numerical integration depends significantly on precision of zeroes \xi_i. There is no easy-to-use analytical expression for \xi_i and usually they are computed numerically by root-finding algorithms.

Surprisingly popular mathematical references/software provide low-precision \xi_i,\,w_i (only 4-6 correct digits) for limited n\le 5. Obviously such situation doesn’t reflect contemporary computational capabilities and accuracy demands of many applications. For instance, commonly used floating-point type double of IEEE 754 standard is capable to store values with 16 digits precision (machine epsilon is about 1e-16 for 64 bits double).

This page aims to provide software libraries for calculation of high-precision abscissas \xi_i and weights w_i for any desired n.

Gauss-Legendre Quadrature for C/C++

This open-source library implements numerical integration based on Gauss-Legendre quadrature of any order.

Pre-calculated high-precision abscissas and weights with 25 correct decimal places are used for specific orders n=2,…, 20, 32, 64, 96, 100, 128, 256, 512, 1024. Although all values can be found in the source code human-friendly tables with abscissas and weights are available too. Nodes for all other n are generated on the fly with 1e-10 precision by fast root-finding algorithm.

Library also includes routine for numerical integration over 2D rectangle using product of two 1D Gaussian quadratures.

If you are looking for numerical integration over the unit disk (2D sphere) you might be interested in this page Cubature formulas for the unit disk. It contains derivation details and source code in C/C++.

This C/C++ library has been included without modifications in few prominent open source packages (e.g. GNU Scientific Library) as well as in several commercial software products (financial, medical, astronomical, etc.).

Multiprecision Computing Toolbox for MATLAB worth mentioning in particular. It provides Gauss-Legendre as well as a full set of Gaussian type quadrature (Jacobi, Chebyshev, Laguerre, Hermit and Gegenbauer) in arbitrary precision. Reference tables on this page can be easily re-computed using this extension.

Gauss-Legendre Quadrature for MATLAB

There are several routines for Gaussian quadrature implemented using Matlab programming language. All of them suffer from low performance and even some of them provide low accuracy results.

To alleviate the situation I am developing quadrature toolbox which has Matlab interface (m-files) and delegates all the heavy lifting to high-performance MEX library written in C/C++ (see previous section). So far it includes only Gauss-Legendre rule with many others to come:

Quadrature toolbox consists of three files:


Toolbox was created and tested using Matlab R2010b, both on 32 and 64 bits versions. Usage example:

  ans =

GaussLegendre() accepts Matlab function (inline or from m-file), integration boundaries [a,b], order of quadrature (1024 in example) and desired precision for nodes as optional parameter. Also it returns vectors of abscissas and weights if such outputs are specified.

Performance comparison with the fastest Matlab-language routine (x-quadrature order, y-execution time, green-quadrature toolbox, red-competitor, lower is better):

Besides being fast library usually provides results with better accuracy thanks to pre-calculated high-precision abscissas and weights.

I am extending quadrature toolbox for Matlab with other rules and adaptive methods. I would really appreciate any feature suggestions and advices on what methods to include first.

Abscissae and Weights of Gauss-Legendre Quadrature

Table below lists Gauss-Legendre quadrature nodes for n=2,…, 20, 32, 64, 100 with the accuracy of 25 decimal digits. Data for higher order formulae n = 128, 256, 512, 1024 can be found in C/C++ library source code.

High-precision Abscissae and Weights of Gaussian Quadrature.
Correctly rounded to 25 decimal digits to the nearest.
nAbscissae \xi_iWeights w_i

1 Star2 Stars3 Stars4 Stars5 Stars (21 votes, average: 5.00)


  1. Ben
    Posted June 17, 2011 at 6:47 pm | #

    It’s worth looking at eigenvalue methods to calculate the nodes and weights. Basically the 3 term polynomial recurrence relation allows the roots to be expressed as eigenvalues of a tridiagonal matrix; with appropriate rescaling the matrix is transformed to symmetric tridiagonal and hence the roots can be determined to arbitrary precision e.g. with your MPFR library. Similarly the weights can be expressed in terms of the eigenvectors.

    • Posted June 18, 2011 at 12:59 pm | #

      Derivation of quadrature nodes using eigenvalues/vectors of Jacobi matrix is universal approach.

      Although it is very powerful I find it inefficient for particular well studied case – Gauss-Legendre quadrature.
      We have nice analytic formula for close approximation of Legendre polynomial roots. Starting from there Newton iterations converge to “real” root very quickly with any desired precision. I think even Golub-Welsch eigen solver optimized for symmetric 3-diagonal matrices cannot beat that.

      BTW I’ve used eigen method to derive quadrature for less-known weight function \omega(x)=|x| here: Cubature formulas for the unit disk.

      Thanks for reading pages on my site ;-).

  2. Gowri Priya
    Posted July 14, 2011 at 4:36 pm | #

    Thank u very much for giving useful information

  3. vladan
    Posted August 4, 2011 at 10:12 pm | #

    Do you know where I can find “Abscissae and Weights of Gauss-Legendre Quadrature” for “Associated Legendre”?

    • Posted August 5, 2011 at 9:27 am | #

      Do you mean quadrature based on “Associated Legendre” polynomials? I don’t know this kind of quadrature.
      Could you be more specific – what task do you solve ?

  4. ankit
    Posted October 12, 2011 at 10:33 pm | #

    excellent work….helped me alot..

  5. Posted October 24, 2011 at 5:38 pm | #

    show the difference between langrange polynomial and newton interpolating polynomial

  6. Sergey
    Posted November 15, 2011 at 5:26 pm | #

    Brilliant job!!! Thank you

  7. Ali
    Posted January 14, 2012 at 12:35 am | #

    nice job,but I really wants to know how to use it for double integration,if some how it’s doable please let me know?

  8. Stephen
    Posted January 17, 2012 at 3:32 pm | #

    I came across your site while searching for Gauss-Laguerre quadrature. Do you plan to include this in your toolbox? I am currently using the eigen method, but I expect Newton’s method with the initial approximations from Stroud and Secrest would be faster.

    • Posted January 17, 2012 at 4:47 pm | #

      I completely agree with you, direct root search with good starting point is much faster than eigen-solving Jacobi matrix.

      I have plans to cover all commonly used quadrature, plus to port QuadPack adaptive integration rules to Matlab. I will release them as a part of Multiprecision Computing Toolbox for MATLAB.

      It already has Gauss-Legendre and Gauss-Kronrod rules which is the most important for practical applications since it allows re-use of lower order quadrature nodes in higher order rules. Thus saving on function evaluations.

  9. Pablo
    Posted February 17, 2012 at 2:37 am | #

    Hi Pavel!

    First of all congrats for your excellent job.

    I’m currently using your code to integrate the elements of the stiffness matrix in a solver for Poisson Eq using weak form of Finite Elements Method.

    The point is that the code works well when I try to evaluate numerically the integral of a function in ‘x’, but I don’t manage to integrate a constant (i.e. imagine we want to do gausslegendre(@(x) 1,-1,1). The result of this operation is, as far as I’m concerned, random.

    I could neither do this with the original version of quad.m, but after performing some changes to the code now I’m able to do it with a modified ver called myquad.m

    Could it be possible to perform the integration of a constant using your code?

    Many thanks and congrats again for your work!

    Best regards from El Salvador


  10. Posted February 28, 2012 at 6:47 pm | #

    Hi Pavel,

    I am looking for a way to get forces from pressure distribution on a plate and this looks like exactly what i need.

    Many Thanks,

    Balaji Sankar.

  11. isen1988
    Posted October 26, 2012 at 3:18 am | #

    Thank god! I find you!My project takes much time doing integration. So, I want to write a integration mex file to minus time,
    But I just can’t find out a way to pass in the function pointer to the mex file. That has really puzzled me a lot. How do you solve this ?

  12. Grunde Waag
    Posted November 30, 2012 at 12:22 am | #

    Nice pages.

    I’m integrating some very spiky functions. They are almost zero most of the time and contains one or more spikes. What is the best scheme to use for such a situation? With the trapezoidal rule I have to sample very dense for the integration to converge. Higher order schemes doesn’t seem to improve the situation. Any advice?

    • Posted November 30, 2012 at 11:49 am | #

      You could try adaptive integration methods based on Gauss-Kronrod rules – it will isolate “bad” regions and sample it with finer grid.
      Or look up methods for integrating “oscillatory functions”.

  13. Amao
    Posted December 5, 2012 at 10:26 pm | #

    Hi Pavel,

    many thanks for the ultrafast code. However, i found the matlab code can not treat complex function. Do you have a solution to it? Thanks in advance.


    • Posted December 6, 2012 at 1:34 pm | #

      Hi, you can re-formulate integration of complex function in real domain.

  14. Tejas
    Posted January 10, 2013 at 9:33 am | #


    thanks for the neat explanation of the Gauss-Legendre quadrature rule!
    I have 2 questions for you and would really appreciate if you could help in any way

    1) Say we have two functions f(x) and h(x) defined over the same interval (a,b) and we wish to integrate the product f(x)*h(x), does the coefficient term (Jacobi) appear only once in this case?

    2) now if we have f(x) over (a,b) and h(x) over (a,c) [ also lets assume ab > ac ], can we integrate the product f(x)*h(x) over (a,c) [ since ac is within ab] using Gauss-Legendre quadrature rule? if yes, would the Jacobi be (ac -> zeta mapping) or some combination of (ac -> zeta; ab-> zeta mapping)?


  15. Robert
    Posted March 7, 2013 at 5:12 pm | #


    Thanks for making these procedures available.

    I need the mex files for Linux (mexa64), so I haven’t been able to them. Is it possible for you to share the source code so I can compile the files myself or add this type of mex file?


    • Lionel
      Posted August 7, 2013 at 1:28 am | #

      Hi all,

      I’m also working on linux 64 and i facing exactly the same problem. Would it be possible for you to share them? (mexa64 or source code to compile it myself).

      Thank you for publishing these codes,


  16. David L
    Posted April 29, 2013 at 5:05 pm | #

    Thanks for the mex implementation.. its enormous fast!
    Nevertheless I suspect a memory leak (at least in) quadlab.mexw64.
    When you call the function several times in Matlab, it seems that the used memory increases.
    Try the following:
    for i =1:100000, quadgk(@sinc,0,1);end
    for i =1:100000, GaussLegendre(@sinc,0,1,1024);end
    for i =1:100000, quadgk(@sinc,0,1);end
    for i =1:100000, GaussLegendre(@sinc,0,1,1024);end
    % check memory in task manager
    for i =1:100000, quadgk(@sinc,0,1);end
    % check memory in task manager
    for i =1:100000, GaussLegendre(@sinc,0,1,1024);end
    % check memory in task manager

    After the first 2 calls of the for loops the memory is increasing for both integration methods. Later the memory stays constant for the quadgk command but keeps increasing with the GaussLegendre command.


    • Posted June 22, 2014 at 9:31 pm | #

      Thank you for the report. Memory leak is now fixed – please re-download the library.

  17. annel
    Posted September 10, 2013 at 3:50 pm | #

    Great code, thanks !


  18. Posted June 22, 2014 at 10:30 am | #

    Your Gauss-Legendere MATLAB code has speeded up a computational economics project I am working on: really nice work! But, I have two questions that are similar to other posters:

    i. would it be possible to get the source code to create .mexa64 files for LINUX/UNIX?
    ii. is it expected that the memory usage can get so high when quadlab is called? I have tried doing some parallel implementation on a 96GB windows workstation and all of the memory is being used up.



    • Posted June 22, 2014 at 9:49 pm | #

      Thank you for your feedback!

      Memory leak is now fixed. I have just released the new version – please re-download the lib (using the same link above).
      It still has only Windows binaries though.

      I am not ready to go open source, but Linux version is in my todo list.
      (It requires more time – I have to remove some Windows-specifics first).

      Let me know how is it working now,

      Will keep you updated on Linux version.

      • Andrew
        Posted June 23, 2014 at 10:43 pm | #

        I tried the new mex binaries. Memory usage is much more stable, but it runs much more slowly (by a factor of 4). Answers are identical up to 1e-15 on my example, so I may switch between the old and new ones depending on how I am using them.


        • Posted June 24, 2014 at 12:24 am | #

          Now fixed and should run even faster than the very first version (thanks to Intel C++). Please re-download & test the library.

  19. Posted July 9, 2014 at 5:02 am | #

    I think the scientific way to write a code, which is not filling it with numbers ,a much better way if you write a short piece of code to regenerate those weights and roots instead of putting them like this. Moreover, you will reduce the error that comes from human

    • Posted July 9, 2014 at 6:10 pm | #

      Pre-computed coefficients are needed for high speed. Re-generation of high-precision coefficients every time would take forever.

      • Posted July 10, 2014 at 2:03 pm | #

        Thanks for this useful comment, you are absolutely right. I would say if you put a commented code, since you already have it, for those people who don’t worry about the speed but the accuracy. Such a code make your post wider and much helpful popular .

        • Posted July 10, 2014 at 2:20 pm | #

          Library includes the code to generate nodes & weights for any required order (in double precision).
          Only few particular cases are pre-computed, all others are generated on the fly.

          If you are interested in generating coefficients in arbitrary precision, then check out this MATLAB code:
          Arbitrary Precision Gauss-Legendre Quadrature

          It requires multiprecision toolbox for MATLAB, which includes computation of nodes & weights of all common Gauss-type quadrature in arbitrary precision. See here:
          Abscissas and Weights of Classical Gaussian Quadrature Rules

          Btw, I am author of the toolbox :).

          • Posted July 11, 2014 at 2:16 pm | #

            Thanks Pavel . It is good to know that I am dealing with professional people 🙂

            I already coded Gauss quadrature rule using C++, and every time I regenerate the weights and the roots using a little piece of code, but when I need to find the weights and the roots where N=64 then it takes so long time. does that what happen with you when you generate such roots and weights.

          • Posted July 11, 2014 at 3:48 pm | #

            It as always better to do a little research before coding. Just count number of operation (algorithmic complexity). Or learn how it is done in other libraries.

  20. Simone
    Posted April 4, 2015 at 11:54 pm | #

    Thanks for the code, it works perfectly! However, is it possible to use the Gauss-Legendre integration also for double integrals?

    • labd1
      Posted October 2, 2015 at 11:31 pm | #

      How is it possible to implement this code for 2D integration please?

  21. Ali
    Posted October 14, 2015 at 3:35 am | #

    First, I thank you for your useful mex function. It helped me a lot with my thesis.
    One problem though: I need to integrate complex functions. Since your function does not cover complex variables, I integrate once imag(f(x)) and once real(f(x)). This doubles the number of integrals I have to compute and integration is a serious bottleneck in my programs.
    I would appreciate any suggestions.

  22. Jay
    Posted April 12, 2016 at 10:28 pm | #

    Hi. I was really impressed by your Gauss-Legendre Quadrature for C/C++.
    I’ve intersted, is the ability to estimate of the absolute error in the result. like it implements here (abserr) ? Thanks!

    • Posted April 13, 2016 at 1:44 pm | #

      You can run two quadrature of different orders and difference in their results will give you the absolute error.

      • Jay
        Posted April 18, 2016 at 9:38 pm | #

        I guess that for such approach there is need
        I think that different orders is not enough, different grid is also needed.

        • Posted April 18, 2016 at 9:46 pm | #

          Yes, of course, Gauss-Kronrod or Curtis-Clenshaw are the best methods for such purposes.
          I thought you are restricted to use Gauss-Legendre.

          • Jay
            Posted April 19, 2016 at 7:50 pm | #

            I mean that that could be good addon to your lib;)

          • Posted April 20, 2016 at 11:00 am | #

            Absolutely – been thinking about it for some time. Also all Gauss-type quadrature, multi-dimensional cubature for various types of domains (disk, hypercube, simplex, etc.), Monte-Carlo integration, etc. The only question is where to get the time & $ :).

  23. Justin
    Posted August 20, 2016 at 11:09 pm | #

    Hi Pavel, thanks for hosting such an informative thread. In situations where we integrate over the unit sphere, the basis functions become the Associated Legendre polynomials P_{n}^{m}(\theta), which certainly have unique \xi_i and w_i compared to P_n^0. In this case, it is unclear to me how the algorithm should be changed. And advice is appreciated!

  24. Posted April 28, 2024 at 8:57 pm | #

    Hi Pavel, thank you for this amazing work. Why do you not load this code in a GitHub Repository? I’m using your lib for my robotic application and it could be very useful to share this code in GitHub. If you don’t want to do that, I can publish it as a public repository, referencing this page.
    Thank you again for this library!

Post a Comment

Your email is never published nor shared.

Use native LaTeX syntax to include formulas: $ ... $, \[ ... \], etc. Do not forget to preview comment before posting.

Also you may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

Subscribe without commenting