Adaptive 2D quadrature using iterated 1D 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 Gaussian quadrature rules and this is the way we go. We use the 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 (although in special cases it may be much more efficient tp use different orders, if for example f is polynomial in one or both variables.)
• For the theoretical background of Gaussian quadrature rules see there.
• This is not the most efficient way to do 2D quadrature, 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 wether you know the exact result of this case or not.
• Then you specify the desired precision via the variable epsrel (ideally (I - I_approx)/I and epsabs , the desired absolute error (ideally I - I_approx)n. Attention: these parameters are applied in the ''inner''' integration (i.e. evaluation of f(x,y) as well as in the ''outer'' integration (integration of F). There are internal restrictions: epsrel must be in [0.01,1.0^-12] and epsabs in [1.0^-6,1.0^-6,1.0^-16]. Wrong input will be corrected without notice. li> You also select the number of Gauss nodes to be used on each ''panel'' like in the 1D case using the input ''key''. This takes the values 1 to 6. and has the meaning that (n,2n+1) Gauss-Kronrod pairs are used reaching order 3n+1 (exact for polynomials of degree <=3n) for n=7,10,15,20,25,30 respectively.

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?
• Despite the boundary singularity Testcase 3 makes no much 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