Pavel Holoborodko Numerical Integration

Numerical Integration

Posted By Pavel Holoborodko On October 29, 2008 @ 5:46 pm In | 12 Comments

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 to mention separately. It provides multiple precision version of the library which allows calculation of Gauss-Legendre quadrature with arbitrary accuracy for any order. In particular, 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

12 Comments (Open | Close)

12 Comments To "Numerical Integration"

#1 Comment By Ben On June 17, 2011 @ 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.

#2 Comment By Pavel Holoborodko On June 18, 2011 @ 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 .

#3 Comment By Gowri Priya On July 14, 2011 @ 4:36 pm

Thank u very much for giving useful information

#4 Comment By vladan On August 4, 2011 @ 10:12 pm

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

#5 Comment By Pavel Holoborodko On August 5, 2011 @ 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 ?

#6 Comment By ankit On October 12, 2011 @ 10:33 pm

excellent work….helped me alot..

#7 Comment By edu On October 24, 2011 @ 5:38 pm

show the difference between langrange polynomial and newton interpolating polynomial

#8 Comment By Sergey On November 15, 2011 @ 5:26 pm

Brilliant job!!! Thank you

#9 Comment By Ali On January 14, 2012 @ 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?

#10 Comment By Stephen On January 17, 2012 @ 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.

#11 Comment By Pavel Holoborodko On January 17, 2012 @ 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-Krondrod 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.

#12 Comment By Pablo On February 17, 2012 @ 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

Article printed from Pavel Holoborodko: http://www.holoborodko.com/pavel

URL to article: http://www.holoborodko.com/pavel/numerical-methods/numerical-integration/

Click here to print.