Adaptive 2D quadrature using iterated 1D Simpson quadrature

Directly to the input form

• This application deals with the computation of integrals of the type
  I = ∫_B f(x,y) dxdy
  where B is a region in the plane. In most cases (and possibly after changing the role of the two free variables x and y) it will be possible to cut B into finitely many regiones of the special form
  B = { (x,y): a ≤ x ≤ b , ψ(x) ≤ y ≤ φ(x) }
  with two smooth functions ψ and φ , a so called ''normal area w.r.t. the x-axis''.
  This allows us to cut down the job to a nested 1D quadrature using
  I = ∫_a^b F(x) dx
  where
  F(x) = ∫_ψ(x)^φ(x)f(x,y)dy .
  Hence every evaluation of F involves a 1D quadrature. With f, ψ and φ smooth functions w.r.t. to x and f also continuous smooth w.r.t. y F is smooth w.r.t. x and this opens the job for 1D quadrature rules and this is the way we go. We use the adaptive Simpson code presented in this section for the 1D case, but not, as you might suspect, in a recursive way but simply by using a copy of the code with a ''2'' appended to each subroutine name for realizing F. We use the same order of accuracy for both variables and the same adaptive scheme.
• For the theoretical background of the adaptive Simpson code see there.
• This is not the most efficient way to do the job, but specialized (and much more efficient) rules require often approximating also the region B in a manner which fits in precision with these special rules.

Input

• You first specify the integrand. There are five test cases predefined but you also have the option to define f yourself, using a piece of code written following FORTRAN rules. The latter is case 6. Don't forget to click it active!
• In case 6 You also specify the interval [a,b] and have the choice to specify whether you know the exact result of this case or not.
• Then you specify the desired precision via the variable
      genau, (ideally this means (I - I_approx)/I).
  The bound for absolute error is fixed here to machine precision squared.
  Attention: these parameters are applied in the ''inner''' integration (i.e. evaluation of f(x,y), as well as in the ''outer'' integration (the integration of F). There are internal restrictions: genau must be in [1.0^-12,0.01], hmax ≤ b-a and 1.0^-6 ≤ hmin ≤ hmax . Wrong input will be corrected without notice to fit these bounds. You may also select whether the integrand should be plotted over B.

Output

• You get a short protocol about your input and the obtained results.
• A plot shows the region B, the x-grid and along the lines x constant and ψ(x) ≤ y ≤ φ(x) the evaluation points of f. This may not cover the whole region in case of premature termination of the code. You also may see a plot of the integrand over the region.

Questions?!

• If f is polynomial in y and ψ and φ polynomial in x how many nodes will you need at least to get the integral exactly (up to roundoff)? How many evaluations reports the code and why? (Hint: think about the acceptance rule of the adaptivity scheme and the order of the Simpson rule. Hence: where no adaptivity is used at all?
• Testcase one has f=1 Why uses the code so much function values? And why the progressive narrowing of the grid near 1?
• Testcase 3 makes trouble here , why? Compare with the Gauss-Kronrod application!
• explain the result of testcase 5!
• What can be observed if you use discontinuous or nonsmooth functions (for example sign(...))?

To the input form