Improper integrals

Directly to the input form

 
  • Gauss-Kronrod rules are quadrature formulas of very high order (3n+1). They are obtained by adding to a classical Gauss rule with n nodes and order 2n, (that is exactness for polynomials of degree <= 2n-1) n+1 further nodes in proper position. These rules integrate with only 2n+1 function evaluations any polynomial of degree 3n exactly. With these function evaluations one can compute now two approximate integral values, one of order 2n and error of the type
    C2nf(2n)(x~)(b-a)2n+1
    and one of order 3n+1 with error
    C~3n+1f(3n+1)(x#)(b-a)3n+2.
    For small b-a the difference of these two values will give a good error estimation for the value of smaller order. The nodes of these formulas lie always in the interior of the interval. Hence they are usable as long as the integrand has no singularity inside the interval of integration even if it has an integrable singularity at the boundary.
  • The nodes of these formulas, originally located in [-1,1], can be transformed linearly to any finite interval [xl,,xr].
  • An improper integral over an infinite interval is transformed into an improper integral over a finite interval first, therefore. For an integral b f(x) dx one uses the transformation z = 1/(x-b+1) yielding the interval (0,1] with the singularity now usually at z=0. If originally a singularity had been at x=b too this is now located at z=1.
  • Now ]0,1[ is cut into smaller pieces ]zi,zi+1[ with the aim to satisfy the requirement
    |erri| <= (zi+1-zi)/(b-a) *( epsrel |Ii| +epsabs)
    for the integration error, where erri is the difference described above. This is achieved by halving the interval length where necessary. Ii is the estimated true value of the integral over the subinterval. Would this hold true, then the total relative error would be approximately epsrel.
  • Improper integrals over (-∞,b] are treated by reversing the interval (and the sign) and such over R as a whole are decomposed into two using 0 for the cut. In this latter case you finally will see the transformed integrand on [-1,1] instead of [0,1].
  • the numerical code behind this is dqagi from quadpack.
  • the number of grid refinements is limited to 2000 here.
 

Input

 
  • You first specify the integrand. There are some predefined cases where the true value of the integral is known. Here the interval is always [0,∞[.
  • You also have the possibility to define an integrand yourself using FORTRAN- conventions. In that case you also specify the interval resp. the type of integral.
  • You specify the two parameters epsrel and epsabs.
  • Here we have n fixed to n=7, that is order 21.
 

Output

 
  • You get a short protocol specifying the integrand, giving the approximate integral, the error estimate and in the predefined cases the true error as well as the number of function evaluations used.
  • A plot shows the integrand (after transformation to ]0,1]) (more precisely the piecewise linear interpolant of the function values used) shifted above the z-axis (i.e. after subtraction of the minimal function value) and as negative impulses the grid of the subintervals.
  • If the required precision is unreachable with the internal limits for memory, then you get an appropriate error message
  • A warning is at place: if you use trigonometric functions here and the integrand decays slowly, then excessive large arguments may be used in evaluating these, which, due to the internally used reduction modulo π/4 might result in nonsensical function values.
  • A further warning is at place: think about the limitations of such a plot: if, say, an integrand varies over several powers of 10 and this at the boundary of the interval and in the remaining part it is small, then you might see little more than the axis!
 

Questions ?!

 
  • If the interval has infinite length, then the integrand must decay to zero at infinity. Which role plays the type of decay here ?
  • What takes place if there is a point of discontinuity or even singularity inside the interval?
  • What takes place if the integral doesn't exist in the Riemann sense?
 

To the input form

 
 
Back to the top!

18.02.2015