The Gauss–Kronrod quadrature[1] includes two distinct sets of abscissae – Gauss nodes interleaved with Kronrod nodes. The special geometry of this system enables very efficient error estimation. In the beginning, the integral is approximated using Gauss quadrature. Then, the already evaluated function values are reused to produce higher order approximation, in combination with new values sampled at the Kronrod nodes. The difference between the two rules is a good indicator of how close the estimation is to the actual behavior of the integral. This technique is normally used as the stopping criteria in many adaptive numerical integration methods.

Here we focus on the Kronrod nodes computing in arbitrary precision, while using an efficient algorithm previously proposed by Laurie[2]. (See [3], [4] for other methods).

Existing code porting to arbitrary precision

The MATLAB implementation of the Laurie method can be found in OPQ[5] – a companion set of routines used in the outstanding monograph “Orthogonal Polynomials: Computation and Approximation” by Walter Gautschi[6]. Our goal is to enable the existing MATLAB program to conduct these calculations in arbitrary precision.

The Laurie algorithm is implemented in the following four files: r_kronrod.m, kronrod.m, r_jacobi01.m and r_jacobi.m. Examples of the usage of these functions can be found here [6, p. 171]:

ab=r_jacobi01(61);
for N=[5 10 25 40]
   xw=kronrod(N,ab)      % Kronrod abscissae & weights for [0, 1]
end

After a thorough examination of the r_jacobi01.m and r_jacobi.m, we conclude that they are prepared for multiprecision support without any additional changes – we only need to supply the multiprecision numbers as parameters of the functions. After these are input, the toolbox takes care of the rest. For example, the r_jacobi.m code is shown below:

function ab=r_jacobi(N,a,b)
if nargin<2, a=0; end;  if nargin<3, b=a; end
if((N<=0)|(a<=-1)|(b<=-1)) error('parameter(s) out of range'), end
nu=(b-a)/(a+b+2);
mu=2^(a+b+1)*gamma(a+1)*gamma(b+1)/gamma(a+b+2);
if N==1, ab=[nu mu]; return, end 
N=N-1; n=1:N; nab=2*n+a+b;
A=[nu (b^2-a^2)*ones(1,N)./(nab.*(nab+2))];
n=2:N; nab=nab(n);
B1=4*(a+1)*(b+1)/((a+b+2)^2*(a+b+3));
B=4*(n+a).*(n+b).*n.*(n+a+b)./((nab.^2).*(nab+1).*(nab-1));
ab=[A' [mu; B1; B']];

If the input parameters a and b are of the multiprecision type, then all of the subsequent calculations inside the function will be conducted with high precision, thanks to the overloaded functions provided by the toolbox.

The core functions require only small cosmetic changes:

% r_kronrod.m:
a=mp(zeros(2*N+1,1)); b=a;						
...
s=mp(zeros(floor(N/2)+2,1)); t=s; t(2)=b(N+2);
 
% kronrod.m:
J=mp(zeros(2*N+1));

We simply wrap zeros with the mp() constructor to convert the zero-matrices into arbitrary precision. This completes the process of porting existing code to multiprecision.

Numerical Results

With the Multiprecision Computing Toolbox, we wish to accurately compute the Kronrod nodes of any order N on the interval [-1,1]. Below is an additional driver routine that helps us accomplish this task:

function xw=mpkronrod(N)
 ab=r_jacobi01(2*N,mp(0),mp(0));
 t=kronrod(N,ab);
 xw(:,1) = 2*t(:,1)-1;      % translate abscissae to [-1, 1]
 xw(:,2) = 2*t(:,2);        % translate weights to [-1, 1]
end

Now, the Kronrod nodes can be computed with any required degree of accuracy (e.g. 25 digits):

>> mp.Digits(25);
>> mp.FollowMatlabNumericFormat(true);     % Abide the MATLAB output format specification
 
>> format longE
 
>> xw=mpkronrod(10)
xw = 
-0.9956571630258080807355273     0.0116946388673718742780644    
-0.9739065285171717200779640     0.0325581623079647274788190    % *
-0.9301574913557082260012072     0.0547558965743519960313813    
-0.8650633666889845107320967     0.0750396748109199527670431    % *
-0.7808177265864168970637176     0.0931254545836976055350655    
-0.6794095682990244062343274     0.1093871588022976418992106    % *
-0.5627571346686046833390001     0.1234919762620658510779581    
-0.4333953941292471907992659     0.1347092173114733259280540    % *
-0.2943928627014601981311266     0.1427759385770600807970943    
-0.1488743389816312108848260     0.1477391049013384913748415    % *
 0.0000000000000000000000000     0.1494455540029169056649365    
 0.1488743389816312108848260     0.1477391049013384913748415    % *
 0.2943928627014601981311266     0.1427759385770600807970943    
 0.4333953941292471907992659     0.1347092173114733259280540    % *
 0.5627571346686046833390001     0.1234919762620658510779581    
 0.6794095682990244062343274     0.1093871588022976418992106    % *
 0.7808177265864168970637176     0.0931254545836976055350655    
 0.8650633666889845107320967     0.0750396748109199527670431    % *
 0.9301574913557082260012072     0.0547558965743519960313813    
 0.9739065285171717200779640     0.0325581623079647274788190    % *
 0.9956571630258080807355273     0.0116946388673718742780644

The newly computed values by toolbox coincide with publicly available high precision tables listed in QUADPACK[7], and with original double-Laurie algorithm as well. The calculated Gauss nodes (marked with *) are correct in all 25 digits (e.g. compare with the High precision abscissae and weights of Gauss-Legendre quadrature).

Required accuracy can be (reasonably) high:

>> mp.Digits(300);
>> tic; xw300=mpkronrod(10); toc;
Elapsed time is 0.436994 seconds.
 
>> mp.Digits(350);
>> tic; xw350=mpkronrod(10); toc;
Elapsed time is 0.498385 seconds.

Also, it is important that the first digits of the higher precision values coincide with all of the digits of the lower precision values (up to mp-machine epsilon):

>> norm(xw300-xw350,1)
ans = 
     1.39157....5665e-303

This demonstrates that the computation is stable and uses only the multiprecision arithmetic (otherwise the digits would be random after 15th decimal place).

The proven algorithm, our verified implementation (without significant changes for multiprecision), stability, and validation with existing tables ensures the accuracy of our newly calculated Kronrod nodes.

Download

October 16, 2013, Ver. 3.5.5.4731. Toolbox now provides full-featured quadgk routine in arbitrary precision.

January 16, 2016, Ver. 3.9.0.9978. The mp.GaussKronrod function has been added to toolbox to match all other functions related to Gaussian quadrature: mp.GaussLegendre, mp.GaussJacobi, mp.GaussLaguerre, mp.GaussHermite, mp.GaussChebyshevType1, mp.GaussChebyshevType2 and mp.GaussGegenbauer.

The calling syntax is: [x,w] = mp.GaussKronrod(order). The function returns nodes and weights for quadrature of specified order. Run help mp.GaussKronrod for more details.

(If you are still interested in source code of the mpkronrod and related routines, please use following link for download: Multiprecision Kronrod Nodes in MATLAB).

Tabulated Gauss-Kronrod weights and abscissae

The tables provided below display the most popular pairs of the Gauss-Kronrod quadrature nodes and weights calculated in arbitrary precision by the Multiprecision Computing Toolbox. The Gauss points are calculated using the built-in mp.GaussLegendre() function from the toolbox, and the Kronrod points are calculated using the MATLAB program converted to arbitrary precision.

All of the values are accurate up to quadruple precision (34 decimal places). These computations can easily be reproduced by running the following code:

>> mp.Digits(50);                         % Use extra precision
 
>> N = 7;                                 % Setup Gauss quadrature order
>> [xg,wg] = mp.GaussLegendre(N);         % Gauss nodes & weights for [-1,1]
>> num2str([xg wg],'%45.34e\t\t%40.34e')
ans =
    -9.4910791234275852452618968404785126e-01		1.2948496616886969327061143267908202e-01
    -7.4153118559939443986386477328078841e-01		2.7970539148927666790146777142377958e-01
    -4.0584515137739716690660641207696146e-01		3.8183005050511894495036977548897513e-01
     0.0000000000000000000000000000000000e+00		4.1795918367346938775510204081632653e-01
     4.0584515137739716690660641207696146e-01		3.8183005050511894495036977548897513e-01
     7.4153118559939443986386477328078841e-01		2.7970539148927666790146777142377958e-01
     9.4910791234275852452618968404785126e-01		1.2948496616886969327061143267908202e-01
 
>> [xk,wk] = mp.GaussKronrod(N);          % Corresponding Kronrod nodes & weights for [-1,1]
>> num2str([xk wk],'%45.34e\t\t%40.34e')
ans =
    -9.9145537112081263920685469752632852e-01		2.2935322010529224963732008058969592e-02
    -9.4910791234275852452618968404785126e-01		6.3092092629978553290700663189204287e-02
    -8.6486442335976907278971278864092620e-01		1.0479001032225018383987632254151802e-01
    -7.4153118559939443986386477328078841e-01		1.4065325971552591874518959051023792e-01
    -5.8608723546769113029414483825872960e-01		1.6900472663926790282658342659855028e-01
    -4.0584515137739716690660641207696146e-01		1.9035057806478540991325640242101368e-01
    -2.0778495500789846760068940377324491e-01		2.0443294007529889241416199923464908e-01
     1.8367099231598242312011508394097589e-40		2.0948214108472782801299917489171426e-01
     2.0778495500789846760068940377324491e-01		2.0443294007529889241416199923464908e-01
     4.0584515137739716690660641207696146e-01		1.9035057806478540991325640242101368e-01
     5.8608723546769113029414483825872960e-01		1.6900472663926790282658342659855028e-01
     7.4153118559939443986386477328078841e-01		1.4065325971552591874518959051023792e-01
     8.6486442335976907278971278864092620e-01		1.0479001032225018383987632254151802e-01
     9.4910791234275852452618968404785126e-01		6.3092092629978553290700663189204287e-02
     9.9145537112081263920685469752632852e-01		2.2935322010529224963732008058969592e-02

References

1. A.S. Kronrod (1965). Nodes and weights of quadrature formulas. Sixteen-place tables. New York: Consultants Bureau (Authorized translation from the Russian). ^
2. D. P. Laurie (1997). Calculation of Gauss-Kronrod Quadrature Rules. Mathematics of Computation, 66(219), 1133-1145. ^
3. R. Piessens, M. Branders (1974). A Note on the Optimal Addition of Abscissas to Quadrature Formulas of Gauss and Lobatto. Mathematics of Computation, 28(125), 135-139. ^
4. G. Monegato (1978). Some remarks on the construction of extended Gaussian quadrature rules. Mathematics of Computation, 32(141), 247-252. ^
5. W. Gautschi (2002). OPQ: a Matlab suite of programs for generating orthogonal polynomials and related quadrature rules. ^
6. W. Gautschi (2004). Orthogonal Polynomials: Computation and Approximation. USA: Oxford University Press. ^
7. R. Piessens, E. de Doncker–Kapenga, C. Uberhuber, D. Kahaner (1983). QUADPACK, A Subroutine Package for Automatic Integration Springer-Verlag. ^