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.

Gauss-Legendre Quadrature.zip

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:

Gauss-Legendre Quadrature for Matlab.zip

Quadrature toolbox consists of three files:

GaussLegendre.m 
quadlab.mexw32 
quadlab.mexw64

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

>>GaussLegendre(@cos,-pi/2,pi/2,1024)    
  ans =
       2.000000000000000

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.
$n$	Abscissae $\xi_i$	Weights $w_i$
2	±0.5773502691896257645091488	1.0000000000000000000000000
3	0	0.8888888888888888888888889
3	±0.7745966692414833770358531	0.5555555555555555555555556
4	±0.3399810435848562648026658	0.6521451548625461426269361
4	±0.8611363115940525752239465	0.3478548451374538573730639
5	0	0.5688888888888888888888889
	±0.5384693101056830910363144	0.4786286704993664680412915
	±0.9061798459386639927976269	0.2369268850561890875142640
6	±0.2386191860831969086305017	0.4679139345726910473898703
	±0.6612093864662645136613996	0.3607615730481386075698335
	±0.9324695142031520278123016	0.1713244923791703450402961
7	0	0.4179591836734693877551020
	±0.4058451513773971669066064	0.3818300505051189449503698
	±0.7415311855993944398638648	0.2797053914892766679014678
	±0.9491079123427585245261897	0.1294849661688696932706114
8	±0.1834346424956498049394761	0.3626837833783619829651504
	±0.5255324099163289858177390	0.3137066458778872873379622
	±0.7966664774136267395915539	0.2223810344533744705443560
	±0.9602898564975362316835609	0.1012285362903762591525314
9	0	0.3302393550012597631645251
	±0.3242534234038089290385380	0.3123470770400028400686304
	±0.6133714327005903973087020	0.2606106964029354623187429
	±0.8360311073266357942994298	0.1806481606948574040584720
	±0.9681602395076260898355762	0.0812743883615744119718922
10	±0.1488743389816312108848260	0.2955242247147528701738930
	±0.4333953941292471907992659	0.2692667193099963550912269
	±0.6794095682990244062343274	0.2190863625159820439955349
	±0.8650633666889845107320967	0.1494513491505805931457763
	±0.9739065285171717200779640	0.0666713443086881375935688
11	0	0.2729250867779006307144835
	±0.2695431559523449723315320	0.2628045445102466621806889
	±0.5190961292068118159257257	0.2331937645919904799185237
	±0.7301520055740493240934163	0.1862902109277342514260976
	±0.8870625997680952990751578	0.1255803694649046246346943
	±0.9782286581460569928039380	0.0556685671161736664827537
12	±0.1252334085114689154724414	0.2491470458134027850005624
	±0.3678314989981801937526915	0.2334925365383548087608499
	±0.5873179542866174472967024	0.2031674267230659217490645
	±0.7699026741943046870368938	0.1600783285433462263346525
	±0.9041172563704748566784659	0.1069393259953184309602547
	±0.9815606342467192506905491	0.0471753363865118271946160
13	0	0.2325515532308739101945895
	±0.2304583159551347940655281	0.2262831802628972384120902
	±0.4484927510364468528779129	0.2078160475368885023125232
	±0.6423493394403402206439846	0.1781459807619457382800467
	±0.8015780907333099127942065	0.1388735102197872384636018
	±0.9175983992229779652065478	0.0921214998377284479144218
	±0.9841830547185881494728294	0.0404840047653158795200216
14	±0.1080549487073436620662447	0.2152638534631577901958764
	±0.3191123689278897604356718	0.2051984637212956039659241
	±0.5152486363581540919652907	0.1855383974779378137417166
	±0.6872929048116854701480198	0.1572031671581935345696019
	±0.8272013150697649931897947	0.1215185706879031846894148
	±0.9284348836635735173363911	0.0801580871597602098056333
	±0.9862838086968123388415973	0.0351194603317518630318329
15	0	0.2025782419255612728806202
	±0.2011940939974345223006283	0.1984314853271115764561183
	±0.3941513470775633698972074	0.1861610000155622110268006
	±0.5709721726085388475372267	0.1662692058169939335532009
	±0.7244177313601700474161861	0.1395706779261543144478048
	±0.8482065834104272162006483	0.1071592204671719350118695
	±0.9372733924007059043077589	0.0703660474881081247092674
	±0.9879925180204854284895657	0.0307532419961172683546284
16	±0.0950125098376374401853193	0.1894506104550684962853967
	±0.2816035507792589132304605	0.1826034150449235888667637
	±0.4580167776572273863424194	0.1691565193950025381893121
	±0.6178762444026437484466718	0.1495959888165767320815017
	±0.7554044083550030338951012	0.1246289712555338720524763
	±0.8656312023878317438804679	0.0951585116824927848099251
	±0.9445750230732325760779884	0.0622535239386478928628438
	±0.9894009349916499325961542	0.0271524594117540948517806
17	0	0.1794464703562065254582656
	±0.1784841814958478558506775	0.1765627053669926463252710
	±0.3512317634538763152971855	0.1680041021564500445099707
	±0.5126905370864769678862466	0.1540457610768102880814316
	±0.6576711592166907658503022	0.1351363684685254732863200
	±0.7815140038968014069252301	0.1118838471934039710947884
	±0.8802391537269859021229557	0.0850361483171791808835354
	±0.9506755217687677612227170	0.0554595293739872011294402
	±0.9905754753144173356754340	0.0241483028685479319601100
18	±0.0847750130417353012422619	0.1691423829631435918406565
	±0.2518862256915055095889729	0.1642764837458327229860538
	±0.4117511614628426460359318	0.1546846751262652449254180
	±0.5597708310739475346078715	0.1406429146706506512047313
	±0.6916870430603532078748911	0.1225552067114784601845191
	±0.8037049589725231156824175	0.1009420441062871655628140
	±0.8926024664975557392060606	0.0764257302548890565291297
	±0.9558239495713977551811959	0.0497145488949697964533349
	±0.9915651684209309467300160	0.0216160135264833103133427
19	0	0.1610544498487836959791636
	±0.1603586456402253758680961	0.1589688433939543476499564
	±0.3165640999636298319901173	0.1527660420658596667788554
	±0.4645707413759609457172671	0.1426067021736066117757461
	±0.6005453046616810234696382	0.1287539625393362276755158
	±0.7209661773352293786170959	0.1115666455473339947160239
	±0.8227146565371428249789225	0.0914900216224499994644621
	±0.9031559036148179016426609	0.0690445427376412265807083
	±0.9602081521348300308527788	0.0448142267656996003328382
	±0.9924068438435844031890177	0.0194617882297264770363120
20	±0.0765265211334973337546404	0.1527533871307258506980843
	±0.2277858511416450780804962	0.1491729864726037467878287
	±0.3737060887154195606725482	0.1420961093183820513292983
	±0.5108670019508270980043641	0.1316886384491766268984945
	±0.6360536807265150254528367	0.1181945319615184173123774
	±0.7463319064601507926143051	0.1019301198172404350367501
	±0.8391169718222188233945291	0.0832767415767047487247581
	±0.9122344282513259058677524	0.0626720483341090635695065
	±0.9639719272779137912676661	0.0406014298003869413310400
	±0.9931285991850949247861224	0.0176140071391521183118620
32	±0.0483076656877383162348126	0.0965400885147278005667648
	±0.1444719615827964934851864	0.0956387200792748594190820
	±0.2392873622521370745446032	0.0938443990808045656391802
	±0.3318686022821276497799168	0.0911738786957638847128686
	±0.4213512761306353453641194	0.0876520930044038111427715
	±0.5068999089322293900237475	0.0833119242269467552221991
	±0.5877157572407623290407455	0.0781938957870703064717409
	±0.6630442669302152009751152	0.0723457941088485062253994
	±0.7321821187402896803874267	0.0658222227763618468376501
	±0.7944837959679424069630973	0.0586840934785355471452836
	±0.8493676137325699701336930	0.0509980592623761761961632
	±0.8963211557660521239653072	0.0428358980222266806568786
	±0.9349060759377396891709191	0.0342738629130214331026877
	±0.9647622555875064307738119	0.0253920653092620594557526
	±0.9856115115452683354001750	0.0162743947309056706051706
	±0.9972638618494815635449811	0.0070186100094700966004071
64	±0.0243502926634244325089558	0.0486909570091397203833654
	±0.0729931217877990394495429	0.0485754674415034269347991
	±0.1214628192961205544703765	0.0483447622348029571697695
	±0.1696444204239928180373136	0.0479993885964583077281262
	±0.2174236437400070841496487	0.0475401657148303086622822
	±0.2646871622087674163739642	0.0469681828162100173253263
	±0.3113228719902109561575127	0.0462847965813144172959532
	±0.3572201583376681159504426	0.0454916279274181444797710
	±0.4022701579639916036957668	0.0445905581637565630601347
	±0.4463660172534640879849477	0.0435837245293234533768279
	±0.4894031457070529574785263	0.0424735151236535890073398
	±0.5312794640198945456580139	0.0412625632426235286101563
	±0.5718956462026340342838781	0.0399537411327203413866569
	±0.6111553551723932502488530	0.0385501531786156291289625
	±0.6489654712546573398577612	0.0370551285402400460404151
	±0.6852363130542332425635584	0.0354722132568823838106931
	±0.7198818501716108268489402	0.0338051618371416093915655
	±0.7528199072605318966118638	0.0320579283548515535854675
	±0.7839723589433414076102205	0.0302346570724024788679741
	±0.8132653151227975597419233	0.0283396726142594832275113
	±0.8406292962525803627516915	0.0263774697150546586716918
	±0.8659993981540928197607834	0.0243527025687108733381776
	±0.8893154459951141058534040	0.0222701738083832541592983
	±0.9105221370785028057563807	0.0201348231535302093723403
	±0.9295691721319395758214902	0.0179517157756973430850453
	±0.9464113748584028160624815	0.0157260304760247193219660
	±0.9610087996520537189186141	0.0134630478967186425980608
	±0.9733268277899109637418535	0.0111681394601311288185905
	±0.9833362538846259569312993	0.0088467598263639477230309
	±0.9910133714767443207393824	0.0065044579689783628561174
	±0.9963401167719552793469245	0.0041470332605624676352875
	±0.9993050417357721394569056	0.0017832807216964329472961
100	±0.0156289844215430828722167	0.0312554234538633569476425
	±0.0468716824215916316149239	0.0312248842548493577323765
	±0.0780685828134366366948174	0.0311638356962099067838183
	±0.1091892035800611150034260	0.0310723374275665165878102
	±0.1402031372361139732075146	0.0309504788504909882340635
	±0.1710800805386032748875324	0.0307983790311525904277139
	±0.2017898640957359972360489	0.0306161865839804484964594
	±0.2323024818449739696495100	0.0304040795264548200165079
	±0.2625881203715034791689293	0.0301622651051691449190687
	±0.2926171880384719647375559	0.0298909795933328309168368
	±0.3223603439005291517224766	0.0295904880599126425117545
	±0.3517885263724217209723438	0.0292610841106382766201190
	±0.3808729816246299567633625	0.0289030896011252031348762
	±0.4095852916783015425288684	0.0285168543223950979909368
	±0.4378974021720315131089780	0.0281027556591011733176483
	±0.4657816497733580422492166	0.0276611982207923882942042
	±0.4932107892081909335693088	0.0271926134465768801364916
	±0.5201580198817630566468157	0.0266974591835709626603847
	±0.5465970120650941674679943	0.0261762192395456763423087
	±0.5725019326213811913168704	0.0256294029102081160756420
	±0.5978474702471787212648065	0.0250575444815795897037642
	±0.6226088602037077716041908	0.0244612027079570527199750
	±0.6467619085141292798326303	0.0238409602659682059625604
	±0.6702830156031410158025870	0.0231974231852541216224889
	±0.6931491993558019659486479	0.0225312202563362727017970
	±0.7153381175730564464599671	0.0218430024162473863139537
	±0.7368280898020207055124277	0.0211334421125276415426723
	±0.7575981185197071760356680	0.0204032326462094327668389
	±0.7776279096494954756275514	0.0196530874944353058653815
	±0.7968978923903144763895729	0.0188837396133749045529412
	±0.8153892383391762543939888	0.0180959407221281166643908
	±0.8330838798884008235429158	0.0172904605683235824393442
	±0.8499645278795912842933626	0.0164680861761452126431050
	±0.8660146884971646234107400	0.0156296210775460027239369
	±0.8812186793850184155733168	0.0147758845274413017688800
	±0.8955616449707269866985210	0.0139077107037187726879541
	±0.9090295709825296904671263	0.0130259478929715422855586
	±0.9216092981453339526669513	0.0121314576629794974077448
	±0.9332885350430795459243337	0.0112251140231859771172216
	±0.9440558701362559779627747	0.0103078025748689695857821
	±0.9539007829254917428493369	0.0093804196536944579514182
	±0.9628136542558155272936593	0.0084438714696689714026208
	±0.9707857757637063319308979	0.0074990732554647115788287
	±0.9778093584869182885537811	0.0065469484508453227641521
	±0.9838775407060570154961002	0.0055884280038655151572119
	±0.9889843952429917480044187	0.0046244500634221193510958
	±0.9931249370374434596520099	0.0036559612013263751823425
	±0.9962951347331251491861317	0.0026839253715534824194396
	±0.9984919506395958184001634	0.0017093926535181052395294
	±0.9997137267734412336782285	0.0007346344905056717304063

(21 votes, average: 5.00)

46 Comments

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.
Reply
- Pavel Holoborodko
  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.
  P.S.
  BTW I’ve used eigen method to derive quadrature for less-known weight function $\omega(x)=|x|$ here: Cubature formulas for the unit disk.
  P.P.S.
  Thanks for reading pages on my site ;-).
  Reply
Gowri Priya
Posted July 14, 2011 at 4:36 pm | #
Thank u very much for giving useful information
Reply
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”?
Reply
- Pavel Holoborodko
  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 ?
  Reply
ankit
Posted October 12, 2011 at 10:33 pm | #
excellent work….helped me alot..
Reply
edu
Posted October 24, 2011 at 5:38 pm | #
show the difference between langrange polynomial and newton interpolating polynomial
Reply
Sergey
Posted November 15, 2011 at 5:26 pm | #
Brilliant job!!! Thank you
Reply
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?
Reply
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.
Reply
- Pavel Holoborodko
  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.
  Reply
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
Pablo
Reply
Balaji Sankar
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.
Reply
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 ?
Reply
Grunde Waag
Posted November 30, 2012 at 12:22 am | #
Hi,
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?
Reply
- Pavel Holoborodko
  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”.
  Reply
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.
Amao
Reply
- Pavel Holoborodko
  Posted December 6, 2012 at 1:34 pm | #
  Hi, you can re-formulate integration of complex function in real domain.
  Reply
Tejas
Posted January 10, 2013 at 9:33 am | #
Hey,
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)?
thanks!
Reply
Robert
Posted March 7, 2013 at 5:12 pm | #
Hi
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?
Thanks!
Reply
- 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,
  Lionel
  Reply
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.
Greets,
David
Reply
- Pavel Holoborodko
  Posted June 22, 2014 at 9:31 pm | #
  Thank you for the report. Memory leak is now fixed – please re-download the library.
  Reply
annel
Posted September 10, 2013 at 3:50 pm | #
Great code, thanks !
Anne
Reply
Andrew
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.
Thanks
Andrew
Reply
- Pavel Holoborodko
  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.
  Reply
  - 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.
    Thanks
    Andrew
    Reply
    - Pavel Holoborodko
      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.
      Reply
Ahmad
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
Cheers.
Reply
- Pavel Holoborodko
  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.
  Reply
  - Ahmad
    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 .
    Reply
    - Pavel Holoborodko
      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 :).
      Reply
      - Ahmad
        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.
        Thanks.
      - Pavel Holoborodko
        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.
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?
Reply
- labd1
  Posted October 2, 2015 at 11:31 pm | #
  How is it possible to implement this code for 2D integration please?
  Reply
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.
Reply
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 http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html (abserr) ? Thanks!
Reply
- Pavel Holoborodko
  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.
  Reply
  - Jay
    Posted April 18, 2016 at 9:38 pm | #
    I guess that for such approach there is need https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Kronrod_rules.
    I think that different orders is not enough, different grid is also needed.
    Reply
    - Pavel Holoborodko
      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.
      Reply
      - Jay
        Posted April 19, 2016 at 7:50 pm | #
        I mean that that could be good addon to your lib;)
      - Pavel Holoborodko
        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 & $ :).
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!
Reply
- Pavel Holoborodko
  Posted August 20, 2016 at 11:32 pm | #
  I don’t know much about the sphere, except that there is Lebedev quadrature: https://en.wikipedia.org/wiki/Lebedev_quadrature
  And some code for it: https://people.sc.fsu.edu/~jburkardt/f_src/sphere_lebedev_rule/sphere_lebedev_rule.html
  Reply
Daniele Caradonna
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!
Reply

Numerical Integration

Gauss-Legendre Quadrature for C/C++

Gauss-Legendre Quadrature for MATLAB

Abscissae and Weights of Gauss-Legendre Quadrature

46 Comments

Post a Comment