|
Improper integrals |
|
|
|
- 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?
|
|
|
|
|