Gauss - Kronrod quadrature

Directly to the input form

• The general theory of quadrature tells that for arbitrary but sufficiently smooth integrand the quadrature error of a formula (sometimes also called a "rule") with only nonnegative weights will be of the order of magnitude (b-a)^M+1 (and not better), where b-a is the width of the interval of integration and M is defined as the smallest integer such that the rule will integrate at least one polynomial of degree M inexact. This integer is called the order of the rule.
• The Gauss quadrature formulas are optimal with respect to the relation of attainable order to number of function evaluations (nodes) used (order = 2*number of nodes). The nodes are the zeroes of the Legendre polynomials (an orthogonal system on [-1,1]) linearly transformed to [a,b]. But in order to get computable error bounds say by comparing rules of different orders this optimality no longer holds, since the nodes of these rules have at most one node in common, namely the middle of the interval if the number of nodes is odd. Kronrod invented formulas which increase the order by adding nodes to a given Gauss formula. This way one gets order 3n+1, by adding to a classical Gauss rule with n nodes (and order 2n) n+1 more nodes in proper position. That means that with only 2n+1 function evaluations each polynomial of degree 3n is integrated exactly. Moreover one gets an approximate value of order 2n (hence a larger error for small b-a) at practically no cost by using the original Gauss nodes alone. Then one has two approximate values for the integral and their difference serves as an estimate of the error of the rule of lower order (the Gauss rule here). This idea of error estimation is used in different situations.
• Based on this one can construct adaptive integrators which work on arbitrarily large intervals. The codes used here comes from quadpack. The interval [a,b] is cut into subintervals [x_i,x_i+1] which will be halved or doubled until for each subinterval there holds
  |err_i| <= (x_i+1-x_i)/(b-a) *( epsrel |I_i| + epsabs )
  where err_i is the error estimation for this subinterval obtained as described above. The hope behind this is that the total relative error will be epsrel with an additional absolute error of epsabs, which normally is chosen much smaller than epsrel. epsrel and epsabs are user defined. A large interval is cut into pieces by default first. As approximate value for the integral one takes the sum of the values obtained from the higher order rule, of course.
• Here a large selection of orders is available, up to n=30 i.e. order 91.
• The error estimation technique assumes high order smoothness of the integrand, but indeed the Gauss rules are convergent for every continuous integrand.

Input

• You first specify the integrand. There is a large selection of predefined functions with increasing regularity. The degree of differentiability is specified in the table. The interval of integration is [0,1] for all these predefined cases.
• You also have the possibility to define the integrand yourself, using FORTRAN conventions. In this case you must specify the interval of integration too.
• You specify the two parameters epsrel and epsabs. The code tries to obtain
  |total error| <= epsabs+epsrel*|I|
  where I is the true value.
• You choose n (getting order 3*n+1) from a list.

Output

• You get a short protocol specifying the integrand, the approximated value of the integral, the error estimate (sum of the estimates for the subintervals) , for the predefined cases the true error and the number of function evaluations needed.
• A plot shows the integrand (shifted above the x-axis) and as negative impulses the positions of the subintervals (not the nodes inside) used.

Questions ?!

• Compare the efficiency of using different orders !
• Which role plays the smoothness of the integrand here?

To the input form