Numerical Integration
*********************

Octave comes with several built-in functions for computing the integral of a function numerically (termed quadrature).  These functions all solve 1-dimensional integration problems.

Functions of One Variable
=========================

Octave supports five different adaptive quadrature algorithms for computing the integral of a function f over the interval from a to b.  These are

quad
     Numerical integration based on Gaussian quadrature.

quadv
     Numerical integration using an adaptive vectorized Simpson's rule.

quadl
     Numerical integration using an adaptive Lobatto rule.

quadgk
     Numerical integration using an adaptive Gauss-Konrod rule.

quadcc
     Numerical integration using adaptive Clenshaw-Curtis rules.

     In addition, the following functions are also provided:

integral
     A compatibility wrapper function that will choose between quadv and quadgk depending on the integrand and options chosen.

trapz, cumtrapz
     Numerical integration of data using the trapezoidal method.

The best quadrature algorithm to use depends on the integrand.  If you have empirical data, rather than a function, the choice is trapz or cumtrapz.  If you are uncertain about the characteristics of the integrand, quadcc will be the most robust as it can handle discontinuities, singularities, oscillatory functions, and infinite intervals. When the integrand is smooth quadgk may be the fastest of the algorithms.

  Function   Characteristics
  quad     Low accuracy with nonsmooth integrands
  quadv    Medium accuracy with smooth integrands
  quadl    Medium accuracy with smooth integrands.  Slower than quadgk.
  quadgk   Medium accuracy (1e-6 -- 1e-9) with smooth integrands.
           Handles oscillatory functions and infinite bounds
  quadcc   Low to High accuracy with nonsmooth/smooth integrands
           Handles oscillatory functions, singularities, and infinite bounds

Here is an example of using quad to integrate the function

       F(X) = X * sin (1/X) * sqrt (abs (1 - X))

from X = 0 to X = 3.

This is a fairly difficult integration (plot the function over the range of integration to see why).

The first step is to define the function:

     function y = f (x)
       y = x .* sin (1./x) .* sqrt (abs (1 - x));
     endfunction

Note the use of the `dot' forms of the operators.  This is not necessary for the quad integrator, but is required by the other integrators.  In any case, it makes it much easier to generate a set of points for plotting because it is possible to call the function with a vector argument to produce a vector result.

The second step is to call quad with the limits of integration:

     [q, ier, nfun, err] = quad ("f", 0, 3)
          => 1.9819
          => 1
          => 5061
          => 1.1522e-07

Although quad returns a nonzero value for IER, the result is reasonably accurate (to see why, examine what happens to the result if you move the lower bound to 0.1, then 0.01, then 0.001, etc.).

The function "f" can be the string name of a function or a function handle.  These options make it quite easy to do integration without having to fully define a function in an m-file.  For example:

     # Verify gamma function = (n-1)! for n = 4
     f = @(x) x.^3 .* exp (-x);
     quadcc (f, 0, Inf)
          => 6.0000

 -- Q =: quad (F, A, B)
 -- Q =: quad (F, A, B, TOL)
 -- Q =: quad (F, A, B, TOL, SING)
 -- [Q, IER, NFEV, ERR] =: quad (...)
     Numerically evaluate the integral of F from A to B using Fortran routines from QUADPACK.

     F is a function handle, inline function, or a string containing the name of the function to evaluate.  The function must have the form y = f (x) where Y and X are scalars.

     A and B are the lower and upper limits of integration.  Either or both may be infinite.

     The optional argument TOL is a vector that specifies the desired accuracy of the result.  The first element of the vector is the desired absolute tolerance, and the second element is the desired relative tolerance.  To choose a relative test only, set the absolute tolerance to zero.  To choose an absolute test only, set the relative tolerance to zero.  Both tolerances default to sqrt (eps) or approximately 1.5e-8.

     The optional argument SING is a vector of values at which the integrand is known to be singular.

     The result of the integration is returned in Q.

     IER contains an integer error code (0 indicates a successful integration).

     NFEV indicates the number of function evaluations that were made.

     ERR contains an estimate of the error in the solution.

     The function quad_options can set other optional parameters for quad.

     Note: because quad is written in Fortran it cannot be called recursively.  This prevents its use in integrating over more than one variable by routines dblquad and triplequad. quad_options, quadv, quadl, quadgk, quadcc, trapz, dblquad, triplequad

 -- quad_options: ()
 -- val =: quad_options (OPT)
 -- quad_options: (OPT, VAL)
     Query or set options for the function quad.

     When called with no arguments, the names of all available options and their current values are displayed.

     Given one argument, return the value of the option OPT.

     When called with two arguments, quad_options sets the option OPT to value VAL.

     Options include

     "absolute tolerance"
          Absolute tolerance; may be zero for pure relative error test.

     "relative tolerance"
          Non-negative relative tolerance.  If the absolute tolerance is zero, the relative tolerance must be greater than or equal to max (50*eps, 0.5e-28).

     "single precision absolute tolerance"
          Absolute tolerance for single precision; may be zero for pure relative error test.

     "single precision relative tolerance"
          Non-negative relative tolerance for single precision.  If the absolute tolerance is zero, the relative tolerance must be greater than or equal to max (50*eps, 0.5e-28).

 -- Q =: quadv (F, A, B)
 -- Q =: quadv (F, A, B, TOL)
 -- Q =: quadv (F, A, B, TOL, TRACE)
 -- Q =: quadv (F, A, B, TOL, TRACE, P1, P2, ...)
 -- [Q, NFEV] =: quadv (...)

     Numerically evaluate the integral of F from A to B using an adaptive Simpson's rule.

     F is a function handle, inline function, or string containing the name of the function to evaluate.  quadv is a vectorized version of quad and the function defined by F must accept a scalar or vector as input and return a scalar, vector, or array as output.

     A and B are the lower and upper limits of integration.  Both limits must be finite.

     The optional argument TOL defines the absolute tolerance used to stop the adaptation procedure.  The default value is 1e-6.

     The algorithm used by quadv involves recursively subdividing the integration interval and applying Simpson's rule on each subinterval. If TRACE is true then after computing each of these partial integrals display: (1) the total number of function evaluations, (2) the left end of the subinterval, (3) the length of the subinterval, (4) the approximation of the integral over the subinterval.

     Additional arguments P1, etc., are passed directly to the function F.  To use default values for TOL and TRACE, one may pass empty matrices ([]).

     The result of the integration is returned in Q.

     The optional output NFEV indicates the total number of function evaluations performed.

     Note: quadv is written in Octave's scripting language and can be used recursively in dblquad and triplequad, unlike the quad function. quad, quadl, quadgk, quadcc, trapz, dblquad, triplequad, integral, integral2, integral3

 -- Q =: quadl (F, A, B)
 -- Q =: quadl (F, A, B, TOL)
 -- Q =: quadl (F, A, B, TOL, TRACE)
 -- Q =: quadl (F, A, B, TOL, TRACE, P1, P2, ...)
 -- [Q, NFEV] =: quadl (...)

     Numerically evaluate the integral of F from A to B using an adaptive Lobatto rule.

     F is a function handle, inline function, or string containing the name of the function to evaluate.  The function F must be vectorized and return a vector of output values when given a vector of input values.

     A and B are the lower and upper limits of integration.  Both limits must be finite.

     The optional argument TOL defines the absolute tolerance with which to perform the integration.  The default value is 1e-6.

     The algorithm used by quadl involves recursively subdividing the integration interval.  If TRACE is defined then for each subinterval display: (1) the total number of function evaluations, (2) the left end of the subinterval, (3) the length of the subinterval, (4) the approximation of the integral over the subinterval.

     Additional arguments P1, etc., are passed directly to the function F.  To use default values for TOL and TRACE, one may pass empty matrices ([]).

     The result of the integration is returned in Q.

     The optional output NFEV indicates the total number of function evaluations performed.

     Reference: W. Gander and W. Gautschi, Adaptive Quadrature - Revisited, BIT Vol. 40, No. 1, March 2000, pp. 84--101. https://www.inf.ethz.ch/personal/gander/ quad, quadv, quadgk, quadcc, trapz, dblquad, triplequad, integral, integral2, integral3

 -- Q =: quadgk (F, A, B)
 -- Q =: quadgk (F, A, B, ABSTOL)
 -- Q =: quadgk (F, A, B, ABSTOL, TRACE)
 -- Q =: quadgk (F, A, B, "PROP", VAL, ...)
 -- [Q, ERR] =: quadgk (...)

     Numerically evaluate the integral of F from A to B using adaptive Gauss-Kronrod quadrature.

     F is a function handle, inline function, or string containing the name of the function to evaluate.  The function F must be vectorized and return a vector of output values when given a vector of input values (See property "ArrayValued" for an exception to this rule).

     A and B are the lower and upper limits of integration.  Either or both limits may be infinite or contain weak end singularities.  Variable transformation will be used to treat any infinite intervals and weaken the singularities.  For example:

          quadgk (@(x) 1 ./ (sqrt (x) .* (x + 1)), 0, Inf)

     Note that the formulation of the integrand uses the element-by-element operator ./ and all user functions to quadgk should do the same.

     The optional argument ABSTOL defines the absolute tolerance used to stop the integration procedure.  The default value is 1e-10 (1e-5 for single).

     The algorithm used by quadgk involves subdividing the integration interval and evaluating each subinterval.  If TRACE is true then after computing each of these partial integrals display: (1) the number of subintervals at this step, (2) the current estimate of the error ERR, (3) the current estimate for the integral Q.

     The behavior of the algorithm can be configured by passing arguments to quadgk as pairs "PROP", VAL.  Valid properties are

     AbsTol
          Define the absolute error tolerance for the quadrature.  The default absolute tolerance is 1e-10 (1e-5 for single).

     RelTol
          Define the relative error tolerance for the quadrature.  The default relative tolerance is 1e-6 (1e-4 for single).

     ArrayValued
          When set to true, the function F produces an array output for a scalar input.  The default is false which requires that F produce an output that is the same size as the input.  For example,

               quadgk (@(x) x .^ (1:5), 0, 2, "ArrayValued", 1)

          will integrate [x.^1, x.^2, x.^3, x.^4, x.^5] in one function call rather than having to repeatedly define a single anonymous function and use a normal invocation of quadgk.

     WayPoints
          Specify points which will become endpoints for subintervals in the algorithm which can result in significantly improved estimation of the error in the integral, faster computation, or both.  It can be useful to specify more subintervals around a region where the integrand is rapidly changing or to flag locations where there is a discontinuity in the first derivative of the function.  For example, the signum function has a discontinuity at x == 0 and by specifying a waypoint

               quadgk (@(x) sign (x), -0.5, 1, "Waypoints", [0])

          the error bound is reduced from 4e-7 to 1e-13.

          If the function has singularities within the region of integration those should not be addressed with waypoints.  Instead, the overall integral should be decomposed into a sum of several smaller integrals such that the singularity occurs as one of the bounds of integration in the call to quadgk.

          If any of the waypoints are complex then contour integration is performed as documented below.

     MaxIntervalCount
          quadgk initially subdivides the interval on which to perform the quadrature into 10 intervals or, if WayPoints are given, at each waypoint. Subintervals that have an unacceptable error are subdivided and re-evaluated.  If the number of subintervals exceeds 650 subintervals at any point then a poor convergence is signaled and the current estimate of the integral is returned.  The property "MaxIntervalCount" can be used to alter the number of subintervals that can exist before exiting.

     Trace
          If logically true quadgk prints information on the convergence of the quadrature at each iteration.

     If any of A, B, or WAYPOINTS is complex then the quadrature is treated as a contour integral along a piecewise linear path defined by [A, WAYPOINTS(1), WAYPOINTS(2), ..., B]. In this case the integral is assumed to have no edge singularities.  For example,

          quadgk (@(z) log (z), 1+1i, 1+1i, "WayPoints",
                  [-1+1i, -1-1i, +1-1i])

     integrates log (z) along the square defined by [1+1i, -1+1i, -1-1i, +1-1i].

     The result of the integration is returned in Q.

     ERR is an approximate bound on the error in the integral abs (Q - I), where I is the exact value of the integral.  If the adaptive integration did not converge, the value of ERR will be larger than the requested tolerance.  If only a single output is requested then a warning will be emitted when the requested tolerance is not met.  If the second output ERR is requested then no warning is issued and it is the responsibility of the programmer to inspect and determine whether the results are satisfactory.

     Reference: L.F. Shampine, "Vectorized adaptive quadrature in MATLAB", Journal of Computational and Applied Mathematics, pp. 131--140, Vol 211, Issue 2, Feb 2008.

     quad, quadv, quadl, quadcc, trapz, dblquad, triplequad, integral, integral2, integral3

 -- Q =: quadcc (F, A, B)
 -- Q =: quadcc (F, A, B, TOL)
 -- Q =: quadcc (F, A, B, TOL, SING)
 -- [Q, ERR, NR_POINTS] =: quadcc (...)
     Numerically evaluate the integral of F from A to B using doubly-adaptive Clenshaw-Curtis quadrature.

     F is a function handle, inline function, or string containing the name of the function to evaluate.  The function F must be vectorized and must return a vector of output values if given a vector of input values. For example,

          f = @(x) x .* sin (1./x) .* sqrt (abs (1 - x));

     which uses the element-by-element ``dot'' form for all operators.

     A and B are the lower and upper limits of integration.  Either or both limits may be infinite.  quadcc handles an infinite limit by substituting the variable of integration with x = tan (pi/2*u).

     The optional argument TOL is a 1- or 2-element vector that specifies the desired accuracy of the result.  The first element of the vector is the desired absolute tolerance, and the second element is the desired relative tolerance. To choose a relative test only, set the absolute tolerance to zero.  To choose an absolute test only, set the relative tolerance to zero.  The default absolute tolerance is 1e-10 (1e-5 for single), and the default relative tolerance is 1e-6 (1e-4 for single).

     The optional argument SING contains a list of points where the integrand has known singularities, or discontinuities in any of its derivatives, inside the integration interval.  For the example above, which has a discontinuity at x=1, the call to quadcc would be as follows

          int = quadcc (f, a, b, [], [ 1 ]);

     The result of the integration is returned in Q.

     ERR is an estimate of the absolute integration error.

     NR_POINTS is the number of points at which the integrand was evaluated.

     If the adaptive integration did not converge, the value of ERR will be larger than the requested tolerance.  If only a single output is requested then a warning will be emitted when the requested tolerance is not met.  If the second output ERR is requested then no warning is issued and it is the responsibility of the programmer to inspect and determine whether the results are satisfactory.

     quadcc is capable of dealing with non-numeric values of the integrand such as NaN or Inf.  If the integral diverges, and quadcc detects this, then a warning is issued and Inf or -Inf is returned.

     Note: quadcc is a general purpose quadrature algorithm and, as such, may be less efficient for a smooth or otherwise well-behaved integrand than other methods such as quadgk.

     The algorithm uses Clenshaw-Curtis quadrature rules of increasing degree in each interval and bisects the interval if either the function does not appear to be smooth or a rule of maximum degree has been reached.  The error estimate is computed from the L2-norm of the difference between two successive interpolations of the integrand over the nodes of the respective quadrature rules.

     Reference: P. Gonnet, Increasing the Reliability of Adaptive Quadrature Using Explicit Interpolants, ACM Transactions on Mathematical Software, Vol. 37, Issue 3, Article No. 3, 2010. quad, quadv, quadl, quadgk, trapz, dblquad, triplequad

 -- Q =: integral (F, A, B)
 -- Q =: integral (F, A, B, PROP, VAL, ...)
 -- [Q, ERR] =: integral (...)

     Numerically evaluate the integral of F from A to B using adaptive quadrature.

     integral is a wrapper for quadcc (general real-valued, scalar integrands and limits), and quadgk (integrals with specified integration paths and array-valued integrands) that is intended to provide MATLAB compatibility.  More control of the numerical integration may be achievable by calling the various quadrature functions directly.

     F is a function handle, inline function, or string containing the name of the function to evaluate.  The function F must be vectorized and return a vector of output values when given a vector of input values.

     A and B are the lower and upper limits of integration.  Either or both limits may be infinite or contain weak end singularities.  If either or both limits are complex, integral will perform a straight line path integral.  Alternatively, a complex domain path can be specified using the "Waypoints" option (see below).

     Additional optional parameters can be specified using "PROPERTY", VALUE pairs.  Valid properties are:

     Waypoints
          Specifies points to be used in defining subintervals of the quadrature algorithm, or if A, B, or WAYPOINTS are complex then the quadrature is calculated as a contour integral along a piecewise continuous path.  For more detail, see quadgk.

     ArrayValued
          integral expects F to return a scalar value unless ARRAYVALUED is specified as true.  This option will cause integral to perform the integration over the entire array and return Q with the same dimensions as returned by F.  For more detail see quadgk.

     AbsTol
          Define the absolute error tolerance for the quadrature.  The default absolute tolerance is 1e-10 (1e-5 for single).

     RelTol
          Define the relative error tolerance for the quadrature.  The default relative tolerance is 1e-6 (1e-4 for single).

     The optional output ERR contains the absolute error estimate used by the called integrator.

     Adaptive quadrature is used to minimize the estimate of error until the following is satisfied:

            ERROR <= max (ABSTOL, RELTOL*|Q|).

     Known MATLAB incompatibilities:

      1. If tolerances are left unspecified, and any integration limits or waypoints are of type single, then Octave's integral functions automatically reduce the default absolute and relative error tolerances as specified above.  If tighter tolerances are desired they must be specified. MATLAB leaves the tighter tolerances appropriate for double inputs in place regardless of the class of the integration limits.

     integral2, integral3, quad, quadgk, quadv, quadl, quadcc, trapz, dblquad, triplequad

Sometimes one does not have the function, but only the raw (x, y) points from which to perform an integration.  This can occur when collecting data in an experiment.  The trapz function can integrate these values as shown in the following example where "data" has been collected on the cosine function over the range [0, pi/2).

     x = 0:0.1:pi/2;  # Uniformly spaced points
     y = cos (x);
     trapz (x, y)
          => 0.99666

The answer is reasonably close to the exact value of 1.  Ordinary quadrature is sensitive to the characteristics of the integrand.  Empirical integration depends not just on the integrand, but also on the particular points chosen to represent the function.  Repeating the example above with the sine function over the range [0, pi/2) produces far inferior results.

     x = 0:0.1:pi/2;  # Uniformly spaced points
     y = sin (x);
     trapz (x, y)
          => 0.92849

However, a slightly different choice of data points can change the result significantly.  The same integration, with the same number of points, but spaced differently produces a more accurate answer.

     x = linspace (0, pi/2, 16);  # Uniformly spaced, but including endpoint
     y = sin (x);
     trapz (x, y)
          => 0.99909

In general there may be no way of knowing the best distribution of points ahead of time.  Or the points may come from an experiment where there is no freedom to select the best distribution.  In any case, one must remain aware of this issue when using trapz.

 -- Q =: trapz (Y)
 -- Q =: trapz (X, Y)
 -- Q =: trapz (..., DIM)

     Numerically evaluate the integral of points Y using the trapezoidal method.

     trapz (Y) computes the integral of Y along the first non-singleton dimension.  When the argument X is omitted an equally spaced X vector with unit spacing (1) is assumed. trapz (X, Y) evaluates the integral with respect to the spacing in X and the values in Y.  This is useful if the points in Y have been sampled unevenly.

     If the optional DIM argument is given, operate along this dimension.

     Application Note: If X is not specified then unit spacing will be used.  To scale the integral to the correct value you must multiply by the actual spacing value (deltaX).  As an example, the integral of x^3 over the range [0, 1] is x^4/4 or 0.25.  The following code uses trapz to calculate the integral in three different ways.

          x = 0:0.1:1;
          y = x.^3;
          ## No scaling
          q = trapz (y)
            => q = 2.5250
          ## Approximation to integral by scaling
          q * 0.1
            => 0.25250
          ## Same result by specifying X
          trapz (x, y)
            => 0.25250

     cumtrapz

 -- Q =: cumtrapz (Y)
 -- Q =: cumtrapz (X, Y)
 -- Q =: cumtrapz (..., DIM)
     Cumulative numerical integration of points Y using the trapezoidal method.

     cumtrapz (Y) computes the cumulative integral of Y along the first non-singleton dimension.  Where trapz reports only the overall integral sum, cumtrapz reports the current partial sum value at each point of Y.

     When the argument X is omitted an equally spaced X vector with unit spacing (1) is assumed.  cumtrapz (X, Y) evaluates the integral with respect to the spacing in X and the values in Y.  This is useful if the points in Y have been sampled unevenly.

     If the optional DIM argument is given, operate along this dimension.

     Application Note: If X is not specified then unit spacing will be used.  To scale the integral to the correct value you must multiply by the actual spacing value (deltaX). trapz, cumsum

Orthogonal Collocation
======================

 -- [R, AMAT, BMAT, Q] =: colloc (N, "left", "right")
     Compute derivative and integral weight matrices for orthogonal collocation.

     Reference: J. Villadsen, M. L. Michelsen, Solution of Differential Equation Models by Polynomial Approximation.

Here is an example of using colloc to generate weight matrices for solving the second order differential equation U' - ALPHA * U'' = 0 with the boundary conditions U(0) = 0 and U(1) = 1.

First, we can generate the weight matrices for N points (including the endpoints of the interval), and incorporate the boundary conditions in the right hand side (for a specific value of ALPHA).

     n = 7;
     alpha = 0.1;
     [r, a, b] = colloc (n-2, "left", "right");
     at = a(2:n-1,2:n-1);
     bt = b(2:n-1,2:n-1);
     rhs = alpha * b(2:n-1,n) - a(2:n-1,n);

Then the solution at the roots R is

     u = [ 0; (at - alpha * bt) \ rhs; 1]
          => [ 0.00; 0.004; 0.01 0.00; 0.12; 0.62; 1.00 ]

Functions of Multiple Variables
===============================

Octave includes several functions for computing the integral of functions of multiple variables.  This procedure can generally be performed by creating a function that integrates f with respect to x, and then integrates that function with respect to y.  This procedure can be performed manually using the following example which integrates the function:

     f(x, y) = sin(pi*x*y) * sqrt(x*y)

for x and y between 0 and 1.

Using quadgk in the example below, a double integration can be performed.  (Note that any of the 1-D quadrature functions can be used in this fashion except for quad since it is written in Fortran and cannot be called recursively.)

     function q = g(y)
       q = ones (size (y));
       for i = 1:length (y)
         f = @(x) sin (pi*x.*y(i)) .* sqrt (x.*y(i));
         q(i) = quadgk (f, 0, 1);
       endfor
     endfunction

     I = quadgk ("g", 0, 1)
           => 0.30022

The algorithm above is implemented in the function dblquad for integrals over two variables.  The 3-D equivalent of this process is implemented in triplequad for integrals over three variables.  As an example, the result above can be replicated with a call to dblquad as shown below.

     I = dblquad (@(x, y) sin (pi*x.*y) .* sqrt (x.*y), 0, 1, 0, 1)
           => 0.30022

 -- Q =: dblquad (F, XA, XB, YA, YB)
 -- Q =: dblquad (F, XA, XB, YA, YB, TOL)
 -- Q =: dblquad (F, XA, XB, YA, YB, TOL, QUADF)
 -- Q =: dblquad (F, XA, XB, YA, YB, TOL, QUADF, ...)
     Numerically evaluate the double integral of F.

     F is a function handle, inline function, or string containing the name of the function to evaluate.  The function F must have the form z = f(x,y) where X is a vector and Y is a scalar.  It should return a vector of the same length and orientation as X.

     XA, YA and XB, YB are the lower and upper limits of integration for x and y respectively.  The underlying integrator determines whether infinite bounds are accepted.

     The optional argument TOL defines the absolute tolerance used to integrate each sub-integral.  The default value is 1e-6.

     The optional argument QUADF specifies which underlying integrator function to use.  Any choice but quad is available and the default is quadcc.

     Additional arguments, are passed directly to F.  To use the default value for TOL or QUADF one may pass ':' or an empty matrix ([]). integral2, integral3, triplequad, quad, quadv, quadl, quadgk, quadcc, trapz

 -- Q =: triplequad (F, XA, XB, YA, YB, ZA, ZB)
 -- Q =: triplequad (F, XA, XB, YA, YB, ZA, ZB, TOL)
 -- Q =: triplequad (F, XA, XB, YA, YB, ZA, ZB, TOL, QUADF)
 -- Q =: triplequad (F, XA, XB, YA, YB, ZA, ZB, TOL, QUADF, ...)
     Numerically evaluate the triple integral of F.

     F is a function handle, inline function, or string containing the name of the function to evaluate.  The function F must have the form w = f(x,y,z) where either X or Y is a vector and the remaining inputs are scalars.  It should return a vector of the same length and orientation as X or Y.

     XA, YA, ZA and XB, YB, ZB are the lower and upper limits of integration for x, y, and z respectively.  The underlying integrator determines whether infinite bounds are accepted.

     The optional argument TOL defines the absolute tolerance used to integrate each sub-integral.  The default value is 1e-6.

     The optional argument QUADF specifies which underlying integrator function to use.  Any choice but quad is available and the default is quadcc.

     Additional arguments, are passed directly to F.  To use the default value for TOL or QUADF one may pass ':' or an empty matrix ([]). integral3, integral2, dblquad, quad, quadv, quadl, quadgk, quadcc, trapz

The recursive algorithm for quadrature presented above is referred to as "iterated".  A separate 2-D integration method is implemented in the function quad2d.  This function performs a "tiled" integration by subdividing the integration domain into rectangular regions and performing separate integrations over those domains.  The domains are further subdivided in areas requiring refinement to reach the desired numerical accuracy.  For certain functions this method can be faster than the 2-D iteration used in the other functions above.

 -- Q =: quad2d (F, XA, XB, YA, YB)
 -- Q =: quad2d (F, XA, XB, YA, YB, PROP, VAL, ...)
 -- [Q, ERR, ITER] =: quad2d (...)

     Numerically evaluate the two-dimensional integral of F using adaptive quadrature over the two-dimensional domain defined by XA, XB, YA, YB using tiled integration.  Additionally, YA and YB may be scalar functions of X, allowing for the integration over non-rectangular domains.

     F is a function handle, inline function, or string containing the name of the function to evaluate.  The function F must be of the form z = f(x,y), and all operations must be vectorized such that X and Y accept array inputs and return array outputs of the same size. (It can be assumed that X and Y will either be same-size arrays or one will be a scalar.)  The underlying integrators will input arrays of integration points into F and/or use internal vector expansions to speed computation that can produce unpredictable results if F is not restricted to elementwise operations.  For integrands where this is unavoidable, the ("Vectorized") option described below may produce more reliable results.

     Additional optional parameters can be specified using "PROPERTY", VALUE pairs.  Valid properties are:

     AbsTol
          Define the absolute error tolerance for the quadrature.  The default value is 1e-10 (1e-5 for single).

     RelTol
          Define the relative error tolerance for the quadrature.  The default value is 1e-6 (1e-4 for single).

     MaxFunEvals
          The maximum number of function calls to the vectorized function F. The default value is 5000.

     Singular
          Enable/disable transforms to weaken singularities on the edge of the integration domain.  The default value is TRUE.

     Vectorized
          Enable or disable vectorized integration.  A value of false forces Octave to use only scalar inputs when calling the integrand, which enables integrands f(x,y) that have not been vectorized or only accept scalar values of X or Y.  The default value is true. Note that this is achieved by wrapping f(x,y) with the function arrayfun, which may significantly decrease computation speed.

     FailurePlot
          If quad2d fails to converge to the desired error tolerance before MaxFunEvals is reached, a plot of the areas that still need refinement is created.  The default value is FALSE.

     Adaptive quadrature is used to minimize the estimate of error until the following is satisfied:

                  ERROR <= max (ABSTOL, RELTOL*|Q|)

     The optional output ERR is an approximate bound on the error in the integral abs (Q - I), where I is the exact value of the integral.  The optional output ITER is the number of vectorized function calls to the function F that were used.

     Example 1 : integrate a rectangular region in x-y plane

          F = @(X,Y) 2*ones (size (X));
          Q = quad2d (F, 0, 1, 0, 1)
            => Q =  2

     The result is a volume, which for this constant-value integrand, is just LENGTH * WIDTH * HEIGHT.

     Example 2 : integrate a triangular region in x-y plane

          F = @(X,Y) 2*ones (size (X));
          YMAX = @(X) 1 - X;
          Q = quad2d (F, 0, 1, 0, YMAX)
            => Q =  1

     The result is a volume, which for this constant-value integrand F = 2, is the Triangle Area x Height or 1/2 * BASE * WIDTH * HEIGHT.

     Example 3 : integrate a non-vectorized function over a square region

          F = @(X,Y) sinc (X) * sinc (Y));
          Q = quad2d (F, -1, 1, -1, 1)
            => Q =  12.328  (incorrect)
          Q = quad2d (F, -1, 1, -1, 1, "Vectorized", false)
            => Q =  1.390 (correct)
          F = @(X,Y) sinc (X) .* sinc (Y);
          Q = quad2d (F, -1, 1, -1, 1)
            => Q =  1.390  (correct)

     The first result is incorrect as the non-elementwise operator between the sinc functions in F create unintended matrix multiplications between the internal integration arrays used by quad2d.  In the second result, setting "Vectorized" to false forces quad2d to perform scalar internal operations to compute the integral, resulting in the correct numerical result at the cost of about a 20x increase in computation time. In the third result, vectorizing the integrand F using the elementwise multiplication operator gets the correct result without increasing computation time.

     Programming Notes: If there are singularities within the integration region it is best to split the integral and place the singularities on the boundary.

     Known MATLAB incompatibility: If tolerances are left unspecified, and any integration limits are of type single, then Octave's integral functions automatically reduce the default absolute and relative error tolerances as specified above.  If tighter tolerances are desired they must be specified.  MATLAB leaves the tighter tolerances appropriate for double inputs in place regardless of the class of the integration limits.

     Reference: L.F. Shampine, MATLAB program for quadrature in 2D, Applied Mathematics and Computation, pp. 266--274, Vol 1, 2008.

     integral2, dblquad, integral, quad, quadgk, quadv, quadl, quadcc, trapz, integral3, triplequad

Finally, the functions integral2 and integral3 are provided as general 2-D and 3-D integration functions.  They will auto-select between iterated and tiled integration methods and, unlike dblquad and triplequad, will work with non-rectangular integration domains.

 -- Q =: integral2 (F, XA, XB, YA, YB)
 -- Q =: integral2 (F, XA, XB, YA, YB, PROP, VAL, ...)
 -- [Q, ERR] =: integral2 (...)

     Numerically evaluate the two-dimensional integral of F using adaptive quadrature over the two-dimensional domain defined by XA, XB, YA, YB (scalars may be finite or infinite).  Additionally, YA and YB may be scalar functions of X, allowing for integration over non-rectangular domains.

     F is a function handle, inline function, or string containing the name of the function to evaluate.  The function F must be of the form z = f(x,y), and all operations must be vectorized such that X and Y accept array inputs and return array outputs of the same size. (It can be assumed that X and Y will either be same-size arrays or one will be a scalar.)  The underlying integrators will input arrays of integration points into F and/or use internal vector expansions to speed computation that can produce unpredictable results if F is not restricted to elementwise operations.  For integrands where this is unavoidable, the ("Vectorized") option described below may produce more reliable results.

     Additional optional parameters can be specified using "PROPERTY", VALUE pairs.  Valid properties are:

     AbsTol
          Define the absolute error tolerance for the quadrature.  The default value is 1e-10 (1e-5 for single).

     RelTol
          Define the relative error tolerance for the quadrature.  The default value is 1e-6 (1e-4 for single).

     Method
          Specify the two-dimensional integration method to be used, with valid options being "auto" (default), "tiled", or "iterated".  When using "auto", Octave will choose the "tiled" method unless any of the integration limits are infinite.

     Vectorized
          Enable or disable vectorized integration.  A value of false forces Octave to use only scalar inputs when calling the integrand, which enables integrands f(x,y) that have not been vectorized or only accept scalar values of X or Y.  The default value is true. Note that this is achieved by wrapping f(x,y) with the function arrayfun, which may significantly decrease computation speed.

     Adaptive quadrature is used to minimize the estimate of error until the following is satisfied:

                  ERROR <= max (ABSTOL, RELTOL*|Q|)

     ERR is an approximate bound on the error in the integral abs (Q - I), where I is the exact value of the integral.

     Example 1 : integrate a rectangular region in x-y plane

          F = @(X,Y) 2*ones (size (X));
          Q = integral2 (F, 0, 1, 0, 1)
            => Q =  2

     The result is a volume, which for this constant-value integrand, is just LENGTH * WIDTH * HEIGHT.

     Example 2 : integrate a triangular region in x-y plane

          F = @(X,Y) 2*ones (size (X));
          YMAX = @(X) 1 - X;
          Q = integral2 (F, 0, 1, 0, YMAX)
            => Q =  1

     The result is a volume, which for this constant-value integrand F = 2, is the Triangle Area x Height or 1/2 * BASE * WIDTH * HEIGHT.

     Example 3 : integrate a non-vectorized function over a square region

          F = @(X,Y) sinc (X) * sinc (Y));
          Q = integral2 (F, -1, 1, -1, 1)
            => Q =  12.328  (incorrect)
          Q = integral2 (F, -1, 1, -1, 1, "Vectorized", false)
            => Q =  1.390 (correct)
          F = @(X,Y) sinc (X) .* sinc (Y);
          Q = integral2 (F, -1, 1, -1, 1)
            => Q =  1.390  (correct)

     The first result is incorrect as the non-elementwise operator between the sinc functions in F create unintended matrix multiplications between the internal integration arrays used by integral2.  In the second result, setting "Vectorized" to false forces integral2 to perform scalar internal operations to compute the integral, resulting in the correct numerical result at the cost of about a 20x increase in computation time. In the third result, vectorizing the integrand F using the elementwise multiplication operator gets the correct result without increasing computation time.

     Programming Notes: If there are singularities within the integration region it is best to split the integral and place the singularities on the boundary.

     Known MATLAB incompatibility: If tolerances are left unspecified, and any integration limits are of type single, then Octave's integral functions automatically reduce the default absolute and relative error tolerances as specified above.  If tighter tolerances are desired they must be specified.  MATLAB leaves the tighter tolerances appropriate for double inputs in place regardless of the class of the integration limits.

     Reference: L.F. Shampine, MATLAB program for quadrature in 2D, Applied Mathematics and Computation, pp. 266--274, Vol 1, 2008.

     quad2d, dblquad, integral, quad, quadgk, quadv, quadl, quadcc, trapz, integral3, triplequad

 -- Q =: integral3 (F, XA, XB, YA, YB, ZA, ZB)
 -- Q =: integral3 (F, XA, XB, YA, YB, ZA, ZB, PROP, VAL, ...)

     Numerically evaluate the three-dimensional integral of F using adaptive quadrature over the three-dimensional domain defined by XA, XB, YA, YB, ZA, ZB (scalars may be finite or infinite).  Additionally, YA and YB may be scalar functions of X and ZA, and ZB maybe be scalar functions of X and Y, allowing for integration over non-rectangular domains.

     F is a function handle, inline function, or string containing the name of the function to evaluate.  The function F must be of the form z = f(x,y,z), and all operations must be vectorized such that X, Y, and Z accept array inputs and return array outputs of the same size.  (It can be assumed that X, Y, and Z will either be same-size arrays or scalars.)  The underlying integrators will input arrays of integration points into F and/or use internal vector expansions to speed computation that can produce unpredictable results if F is not restricted to elementwise operations.  For integrands where this is unavoidable, the ("Vectorized") option described below may produce more reliable results.

     Additional optional parameters can be specified using "PROPERTY", VALUE pairs.  Valid properties are:

     AbsTol
          Define the absolute error tolerance for the quadrature.  The default value is 1e-10 (1e-5 for single).

     RelTol
          Define the relative error tolerance for the quadrature.  The default value is 1e-6 (1e-4 for single).

     Method
          Specify the two-dimensional integration method to be used, with valid options being "auto" (default), "tiled", or "iterated".  When using "auto", Octave will choose the "tiled" method unless any of the integration limits are infinite.

     Vectorized
          Enable or disable vectorized integration.  A value of false forces Octave to use only scalar inputs when calling the integrand, which enables integrands f(x,y,z) that have not been vectorized or only accept scalar values of X, Y, or Z.  The default value is true.  Note that this is achieved by wrapping f(x,y,z) with the function arrayfun, which may significantly decrease computation speed.

     Adaptive quadrature is used to minimize the estimate of error until the following is satisfied:

                  ERROR <= max (ABSTOL, RELTOL*|Q|)

     ERR is an approximate bound on the error in the integral abs (Q - I), where I is the exact value of the integral.

     Example 1 : integrate over a rectangular volume

          F = @(X,Y,Z) ones (size (X));
          Q = integral3 (F, 0, 1, 0, 1, 0, 1)
            => Q =  1.00000

     For this constant-value integrand, the result is a volume which is just LENGTH * WIDTH * HEIGHT.

     Example 2 : integrate over a spherical volume

          F = @(X,Y) ones (size (X));
          YMAX = @(X) sqrt (1 - X.^2);
          ZMAX = @(X,Y) sqrt (1 - X.^2 - Y.^2);
          Q = integral3 (F, 0, 1, 0, YMAX, 0, ZMAX)
            => Q =  0.52360

     For this constant-value integrand, the result is a volume which is 1/8th of a unit sphere or 1/8 * 4/3 * pi.

     Example 3 : integrate a non-vectorized function over a cubic volume

          F = @(X,Y) sinc (X) * sinc (Y), * sinc (Z);
          Q = integral3 (F, -1, 1, -1, 1, -1, 1)
            => Q =  14.535  (incorrect)
          Q = integral3 (F, -1, 1, -1, 1, -1, 1, "Vectorized", false)
            => Q =  1.6388 (correct)
          F = @(X,Y,Z) sinc (X) .* sinc (Y), .* sinc (Z);
          Q = integral3 (F, -1, 1, -1, 1, -1, 1)
            => Q =  1.6388  (correct)

     The first result is incorrect as the non-elementwise operator between the sinc functions in F create unintended matrix multiplications between the internal integration arrays used by integral3.  In the second result, setting "Vectorized" to false forces integral3 to perform scalar internal operations to compute the integral, resulting in the correct numerical result at the cost of about a 30x increase in computation time. In the third result, vectorizing the integrand F using the elementwise multiplication operator gets the correct result without increasing computation time.

     Programming Notes: If there are singularities within the integration region it is best to split the integral and place the singularities on the boundary.

     Known MATLAB incompatibility: If tolerances are left unspecified, and any integration limits are of type single, then Octave's integral functions automatically reduce the default absolute and relative error tolerances as specified above.  If tighter tolerances are desired they must be specified.  MATLAB leaves the tighter tolerances appropriate for double inputs in place regardless of the class of the integration limits.

     Reference: L.F. Shampine, MATLAB program for quadrature in 2D, Applied Mathematics and Computation, pp. 266--274, Vol 1, 2008.

     triplequad, integral, quad, quadgk, quadv, quadl, quadcc, trapz, integral2, quad2d, dblquad

The above integrations can be fairly slow, and that problem increases exponentially with the dimensionality of the integral.  Another possible solution for 2-D integration is to use Orthogonal Collocation as described in the previous section (see Orthogonal Collocation).  The integral of a function f(x,y) for x and y between 0 and 1 can be approximated using n points by the sum over i=1:n and j=1:n of q(i)*q(j)*f(r(i),r(j)), where q and r is as returned by colloc (n).  The generalization to more than two variables is straight forward.  The following code computes the studied integral using n=8 points.

     f = @(x,y) sin (pi*x*y') .* sqrt (x*y');
     n = 8;
     [t, ~, ~, q] = colloc (n);
     I = q'*f(t,t)*q;
           => 0.30022

It should be noted that the number of points determines the quality of the approximation.  If the integration needs to be performed between a and b, instead of 0 and 1, then a change of variables is needed.