Adaptive quadrature using Simpson's rule

Directly to the input form

• The numerical approximation of definite integrals can be based on integrating interpolation polynomials exactly. One replaces the integrand f by an interpolation polynomial, be it on the whole interval [a,b] or piecewise on small subintervals of this. These polynomials can be integrated analytically and the obtained value resp. their sum replaces the true value. The error of this process hence simply is the integral of the interpolation remainder.
• Here we use Simpson's rule in its composed form, but not on an equidistant grid but rather on a locally variable grid. Simpson's rule is obtained by integrating a second degree polynomial, interpolating f at three points in a symmetric position: xl=x,x+h,xr=x+2*h. This yields the value
  I₁ = (1/3)*h*(f(xl)+4*f(xl+h)+f(xr)) , xr=xl+2h .
  The abbreviation I₁ becomes clear below.
• One would originally guess that this is exact for second order polynomials only, but indeed it turns out that it is exact for third order polynomials too. This is due to its symmetry. The length of the subintervals [xl,xr] is computed adaptively. The total error of this process can be bounded from above by
  (1/90) × max{ |f⁽⁴⁾(x)| : x in [a,b] } × ∑_i h_i⁵.
  This bound however can be quite crude.
• The method however works if the integrand is of lower differentiability, the order of convergence is at least like O(h_max^k) if f is k, k=1,2,3 times differentiable and even if it is only continuous, convergence is guaranteed for h_max → 0 since the formula is a special form of a Riemannian sum.
• Here we determine the length of the subintervals by estimating the quadrature error locally: a more precise analysis of the remainder of Simpson's rule yields as error the expression
  E₁ = (1/180)*h⁴*(f⁽³⁾(xr)-f⁽³⁾(xl))+O(h⁶)
  provided f is five times continuously differentiable. Hence, if we use also the second estimate for the integral, obtained from
  I₂ = (1/6)*h*(f(xl)+4*f(xl+h/2)+2*f(xl+h)+4*f(xr-h/2)+f(xr))
  (i.e. using Simpson's rule twice with the stepsize halved), then we get as error of this
  E₂ = (1/180)*h⁴/16*(f⁽³⁾(xr)-f⁽³⁾(xl))+O(h⁶)
  Hence, up to an error O(h⁶) the error E₁ can be obtained using
  I_true - I₁ = E₁
  I_true - I₂ = E₂.
  Hence
  I₂ - I₁ = E₁ - E₂ = (15/16) E₁
  and hence
  E₁ = (16/15)*(I₂ - I₁) = errest(h).
  But errest(h) is computable!
• Given a desired bound of the total error eps for the integral we require
  |errest(h)| <= eps*(xr-xl)/(b-a)
  as precision for the subinterval. Now let c*h_cur be the largest stepsize compatible with this and errest be obtained with the stepsize h_cur . Inserting this into the expression for E₁ yields
  c = (eps*h_cur*15/(16*(b-a)*(|I₂ - I₁|)) ^0.25.
  If it turns out that c >= 1 then we accept I₂ for this subinterval and proceed by replacing xl with xr and h_cur by c*h_cur. Otherwise the current step is rejected and the step repeated with a reduced h, in this case we use here 0.9*c*h_cur.
• This adaptive control of the stepsize is enhanced by providing two ''safety'' parameters h_min, h_max which never will be exceeded.
• Summing up all values errest we get an estimate of the total error.
• One must be aware of the fact that this is not a strict error bound, we neglected higher order terms in our arguments. Clearly this can fail.

Input

• You first specify the integrand. There are two test cases predefined but you also have the option to define f yourself, using a piece of code written following FORTRAN rules.
• You also specify the interval and have the choice to specify eps. Attention: eps is an absolute error, if f can be quite large this never might be reachable. We provide a minimal stepsize 1.0e-6*(b-a) (that means you may use about one million function values), but we terminate prematurely if this is attained.
• Because of the neglection of the higher order terms, we also require a maximal stepsize h_max to be specified. You might play a bit with h_max which however is internally bounded from above by (b-a)/10

Output

• You get a short protocol about your input and the obtained result.
• A plot shows the interand (shifted above the x-axis, you see f-min{f}) and below the x-axis a grid showing the subintervals actually used in the process. This can be shorter than [a,b] in case of premature termination, in which case the results are for the indicated subinterval only.

Questions?!

• In which cases you will get trouble ?
• What can be said about the effort, measured in function evaluations ?
• How should h_max be restricted at least?
• Which effect has the choice of h_max on the total effort ?
• What can be observed if you use discontinuous or nonsmooth functions (for example sign(...))?

To the input form