Double Exponential Quadrature

Directly to the input form

• In another application we compute improper integrals by transformation of the interval to [0,1] and applying a high order Gauss-Kronrod rule there. Here we go the other way around. We assume a possibly improper Riemann integrable function on an interval [a,b], transform this interval to the whole real axis by applying first a transformation of [a,b] to [-1,1]
  ξ = (2/(b-a))*(x - (a+b)/2)
  and afterwards the transformation
  ζ = arsinh( (2/π) artanh ( ξ ) )
• The background for doing this is the Poisson summation formula, which allows for approximation of integrals over the real axis by the composed equidistant trapezoidal sum:
  Theorem: assume that F is without singularities on the real axis , that
  ∫_R |F(x)| dx exists and that ω > 0 is the smallest distance of any possible complex singularity of F from the real line. Then
  ∫_R F(x)dx = lim _{k to ∞} h ∑ _{i=-k to k} F(x_i) + O(exp(-ω π/h))
  with h the distance on the grid of the x_i.
• This means that the error in this approximation decreases exponentially with decreasing h and that halving h roughly will double the number of correct digits. In order to be useful F must decay fast with increasing magnitude of x.
• Translating this back into the variable range of ξ one gets
  ∫ _{-1 to 1} f(ξ) d ξ = lim_{k to ∞} ∑_{i = -k to k} w_i f(ξ_i) + O(exp(- 2 ω π/h))
  with
  ξ_i = tanh( π/2 sinh( i h ) )
  w_i = (h π/2 cosh( i h ) )/( cosh²(π/2 sinh( i h ) )
  Hence the nodes of the integration become more and more dense at the boundaries of the interval [-1,1] and the weights will decay extremely fast there.
• Since we allow singularities at -1 and 1 we avoid numerical trouble by staying away from the boundaries by 10 ε with ε the machine precision. In the output you will get the information of the true evaluation interval, which is also the plot-interval. Moreover we also avoid trouble with overflow or underflow of the hyperbolic cos and sin by limiting the product i h to 3.135.
• Our first approximation works with 2k₀+1 function values with k₀ user defined. Afterwards the grid is refined by halving h and this at most 20 times. We finish earlier if two successive approximations of the integral differ by less than
  10^-14(∫|F(.)|dx + 10^-14)
  with the integral for the absvalue function approximated by the same method simultaneously.

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,1[.
• You also have the possibility to define an integrand yourself using FORTRAN- conventions. In that case you also specify the interval. You may also specify whether you know the true integral value in which case you type this in the specified field.
• Finally you fix k₀. This must be 3 at least, i.e. 7 function values are required.

Output

• You get a short protocol specifying the integrand, giving the approximate integral, and in the case of the exact value known the true error or otherwise the error estimate. You also get the number of function evaluations used.
• A plot shows the integrand (more precisely the piecewise linear interpolant of the function values used) shifted above the x-axis (i.e. after subtraction of the minimal function value) and as negative impulses the grid of the x-values used.
• A 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!
• A second plot shows you in a double logarithmic scale the development of the error over the h-axis if the true value of the integral is given.

Questions ?!

• Can the exponential decay of the error be verified?
• 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